PlasmoBenders.jl: Graph-Based Benders Decomposition

• PlasmoBenders.jl is an open-source tool for graph-native Benders decomposition that leverages hierarchical OptiGraph models for scalable optimization.
• It integrates multi-parametric programming surrogates to generate exact Benders cuts, significantly reducing subproblem solution times.
• The framework supports distributed-memory execution, allowing use in shared-memory, HPC, and cloud environments via RemoteOptiGraph.

PlasmoBenders.jl is an open-source, graph-native Benders decomposition capability in the Plasmo.jl ecosystem. It operates directly on graph-structured optimization models defined with the OptiGraph abstraction, where nodes encode modular subproblems and edges encode algebraic couplings, and it has been presented both as a Julia implementation of graph-based Benders decomposition for hierarchical optimization and as a framework for accelerating Benders by embedding multi-parametric programming surrogates for optimization subproblems. In later work, the same decomposition logic is extended to distributed-memory execution through RemoteOptiGraph, so that the Benders meta-algorithm can run on shared-memory, HPC, and cloud settings within a unified graph-based API (Cole et al., 3 Jan 2025, Brahmbhatt et al., 1 Aug 2025, Cole et al., 18 Nov 2025).

1. Position within the Plasmo.jl modeling framework

The immediate conceptual background for PlasmoBenders.jl is Plasmo.jl’s graph-based view of optimization. In the original algebraic-graph formulation, an optimization model is represented as a hypergraph whose nodes carry component models and whose hyperedges impose linking constraints. This makes separability explicit: removing link constraints yields independent node subproblems, which is precisely the structural condition exploited by decomposition schemes such as Benders and Lagrangian methods. The same body of work also introduced a computing-graph abstraction in which nodes execute tasks over attributes and edges communicate attributes with delays, thereby supplying a stateful runtime for orchestrating centralized, distributed, or hierarchical algorithms (Jalving et al., 2018).

Within that ecosystem, PlasmoBenders.jl specializes the graph abstraction to Benders decomposition. The hierarchical-optimization formulation describes it as a Julia implementation of graph-based Benders decomposition, or gBD, for nested optimization problems represented by OptiGraph. The distributed-memory formulation describes it as the graph-native Benders decomposition capability in the Plasmo.jl ecosystem. The multi-parametric programming formulation describes it as an open-source Benders decomposition engine that operates on graph-structured optimization models defined with Plasmo.jl’s OptiGraph abstraction. Across these formulations, the common design objective is to expose and exploit problem structure by preserving graph-level information rather than flattening the model into a monolithic algebraic program (Cole et al., 3 Jan 2025, Brahmbhatt et al., 1 Aug 2025, Cole et al., 18 Nov 2025).

A recurrent technical point is that OptiGraph supports hierarchical nesting. Any OptiNode can embed an OptiGraph, so subgraphs may represent planning, scheduling, operations, scenario blocks, or time blocks at different granularities. This is significant because classical nested Benders methods were historically associated with time-line structures, whereas the graph-based formulation generalizes the same logic to arbitrary hierarchical graphs. This suggests that PlasmoBenders.jl is best understood not merely as a solver wrapper, but as a decomposition layer that is structurally coupled to graph-based model design (Cole et al., 3 Jan 2025).

2. Mathematical structure and graph-based Benders formulation

The generic OptiGraph model used by PlasmoBenders.jl is written as

$\min f\left(\{x_n\}_{n\in N(G)}\right)$

subject to local node constraints

$x_n \in \mathcal{X}_n,\quad n\in N(G),$

and linking constraints

$g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$

In an equivalent descriptive form, nodes store local variables, constraints, and objectives, while edges store linking constraints that couple variables living on two or more nodes. In the hypergraph language used by Plasmo.jl, these are undirected hyperedges because a single linking constraint may involve more than two nodes (Brahmbhatt et al., 1 Aug 2025, Cole et al., 18 Nov 2025).

For classical Benders decomposition, the package adopts the standard separation between complicating variables and recourse variables. In the two-stage LP form used in the multi-parametric paper, the master minimizes $c^\top x + \theta$ subject to $Ax \le b$, $x \in X$, and iteratively added Benders cuts. For a fixed master decision $x$, the subproblem is

$$Q(x) = \min_y q^\top y \quad \text{s.t. } Wy \ge h - Tx,\; y \in Y,$$

with dual

$$Q(x) = \max_\pi \pi^\top(h-Tx) \quad \text{s.t. } W^\top \pi \le q,\; \pi \ge 0.$$

If the subproblem is feasible and bounded, an optimal dual solution yields an optimality cut $\theta \ge \pi^\top(h-Tx)$; if it is infeasible, an extreme ray yields the feasibility cut $x_n \in \mathcal{X}_n,\quad n\in N(G),$0 (Brahmbhatt et al., 1 Aug 2025).

