How Benders Decomposition Works Part I: Optimality Cuts

military planners with coordination problems on a scale that had rarely been imagined. Aircraft, personnel, fuel, equipment, training, and supplies had to be allocated across an immense system, often under severe time pressure and uncertainty. This environment accelerated the development of operations research, a discipline founded on the idea that complex decisions could be studied systematically through data, scientific reasoning, and mathematics rather than addressed solely through experience and intuition [1]. George Dantzig was among the researchers immersed in this new world of large-scale planning. During the war, he worked for the United States Army Air Forces, where he contributed to the analysis and planning of military operations. After the war, the Air Force challenged him to help mechanize planning processes that had previously been performed largely by hand. In 1947, Dantzig formulated the general linear programming problem and developed the simplex method, providing a systematic procedure for finding the best allocation of limited resources among competing activities [2, 3]. This research continued through Project SCOOP, the Scientific Computation of Optimal Programs, and later at the RAND Corporation, which became one of the principal centers for the early development of linear programming and mathematical optimization [2]. This mathematical revolution eventually reached the Shell Laboratories in Amsterdam, where Jacques Benders began working in 1955. Together with his office mate, Guus Zoutendijk, he encountered difficult refinery-planning problems involving production decisions, material flows, processing resources, and costs. After visiting Shell’s Californian branch and learning about the potential of linear programming, they returned to the Netherlands and began implementing the simplex method for production planning in the refinery and petrochemical industries. Mathematical optimization promised to help coordinate industrial systems whose many interconnected decisions could easily overwhelm unaided human reasoning [4]. However, the promise of linear programming came with a frustrating limitation. It worked beautifully when decisions could vary continuously, such as how much material to process, produce, or transport. Many important industrial decisions were not continuous. They required choosing an option or rejecting it, activating an alternative or leaving it unused, and answering yes or no. Introducing these discrete variables transformed an otherwise manageable linear program into a much harder mixed-variable optimization problem. Benders focused precisely on this complication. In 1959, he and two Shell colleagues presented an algorithm for linear programs containing binary variables at a RAND symposium, reporting computational experiments performed on a Ferranti Mark I computer whose addition time was approximately 1.2 milliseconds [4]. Benders could have searched only for a faster way to solve the complete model, but instead he questioned whether the complete model needed to be solved all at once. The discrete variables created much of the difficulty, yet once their values were fixed, the remaining decisions often formed an ordinary linear program. This suggested a different strategy: isolate the complicating decisions, allow one problem to choose them, and let a second problem determine their operational consequences. Instead of treating the optimization model as an indivisible object, Benders divided it into parts that could communicate with and progressively improve one another [5, 6]. This insight produced something more subtle than simply solving two smaller problems. The first problem, which we now call the master problem, proposes a high-level decision while initially possessing only an incomplete understanding of its true consequences. The subproblem evaluates that proposal and returns mathematical information showing where the master has been too optimistic. This information is converted into a new constraint, known as a Benders cut, and added to the master before it tries again. With each iteration, the master develops a more accurate representation of the operational reality hidden behind its decisions. The idea became the foundation of Benders’s doctoral thesis, Partitioning in Mathematical Programming, completed at Utrecht University in 1960 under the supervision of Hans Freudenthal [4, 5]. In 1962, he published the method in the now-classic paper “Partitioning Procedures for Solving Mixed-Variables Programming Problems” [6]. More than six decades later, Benders decomposition remains one of the most influential decomposition techniques in mathematical optimization. It has been extended, strengthened, and successfully applied to facility location, transportation, supply-chain design, energy systems, stochastic programming, scheduling, and many other difficult optimization problems [7]. Benders decomposition has also become a staple in my own operations research toolkit, but I must admit that it took me considerably longer than it should have to understand it properly. Its central insight is elegant and intuitive, yet many learning materials introduce it through dense notation, extreme points, dual rays, and lengthy derivations before clearly explaining what the algorithm is trying to accomplish. As a result, a technique built around a remarkably simple conversation between two optimization problems can initially appear far more mysterious than it really is. This series of three articles is my attempt to provide the Benders decomposition tutorial that I would have loved to have during the first years of my PhD. My objective is to make the method accessible to anyone fascinated by operations research and mathematical optimization, without sacrificing the mathematical reasoning required to understand why it works. We will progress from classical optimality cuts, to infeasible subproblems and feasibility cuts, and finally to logic-based Benders decomposition for combinatorial problems. In this first article, we begin with the most accessible possible setting: a problem in which obtaining a feasible solution is straightforward, but identifying the optimal one remains challenging. We will first explain the general mechanics of Benders decomposition and the roles played by the master problem, the subproblem, and the optimality cuts. We will then work through a small toy problem by hand, deriving each cut step by step until the lower and upper bounds converge to the optimal solution. Finally, we will apply the same machinery to the uncapacitated facility location problem and implement the complete algorithm in Python using Pyomo and the open-source HiGHS solver. 1. The basic idea behind Benders decomposition Many optimization problems combine two fundamentally different types of decisions. The first are strategic decisions that determine the structure of the solution and are often binary or integer. The second are operational decisions that become much easier to optimize once the strategic choices are known. Facility location provides a natural example. Deciding which facilities to open is a strategic decision, while assigning customers to the selected facilities is an operational decision. Benders decomposition separates these two layers instead of solving them simultaneously. Consider an optimization problem of the following form: minx,y⁡c⊤x+d⊤y\min_{x,y} \quad c^\top x + d^\top y subject to Ax+By≥b,x∈X,y≥0.Ax + By \geq b, \qquad x \in X, \qquad y \geq 0. The variables (x)(x) represent the difficult strategic decisions, while (y)(y)contains the operational decisions. If we temporarily fix the strategic variables at some value (x‾\bar{x}), the remaining optimization problem becomes: Q(x‾)=miny⁡d⊤yQ(\bar{x}) = \min_y \quad d^\top y Subject to By≥b−Ax‾,y≥0By \geq b – A\bar{x}, \qquad y \geq 0 The function Q(x)Q(x) represents the best operational cost associated with a strategic decision (xx). Therefore, the original problem can be viewed as: minx∈X⁡c⊤x+Q(x).\min_{x \in X} \quad c^\top x + Q(x). The difficulty is that the master problem does not initially know the function Q(x)Q(x). Benders decomposition introduces a variable (θ\theta) to approximate it: minx,θ⁡c⊤x+θ.\min_{x,\theta} \quad c^\top x + \theta. The master and subproblem are then solved iteratively, with each subproblem solution providing new information about Q(x)Q(x). The next question is where this information comes from and how it can be expressed as a valid constraint. 2. Where do Benders cuts come from? For a fixed strategic decision (xx), the operational subproblem is: Q(x)=miny⁡d⊤yQ(x)=\min_y \quad d^\top y subject to By≥b−Ax,y≥0.By\geq b-Ax, \qquad y\geq 0. Because this is a linear program, it has the following dual: maxu⁡(b−Ax)⊤u\max_u \quad (b-Ax)^\top u subject to B⊤u≤d,u≥0.B^\top u\leq d, \qquad u\geq 0. Every feasible dual solution (uku^k) provides a lower bound on the operational cost: Q(x)≥(b−Ax)⊤uk.Q(x)\geq (b-Ax)^\top u^k. Since (θ\theta) represents the master problem’s approximation of Q(x)Q(x), we can impose: θ≥(b−Ax)⊤uk.\theta\geq (b-Ax)^\top u^k. This inequality is a Benders optimality cut. Its key advantage is that it does not only describe the operational cost of the current master solution. It remains valid for other values of (xx), allowing one subproblem solve to improve the master’s understanding of several possible strategic decisions. 3. A small example, solved step by step The previous sections introduced the formal structure of Benders decomposition and explained where optimality cuts come from. However, these ideas can still feel abstract when they are presented only through general notation. To make the mechanism more concrete, we will now work through a very small example by hand and follow the algorithm one iteration at a time. The purpose of this example is not to solve a realistic optimization problem. It is to make the interaction between the master problem and the subproblem completely visible. We will see how the master begins with an overly optimistic estimate, how the subproblem reveals the true operational cost of that decision, and how the dual solution generates a cut that makes the master slightly wiser in the next iteration. 3.1 The toy optimization problem So, consider the following small optimization problem: min⁡6x1+4x2+y1+y2\min \quad 6x_1 + 4x_2 + y_1 + y_2 subject to y1≥10−8x1−4x2,y_1 \geq 10 – 8x_1 – 4x_2, y2≥6−2x1−6x2,y_2 \geq 6 – 2x_1 – 6x_2, x1,x2∈0,1,y1,y2≥0.x_1,x_2 \in {0,1}, \qquad y_1,y_2 \geq 0. We can interpret (x1x_1) and (x2x_2) as two strategic investment decisions. Activating the first option costs 6 units, while activating the second costs 4 units. The variables (y1y_1) and (y2y_2) represent the remaining operational requirements after those strategic decisions have been made, with each unit carrying a cost of 1. 3.2 The initial master problem Benders decomposition begins by placing the strategic variables (x1x_1) and (x2x_2) in the master problem. The operational variables (y1y_1) and (y2y_2) are temporarily removed and replaced by the variable (θ\theta), which represents an estimate of their combined cost. The initial master problem is therefore: min⁡6x1+4x2+θ\min \quad 6x_1 + 4x_2 + \theta subject to θ≥0,\theta \geq 0, x1,x2∈0,1.x_1,x_2 \in {0,1}. At this stage, the master has not yet received any Benders cuts. It only knows that the operational cost cannot be negative. Consequently, its optimal solution is: x‾1=0,x‾2=0,θ=0.\bar{x}_1=0, \qquad \bar{x}_2=0, \qquad \theta=0. The master therefore proposes making no strategic investment and assumes that the operational cost will also be zero. Its objective value is LB=0.LB=0. This value provides a lower bound on the optimal objective because the master is currently underestimating the operational cost. The solution is not yet valid for the original problem, since the operational requirements represented by (y1y_1) and (y2y_2) have not been evaluated. The next step is therefore to send the decision (x‾=(0,0)\bar{x}=(0,0)) to the subproblem and determine its true operational consequences. 3.3 Evaluating the master decision The master problem has proposed the strategic decision: x‾1=0,x‾2=0.\bar{x}_1=0, \qquad \bar{x}_2=0. The subproblem is obtained by fixing the strategic variables (x1x_1) and (x2x_2) at these values and optimizing only the operational variables (y1y_1) and (y2y_2). Recall the objective function of the original problem: min⁡6x1+4x2+y1+y2.\min \quad 6x_1+4x_2+y_1+y_2. The terms (6x1+4x26x_1+4x_2) represent the strategic cost and are already accounted for in the master problem. Therefore, the subproblem only needs to minimize the operational cost: Q(x‾)=min⁡y1+y2.Q(\bar{x})=\min \quad y_1+y_2. We now substitute the master solution (x‾1=0\bar{x}_1=0) and (x‾2=0\bar{x}_2=0) into each operational constraint. The first constraint becomes: y1≥10−8(0)−4(0)=10.y_1\geq 10-8(0)-4(0)=10. The second constraint becomes: y2≥6−2(0)−6(0)=6.y_2\geq 6-2(0)-6(0)=6. The complete subproblem is therefore: Q(0,0)=min⁡y1+y2Q(0,0)=\min \quad y_1+y_2 subject to y1≥10,y_1\geq 10, y2≥6,y_2\geq 6, y1,y2≥0.y_1,y_2\geq 0. Because both variables have positive coefficients in the objective function, the subproblem chooses the smallest feasible values: y1=10,y2=6.y_1=10, \qquad y_2=6. The operational cost associated with the master decision is therefore: Q(0,0)=10+6=16.Q(0,0)=10+6=16. The strategic cost of choosing (x1=x2=0(x_1=x_2=0) is zero, so the total cost of this feasible solution to the original problem is: 6(0)+4(0)+16=16.6(0)+4(0)+16=16. This gives us an upper bound of (16), while the master objective provides a lower bound of (0): LB=0,UB=16.LB=0, \qquad UB=16. 3.4 Constructing the dual subproblem The subproblem tells us the operational cost associated with the current master decision. However, to generate a Benders cut, we need more than the value (Q(0,0)=16Q(0,0)=16). We need an expression that remains valid when the master later considers different values of (x1x_1) and (x2x_2). This expression comes from the dual of the subproblem. Before substituting the current master solution, let us write the subproblem using the generic strategic variables (x1x_1) and (x2x_2): Q(x)=min⁡y1+y2Q(x)=\min \quad y_1+y_2 subject to y1≥10−8x1−4x2,y_1\geq 10-8x_1-4x_2, y2≥6−2x1−6x2,y_2\geq 6-2x_1-6x_2, y1,y2≥0.y_1,y_2\geq 0. To construct the dual, we use two basic rules: Each constraint in the primal subproblem creates one dual variable. Each variable in the primal subproblem creates one constraint in the dual. The first primal constraint receives the dual variable (u1u_1), while the second receives (u2u_2): u1⟷y1≥10−8x1−4x2,u_1 \longleftrightarrow y_1\geq 10-8x_1-4x_2, u2⟷y2≥6−2x1−6x2.u_2 \longleftrightarrow y_2\geq 6-2x_1-6x_2. Because the primal is a minimization problem and its constraints use the greater-than-or-equal-to direction, both dual variables must be nonnegative: u1,u2≥0.u_1,u_2\geq 0. The dual objective is constructed by multiplying the right-hand side of each primal constraint by its corresponding dual variable. Therefore, the dual objective is: max⁡(10−8x1−4x2)u1+(6−2x1−6x2)u2.\max \quad (10-8x_1-4x_2)u_1 + (6-2x_1-6x_2)u_2. We now construct the dual constraints. There is one dual constraint for each primal variable. The variable (y1y_1) appears only in the first primal constraint, with a coefficient of 1. Its coefficient in the primal objective is also 1. Therefore, the corresponding dual constraint is: u1≤1.u_1\leq 1. Similarly, (y2y_2) appears only in the second primal constraint, also with a coefficient of 1, and its objective coefficient is 1. Its dual constraint is therefore: u2≤1.u_2\leq 1. Combining these elements, the dual of the general subproblem is: max⁡(10−8x1−4x2)u1+(6−2x1−6x2)u2\max \quad (10-8x_1-4x_2)u_1 + (6-2x_1-6x_2)u_2 subject to 0≤u1≤1,0\leq u_1\leq 1, 0≤u2≤1.0\leq u_2\leq 1. We can now evaluate this dual at the current master solution (x‾1=0\bar{x}_1=0) and (x‾2=0\bar{x}_2=0). Substituting these values into the objective gives: max⁡10u1+6u2\max \quad 10u_1+6u_2 subject to 0≤u1≤1,0\leq u_1\leq 1, 0≤u2≤1.0\leq u_2\leq 1. Both coefficients in the objective are positive, so the objective is maximized by setting both dual variables to their largest feasible values: u1=1,u2=1.u_1^=1, \qquad u_2^=1. The optimal dual objective value is therefore: 10(1)+6(1)=16.10(1)+6(1)=16. This matches the optimal value of the primal subproblem, Q(0,0)=16. This equality is a consequence of strong duality. Because the primal and dual subproblems are both feasible linear programs, their optimal objective values are equal. The dual solution (u1=u2=1u_1 =u_2 =1) will now allow us to transform the operational information obtained at (x=(0,0)x=(0,0)) into a Benders cut that remains valid for other strategic decisions. 3.5 Generating the first Benders cut The optimal dual solution obtained for the current master decision is, u1=1,u2=1.u_1=1, u_2=1.Recall that the objective of the general dual subproblem is: (10−8x1−4x2)u1+(6−2x1−6x2)u2.(10-8x_1-4x_2)u_1+(6-2x_1-6x_2)u_2. To generate the Benders cut, we substitute the optimal dual values into this expression: θ≥(10−8x1−4x2)(1)+(6−2x1−6x2)(1).\theta \geq (10-8x_1-4x_2)(1)+(6-2x_1-6x_2)(1). After collecting terms, we obtain: θ≥16−10x1−10x2.\theta \geq 16-10x_1-10x_2. This inequality is the first Benders optimality cut. At the current solution (x1=x2=0x_1=x_2=0), it forces (θ\theta) to be at least 16, which is exactly the operational cost found by the subproblem. However, the cut also describes how this lower bound changes when either strategic option is selected. We now add the cut to the master problem: min⁡6x1+4x2+θ\min \quad 6x_1+4x_2+\theta subject to θ≥0,\theta\geq 0, θ≥16−10x1−10x2,\theta\geq 16-10x_1-10x_2, x1,x2∈0,1.x_1,x_2\in{0,1}. The master can no longer choose (x1=x2=0x_1=x_2=0) while pretending that the operational cost is zero. If it makes no investment, the new cut requires (θ≥16\theta\geq16), giving a total estimated cost of 16. If it selects only the second option, the cut gives: θ≥16−10(0)−10(1)=6.\theta\geq 16-10(0)-10(1)=6. The corresponding master objective value is: 6(0)+4(1)+6=10.6(0)+4(1)+6=10. The updated master therefore moves away from the original decision and proposes: x‾1=0,x‾2=1,θ=6.\bar{x}_1=0, \qquad \bar{x}_2=1, \qquad \theta=6. Its new objective value is: LB=10.LB=10. The master has now learned that avoiding strategic investment creates a substantial operational burden. The next step is to send the new decision (x‾=(0,1)\bar{x}=(0,1)) to the subproblem and verify whether the estimated operational cost (θ=6\theta=6) is accurate. 3.6 The second iteration and convergence After adding the first optimality cut, the master proposes the solution, with an objective value of (10), which becomes the new lower bound LB=10.LB=10. It is important to highlight, that the updated master is actually indifferent between (x=(0,1)x=(0,1) with (θ=6\theta = 6)) and (x=(1,1)x=(1,1) with (θ=0\theta = 0)), since both produce an objective value of (10). The solution selected may therefore depend on the solver’s tie-breaking rules. For this explanation, we will continue with (x=(0,1)x=(0,1)). To evaluate this decision, we once again fix the strategic variables in the original operational constraints of the subproblem. The first constraint becomes: y1≥10−8(0)−4(1)=6,y_1\geq 10-8(0)-4(1)=6, while the second becomes: y2≥6−2(0)−6(1)=0.y_2\geq 6-2(0)-6(1)=0. The resulting subproblem is then: Q(0,1)=min⁡y1+y2Q(0,1)=\min \quad y_1+y_2 subject to y1≥6,y_1\geq 6, y2≥0,y_2\geq 0, y1,y2≥0.y_1,y_2\geq 0. The optimal operational solution is therefore: y1=6,y2=0,y_1=6, \qquad y_2=0, with an operational cost of, Q(0,1)=6.Q(0,1)=6. This is exactly the value predicted by the master through (θ=6\theta=6). Adding the strategic cost gives a total feasible cost of: 6(0)+4(1)+6=10.6(0)+4(1)+6=10. This solution improves the previous upper bound of (16), so, UB=10.UB=10.We now have: LB=10,UB=10.LB=10, \qquad UB=10. The lower and upper bounds coincide, which proves that the current solution (x1=1,x2=0,y1=6,y2=0x_1 = 1, x_2=0,y_1=6,y_2=0) is optimal. The master cannot find a solution costing less than (10), while the subproblem has confirmed that a feasible solution with cost (10) exists. Benders decomposition has therefore converged. For the sake of completeness, the dual subproblem at (x=(0,1)x=(0,1)) is: max⁡6u1\max \quad 6u_1 subject to 0≤u1≤1,0≤u2≤1.0\leq u_1\leq 1, \qquad 0\leq u_2\leq 1. One optimal dual solution is (u1=1u_1=1) and (u2=0u_2=0), which would generate the additional cut: θ≥10−8x1−4x2.\theta\geq 10-8x_1-4x_2. However, this cut is not required to prove optimality because the lower and upper bounds have already converged. 4. Benders decomposition for the uncapacitated facility location problem Now that we understand how Benders optimality cuts are generated, we can apply the machinery to one of the quintessential problems in operations research: the uncapacitated facility location problem, or UFL. Facility location problems appear whenever an organization must decide where to establish resources that will serve geographically distributed demand. Typical applications include locating warehouses, distribution centers, hospitals, emergency services, data centers, and charging stations. The strategic decision is where to open the facilities, while the operational decision is how customers should be assigned to the selected locations. The uncapacitated version assumes that any open facility can serve any number of customers. This makes it an especially convenient introduction to Benders decomposition because, provided that at least one facility is opened and every customer can be served by every facility, the assignment subproblem is always feasible. 4.1 The uncapacitated facility location model Let (II) denote the set of candidate facilities and (JJ) the set of customers. Opening facility (ii) incurs a fixed cost (fif_i), while assigning customer (jj) to facility (ii) incurs a service cost (qijq_{ij}). We define the strategic variables: xi=1 if facility i is opened, and 0 otherwise.x_i = 1 \text{ if facility } i \text{ is opened, and } 0 \text{ otherwise}. The operational variables are: yij=the fraction of customer j served by facility i.y_{ij} = \text{the fraction of customer } j \text{ served by facility } i. The complete, or monolithic, formulation is then: min⁡∑ifixi+∑i∑jqijyij\min \quad \sum_i f_i x_i + \sum_i \sum_j q_{ij}y_{ij} subject to ∑iyij=1∀j,\sum_i y_{ij}=1 \qquad \forall j, yij≤xi∀i,j,y_{ij}\leq x_i \qquad \forall i,j, ∑ixi≥1,\sum_i x_i\geq 1, xi∈0,1,yij≥0.x_i\in{0,1}, \qquad y_{ij}\geq 0. The first constraint requires every customer to be fully assigned. The second prevents customers from being assigned to closed facilities, the third constraints, forces to open at least one facility. Although the assignment variables are continuous, the structure of the problem ensures that an optimal solution assigns each customer to one of the cheapest open facilities. 4.2 The master problem In the Benders decomposition, the facility-opening variables remain in the master problem. The assignment variables are removed and their cost is represented by (θ\theta): min⁡∑ifixi+θ\min \quad \sum_i f_i x_i+\theta subject to ∑ixi≥1,\sum_i x_i\geq 1, xi∈0,1,θ≥0,x_i\in{0,1}, \qquad \theta\geq 0, Please notice that the constraints above, would be joined together with the Benders cuts generated during the different iterations of the Benders algorithm. The master therefore decides which facilities to open while gradually learning the resulting assignment cost. 4.3 The assignment subproblem Once the master produces a facility-opening decision (x‾\bar{x}), we fix those values and solve the operational assignment problem: min⁡∑i∑jqijyij\min \quad \sum_i \sum_j q_{ij}y_{ij} subject to ∑iyij=1∀j,\sum_i y_{ij}=1 \qquad \forall j, yij≤x‾i∀i,j,y_{ij}\leq \bar{x}_i \qquad \forall i,j, yij≥0.y_{ij}\geq 0. 4.4 The dual subproblem and optimality cut To generate a Benders cut, we construct the dual of the assignment subproblem. Let (αj\alpha_j) be the dual variable associated with the assignment equality for customer (jj). Since the corresponding primal constraint is an equality, (αj\alpha_j) is unrestricted in sign. We rewrite the linking constraint as: −yij≥−x‾i-y_{ij}\geq-\bar{x}_i and associate it with a nonnegative dual variable (βij\beta_{ij}). The dual subproblem is then: max⁡∑jαj−∑i∑jx‾iβij\max \quad \sum_j \alpha_j – \sum_i \sum_j \bar{x}_i\beta_{ij} subject to αj−βij≤qij∀i,j,\alpha_j-\beta_{ij}\leq q_{ij} \qquad \forall i,j, βij≥0,αj unrestricted.\beta_{ij}\geq 0, \qquad \alpha_j \text{ unrestricted}. Suppose that solving the dual produces the optimal values (αjk\alpha_j^k) and (βijk\beta_{ij}^k). The corresponding Benders optimality cut is then: θ≥∑jαjk−∑i∑jβijkxi\theta \geq \sum_j \alpha_j^k – \sum_i \sum_j \beta_{ij}^k x_i This cut gives the master a new lower bound on the assignment cost. The master is then solved again, a new facility configuration is evaluated, and additional cuts are generated until the lower and upper bounds converge. 5. Implementing Benders decomposition with Pyomo and HiGHS We now have all the mathematical machinery required to implement the Benders decomposition algorithm discussed in the previous sections. The master problem will select the facilities to open, the dual subproblem will evaluate the resulting customer-assignment cost, and the generated optimality cuts will progressively improve the master problem. To illustrate the complete procedure, we will solve a small uncapacitated facility location instance containing five candidate facilities and twenty customers. The instance is stored in JSON format so that the optimization model remains separate from the data and can easily be reused or modified. Both the instance and the complete notebook can be downloaded from the accompanying GitHub repository. The implementation uses Pyomo to formulate the optimization models and the open-source HiGHS solver to solve the master and subproblems. We will first load and prepare the data, then construct each component of the decomposition algorithm before combining them into the complete iterative procedure. 5.1 Loading and preparing the instance We begin by installing the required packages and importing the libraries used throughout the implementation. The JSON instance is then downloaded directly from the GitHub repository and converted into the sets and dictionaries required by Pyomo. !pip -q install pyomo highspy import requests import pandas as pd import matplotlib.pyplot as plt import pyomo.environ as pyo from pyomo.opt import SolverFactory, TerminationCondition # Load the instance from GitHub JSON_URL = ( "https://raw.githubusercontent.com/ceche1212/" "Benders_Tutorials_TDS/refs/heads/main/Data/" "ufl_benders_5x20_instance.json" ) response = requests.get(JSON_URL, timeout=30) response.raise_for_status() instance = response.json() # Sets facilities = [facility["id"] for facility in instance["facilities"]] customers = [customer["id"] for customer in instance["customers"]] # Facility parameters fixed_cost = { facility["id"]: facility["fixed_cost"] for facility in instance["facilities"] } facility_coordinates = { facility["id"]: (facility["x"], facility["y"]) for facility in instance["facilities"] } # Customer parameters demand = { customer["id"]: customer["demand"] for customer in instance["customers"] } customer_coordinates = { customer["id"]: (customer["x"], customer["y"]) for customer in instance["customers"] } # Assignment costs assignment_cost = { (record["facility"], record["customer"]): record["cost"] for record in instance["assignment_costs"] } Before building the optimization models, we perform a small validation to confirm that all identifiers are unique and that every facility-customer combination has an assignment cost. if len(facilities) != len(set(facilities)): raise ValueError("Duplicate facility IDs were found.") if len(customers) != len(set(customers)): raise ValueError("Duplicate customer IDs were found.") expected_pairs = { (facility, customer) for facility in facilities for customer in customers } available_pairs = set(assignment_cost) if expected_pairs != available_pairs: missing_pairs = expected_pairs - available_pairs extra_pairs = available_pairs - expected_pairs raise ValueError( f"Missing assignment costs: {sorted(missing_pairs)}\n" f"Unexpected assignment costs: {sorted(extra_pairs)}" ) print("Instance loaded successfully") print(f"Name: {instance['instance_name']}") print(f"Facilities: {len(facilities)}") print(f"Customers: {len(customers)}") print(f"Assignment costs: {len(assignment_cost)}") The JSON data is transformed into ordinary Python lists and dictionaries. The facility and customer coordinates will later be used to visualize the final solution, while fixed_cost and assignment_cost contain the parameters required by the master and subproblems. 5.2 Building the initial master problem We first construct the master problem, which contains only the strategic facility-opening decisions and the variable (θ\theta), representing the current approximation of the customer-assignment cost. The ConstraintList is initially empty because the Benders optimality cuts will be generated and added during the iterative algorithm. # ============================================================ # Build the initial master problem # ============================================================ master = pyo.ConcreteModel() master.FACILITIES = pyo.Set( initialize=facilities, ordered=True ) master.fixed_cost = pyo.Param( master.FACILITIES, initialize=fixed_cost ) # Facility-opening decisions and assignment-cost approximation master.x = pyo.Var( master.FACILITIES, domain=pyo.Binary ) master.theta = pyo.Var( domain=pyo.NonNegativeReals ) # At least one facility must be opened master.OpenAtLeastOne = pyo.Constraint( expr=sum(master.x[i] for i in master.FACILITIES) >= 1 ) # Optimality cuts will be added during the algorithm master.BendersCuts = pyo.ConstraintList() master.TotalCost = pyo.Objective( expr=( sum( master.fixed_cost[i] * master.x[i] for i in master.FACILITIES ) + master.theta ), sense=pyo.minimize ) # ============================================================ # Configure and solve with HiGHS # ============================================================ solver = SolverFactory("appsi_highs") solver.highs_options.update({ "output_flag": True, "mip_rel_gap": 0.0, "time_limit": 360 }) results = solver.solve(master, tee=True) if results.solver.termination_condition != TerminationCondition.optimal: raise RuntimeError( "The initial master problem was not solved to optimality." ) # ============================================================ # Extract the initial solution # ============================================================ x_solution = { i: int(round(pyo.value(master.x[i]))) for i in master.FACILITIES } theta_value = pyo.value(master.theta) fixed_cost_value = sum( fixed_cost[i] * x_solution[i] for i in master.FACILITIES ) master_objective_value = pyo.value(master.TotalCost) print("\nInitial master solution") for i, value in x_solution.items(): status = "Open" if value else "Closed" print(f"Facility {i}: x = {value} ({status})") print(f"\nFixed opening cost: {fixed_cost_value:,.2f}") print(f"Estimated assignment cost: {theta_value:,.2f}") print(f"Master objective value: {master_objective_value:,.2f}") Initial master solution Facility F1: x = 0 (Closed) Facility F2: x = 0 (Closed) Facility F3: x = 0 (Closed) Facility F4: x = 1 (Open) Facility F5: x = 0 (Closed) Fixed opening cost: 2,100.00 Estimated assignment cost: 0.00 Master objective value: 2,100.00 The initial master has no information about the actual assignment cost beyond its nonnegativity, so it sets (θ=0\theta=0). It therefore opens only the facility with the lowest fixed cost, producing an optimistic lower bound for the original problem. The dual subproblem will now evaluate this facility configuration and generate the first optimality cut. 5.3 Evaluating the initial decision and generating the first cut The initial master solution opens the facility with the lowest fixed cost while assuming that customer assignments cost nothing. We now evaluate this decision using the dual assignment subproblem derived earlier. The function below receives the current facility-opening vector (x‾\bar{x}), solves the dual problem, and returns the information needed to construct an optimality cut. # ============================================================ # Build the dual assignment subproblem # ============================================================ def build_dual_subproblem(x_bar): dual = pyo.ConcreteModel() dual.FACILITIES = pyo.Set( initialize=facilities, ordered=True ) dual.CUSTOMERS = pyo.Set( initialize=customers, ordered=True ) dual.assignment_cost = pyo.Param( dual.FACILITIES, dual.CUSTOMERS, initialize=assignment_cost ) dual.x_bar = pyo.Param( dual.FACILITIES, initialize=x_bar ) # alpha_j is unrestricted; beta_ij is nonnegative dual.alpha = pyo.Var( dual.CUSTOMERS, domain=pyo.Reals ) dual.beta = pyo.Var( dual.FACILITIES, dual.CUSTOMERS, domain=pyo.NonNegativeReals ) dual.DualFeasibility = pyo.Constraint( dual.FACILITIES, dual.CUSTOMERS, rule=lambda dual, i, j: ( dual.alpha[j] - dual.beta[i, j] = cut_expression ) print("First Benders optimality cut added successfully.") The parameters x_bar fix the facility decisions proposed by the master. The variables (αj\alpha_j) correspond to the customer-assignment equalities, while (βij\beta_{ij}) correspond to the constraints preventing assignments to closed facilities. Solving the dual provides both the true assignment cost of the current configuration and the coefficients required for the cut. The cut is then added to master.BendersCuts. When the master is solved again, it can no longer assign an unrealistically low value to (θ\theta) for the same facility configuration. 5.4 Automating the Benders iterations The previous cells performed the first Benders iteration explicitly. We can now automate the same sequence: solve the master problem, evaluate its facility configuration with the dual subproblem, update the lower and upper bounds, and add a new optimality cut. The procedure stops when both bounds coincide, proving that the current solution is optimal. # ============================================================ # Benders decomposition algorithm # ============================================================ maximum_iterations = 100 tolerance = 1e-6 lower_bound = -float("inf") upper_bound = float("inf") best_x_solution = None best_fixed_cost = None best_assignment_cost = None best_total_cost = None iteration_results = [] # Suppress the HiGHS log during the iterative procedure solver.highs_options["output_flag"] = False solver.highs_options["log_to_console"] = False for iteration in range(1, maximum_iterations + 1): # -------------------------------------------------------- # 1. Solve the master problem # -------------------------------------------------------- master_results = solver.solve(master, tee=False) if master_results.solver.termination_condition != TerminationCondition.optimal: raise RuntimeError( f"Master problem failed at iteration {iteration}." ) x_bar = { i: int(round(pyo.value(master.x[i]))) for i in master.FACILITIES } theta_value = pyo.value(master.theta) fixed_cost_value = sum( fixed_cost[i] * x_bar[i] for i in facilities ) lower_bound = pyo.value(master.TotalCost) # -------------------------------------------------------- # 2. Solve the dual subproblem # -------------------------------------------------------- dual = build_dual_subproblem(x_bar) dual_results = solver.solve(dual, tee=False) if dual_results.solver.termination_condition != TerminationCondition.optimal: raise RuntimeError( f"Dual subproblem failed at iteration {iteration}." ) assignment_cost_value = pyo.value(dual.AssignmentCost) total_cost = fixed_cost_value + assignment_cost_value # -------------------------------------------------------- # 3. Update the best feasible solution # -------------------------------------------------------- if total_cost = cut_expression ) else: raise RuntimeError( "Maximum number of iterations reached without convergence." ) After nine iterations, the Benders algorithm converges to the optimal solution. At that point, the lower bound provided by the master problem and the upper bound obtained from the evaluated subproblem reach the same value, which proves optimality. ---------------------------------------------------------------------- Iteration: 1 Open facilities: ['F2', 'F3', 'F5'] Fixed cost: 7,000.00 Theta: 56.00 Actual assignment cost: 8,955.00 Lower bound: 7,056.00 Upper bound: 15,955.00 Relative gap: 55.775619% ---------------------------------------------------------------------- Iteration: 2 Open facilities: ['F1', 'F2'] Fixed cost: 4,600.00 Theta: 7,306.00 Actual assignment cost: 13,781.00 Lower bound: 11,906.00 Upper bound: 15,955.00 Relative gap: 25.377625% ---------------------------------------------------------------------- Iteration: 3 Open facilities: ['F2', 'F4', 'F5'] Fixed cost: 6,600.00 Theta: 6,403.00 Actual assignment cost: 8,143.00 Lower bound: 13,003.00 Upper bound: 14,743.00 Relative gap: 11.802211% ---------------------------------------------------------------------- Iteration: 4 Open facilities: ['F1', 'F3', 'F4', 'F5'] Fixed cost: 9,300.00 Theta: 4,206.00 Actual assignment cost: 5,455.00 Lower bound: 13,506.00 Upper bound: 14,743.00 Relative gap: 8.390423% ---------------------------------------------------------------------- Iteration: 5 Open facilities: ['F3', 'F5'] Fixed cost: 4,800.00 Theta: 8,955.00 Actual assignment cost: 13,792.00 Lower bound: 13,755.00 Upper bound: 14,743.00 Relative gap: 6.701485% ---------------------------------------------------------------------- Iteration: 6 Open facilities: ['F2', 'F3', 'F4'] Fixed cost: 6,800.00 Theta: 7,057.00 Actual assignment cost: 8,750.00 Lower bound: 13,857.00 Upper bound: 14,743.00 Relative gap: 6.009632% ---------------------------------------------------------------------- Iteration: 7 Open facilities: ['F2', 'F5'] Fixed cost: 4,500.00 Theta: 9,594.00 Actual assignment cost: 11,243.00 Lower bound: 14,094.00 Upper bound: 14,743.00 Relative gap: 4.402089% ---------------------------------------------------------------------- Iteration: 8 Open facilities: ['F1', 'F5'] Fixed cost: 4,700.00 Theta: 9,594.00 Actual assignment cost: 12,241.00 Lower bound: 14,294.00 Upper bound: 14,743.00 Relative gap: 3.045513% ---------------------------------------------------------------------- Iteration: 9 Open facilities: ['F2', 'F4', 'F5'] Fixed cost: 6,600.00 Theta: 8,143.00 Actual assignment cost: 8,143.00 Lower bound: 14,743.00 Upper bound: 14,743.00 Relative gap: 0.000000% Benders decomposition converged. 5.5 Recovering the final facility and customer decisions Once the bounds converge, we can report the optimal facility configuration and its cost. The dual subproblem provides the assignment cost and the coefficients required to generate cuts, but it does not directly return the primal assignment variables (yijy_{ij}). We therefore solve the primal assignment subproblem one final time, fixing the facilities at the best solution found by Benders. # ============================================================ # Display the final Benders solution # ============================================================ print("\nFinal Benders solution") print("=" * 50) for i in facilities: status = "Open" if best_x_solution[i] else "Closed" print(f"Facility {i}: {status}") print(f"\nFixed opening cost: {best_fixed_cost:,.2f}") print(f"Assignment cost: {best_assignment_cost:,.2f}") print(f"Total cost: {best_total_cost:,.2f}") print(f"Iterations: {iteration}") # ============================================================ # Build the final primal assignment subproblem # ============================================================ assignment_model = pyo.ConcreteModel() assignment_model.FACILITIES = pyo.Set( initialize=facilities, ordered=True ) assignment_model.CUSTOMERS = pyo.Set( initialize=customers, ordered=True ) assignment_model.assignment_cost = pyo.Param( assignment_model.FACILITIES, assignment_model.CUSTOMERS, initialize=assignment_cost ) assignment_model.x_bar = pyo.Param( assignment_model.FACILITIES, initialize=best_x_solution ) assignment_model.y = pyo.Var( assignment_model.FACILITIES, assignment_model.CUSTOMERS, domain=pyo.NonNegativeReals ) assignment_model.AssignEachCustomer = pyo.Constraint( assignment_model.CUSTOMERS, rule=lambda model, j: sum( model.y[i, j] for i in model.FACILITIES ) == 1 ) assignment_model.UseOnlyOpenFacilities = pyo.Constraint( assignment_model.FACILITIES, assignment_model.CUSTOMERS, rule=lambda model, i, j: ( model.y[i, j] 1e-6 ] assignment_table = pd.DataFrame(assignment_rows) display( assignment_table.style.format({ "assignment": "{:.2f}", "demand": "{:,.0f}", "assignment cost": "{:,.2f}" }) ) The final output identifies the facilities selected by the master problem and separates the total objective into fixed opening and customer-assignment costs. We then fix these facility decisions in the primal subproblem and recover the corresponding (yijy_{ij}) values. The resulting table shows which open facility serves each customer, together with the associated demand and assignment cost. customerfacilityassignmentdemandassignment cost0C1F21.002537.001C2F21.003562.002C3F21.001215.003C4F21.0041,024.004C5F21.002113.005C6F21.003129.006C7F21.002195.007C8F21.001115.008C9F21.003840.009C10F21.002691.0010C11F21.0041,152.0011C12F21.001437.0012C13F41.003206.0013C14F41.002113.0014C15F41.004525.0015C16F41.001131.0016C17F51.002138.0017C18F51.003170.0018C19F51.002262.0019C20F51.004588.00Table 1. Optimal Solution, open facilities and assingments. The assignment table is useful for verification, but a figure provides a much clearer overview of the solution. In the next cell, we create a simple plot showing which facilities are open, which are closed, and how customers are assigned to the selected facilities. # ============================================================ # Visualize the final facility-location solution # ============================================================ fig, ax = plt.subplots(figsize=(12, 8)) open_facilities = [i for i in facilities if best_x_solution[i] == 1] closed_facilities = [i for i in facilities if best_x_solution[i] == 0] # Plot assignment connections for _, row in assignment_table.iterrows(): facility, customer = row["facility"], row["customer"] fx, fy = facility_coordinates[facility] cx, cy = customer_coordinates[customer] ax.plot( [fx, cx], [fy, cy], linewidth=1.2, alpha=0.55 * row["assignment"], zorder=1 ) # Plot customers ax.scatter( [customer_coordinates[j][0] for j in customers], [customer_coordinates[j][1] for j in customers], marker="o", s=70, label="Customers", zorder=3 ) # Plot open facilities ax.scatter( [facility_coordinates[i][0] for i in open_facilities], [facility_coordinates[i][1] for i in open_facilities], marker="*", s=350, label="Open facilities", zorder=4 ) # Plot closed facilities ax.scatter( [facility_coordinates[i][0] for i in closed_facilities], [facility_coordinates[i][1] for i in closed_facilities], marker="X", s=150, alpha=0.45, label="Closed facilities", zorder=2 ) # Add customer labels for customer in customers: x, y = customer_coordinates[customer] ax.annotate( customer, (x, y), xytext=(5, 5), textcoords="offset points", fontsize=9 ) # Add facility labels for facility in facilities: x, y = facility_coordinates[facility] status = "Open" if best_x_solution[facility] else "Closed" ax.annotate( f"{facility}\n{status}", (x, y), xytext=(7, 7), textcoords="offset points", fontsize=10, fontweight="bold" ) ax.set_title("Optimal Uncapacitated Facility Location Solution", fontsize=15) ax.set_xlabel("X coordinate") ax.set_ylabel("Y coordinate") ax.set_aspect("equal", adjustable="box") ax.grid(alpha=0.25) ax.legend() plt.tight_layout() plt.show() The plot provides a bird’s-eye view of the optimal solution. Open facilities are shown with star markers, closed facilities with X markers, and customers with circles. The connecting segments indicate which facility serves each customer, making it easy to interpret the geographic structure of the solution at a glance. Figure 2. Eye view of the optimal solution with the open facilities and customers assigned to each (Image generated by the author) Conclusions In this article, we explored the central mechanism of classical Benders decomposition. The master problem selected which facilities to open, the subproblem evaluated the resulting assignment cost, and the dual solution generated optimality cuts that gradually corrected the master’s initially optimistic estimate. We then implemented the complete algorithm in Python using Pyomo and the open-source HiGHS solver. Benders decomposition is an extraordinarily powerful technique, but it is not guaranteed to outperform the equivalent monolithic formulation. For a small instance, a modern solver may solve the complete model almost immediately. In that situation, constructing separate master and subproblems, transferring solutions between them, generating cuts, and repeatedly solving the models may introduce more computational overhead than benefit. Decomposition should therefore not be adopted simply because it appears more sophisticated. It is not a silver bullet either. Benders decomposition requires a problem with an exploitable structure, where fixing a set of complicating variables leaves a subproblem that is substantially easier to solve. Even when this structure exists, performance depends heavily on the quality of the generated cuts. Weak cuts may provide very little information, forcing the algorithm to complete many iterations before the master obtains a sufficiently accurate representation of the subproblem. In practice, much of the art of designing an effective Benders algorithm lies in producing stronger cuts, removing redundant ones, generating multiple cuts, or exploiting additional problem-specific knowledge. The implementation presented in this article is commonly known as an outer-loop Benders algorithm. Python explicitly solves the master, solves the subproblem, adds a cut, and then solves the master again. This approach is extremely useful for learning because every step of the algorithm remains visible. However, the master is repeatedly passed back to the solver instead of allowing the entire procedure to operate directly inside a single branch-and-cut search. Commercial solvers such as Gurobi and CPLEX support callback mechanisms through which cuts can be generated while the solver is exploring the branch-and-cut tree. Rather than repeatedly terminating and restarting the master optimization process, a callback can inspect candidate or relaxation solutions and add lazy constraints or user cuts during the search. Depending on the problem and the strength of the cuts, this can produce a substantial improvement in computational performance. Nevertheless, the value of Benders decomposition is not limited to solving a problem faster. Sometimes decomposition makes it possible to solve a problem that cannot even be represented as a monolithic model on the available hardware. I am currently encountering precisely this situation in one of my research projects involving a variation of the facility location problem. On our research machine, which has 128 GB of RAM, attempting to construct the complete monolithic model exhausts the available memory before the solver can even begin optimizing it. The machine simply reports that there is insufficient memory to build the formulation. With Benders decomposition, however, we do not need to create and maintain all master and operational variables simultaneously. By separating the problem and generating only the information required by the master, we can tackle instances that would otherwise remain completely inaccessible. This is another side of the power of decomposition. Sometimes Benders helps us solve a model more efficiently. Sometimes it allows us to solve a model that cannot fit in memory in the first place. To keep this first introduction as friendly as possible, we deliberately considered a setting in which the subproblem was always feasible. By requiring at least one facility to open and assuming that every facility could serve every customer, each master decision could always be translated into a valid assignment. The only question was how expensive that assignment would be. But Benders decomposition is not always this well behaved. What happens when the master proposes a strategic decision for which no feasible operational plan exists? In that situation, the subproblem cannot return an optimal cost because there is no feasible solution to evaluate. Instead, the algorithm requires a different type of feedback: a feasibility cut that blocks the current decision and, ideally, an entire family of similarly impossible decisions. That will be the subject of How Benders Decomposition Works, Part II: Feasibility Cuts. We will study the capacitated facility location problem, where each facility can serve only a limited amount of demand. The master may therefore open insufficient capacity and render the assignment subproblem infeasible. We will examine how feasibility cuts allow Benders decomposition to learn from these failures and continue searching for an optimal solution. I sincerely hope you found this article useful and that it made Benders decomposition feel a little less mysterious. Feel free to leave a comment or show your appreciation with a clap 👏. You can also follow the latest work from Sávila Education and connect with me on LinkedIn. All the code and data used in this article can be found in the accompanying GitHub repository. Thank you for taking the time to read. See you in Part II. References [1] McCloskey, J. F. (1987). U.S. operations research in World War II. Operations Research, 35(6), 910-925. DOI: 10.1287/opre.35.6.910. [2] Dantzig, G. B. (1982). Reminiscences about the origins of linear programming. Operations Research Letters, 1(2), 43-48. DOI: 10.1016/0167-6377(82)90043-8. [3] Dantzig, G. B. (1951). Maximization of a linear function of variables subject to linear inequalities. In T. C. Koopmans (Ed.), Activity analysis of production and allocation (pp. 339-347). John Wiley & Sons. [4] Aardal, K. I., Hurkens, C. A. J., & Lenstra, J. K. (2025). Jacques Benders and his decomposition algorithm. Operations Research Letters, 63, Article 107361. DOI: 10.1016/j.orl.2025.107361. [5] Benders, J. F. (1960). Partitioning in mathematical programming [Doctoral dissertation, Utrecht University]. [6] Benders, J. F. (1962). Partitioning procedures for solving mixed-variables programming problems. Numerische Mathematik, 4, 238-252. DOI: 10.1007/BF01386316. [7] Rahmaniani, R., Crainic, T. G., Gendreau, M., & Rei, W. (2017). The Benders decomposition algorithm: A literature review. European Journal of Operational Research, 259(3), 801-817. DOI: 10.1016/j.ejor.2016.12.005.