Where the Proof Runs Out

TL;DR

• A model can be continuous and still be non-convex. Haverly's pooling problem has no integer variables; products of flow and concentration are enough to break convexity.
• On the same seven-variable pooling model, the same local SLSQP configuration returns feasible plans worth $0, $100 and $400 from three fixed starting points, reporting success each time. Its status reports satisfaction of local stopping tests, not global optimality.
• A global method needs a valid bound that covers territory it has not visited. For bilinear terms, McCormick inequalities provide convex relaxations over bounded boxes; spatial branch-and-bound tightens them by subdividing those boxes.
• The four McCormick inequalities describe the exact convex hull of one bounded bilinear graph. That hull is still a relaxation of the equality, so wide variable bounds can leave a weak proof.
• “Hard” ranges from expensive to undecidable. Jeroslow's impossibility result concerns a general class with unbounded integer variables and mixed-sign quadratic constraints. Explicitly bounded pure-integer models remain finite and decidable.

I can make the same local algorithm report success three times on the same small problem. One run earns nothing. Another earns $100. A third earns $400. Nothing crashed. Every variable is continuous. The only input I changed was the point from which the search began.

This is the fifth and final article in the series, and the one closest to the edge of my own practice. The calculation below uses fixed, recorded inputs and is checked against an independent global benchmark. The wider subject, non-convex global optimisation, is one I am learning rather than claiming as home ground. That distinction matters in an article about the limits of claims.

No integers required

In 1978, C. A. Haverly published a tiny refinery-planning example that became the canonical pooling problem.[1] Two feeds, A and B, enter a common pool. A third feed, C, can go directly to either of two products. The pool destroys the identity of its inputs: once A and B are mixed, every unit leaving has the same sulphur concentration.

The data fit in a paragraph. A contains 3 per cent sulphur and costs $6 per unit. B contains 1 per cent and costs $16. C contains 2 per cent and costs $10. Product X sells for $9, accepts at most 2.5 per cent sulphur and has demand capped at 100 units. Product Y sells for $15, accepts at most 1.5 per cent and has demand capped at 200.

Let a and b be the flows into the pool, p_X and p_Y its outputs, c_X and c_Y the direct flows from C, and q the pool's sulphur concentration. The model's constraints include:

The objective is revenue from X and Y minus the three feed costs. Every decision variable can take fractional values. The difficulty comes from the products q(a + b), qp_X and qp_Y. They multiply one decision by another.

A bilinear equality need not preserve feasibility under averaging: two feasible points can have an infeasible midpoint. The resulting feasible region is non-convex, so a stopping test at one point does not bound the rest.

This is also why convex MINLP is easy to misread. Its integer coordinates remain discrete; convex describes the nonlinear structure after those restrictions are relaxed. Haverly's original example goes the other way: it is a completely continuous, non-convex nonlinear program, or NLP. Integrality is absent and the global proof is already difficult.

What success meant

I solved that model with SciPy's SLSQP method from three fixed initial guesses selected from a nine-point sequence generated with seed 20260820. The three runs illustrate possible terminations; they do not estimate the relative sizes of attraction basins. The solver settings, equations and tolerances were otherwise identical, and each run returned a feasible plan with status success.

At $0, every flow is zero. The returned value q = 2.11234 is mathematically admissible but has no physical meaning because the pool is empty; this empty-pool degeneracy helps explain why the algebraic formulation admits a zero-profit termination. At $100, 50 units of A pass through the pool and meet 50 direct units of C in product X. At $400, 100 units of B pass through the pool and meet 100 direct units of C in product Y. That last blend contains exactly 1.5 per cent sulphur, fills Y's demand and earns 3,000 − 1,600 − 1,000 = $400.

The $400 result is independently documented as the global solution; the GAMS Model Library records the same value as −400 under its minimisation sign convention.[2] That global label does not come from SLSQP's status.

success has a narrower meaning than the chart's shared green badges suggest. It says only that SLSQP met its termination tests from that start.[3] It does not establish even strict local optimality, let alone a valid objective bound over the rest of the feasible region.[4] Finding the $400 plan once, or finding it in ninety-nine restarts, would be evidence. The bound is what makes it a certificate.

A bound that survives the starting point

Suppose a model contains w = xy, with finite bounds x ∈ [ℓ_x, u_x] and y ∈ [ℓ_y, u_y]. Garth P. McCormick showed how four linear inequalities can surround that product:[5]

For this isolated bilinear graph on that box, the resulting polytope is its exact convex hull. “Exact convex hull” does not mean that every point in the polytope satisfies w = xy. It is the tightest convex relaxation of the graph, and therefore still contains points the original equality forbids.

The quality of the bound depends on the box. If x and y range widely, the hull leaves ample room between the bilinear graph and its relaxation. Split either interval into smaller pieces and rebuild the four inequalities; the room contracts.

That is the central move in spatial branch-and-bound, a common framework for deterministic global optimisation. The algorithm builds valid convex relaxations over bounded subdomains, uses them to bound what the unexplored region could contain, and subdivides domains where the bound remains loose. Feasible points provide incumbents. Relaxations provide global bounds. Bound tightening, cuts and better convex underestimators make the same outline substantially more powerful in practice.[4][6]

