Fix Kaibel GDP reformulation failures

[bernalde]

Summary

• Refactors Kaibel vapor-composition equations to avoid variable division and fractional powers that break GDP Hull/GAMS expression evaluation.
• Adds bounded reduced-temperature and critical-temperature-gap terms for vapor-pressure expressions.
• Initializes tray GDP choices to the full-column design so Hull perspective expressions evaluate at safe initial values for the GAMS writer.
• Sets total liquid/vapor flow bounds consistently with initialized total flow values and declares tray disjunctions as XOR.
• Adds focused Kaibel regression tests for bounds, Big-M/Hull transformation, vapor-pressure expression form, and GAMS-style post-Hull expression evaluation.

Tests run

• pixi run pytest tests/test_kaibel.py -v --tb=short -> 5 passed
• pixi run pytest tests/test_module_imports.py -v --tb=short -> 68 passed
• pixi run test -> 289 passed, 1 skipped
• pixi run lint -> exit 0; Black clean, critical flake8 count 0, broader flake8 statistics are existing exit-zero findings, typos clean
• git diff --check and git diff --cached --check -> passed
• pixi run gdplib-benchmark run --cases-file /tmp/kaibel_issue68_hull_eval_fix.csv --run-id issue68_kaibel_hull_eval_fix2_20260512 --no-skip-existing --no-summary -> 1 result row, 0 failures; gdp.hull/gams-local/DICOPT reached GAMS and returned intermediateNonInteger under the 60 s limit with objective/upper bound 151820.17183655952

Notes

• Generated benchmark outputs were left under ignored benchmark result directories and were not committed.
• A Pixi editable-package metadata change in pixi.lock was left unstaged as environment churn.
• The short Hull/DICOPT probe confirms the expression-evaluation blocker is fixed; it is not an optimality certificate.

Closes #68

[bernalde]

5 minute Kaibel MINLP reformulation campaign after the Hull vapor-pressure evaluation fix.

Command:

/home/bernalde/.pixi/bin/pixi run gdplib-benchmark run --cases-file /tmp/kaibel_issue68_minlp_5min.csv --run-id issue68_kaibel_minlp_5min_after_hullfix_20260512 --no-skip-existing --no-summary

Campaign result: 6 result rows, 0 benchmark failures. This is the main improvement from the latest fix: the Hull reformulations no longer fail in Pyomo/GAMS expression evaluation.

Reformulation	Solver profile	GAMS solver	Termination	Solver-log interpretation
gdp.bigm	gams-local	DICOPT	optimal wrapper status after resource-limit message	Best integer solution 91561.861135; lower bound 57862.658306; about 36.8% gap using the wrapper bound/incumbent
gdp.hull	gams-local	DICOPT	intermediateNonInteger	Reached GAMS and ran ~309 s, but DICOPT log says no integer solution; relaxed NLP reported at 118352.385586
gdp.bigm	gams-gurobi	Gurobi	maxTimeLimit	Gurobi log says Solution count 0; best bound 59218.293360. Ignore wrapper upper_bound=92000.0 as an initialized/wrapper artifact
gdp.hull	gams-gurobi	Gurobi	maxTimeLimit	Gurobi log says Solution count 0; best bound 101857.427497. Wrapper reports upper_bound=92000.0, but the solver log has no incumbent and the bound/upper-bound pair is inconsistent
gdp.bigm	gams-baron	BARON	maxTimeLimit	BARON lower bound 49000.0; no feasible solution found
gdp.hull	gams-baron	BARON	infeasible	BARON reported infeasible quickly with lower bound 1000.0; no feasible solution

Model sizes from the result JSON:

• Big-M transformed model: 11,928 constraints, 4,512 variables, 178 integer variables.
• Hull transformed model: 31,788 constraints, 17,299 variables, 178 integer variables.

Bottom line: the Hull expression-evaluation blocker is fixed, so all three Hull solver settings now run through the benchmark path. It does not yet push the best feasible incumbent forward in 5 minutes. The only clear feasible integer incumbent in this campaign is Big-M/DICOPT at 91561.861135, which is worse than the earlier 5 minute Big-M/DICOPT incumbent around 89332.39.

Generated artifacts were left under ignored benchmark output paths:

• benchmark_runs/issue68_kaibel_minlp_5min_after_hullfix_20260512
• gdplib/kaibel/benchmark_result/issue68_kaibel_minlp_5min_after_hullfix_20260512