The graph-based extension replaces the two-block view by a tree of subgraphs. The hierarchical gBD framework first partitions an OptiGraph into subgraphs that form an acyclic interconnection, selects a root subgraph, assigns stages by breadth-first traversal, and places cost-to-go variables $x_n \in \mathcal{X}_n,\quad n\in N(G),$1 on subgraphs with children. Parent–child links encode the analog of $x_n \in \mathcal{X}_n,\quad n\in N(G),$2, and cuts generated by child subgraphs are added to the parent. In this sense, gBD generalizes nested Benders or dual dynamic programming from linear temporal hierarchies to arbitrary hierarchical graphs, provided that the inter-subgraph topology can be shaped into a tree by partitioning, aggregation, or lifting (Cole et al., 3 Jan 2025).

PlasmoBenders.jl supports both multi-cut and single-cut variants. In the stochastic-programming notation used in the multi-parametric paper, multi-cut introduces one auxiliary variable per subproblem or scenario and adds cuts of the form

$x_n \in \mathcal{X}_n,\quad n\in N(G),$3

whereas the single-cut variant replaces $x_n \in \mathcal{X}_n,\quad n\in N(G),$4 with a single $x_n \in \mathcal{X}_n,\quad n\in N(G),$5 and aggregates the right-hand sides across subproblems. In the computational study built on scenario-time subgraphs, multi-cut converges in fewer iterations, while single-cut exhibits a more oscillatory search because aggregation reduces information per iteration (Brahmbhatt et al., 1 Aug 2025).

3. Multi-parametric programming surrogates and exact cut reproduction

A distinctive extension of PlasmoBenders.jl is the use of multi-parametric programming surrogates for subproblems. The motivation is that classical Benders repeatedly solves many subproblems at successive master iterates, and this repeated solution often dominates runtime. The proposed alternative is to solve each subproblem once, parametrically, with respect to the master variables and uncertain data, then replace repeated online solves by critical-region lookup and affine evaluation (Brahmbhatt et al., 1 Aug 2025).

For an mp-LP parameterized by $x_n \in \mathcal{X}_n,\quad n\in N(G),$6,

$x_n \in \mathcal{X}_n,\quad n\in N(G),$7

subject to

$x_n \in \mathcal{X}_n,\quad n\in N(G),$8

the parameter space $x_n \in \mathcal{X}_n,\quad n\in N(G),$9 is partitioned into critical regions $g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$0. Within each region, the optimal primal and dual solutions are affine,

$g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$1

and the value function is piecewise-affine convex,

$g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$2

hence globally

$g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$3

In PlasmoBenders.jl, these region descriptions and affine maps are stored and then queried during the Benders iterations (Brahmbhatt et al., 1 Aug 2025).

The central theoretical claim is not merely that the surrogate is accurate, but that, under the LP assumptions stated in the paper, the mp-derived cuts are identical to the classical Benders cuts on the visited region. Writing $g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$4, the mp solution on a critical region yields an affine optimal value and an affine map for the dual multipliers with respect to $g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$5. By differentiating the value function within the active region and using the KKT system, the paper shows that the mp-based cut

$g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$6

is identical to the classical cut

$g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$7

with $g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$8 on the current region. Consequently, replacing subproblem solves with mp region lookup preserves the convergence and optimality guarantees of Benders for LP subproblems, because the same cuts are added at the same master iterates (Brahmbhatt et al., 1 Aug 2025).

This result also addresses a common misunderstanding. The mp component is not presented as an approximate oracle in the LP setting studied there; rather, it reproduces the exact recourse value function and the same dual sensitivities on each visited critical region. The paper further emphasizes that critical-region tracking provides a geometric view of how master iterates traverse the parameter space, which supports interpretability and sensitivity analysis in addition to computational acceleration (Brahmbhatt et al., 1 Aug 2025).

4. Software architecture, APIs, and execution modes

In its OptiGraph form, PlasmoBenders.jl takes a graph, lets the user select a master subgraph, and automatically identifies downstream subgraphs as Benders subproblems. The engine fixes complicating variables on the subproblem subgraphs, solves subproblems or queries surrogates, forms Benders cuts, and updates the master. In the hierarchical gBD implementation, this functionality is exposed through a constructor of the form BendersAlgorithm(G, Gr; ...), which inspects the subgraph-edge topology, assigns stages, creates $g_e(\{x_n\}_{n\in N(e)}) \ge 0,\quad e\in E(G).$9 variables and linking fixings, and manages cut pools. The execution entry point is run_algorithm!(alg) (Brahmbhatt et al., 1 Aug 2025, Cole et al., 3 Jan 2025).