The tree resembles the one in the MILP article. For this continuous Haverly model, its branches divide variable domains; a non-convex MINLP may branch on both integer choices and continuous domains. Global optimality is established when the incumbent and the valid bound meet within declared tolerances. If the run ends with a non-zero gap, the incumbent may be excellent. The proof remains incomplete.

Hard is not one word

The pooling problem is NP-hard. Haugland showed that this remains true with only one quality constraint at each terminal and, under a separate restriction, with at most two sources and two terminals.[7] This is a statement about how computational work can scale. It leaves open algorithms that solve every finite instance, sometimes at forbidding cost.

Robert Jeroslow proved something logically stronger in 1973: there is no algorithm that computes the optimum for every integer program in a class with a linear objective and quadratic constraints.[8] The scope belongs beside the headline. His construction uses integer domains without finite bounds and diagonal quadratic forms with coefficients of mixed signs. It does not make every quadratic integer model undecidable. If all variables in a pure-integer model have finite explicit bounds, only finitely many assignments exist; enumeration restores decidability, though not practicality.

Local, global, NP-hard and undecidable therefore answer different questions. A local optimiser may find the best point without proving it. A deterministic global optimiser can certify a supported, adequately bounded non-convex model when its valid relaxation bound closes to the incumbent within the chosen tolerances. A time limit with a remaining gap cannot. And in Jeroslow's general class, no universal algorithm exists to wait for.

When a solver reports a result, its brand name is less informative than two details: which algorithmic mode ran, and why it terminated.

The edge of my map

The thermodynamic metabolic formulations I have worked with stay on the MILP side of this boundary: continuous fluxes sit beside binary reaction-use and direction variables, but relaxing those binaries leaves a linear model. Nonlinear kinetic relationships would cross into the territory described here.

The boundary is relevant rather than hypothetical. k-OptForce, for example, incorporates nonlinear kinetic descriptions into strain-design optimisation and reformulates the resulting bilevel construction as a single-level MINLP.[9] Similar products, ratios and nonlinear balances appear in process networks, energy systems and engineering design.

The surveys do not offer one general trick that makes these structures easy.[4][10] They point instead to continuing work on tighter convex relaxations, structure-specific cuts, bound tightening and heuristics whose outputs make no global-optimality claim. The modelling question comes first: which nonlinear relationship is essential enough to justify the harder claim?

Where the word stops

Non-convexity does not erase global optimality. It changes what must be built before the word is earned. A local termination describes one search. A global certificate must also account for the rest of the feasible region.

The word optimal remains precise only for the stated model and assumptions, and only as far as its last valid bound. Beyond either boundary, a solution may still be useful or impressive. The run has not established the larger claim.

The detour is over. Markowitz has waited long enough.

Data and code

The Post 5 companion folder contains the complete Haverly data, fixed starting points, pinned library versions, recorded runtime, generated results, static figure and verification checks. It records SLSQP's local statuses separately from the independently documented global benchmark.

References

1. C. A. Haverly, “Studies of the Behavior of Recursion for the Pooling Problem”, ACM SIGMAP Bulletin 25 (1978): 19–28.
2. GAMS Development Corporation, “A Pooling Problem (POOL)”, GAMS Model Library, instance haverly1, accessed 20 August 2026.
3. SciPy Developers, “minimize(method='SLSQP')”, SciPy v1.18.0 Manual, accessed 20 August 2026.
4. Pietro Belotti, Christian Kirches, Sven Leyffer, Jeff Linderoth, James Luedtke and Ashutosh Mahajan, “Mixed-Integer Nonlinear Optimization”, Acta Numerica 22 (2013): 1–131.
5. Garth P. McCormick, “Computability of Global Solutions to Factorable Nonconvex Programs: Part I—Convex Underestimating Problems”, Mathematical Programming 10, no. 1 (1976): 147–175.
6. E. M. B. Smith and C. C. Pantelides, “A Symbolic Reformulation/Spatial Branch-and-Bound Algorithm for the Global Optimisation of Nonconvex MINLPs”, Computers & Chemical Engineering 23, nos. 4–5 (1999): 457–478.
7. Dag Haugland, “The Computational Complexity of the Pooling Problem”, Journal of Global Optimization 64, no. 2 (2016): 199–215.
8. Robert G. Jeroslow, “There Cannot Be Any Algorithm for Integer Programming with Quadratic Constraints”, Operations Research 21, no. 1 (1973): 221–224.
9. Anupam Chowdhury, Ali R. Zomorrodi and Costas D. Maranas, “k-OptForce: Integrating Kinetics with Flux Balance Analysis for Strain Design”, PLOS Computational Biology 10, no. 2 (2014): e1003487.
10. Samuel Burer and Adam N. Letchford, “Non-Convex Mixed-Integer Nonlinear Programming: A Survey”, Surveys in Operations Research and Management Science 17, no. 2 (2012): 97–106.