[bernalde]

5 minute GDPOpt campaign for Kaibel after the Hull/GAMS expression fixes, excluding complete enumeration.

Command:

/home/bernalde/.pixi/bin/pixi run gdplib-benchmark run --cases-file /tmp/kaibel_issue68_gdpopt_no_enum_5min.csv --run-id issue68_kaibel_gdpopt_no_enum_5min_20260512 --no-skip-existing --no-summary

Case file used the local GDPOpt GAMS role stack:

• NLP: gams/ipopth
• MIP: gams/gurobi
• MINLP: gams/dicopt
• local MINLP: gams/dicopt
• time limit: 300 s per algorithm
• excluded: gdpopt.enumerate

Campaign result: 4 result rows, 0 benchmark failures.

Strategy	Termination	Iterations	Time (s)	Lower bound	Upper bound / incumbent	Notes
gdpopt.loa	maxTimeLimit	16	303.020	50000.0	83622.020955	Best incumbent from this campaign; reported gap 40.207%
gdpopt.gloa	optimal	2	58.828	89411.129883	89411.129883	GDPOpt bounds converged/crossed; interpret conservatively because Kaibel is nonconvex and role solvers are local for NLP/MINLP
gdpopt.lbb	unknown	0	108.014	n/a	n/a	Log says LBB has known issues; no feasible solution retained by GDPOpt
gdpopt.ric	maxTimeLimit	23	300.659	50000.0	83622.020955	Same best incumbent as LOA; reported gap 40.207%

Bottom line: yes, GDPOpt helps this instance. Without enumeration, LOA/RIC pushed the best feasible incumbent to 83622.020955, which is better than the direct-reformulation 5 minute campaign (91561.861135 in the rerun, and also better than the earlier ~89332.39 Big-M/DICOPT incumbent). The GLOA 89411.129883 result is useful evidence but should not be treated as a rigorous global certificate without additional nonconvex/global verification.

Generated artifacts were left under ignored benchmark output paths:

• benchmark_runs/issue68_kaibel_gdpopt_no_enum_5min_20260512
• gdplib/kaibel/benchmark_result/issue68_kaibel_gdpopt_no_enum_5min_20260512

[bernalde]

Follow-up on the GDPOpt benchmark interpretation: the gdpopt.gloa optimal result from the no-enumeration campaign is not a valid certificate for this Kaibel run.

The contradiction is real:

• gdpopt.loa found a feasible incumbent 83622.020955 at iteration 13. The log shows the NLP subproblem objective 8.3622020955490123e+04, EXIT: Optimal Solution Found, constraint violation about 3.07e-09, and then GDPOpt updated the upper bound:

  13        subproblem   50000.00000   83622.02096     40.21%    246.77  *
  

• gdpopt.ric reached the same incumbent 83622.020955.

• gdpopt.gloa reported 89411.129883 as both lower and upper bound after only 2 iterations:

  1        subproblem   49000.00000   89411.12988     45.20%     57.77  *
  GDPopt exiting--bounds have converged or crossed.
  

For a minimization problem, a feasible incumbent at 83622.020955 proves that 89411.129883 cannot be a valid global lower bound or optimum.

What appears to be happening:

• The original GDPOpt log says the Kaibel model has 2,528 nonlinear constraints and 100 disjunctions.
• The benchmark role stack used gams/ipopth for GDPOpt NLP subproblems. Pyomo's GLOA default NLP solver is couenne; this run overrode that with a local NLP solver.
• GLOA's master/OA lower bound crossed the incumbent, and GDPOpt's generic bound-crossing logic then reported optimal. That bound is contradicted by the LOA/RIC feasible point, so it should be treated as an invalid certificate for this nonconvex/local-solver run.
• One of the GLOA NLP subproblem solves also ended with IPOPT EXIT: Error in step computation! at an almost feasible point before GDPOpt accepted the 89411.129883 incumbent.

Corrected interpretation of the no-enumeration GDPOpt run:

• Best feasible incumbent found: 83622.020955 from LOA/RIC.
• Best reliable lower bound from LOA/RIC in that run: 50000.0.
• GLOA's optimal/LB=UB=89411.129883 line should not be used as a guarantee.
• This strengthens the earlier caveat: on this nonconvex Kaibel instance, GDPOpt wrapper summaries must be checked against solver logs and cross-algorithm incumbent evidence.

[bernalde]