The implementation is tightly integrated with JuMP and MathOptInterface. Standard OptiGraph workflows use @optinode, @variable, @constraint, @objective, and @linkconstraint, then aggregate objectives with set_to_node_objectives and assign solvers with set_optimizer. In the mp-surrogate extension, each subproblem node loads a custom Julia “solver” that implements JuMP/MathOptInterface calls by using stored mp data $c^\top x + \theta$0. When optimize! is called on a subproblem node, the custom solver performs point location and affine evaluation to return primal and dual information as if an LP had been solved. mp surrogate generation itself is performed with PPOPT in Python, while the master can use any MILP/LP solver supported by JuMP, such as HiGHS or Gurobi (Brahmbhatt et al., 1 Aug 2025).

The typical mp-enabled workflow is explicitly staged. First, the full model is built as an OptiGraph. Second, for each subproblem class, the corresponding mp-LP is exported and solved once to obtain the critical regions and affine maps. Third, in Julia, a conventional solver is attached to the master and the custom mp “solver” is attached to the subproblems. Fourth, the user chooses single-cut or multi-cut mode, sets tolerances such as $c^\top x + \theta$1, and calls the PlasmoBenders.jl driver. During the Benders iterations, the master is solved, the current parameter vector $c^\top x + \theta$2 is formed for each subproblem, a critical region satisfying $c^\top x + \theta$3 is identified, the affine value and dual sensitivity are evaluated, and the corresponding cut is added to the master (Brahmbhatt et al., 1 Aug 2025).

A later extension introduces distributed-memory execution through RemoteOptiGraph. In that setting, a RemoteOptiGraph holds a lightweight reference to an OptiGraph on a worker, together with remote subgraphs and InterWorkerEdges that manage linking constraints spanning workers. The same Benders meta-algorithm implemented for OptiGraphs can be executed without modification on RemoteOptiGraphs. The execution pattern is synchronous: the planning graph solves the master, current master decisions are broadcast to remote operation graphs, subproblems solve locally on their workers, dual or slack information is returned to the main worker, and cuts are added to the master. The user-facing remote machinery includes RemoteNodeRef, RemoteEdgeRef, and RemoteVariableRef, while proxy references provide minimal serialized identifiers for cross-worker communication. RemoteOptiGraph is available starting in Plasmo.jl 0.7.0, and PlasmoBenders.jl supports RemoteOptiGraph as of version 0.3.0 (Cole et al., 18 Nov 2025).

5. Computational studies and demonstrated use cases

The most detailed performance study for the mp-surrogate extension concerns a two-stage stochastic capacity expansion problem under uncertainty in product demand, prices, and raw-material availability. The first stage chooses capacity expansion decisions $c^\top x + \theta$4, binary expansion indicators $c^\top x + \theta$5, and cumulative capacities $c^\top x + \theta$6; the recourse problem, indexed by scenario $c^\top x + \theta$7 and time $c^\top x + \theta$8, decides production $c^\top x + \theta$9, buy quantities $Ax \le b$0, and sell quantities $Ax \le b$1. The graph representation uses one master subgraph over all times and one subgraph node for each $Ax \le b$2 recourse problem, yielding $Ax \le b$3 independent subproblem nodes that can be processed in parallel. The experimental sizes are $Ax \le b$4 and $Ax \le b$5. In these experiments, mp-Benders consistently reduces subproblem solve time relative to HiGHS-Benders and Gurobi-Benders, with reported speedups for subproblems reaching up to two orders of magnitude in some instances. The bound plots show convergence to the same optimal value for mp-Benders and solver-based Benders, while multi-cut converges in fewer iterations and with shorter runtimes than single-cut for this problem (Brahmbhatt et al., 1 Aug 2025).

The distributed-memory study uses a mixed-integer capacity expansion model for the western United States. The formulation has $Ax \le b$6 subproblems corresponding to week-long operational subperiods, a planning master containing $Ax \le b$7, $Ax \le b$8, and $Ax \le b$9, and operation graphs placed on workers $x \in X$0 through $x \in X$1. The total model size is 12.17 million variables, of which 360 are integer, and 16.03 million constraints. On an AMD EPYC 96-core, 2.4 GHz machine with 53 threads and Gurobi 12.0.0 at a target gap of $x \in X$2, the reported runtimes are 1.67 hours for RemoteOptiGraph + Benders, 1.22 hours for OptiGraph + Benders, and 13.03 hours for the monolithic formulation. These correspond to speedups of $x \in X$3 and $x \in X$4 for the two Benders variants relative to the monolithic solve. The RemoteOptiGraph run is modestly slower than the non-remote OptiGraph run because of communication overhead, but the same code also ran successfully across multiple cluster nodes and cloud VMs with similar times, which demonstrates portability (Cole et al., 18 Nov 2025).

