GDPopt GLOA convex affine cuts use concave MC++ slopes

[bernalde]

Summary

GDPopt's gdpopt.gloa affine cut generation appears to build the convex upper affine cut with the concave MC++ slope vector.

On current Pyomo main and Pyomo 6.10.0, pyomo/contrib/gdpopt/gloa.py computes both slope vectors:

ccSlope = mc_eqn.subcc()
cvSlope = mc_eqn.subcv()
ccStart = mc_eqn.concave()
cvStart = mc_eqn.convex()

but then builds one cut_body with ccSlope and reuses it for both cuts:

cut_body = sum(
    ccSlope[var] * (var - var.value)
    for var in vars_in_constr
    if not var.fixed
)
concave_cut = cut_body + ccStart >= lb_int
convex_cut = cut_body + cvStart <= ub_int

I believe the convex cut should use cvSlope, not ccSlope. The current code can generate a convex cut that excludes feasible points and causes GLOA to terminate with an incorrect incumbent.

This was found while investigating GDPlib's CSTR benchmark:

• GDPlib issue: Refactor cstr to fix benchmark issues SECQUOIA/gdplib#63
• Related GDPlib follow-up: Plan PyPI release workflow and Pixi platform support SECQUOIA/gdplib#104
• GDPlib workaround PR: Fix CSTR GLOA formulation SECQUOIA/gdplib#122
• Diagnostic comments:
  • Fix CSTR GLOA formulation SECQUOIA/gdplib#122 (comment)
  • Fix CSTR GLOA formulation SECQUOIA/gdplib#122 (comment)

Steps to reproduce the issue

This minimal script does not require a solver. It shows that MC++ returns different subcc() and subcv() slopes for a bilinear equality from the CSTR model, and that reusing ccSlope for the convex upper cut rejects a point satisfying the original equality.

$ python gloa_cv_slope_mwe.py

# gloa_cv_slope_mwe.py
import pyomo.environ as pyo
from pyomo.contrib.mcpp.pyomo_mcpp import McCormick

m = pyo.ConcreteModel()
m.P = pyo.Var(bounds=(0, 10), initialize=0.05)
m.Q = pyo.Var(bounds=(1, 10), initialize=1.000002190520329)
m.F = pyo.Var(bounds=(0, 10), initialize=0.050000453446344)
m.QP = pyo.Var(bounds=(0, 10), initialize=1.0)

# This is the CSTR constraint split_add[A]:
# P[A] * Q[1] - F[A, 1] * QP == 0
expr = m.P * m.Q - m.F * m.QP

mc = McCormick(expr)
cc_slope = mc.subcc()
cv_slope = mc.subcv()
cc_start = mc.concave()
cv_start = mc.convex()

incumbent = {v.name: pyo.value(v) for v in (m.P, m.Q, m.F, m.QP)}

print("incumbent expr:", pyo.value(expr))
print("cc slopes:", {v.name: cc_slope[v] for v in (m.P, m.Q, m.F, m.QP)})
print("cv slopes:", {v.name: cv_slope[v] for v in (m.P, m.Q, m.F, m.QP)})
print("ccStart:", cc_start)
print("cvStart:", cv_start)

# A feasible point for P*Q - F*QP == 0 in the same variable bounds.
feasible_point = {
    "P": 0.05,
    "Q": 1.1,
    "F": 0.055,
    "QP": 1.0,
}
for v in (m.P, m.Q, m.F, m.QP):
    v.set_value(feasible_point[v.name])

bad_convex_cut_body = (
    sum(
        cc_slope[v] * (pyo.value(v) - incumbent[v.name])
        for v in (m.P, m.Q, m.F, m.QP)
    )
    + cv_start
)
correct_convex_cut_body = (
    sum(
        cv_slope[v] * (pyo.value(v) - incumbent[v.name])
        for v in (m.P, m.Q, m.F, m.QP)
    )
    + cv_start
)

print("feasible point expr:", pyo.value(expr))
print("current GLOA convex cut body, using cc slopes:", bad_convex_cut_body)
print("convex cut body using cv slopes:", correct_convex_cut_body)

# The current GLOA cut is body <= 0. This assertion demonstrates the bug.
assert abs(pyo.value(expr)) <= 1e-12
assert bad_convex_cut_body > 0
assert correct_convex_cut_body <= 0

Observed output:

incumbent expr: -3.4392032755015123e-07
cc slopes: {'P': 1.0, 'Q': 10.0, 'F': 0.0, 'QP': 0.0}
cv slopes: {'P': 1.0, 'Q': 0.0, 'F': -10.0, 'QP': 0.0}
ccStart: 0.05002190520328931
cvStart: -0.45000453446343996
feasible point expr: 6.938893903907228e-18
current GLOA convex cut body, using cc slopes: 0.5499735603332718
convex cut body using cv slopes: -0.5

Because the convex cut is body <= 0, the current ccSlope-based cut excludes a feasible point. The same cut using cvSlope is satisfied.

Error Message

There is no Python exception. The incorrect behavior is an invalid affine cut and, in the full GDPopt run, incorrect convergence.

Full GDPlib evidence:

• Original CSTR formulation with gdpopt.gloa, GAMS/BARON profile:
  • reports optimal objective 3.1301824793863147;
  • stops after evaluating only a small subset of the logical configurations.
• Direct transformed Big-M CSTR solve with BARON gives global objective about 3.0620145766.
• When GLOA affine cuts are monkeypatched out, no-good cuts enumerate the full CSTR logical space and find the global solution.
• When GLOA is monkeypatched to use separate cut bodies:
  • concave lower cut with ccSlope;
  • convex upper cut with cvSlope;
    the original CSTR formulation returns the global solution:

Optimal objective value 3.0620115146
active CSTRs: 1,2,3,4,5
recycle: 5
YF[5] = True

I also evaluated GLOA-generated affine cuts from the original CSTR run at a BARON-solved globally feasible point. The point violates generated GDPopt_aff cuts immediately. The first source-labeled invalid cut came from the same constraint used in the MWE:

SOURCE split_add[A]
EXPR P[A]*Q[1] - F[A,1]*QP == 0
KIND convex_ub
VIOLATION: upper bound 0.0, body 0.5528443771764309

Information on your system

Pyomo version: 6.10.0
Python version: 3.12.13 | packaged by conda-forge | (main, Mar 5 2026, 16:50:00) [GCC 14.3.0]
Operating system: Linux-6.6.114.1-microsoft-standard-WSL2-x86_64-with-glibc2.39
How Pyomo was installed: conda-forge environment managed by Pixi
Solver, if applicable: MWE does not require a solver; full CSTR diagnosis used the Pyomo GAMS solver interface with BARON.

Additional information

The current Pyomo main source still appears to reuse ccSlope for the convex cut in pyomo/contrib/gdpopt/gloa.py:

ccSlope = mc_eqn.subcc()
cvSlope = mc_eqn.subcv()
...
cut_body = sum(ccSlope[var] * (var - var.value) ...)
...
convex_cut = cut_body + cvStart <= ub_int

This is why I think the issue belongs in Pyomo's GLOA affine cut generation, not BARON or MC++:

• BARON is not needed for the MWE.
• MC++ provides the distinct cvSlope that makes the convex cut valid.
• A local GLOA monkeypatch using cvSlope for the convex cut fixes the incorrect CSTR result.