Multiple Masters: A Faster Decomposition for Stochastic Expansion

Expansion planning under uncertainty has a shape that has been understood for decades. You decide today what to build (which units, which reinforcements, which contracts) and you find out afterwards what demand, hydrology, fuel prices and renewable output actually did. Write every scenario out explicitly and the whole thing is one large linear program. Decomposition is how you avoid solving it that way.

Classical Benders decomposition does the job and can crawl getting there. On the case below, a Brazilian system model with 176 live candidate projects, multi-cut Benders needs 465 iterations. Accelerated Benders Decomposition with Multiple Masters (a-BDMM), introduced by Soares, Street, Andrade and Garcia in IEEE Transactions on Power Systems [1] for this exact problem, needs 11.

This post is about what a-BDMM changes to get there, and about implementing it as LightBDMM.jl with a coding agent. That combination suited the problem, because the parts of this algorithm that are easy to get wrong are not the parts where a re-implementation usually breaks.

Why one cut per iteration is slow

The two-stage stochastic program is

$$\min_{x \in X} \; c^\top x + \sum_{s} p_s \, Q_s(x), \qquad Q_s(x) = \min_{y_s \ge 0} \left\{ q_s^\top y_s \;:\; W_s y_s \ge h_s - T_s x \right\}$$

with $x$ the investment decision, $Q_s$ the cost of operating that decision under scenario $s$, and $p_s$ the scenario probabilities. $Q_s$ is convex in $x$, which is what makes cutting-plane methods work: solve the recourse problems at a trial point $x^k$, collect subgradients, and add inequalities that under-estimate the recourse cost everywhere.

The classical single-cut method adds exactly one of those inequalities per iteration,

$$\theta \;\ge\; \sum_s p_s \left[ Q_s(x^k) + g_s^{k\top} \left( x - x^k \right) \right],$$

averaging every scenario’s subgradient $g_s^k$ into a single constraint. That average is where the information goes. The master learns that the recourse cost is higher than it thought, but not which scenarios drove it, so each iteration buys one hyperplane and the bound creeps.

The multi-cut variant keeps the scenarios separate, carrying one $\theta_s$ per scenario and adding $S$ inequalities per iteration instead of one. Same subproblems, same work per iteration, far more information retained. It is the sensible baseline, and it is what we compare against throughout.

One master per scenario

BDMM pushes the same idea one step further. Rather than a single master fed by every scenario, it maintains one master per scenario, each anchored on a different recourse problem, all sharing a common cut pool.

The point is where the trial points come from. A single master produces one $x^k$ per iteration, so the cut pool grows along one trajectory through $X$. With $S$ masters, each anchored differently, an iteration produces $S$ distinct trial points, and the recourse problems get evaluated at all of them. The shared pool therefore accumulates subgradients sampled across the decision space rather than along a single path, and the recourse approximation tightens much faster per iteration.

The cost is honest and obvious: an iteration now solves $S$ masters instead of one, each embedding a full dispatch problem. Those masters are also the natural unit to distribute: in the runs below, one MPI rank per scenario plus a controller.

The consensus penalty

Multiple masters create a problem of their own. Nothing so far forces the $S$ masters to agree, and an investment plan is a single plan: you cannot build one portfolio per scenario. Left alone, the masters wander apart and the shared cut pool ends up describing regions no consensus solution will ever visit.

The acceleration in a-BDMM is a Progressive Hedging term that pulls them together. Each scenario master solves

$$\min_{x_s \in X, \, \theta_s} \; c^\top x_s + \theta_s + w_s^\top \left( x_s - \bar{x} \right) + \frac{\rho}{2} \left\lVert x_s - \bar{x} \right\rVert^2$$

where $\bar{x} = \sum_s p_s x_s$ is the current consensus point and $w_s$ are multipliers updated as $w_s \leftarrow w_s + \rho\left(x_s - \bar{x}\right)$. The linear term prices disagreement; the quadratic term penalizes it directly.

Two places it is easy to get wrong

Both of the algorithm’s genuine subtleties are consequences of that penalty, and neither is visible in the description above.

The lower bound needs its own LP. Once the quadratic term is active, the master objective is no longer a relaxation of anything. The penalty $\tfrac{\rho}{2}\lVert x_s - \bar{x}\rVert^2$ is an artifact of the augmentation (it appears nowhere in the original stochastic program), so a master carrying it does not under-estimate the true optimum and its objective is not a valid bound. The bound has to come from a separate solve that keeps the Lagrangian term $w_s^\top(x_s - \bar{x})$ and drops the quadratic one. LightBDMM.jl maintains a second pool of lower-bound masters for exactly this.

A degenerate start can converge to the wrong answer. Initialize every component of $x$ at zero and the consensus penalty immediately rewards agreement around that uninformed point. The masters agree, the convergence test is satisfied, and the answer is wrong. Agreement must not outrank information: the fix is a warmup phase in which the decomposition explores freely for several iterations before the penalties switch on.

The case

