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.

[bernalde]

Maintainer review of head 43c235c.

The math is equivalence-preserving. Multiplying the Raoult constraints through by P (with P in [1.05, 1.2] > 0) and expressing the 1.5 power as tc_gap*sqrt(tc_gap) with new bounded auxiliaries Tr and TcGap is algebraically equivalent to the base form, and the bounds are exact interval images of monotone maps over T in [Tlo, Tup], so nothing feasible is cut. Physically the sqrt singularity at 0 is unreachable: min TC = 512.6 K (methanol) against Tup ≈ 415–420 K keeps TcGap's lower bound ≈ 0.17. The new tests/test_kaibel.py is genuine — the hull/GAMS-writer test zeroes disaggregated variables and evaluates every active constraint body, which raises on base (TC/T division at T=0) and passes at head.

Findings:

Question: feasible-region widening — total_flow_max = max(flow_max, max Ltotal0, max Vtotal0) raises the Vtotal/Ltotal upper bound from 1000 to 1600 mol/s so the initialization fits inside bounds. This enlarges the feasible region relative to the published benchmark, so optimal objectives may not be comparable across gdplib versions. Was relaxing the bound intended, versus scaling the initial totals down to respect the original flow_max? If intended, record it in the model README or module docstring.

Nonblocking: stale size report — gdplib/kaibel/model_size_report.md is untouched but now stale: the PR adds 800 continuous variables (Tr, TcGap: 2 x 4 sections x 25 trays x 4 components) and 800 constraints. Regenerate the per-model size report (and the root README table row) with generate_model_size_report.py.

Nonblocking: naming/docs — m.Tr is documented as "Reduced temperature term" but is TC/T, the inverse reduced temperature; and the model-level constraint names _reduced_temperature_term/_critical_temperature_gap carry leading underscores, unconventional for public Pyomo components. Fix the docstrings at minimum; renaming the components is optional.

Nonblocking: xor no-op wording — xor=True is Pyomo's default; the body presents "declares tray disjunctions as XOR" as part of the fix, which overstates a no-op. Keep the explicitness, correct the narrative.

Issue linkage: issue #68's seven-strategy checklist is unchecked and the issue targets the defunct benchmark_zd branch. This PR fixes the hull/GAMS-writer blocker with smoke coverage but does not demonstrate the full matrix; treat "Closes #68" as Refs #68 unless the remaining checklist is explicitly superseded.

Summary — Blocking: 0. Nonblocking: 3. Questions: 1. Tests: static verification of the algebra, bound images, and physical data ranges; CI is green at this head (6/6 checks). Merge-ready once the size report is regenerated and the bound-widening question is answered. Posting as COMMENT (author account); the main ruleset requires one approving review from another maintainer.

[bernalde]

Addressed the review of head 43c235c in commit 73a98f0:

• flow-bound-widening (question): the widening is intentional — the total-flow upper bound rises from 1000 to 1600 mol/s only so the initialization totals fit inside bounds. This is now documented at total_flow_max in the model, including the caveat that optimal objectives are not directly comparable with earlier gdplib versions.
• stale-size-report: regenerated via generate_model_size_report.py — gdplib/kaibel/model_size_report.md and the kaibel row of the root README table now reflect the added Tr/TcGap variables and constraints. Every other model's report reproduced its committed content exactly, confirming the generator is deterministic at this head.
• tr-naming: m.Tr's doc and the defining constraint's doc now read "Inverse reduced temperature TC/T"; component renames were left alone as optional.
• xor-noop-136: recorded here rather than editing the body — the xor=True edit is an explicit no-op (Pyomo's default), kept for clarity; it is not part of the behavioral fix.

Note on the edit hygiene: the model file has mixed line endings on main; the address commit preserves the original per-line endings, so the diff is the 12-line semantic change only.

Issue linkage reminder from the review stands: treat Closes #68 as Refs #68 unless the remaining strategy checklist is explicitly superseded.

Tests: pixi run pytest tests/test_kaibel.py — 5 passed; black --check clean. CI for the new head was starting when this was posted; the previous head was green (6/6).

[bernalde]

Correction to the previous address pass (supersedes the regeneration claim in the comment above): the first regeneration ran on this branch's pre-pandemic checkout, so it reproduced the old README table format and made the PR CONFLICTING against main's rewritten table — the exact hand-merge trap flagged on PR #140. Fixed in commits 7e3373f (append-only merge of origin/main) and 9c62d3f: the table is now regenerated on the merged tree, so the kaibel columns update inside main's current table with the pandemic/grid/multiperiod_blending columns preserved. pixi run pytest tests/test_kaibel.py — 5 passed at the merged head.