The hierarchical-optimization paper presents additional evidence on how graph structure affects Benders performance. In a capacity expansion model with detailed operations derived from GenX, the graph is transformed into a tree by lifting policy hyperedges into budgeting variables on the planning subgraph, after which a two-stage multi-cut gBD solve with parallel subperiod solves is applied. In that instance, regularization reduces the iteration count from 152 to 63 to reach a tolerance of $x \in X$5. In production cost models with long-duration storage, the choice of partition granularity and root subgraph materially changes convergence: on the 5-bus case, 3 partitions converge in 3 iterations, 5 partitions in 7 iterations, and 10 partitions stall near a $x \in X$6 gap after more than 50 iterations, while choosing the middle block as root improves convergence to 2 iterations when slacks are used for recourse (Cole et al., 3 Jan 2025).

Beyond speedup, the mp study reports an interpretability benefit that is unusual in Benders implementations. Because the mp surrogate partitions parameter space into critical regions, the path of master iterates can be tracked geometrically. In the reported experiments, multi-cut traverses fewer regions and approaches the optimum more directly, while single-cut visits more regions and displays more oscillations. Region active sets identify which constraints dominate the recourse policy in different areas, and the affine matrices $x \in X$7 reveal which recourse variables are sensitive to the master and uncertainty parameters. The paper describes this as supporting scenario clustering and “principal recourse” analysis (Brahmbhatt et al., 1 Aug 2025).

6. Assumptions, limitations, and current research directions

Several limitations follow directly from the mathematical assumptions underlying the current implementations. The mp-surrogate framework focuses on mp-LP subproblems, although the paper states that the machinery extends to convex mp-QP in principle. The identity between mp-derived cuts and classical Benders cuts relies on convexity and on the existence of piecewise-affine solution and value maps. For nonconvex or mixed-integer subproblems, explicit critical-region enumeration becomes harder, and the stated remedies are either to preserve exactness only for convex LP/QP substructures or to fall back to classical subproblem solves on nodes that are not amenable to mp while still using mp where feasible (Brahmbhatt et al., 1 Aug 2025).

A separate limitation concerns scalability of the mp precomputation. The number of critical regions can grow combinatorially with the number of uncertain parameters and the problem dimension, increasing offline computation and memory. The future directions proposed in the paper include modular mp surrogates that partition the parameter set into sub-blocks, adaptive partitioning and on-demand region discovery, and hybrid surrogates such as decision trees or ICNNs inserted in a plug-and-play fashion within the same graph framework. These are presented as possible remedies for the offline burden of explicit region enumeration rather than as implemented features of the current package (Brahmbhatt et al., 1 Aug 2025).

Graph structure itself imposes algorithmic constraints. The hierarchical gBD formulation requires a tree among subgraphs; dense cycles must be removed or reshaped through aggregation, partitioning, or lifting. The distributed-memory paper adds a structural criterion for the current graph-based Benders implementation: linking edges should not connect more than two subgraphs. It also states current modeling restrictions of RemoteOptiGraph, including separability of the overall objective across the primary RemoteOptiGraph and its remote subgraphs, the absence of an InterWorkerObjective, and the representation of linking equalities as two inequalities (Cole et al., 3 Jan 2025, Cole et al., 18 Nov 2025).

Feasibility handling is version- and context-dependent. In the hierarchical gBD implementation, feasibility cuts are not implemented in the software at the time of writing; instead, add_slacks = true can be used to guarantee recourse by penalizing slack variables. In the mp-surrogate study, the case study has relatively complete recourse, so subproblems remain feasible, but the paper notes that if feasibility issues arise, mp analysis can characterize infeasible regions and the usual feasibility cuts can be generated from extreme rays of the dual analogously to classical Benders. This indicates that the underlying mathematical framework includes classical feasibility logic even where the software path emphasized in a given paper is slack-based recourse (Cole et al., 3 Jan 2025, Brahmbhatt et al., 1 Aug 2025).

Two misconceptions are explicitly countered by the literature. First, graph-based Benders is not restricted to time-ordered nested decompositions; the gBD framework is intended for arbitrary hierarchical graphs after suitable graph transformations. Second, the mp component is not merely a heuristic acceleration layer in the LP setting analyzed there; on a computed critical region, the package generates the same cuts as classical Benders and therefore preserves the same convergence and optimality guarantees. A plausible implication is that PlasmoBenders.jl’s principal research significance lies in separating modeling, decomposition logic, and surrogate choice while retaining a single graph-centered representation across serial, parallel, and distributed-memory regimes (Cole et al., 3 Jan 2025, Brahmbhatt et al., 1 Aug 2025, Cole et al., 18 Nov 2025).