A database of the Brazilian system in planning year 2050: 176 live candidate projects out of 221, ten hydro–wind–solar scenarios, and roughly 324,000 variables per scenario in the second stage. The horizon is one month at hourly resolution with operating cost scaled to an annual equivalent. All three algorithms come from LightBDMM.jl and differ only in algorithm; each ran under MPI on eleven ranks (one controller and one worker per scenario) to a 0.1% relative gap, with Xpress barrier and crossover disabled.

Each refinement earns its keep. Multi-cut Benders takes 465 iterations. Multiple masters take that to 39. The consensus penalty takes it to 11, a 42-fold reduction against the baseline. And because an a-BDMM iteration is expensive but there are so few of them, the ordering survives the translation into wall-clock: a-BDMM finished in about a third of the time multi-cut Benders took.

The right-hand panel is the more honest view, and it shows something the iteration count hides. Multi-cut Benders spends most of its run barely moving (it was still above a 90% gap after 27 iterations), then grinds down over hundreds of cheap iterations. a-BDMM’s first iteration already costs more than fifty of Benders’, and it is done in eleven.

The methods also do not agree on the answer, and they disagree in a consistent order. Measured against the monolithic deterministic equivalent of the same problem, a-BDMM lands 0.060% away, BDMM 0.184%, and multi-cut Benders 0.343%. The method that converges in the fewest iterations is also the one that gets closest.

That ordering deserves a caveat rather than a victory lap. Multi-cut Benders certified a tighter gap than a-BDMM, 0.033% against 0.048%, while sitting five times further from the reference. Its bound is tight around a worse incumbent, which is exactly what a stalled cutting-plane method looks like from the inside.

Building it with an agent

That an agent can write a Julia package is not interesting on its own any more. What is interesting is where the failures land when the subject matter is an algorithm rather than an application.

The package (structure, implementation, regression tests, documentation, continuous integration) was generated by Claude Code and corrected iteratively through execution and review. Our existing Benders and Progressive Hedging repositories went in as references for architectural patterns and testing conventions, not for the algorithm. The protocol was deliberately low-interaction: high-level objectives at the start of each cycle, review at the end, and no per-step steering. Later extensions, mixed-integer first-stage decisions with quadratic-penalty linearization and distributed execution over MPI, arrived the same way.

The first implementation already had the shape right, including the Benders workflow and a newsvendor benchmark. Then, in order:

• Cut aggregation was wrong: scenario contributions combined incorrectly, so objective values came out wrong. This is the one conventional software bug in the list, and it was found and fixed quickly.
• The degenerate-start trap appeared exactly as the theory says it should: everything initialized at zero, penalties pushing toward consensus around that point, the method settling into “do nothing” and stopping. The warmup phase emerged with little human intervention.
• The lower bound took much longer, and took real interaction. Using the upper bound as a proxy produced artificial convergence in the first iteration. Using the raw master objectives made the bound oscillate once the penalty activated. Subtracting the quadratic penalties directly sometimes produced $LB > UB$, which a correct implementation cannot reach.

That third one stopped being a software problem quickly. It became a question of which quantities remain valid bounds after the augmentation, how convergence should be measured, and which parts of the objective are artifacts rather than components of the original program, and it resolved only when the implementation recovered the procedure the paper actually describes.

Almost none of the failures in this project were syntactic. They were algorithmic, and the ones that mattered most surfaced only once the package met a real model. Terence Tao’s argument about machine assistance in mathematics, that the value lies in the overhead around research rather than in the hard part, describes this well: the agent was most useful where the work was mechanical and verifiable, and least useful exactly where the mathematics was subtle.

The harness is the reason it works

In a low-interaction protocol, validation is not a quality-assurance afterthought. It is the thing that makes the protocol viable, because it replaces human review of intermediate steps with an external check.

Every regression test is solved twice: once by LightBDMM.jl, and once as the monolithic deterministic equivalent, the extensive-form formulation of the same program, sharing none of the decomposition’s implementation. They are required to agree. That is also how anomalies get caught: not by reading code, but by two independent paths disagreeing about a number. The suite runs in continuous integration on every proposed change, and the package currently holds roughly 98% automated code coverage.

Where this leaves things

a-BDMM has been in the literature since 2022, and two-stage stochastic programming is one of the oldest frameworks in operations research. Nothing here is new mathematics, and this is not a demonstration of autonomous mathematical engineering.

What changed is the cost of the translation. A package with three decomposition variants, MPI execution, a validation harness, benchmarking infrastructure and documentation came together in a single development cycle, and then produced a result on a real model rather than a benchmark. That matters because the constraint on how much of the optimization literature actually gets used has never really been the mathematics. It has been the engineering effort to get a method from a paper into something a planning study can rely on.

References

[1] A. Soares, A. Street, T. Andrade, and J. D. Garcia, “An integrated progressive hedging and Benders decomposition with multiple master method to solve the Brazilian generation expansion problem,” IEEE Trans. Power Syst., vol. 37, no. 5, pp. 4017–4027, Sep. 2022, doi: 10.1109/TPWRS.2022.3141993. Preprint: arXiv:2108.03143.