Follow-up on the GDPOpt/GLOA contradiction:

I instrumented the 5-minute Kaibel GLOA run to map the generated affine cut that caused the premature certificate.

Findings:

• Default gdpopt.gloa still reports optimal with LB = UB = 89411.12988260087, but the debug log shows this comes immediately after master preprocessing declares an affine cut infeasible:
  • Detected an infeasible constraint during FBBT: GDPopt_aff.GDPopt_aff_cons[1613]
  • GDPOpt then applies its generic "bounds have converged or crossed" path and force-sets the dual bound to the incumbent.
• The generated row maps back to the global Kaibel constraint _external_boilup_ratio:
  • source body: Btotal - (1 - bu)*Ltotal[1,2]
  • generated cut pair: GDPopt_aff_cons[1612], GDPopt_aff_cons[1613]
  • at cut generation, the NLP point has bu = 0.895921269464842, Btotal = 49.99999999559216, Ltotal[1,2] = 480.4055520133598
  • the discrete master has already tightened bu.ub to 0.8795650501693494, so the affine cut is built from a point outside the master bound box.
• This is not just the GAMS NLP solver role. As a control, I reran GLOA with mip_presolve=False. It no longer fails through the same FBBT row, but it still returns a false certificate:
  • termination: optimal
  • reported bounds: LB = UB = 90441.55867777138
  • log sequence: after three subproblem solves, the MILP discrete problem is declared infeasible and GDPOpt again exits through "bounds have converged or crossed."
• I reran LOA for 300 seconds with a feasible-incumbent callback to capture the better point. It reproduces the previous benchmark:
  • termination: maxTimeLimit
  • bounds: LB = 50000.0, UB = 83622.02095549012
  • final captured incumbent at iteration 13:
    • bu = 0.877737590487368
    • Btotal = 49.99999999500094
    • Ltotal[1,2] = 408.9564421099134
    • Vtotal[1,1] = 358.956442108772
    • _external_boilup_ratio residual: -3.0702764775014657e-09
    • _internal_boilup_ratio residual: 3.0701698960911017e-09

Conclusion: the GLOA optimal result should be treated as invalid for this nonconvex Kaibel instance. The evidence points to the GLOA affine-cut/master logic cutting off feasible regions, not to the fractional vapor-pressure expression fix or to a particular NLP solver choice. The reliable 5-minute GDPOpt evidence remains the LOA/RIC incumbent 83622.02095549012 with dual bound 50000.0; GLOA should not be reported as a global certificate here.

[bernalde]

Follow-up on whether Pyomo/pyomo#3939 explains the Kaibel GLOA failure:

Yes, it appears to address the invalid GLOA certificate we observed here.

The failing Kaibel row was a generated GLOA convex affine cut:

GDPopt_aff.GDPopt_aff_cons[1613]
source_constraint = _external_boilup_ratio
source_body = Btotal - (1 - bu)*Ltotal[1,2]

Default GLOA produced the false certificate:

optimal
LB = UB = 89411.12988260087

That happened after the discrete master declared the generated affine cut infeasible and GDPOpt exited through "bounds have converged or crossed."

I ran a local control monkeypatch matching the fix proposed in Pyomo #3939: build the concave lower cut with ccSlope, but build the convex upper cut with cvSlope. With that change, the 5-minute Kaibel GLOA run no longer gave the false optimal result:

maxTimeLimit
LB = 54270.33951053024
UB = 89332.67304439889

This does not make Kaibel globally solved in 5 minutes. It changes the result from an invalid certificate to a conservative open-gap result. The best 5-minute feasible incumbent from the no-enumeration GDPOpt campaign remains the LOA/RIC point:

UB = 83622.02095549012
LB = 50000.0

So for this PR, the Kaibel model fix is still valid, but the GLOA benchmark should be interpreted as affected by upstream Pyomo GLOA cut generation until Pyomo #3939 is fixed.

[bernalde]

Heads up from PR #102: while wrapping up the pandemic PR after the base update, I ran the repository Black check/fix and it touched gdplib/kaibel/kaibel_solve_gdp.py. The Kaibel change there is formatting/line-ending cleanup only, including the cleaned m.candidate_trays_main doc line; it is not intended as a Kaibel behavior change.

Relevant commits on #102:

• 8b5da8d
• bf1c1bf

This may be relevant when rebasing or reviewing this Kaibel PR so the formatting noise is not mistaken for part of the Kaibel model fix.