How Benders Decomposition Works, Part II: Feasibility Cuts

of this series, we explored Benders decomposition in its friendliest possible setting. The master problem made strategic decisions, the subproblem evaluated their operational consequences, and every master solution produced a feasible operational plan. The master could underestimate the true operational cost or propose a very poor decision, but the subproblem could always respond with a valid objective value. That value allowed us to generate an optimality cut. The subproblem effectively told the master that its decision was feasible, but that operating under it would cost more than expected. The master incorporated this feedback, improved its approximation of the operational cost, and tried again. Through this iterative exchange, the lower and upper bounds gradually converged to the optimal solution. Many optimization problems, however, are not nearly so cooperative. Suppose, we are dealing with a facility location problem with limited production capacities on each facility, and that the master decides which facilities to open, while the subproblem must assign customer demand to them. The master may select the cheapest facility and obtain an apparently attractive fixed cost, only for the subproblem to discover that the selected facility does not have enough capacity to serve all customers. The same situation can arise when certain customers can only be served by specific facilities, when too few workers are hired, when insufficient machines are activated, or when the selected vehicles cannot satisfy all transportation requirements. In these cases, if we only have optimality cuts, then the subproblem cannot tell the master how expensive its decision will be because no feasible operational solution exists. The problem is no longer that the operational cost has been underestimated. The strategic decision itself cannot be implemented. This creates an important difficulty for Benders decomposition. An infeasible subproblem cannot provide the usual optimality cut, and the corresponding master solution cannot be used to update the upper bound. Simply reporting that the subproblem is infeasible is also insufficient. The master needs mathematical information that explains which strategic decisions must be avoided and how it should modify its next proposal. Benders feasibility cuts provide this information. Unlike optimality cuts, which improve the approximation of the operational cost, feasibility cuts restrict the master to decisions that admit at least one feasible operational solution. A strong feasibility cut does more than reject the current master solution. It captures the structural reason for its failure and can exclude many other decisions that would produce the same operational impossibility. The mathematical tool that allows us to construct these cuts is Farkas’ lemma. This fundamental result from linear programming provides a certificate proving that a system of linear constraints has no feasible solution. The certificate can then be transformed into a constraint involving the master variables, converting the failure of the subproblem into useful feedback for the next master iteration. In this second article, we will begin with a small optimization problem in which the initial master decision makes the subproblem infeasible. We will solve the problem manually, introduce Farkas’ lemma as a theorem of alternatives, construct an explicit certificate of infeasibility, and transform that certificate into a Benders feasibility cut. We will then continue the algorithm and observe how feasibility and optimality cuts work together until the lower and upper bounds converge. Once the mechanism is fully visible, we will apply the same reasoning to the capacitated facility location problem. Unlike the uncapacitated model studied in Part I, opening at least one facility no longer guarantees that every customer can be served. Finally, we will implement the complete algorithm that solves the capacitated facility location problem in Python using Pyomo and the open-source HiGHS solver. The resulting Benders loop will generate feasibility cuts whenever the operational subproblem is impossible, and optimality cuts whenever a feasible operational plan exists. Part I showed how Benders decomposition learns that a decision is more expensive than expected. Part II shows how it learns that a decision cannot work at all. 2. A small Toy problem with an infeasible subproblem Before moving to capacitated facility location problem, let us first examine a small optimization toy problem in which the master proposes a decision that cannot support any feasible operational solution. The example is intentionally simple so that we can follow every step manually and focus entirely on the new difficulty introduced by infeasibility. We will first formulate the problem and solve its initial master problem. We will then show why the corresponding subproblem is infeasible, use Farkas’ lemma to construct a certificate of infeasibility, and transform that certificate into a Benders feasibility cut. Finally, we will continue the algorithm and observe how feasibility and optimality cuts work together until the bounds converge. 2.1 The toy optimization problem Consider a small planning problem in which two resources can be activated to satisfy a demand of six units. The first resource can provide up to six units but is relatively expensive to activate and operate. The second resource is cheaper, but it can provide at most four units and therefore cannot satisfy the demand by itself Let (x1x_1) and (x2x_2) indicate whether the two resources are activated. The continuous variables (y1y_1) and (y2y_2) represent the quantities supplied by each resource. The complete optimization problem is: min⁡5x1+x2+2y1+y2\min \quad 5x_1+x_2+2y_1+y_2 subject to y1+y2≥6,y_1+y_2\geq 6, y1≤6x1,y_1\leq 6x_1, y2≤4x2,y_2\leq 4x_2, x1+x2≥1,x_1+x_2\geq 1, x1,x2∈{0,1},x_1,x_2\in\{0,1\}, y1,y2≥0.y_1,y_2\geq 0. Activating the first resource costs five units, while activating the second costs only one. Once activated, each unit supplied by the first resource costs two units, whereas each unit supplied by the second costs one. The first constraint requires the two resources to jointly satisfy the demand of six units. The next two constraints limit the quantity supplied by each resource according to its activation decision and available capacity. Finally, at least one resource must be activated. The second resource initially appears particularly attractive because it has both the lowest activation cost and the lowest operating cost. However, its capacity is only four units. As we will see, this encourages the initial master problem to select a decision that looks inexpensive but cannot satisfy the operational requirements. 2.2 The initial master problem Benders decomposition begins by placing the strategic activation 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 the operational cost. The initial master problem is: min⁡5x1+x2+θ\min \quad 5x_1+x_2+\theta subject to x1+x2≥1,x_1+x_2\geq 1, θ≥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 feasibility or optimality cuts. It only knows that at least one resource must be activated and that the operational cost cannot be negative. Because the second resource has the lowest activation cost, the initial master solution is: x‾1=0,x‾2=1,θ=0.\bar{x}_1=0, \qquad \bar{x}_2=1, \qquad \theta=0. The corresponding objective value is, LB=1.LB=1. This value provides a lower bound on the optimal objective because the master is currently ignoring the true operational consequences of its decision. The master believes that activating only the second resource is sufficient and that the operational cost can be zero. The next step is to send the decision (x‾=(0,1)\bar{x}=(0,1)) to the subproblem and determine whether a feasible operational plan exists. 2.3 The first subproblem is infeasible The master has proposed the strategic decision: x‾1=0,x‾2=1.\bar{x}_1=0, \qquad \bar{x}_2=1. To evaluate this decision, we fix the activation variables at these values and optimize only the operational variables (y1y_1) and (y2y_2). The subproblem objective is: Q(x‾)=min⁡2y1+y2Q(\bar{x})=\min \quad 2y_1+y_2 subject to y1+y2≥6,y_1+y_2\geq 6, y1≤6x‾1,y_1\leq 6\bar{x}_1, y2≤4x‾2,y_2\leq 4\bar{x}_2, y1,y2≥0.y_1,y_2\geq 0. Substituting the current master solution into the capacity constraints gives: y1≤6(0)=0,y_1\leq 6(0)=0, y2≤4(1)=4.y_2\leq 4(1)=4. The complete subproblem is therefore: Q(0,1)=min⁡2y1+y2Q(0,1)=\min \quad 2y_1+y_2 subject to y1+y2≥6,y_1+y_2\geq 6, y1≤0,y_1\leq 0, y2≤4,y_2\leq 4, y1,y2≥0.y_1,y_2\geq 0. The nonnegativity constraint and the upper bound on (y1y_1) force: y1=0.y_1=0. The largest possible value of (y2y_2) is: y2=4.y_2=4. Consequently, the maximum quantity that can be supplied is: y1+y2≤4.y_1+y_2\leq 4. However, the demand constraint requires: y1+y2≥6.y_1+y_2\geq 6. These two conditions cannot be satisfied simultaneously. The subproblem is therefore infeasible. This result is fundamentally different from the situations studied in Part I. The subproblem does not return a large operational cost. It returns no operational cost at all because the master decision cannot support any feasible operational plan. The current solution cannot be used to update the upper bound, and the usual optimality cut cannot be generated. Instead, we need a mathematical certificate that explains why the subproblem is infeasible and transforms that explanation into a new constraint for the master problem. 2.4 Understanding Farkas’ lemma The subproblem generated by the master decision (x‾=(0,1)\bar{x}=(0,1)) is infeasible. A solver can report this status, but the word “infeasible” alone does not help the master choose a better solution. Benders decomposition requires mathematical information that identifies the source of the infeasibility and converts it into a new constraint. Farkas’ lemma provides this information. It is a theorem of alternatives for systems of linear inequalities. In simple terms, it states that either a system of linear constraints has a feasible solution, or another vector exists that proves that no feasible solution is possible. Both alternatives cannot hold simultaneously. 2.4.1 A system of linear inequalities Consider the following system: By≥h,By\geq h, y≥0.y\geq 0. The vector (yy) contains the variables whose feasibility we want to determine. The matrix (BB) contains their coefficients, while (hh) contains the right-hand sides of the constraints. There are two possible outcomes. The first is that we can find a nonnegative vector (yy) satisfying every constraint. In this case, the system is feasible. The second possibility is that no such vector exists. Farkas’ lemma states that, when the system is infeasible, we can find another vector (rr) satisfying: r≥0,r\geq 0, B⊤r≤0,B^\top r\leq 0, h⊤r>0.h^\top r>0. The vector (rr) is called a certificate of infeasibility, a Farkas certificate, or, in some computational contexts, a dual ray. Farkas’ lemma therefore states that exactly one of the following alternatives holds. The first alternative is (the system is feasible): ∃y≥0such thatBy≥h.\exists\, y\geq 0 \quad \text{such that} \quad By\geq h. The second alternative is (the is not feasible): ∃r≥0such thatB⊤r≤0\exists\, r\geq 0 \quad \text{such that} \quad B^\top r\leq 0 and h⊤r>0.h^\top r>0. If a feasible vector (yy) exists, no Farkas certificate can exist. Conversely, if a Farkas certificate exists, no feasible vector (yy) can exist. 2.4.2 Why does the certificate prove infeasibility? The conditions imposed on (rr) may initially appear abstract. Their meaning becomes clearer if we assume, temporarily, that both a feasible solution (yy) and a certificate (rr) exist. Because we assume (yy) is feasible, we have: By≥h.By\geq h. Because every component of (rr) is nonnegative, we can multiply each constraint by its corresponding component of (rr) and add the resulting inequalities. This gives: r⊤By≥r⊤h.r^\top By\geq r^\top h. The left-hand side can also be written as: r⊤By=y⊤B⊤r.r^\top By=y^\top B^\top r. The certificate requires: B⊤r≤0.B^\top r\leq 0. Because (y≥0y\geq0), multiplying the nonpositive vector (B⊤rB^\top r) by the nonnegative vector (yy) gives: y⊤B⊤r≤0.y^\top B^\top r\leq 0. Therefore: r⊤By≤0.r^\top By\leq 0. However, the final condition of the Farkas certificate requires: r⊤h>0.r^\top h>0. Combining these results gives: 0≥r⊤By≥r⊤h>0.0\geq r^\top By\geq r^\top h>0. This is impossible. A quantity cannot be both nonpositive and strictly positive. The contradiction proves that a feasible vector (yy) and a Farkas certificate (rr) cannot exist simultaneously. Therefore, finding a vector (rr) satisfying the certificate conditions proves that the original system is infeasible. 2.4.3 Applying Farkas’ lemma to a Benders subproblem Consider a generic Benders subproblem for a fixed master solution (x‾\bar{x}): Q(x‾)=miny⁡d⊤yQ(\bar{x})=\min_y \quad d^\top y subject to By≥b−Ax‾,By\geq b-A\bar{x}, y≥0.y\geq 0. For the moment, the objective function is not important. Before minimizing the operational cost, we must determine whether there is any vector (yy) satisfying the constraints. The feasibility system is: By≥b−Ax‾,By\geq b-A\bar{x}, y≥0.y\geq 0. In the notation of Farkas’ lemma, the right-hand side is: h(x‾)=b−Ax‾.h(\bar{x})=b-A\bar{x}. If the subproblem is infeasible, Farkas’ lemma guarantees the existence of a vector (rr) satisfying: r≥0,r\geq 0, B⊤r≤0,B^\top r\leq 0, (b−Ax‾)⊤r>0.\left(b-A\bar{x}\right)^\top r>0. The first two conditions depend only on the structure of the subproblem. They ensure that (rr) is a valid combination of its constraints. The third condition depends on the current master solution (x‾\bar{x}). It shows that this particular strategic decision produces a contradiction in the operational constraints. 2.4.4 From a certificate to a feasibility cut Suppose that the vector (rkr^k) provides a certificate for the current master solution (x‾\bar{x}). We therefore have: (b−Ax‾)⊤rk>0.\left(b-A\bar{x}\right)^\top r^k>0. This inequality proves that (x‾\bar{x}) produces an infeasible subproblem. However, the objective of Benders decomposition is not merely to identify that the current decision has failed. We need to construct a constraint that prevents the master from choosing decisions that produce the same contradiction. For any master decision (xx) that admits a feasible subproblem, the Farkas contradiction must not occur. Consequently, every feasible master decision must satisfy: (b−Ax)⊤rk≤0.\left(b-Ax\right)^\top r^k\leq 0. This is a Benders feasibility cut. The cut can also be expanded as: b⊤rk−x⊤A⊤rk≤0.b^\top r^k-x^\top A^\top r^k\leq 0. Equivalently: x⊤A⊤rk≥b⊤rk.x^\top A^\top r^k\geq b^\top r^k. The current solution (x‾\bar{x}) violates this inequality because: (b−Ax‾)⊤rk>0.\left(b-A\bar{x}\right)^\top r^k>0. The cut therefore removes the current infeasible master solution. More importantly, it can also remove other master decisions that would create the same operational contradiction. This is why the Farkas certificate is more useful than a simple infeasibility status. It explains how the master variables are connected to the failure of the subproblem and translates that explanation into a valid master constraint. 2.4.5 How can the certificate be obtained? For a very small problem, we may be able to identify a certificate manually by inspecting and combining the constraints. This is what we will do with our toy example. For larger problems, the certificate can be obtained computationally. One possibility is to request an infeasibility ray directly from the linear programming solver. Some solver interfaces provide this information when an LP is infeasible. Another possibility is to solve an auxiliary optimization problem that searches for a Farkas certificate. For the current master solution (x‾\bar{x}), consider: maxr⁡(b−Ax‾)⊤r\max_r \quad \left(b-A\bar{x}\right)^\top r subject to B⊤r≤0,B^\top r\leq 0, ∑iri≤1,\sum_i r_i\leq 1, r≥0.r\geq 0. The normalization constraint ∑iri≤1\sum_i r_i\leq 1, is necessary because the certificate conditions define a cone. If (rr) is a certificate, then any positive multiple of (rr) is also a certificate. Without normalization, the auxiliary problem could become unbounded. If the optimal objective value of this auxiliary problem is strictly positive, the resulting vector (rr) satisfies: (b−Ax‾)⊤r>0,\left(b-A\bar{x}\right)^\top r>0, and therefore certifies that the subproblem is infeasible. The corresponding feasibility cut is then: (b−Ax)⊤r≤0.\left(b-Ax\right)^\top r\leq 0. If the auxiliary objective cannot become positive, no Farkas certificate exists for the current right-hand side. Under the conditions of Farkas’ lemma, this means that the original feasibility system admits a solution (yy). The complete procedure can therefore be summarized as follows. First, fix the current master solution (x‾\bar{x}). Second, determine whether the operational constraints admit a feasible (yy). If they do not, find a Farkas certificate (rr). Finally, replace (x‾\bar{x}) with the master variables (xx) in the certificate condition and reverse the strict violation into a valid feasibility requirement: (b−Ax‾)⊤r>0\left(b-A\bar{x}\right)^\top r>0 becomes (b−Ax)⊤r≤0.\left(b-Ax\right)^\top r\leq 0. We can now apply each of these steps to the infeasible subproblem of our toy example. 2.5 Applying Farkas’ lemma to the toy subproblem The master solution (x‾=(0,1)\bar{x}=(0,1)) produced the following infeasible subproblem: Q(0,1)=min⁡2y1+y2Q(0,1)=\min \quad 2y_1+y_2 subject to y1+y2≥6,y_1+y_2\geq 6, y1≤0,y_1\leq 0, y2≤4,y_2\leq 4, y1,y2≥0.y_1,y_2\geq 0. To apply Farkas’ lemma, we temporarily ignore the objective function. The question is not yet how much the operational solution costs, but whether any operational solution exists. 2.5.1 Rewriting the subproblem in a common inequality direction Farkas’ lemma was presented for a system of the form: By≥h.By\geq h. The demand constraint already uses the required greater-than-or-equal-to direction: y1+y2≥6.y_1+y_2\geq 6. The capacity constraints must therefore be multiplied by (-1). Before substituting the current master decision, they can be written as: −y1≥−6x1,-y_1\geq -6x_1, −y2≥−4x2.-y_2\geq -4x_2. The feasibility system for a generic master decision (x) is consequently: y1+y2≥6,y_1+y_2\geq 6, −y1≥−6x1,-y_1\geq -6x_1, −y2≥−4x2,-y_2\geq -4x_2, y1,y2≥0.y_1,y_2\geq 0. Using matrix notation, the vector of operational variables and The coefficient matrix is: y=[y1y2]B=[11−100−1].y= \begin{bmatrix} y_1\\ y_2 \end{bmatrix} \qquad B= \begin{bmatrix} 1 & 1\\ -1 & 0\\ 0 & -1 \end{bmatrix}. The right-hand side depends on the master variables: h(x)=[6−6x1−4x2].h(x)= \begin{bmatrix} 6\\ -6x_1\\ -4x_2 \end{bmatrix}. The operational feasibility system can therefore be written as: By≥h(x),By\geq h(x), y≥0.y\geq 0. For the current master solution (x‾=(0,1)\bar{x}=(0,1)), the right-hand side becomes: h(x‾)=[60−4].h(\bar{x})= \begin{bmatrix} 6\\ 0\\ -4 \end{bmatrix}. 2.5.2 Constructing the certificate We associate one Farkas multiplier with each constraint. Let: r=[r1r2r3].r= \begin{bmatrix} r_1\\ r_2\\ r_3 \end{bmatrix}. The multiplier (r1r_1) is associated with the demand constraint, (r2r_2) with the capacity of the first resource, and (r3r_3) with the capacity of the second resource. A valid Farkas certificate must first satisfy r1,r2,r3≥0.r_1,r_2,r_3\geq 0. It must also satisfy, B⊤r≤0.B^\top r\leq 0. Using the matrix (BB), we obtain: B⊤r=[1−1010−1][r1r2r3].B^\top r= \begin{bmatrix} 1 & -1 & 0\\ 1 & 0 & -1 \end{bmatrix} \begin{bmatrix} r_1\\ r_2\\ r_3 \end{bmatrix}. Therefore: B⊤r=[r1−r2r1−r3].B^\top r= \begin{bmatrix} r_1-r_2\\ r_1-r_3 \end{bmatrix}. The condition (B⊤r≤0B^\top r\leq0) becomes: r1−r2≤0,r1−r3≤0.r_1-r_2\leq 0, \qquad r_1-r_3\leq 0. Equivalently: r2≥r1,r3≥r1.r_2\geq r_1, \qquad r_3\geq r_1. We now need the certificate to prove that the current right-hand side is inconsistent. This requires: h(x‾)⊤r>0.h(\bar{x})^\top r>0. Substituting the current master solution gives: h(x‾)⊤r=6r1−4r3.h(\bar{x})^\top r = 6r_1-4r_3. We must therefore find nonnegative multipliers satisfying: r2≥r1,r3≥r1,6r1−4r3>0.r_2\geq r_1, \qquad r_3\geq r_1, \qquad 6r_1-4r_3>0. The simplest choice is: r1=1,r2=1,r3=1.r_1=1, \quad r_2=1, \quad r_3=1. Thus: r=[111].r= \begin{bmatrix} 1\\ 1\\ 1 \end{bmatrix}. 2.5.3 Transforming the certificate into a feasibility cut We first verify the nonnegativity condition. The Farkas certificate (rr), was obtained for the current master solution (x‾=(0,1)\bar{x}=(0,1)). To create a Benders cut, we return to the generic right-hand side: h(x)=[6−6x1−4x2].h(x)= \begin{bmatrix} 6\\ -6x_1\\ -4x_2 \end{bmatrix}. Using the certificate (r=(1,1,1)r=(1,1,1)), we obtain: h(x)⊤r=6−6x1−4x2.h(x)^\top r = 6-6x_1-4x_2. Any master decision that produces a feasible subproblem must prevent this expression from being strictly positive. It must therefore satisfy: 6−6x1−4x2≤0.6-6x_1-4x_2\leq 0. Rearranging the inequality gives: 6x1+4x2≥6.6x_1+4x_2\geq 6. This is the Benders feasibility cut generated by the infeasible subproblem. The current master (x‾=(0,1)\bar{x}=(0,1)) solution violates the cut because: 6(0)+4(1)=40.\sum_{j\in J} d_jr_j – \sum_{i\in I} K_i\bar{x}_i s_i >0. The final inequality certifies that the current facility-opening decision produces an infeasible transportation system. Any master solution that admits a feasible subproblem must instead satisfy: ∑j∈Jdjrj−∑i∈IKixisi≤0.\sum_{j\in J} d_jr_j – \sum_{i\in I} K_i x_i s_i \leq0. Equivalently, the Benders feasibility cut is: ∑i∈IKisixi≥∑j∈Jdjrj.\sum_{i\in I} K_i s_i x_i \geq \sum_{j\in J} d_jr_j. This cut excludes the current infeasible facility-opening decision and potentially many other decisions that would produce the same capacity or connectivity problem. When every facility can serve every customer, one possible certificate may lead to the familiar total-capacity requirement: ∑i∈IKixi≥∑j∈Jdj.\sum_{i\in I} K_i x_i \geq \sum_{j\in J} d_j. With restricted facility-customer connections, however, total capacity may not be sufficient. Capacity must also be available at facilities capable of reaching the customers that require it. Farkas cuts can capture these more specific sources of infeasibility. The complete Benders algorithm therefore has two branches. If the subproblem is infeasible, it generates a feasibility cut and returns to the master. If the subproblem is feasible, it updates the upper bound and generates an optimality cut. We will now implement this complete procedure in Python using Pyomo and HiGHS. 4. Implementing feasibility cuts with Pyomo and HiGHS We now have all the mathematical machinery required to implement the extended Benders decomposition algorithm discussed in the previous sections. The master problem will select the facilities to open, the feasibility subproblem will determine whether the resulting transportation problem can be operated at all, and the dual subproblem will evaluate its transportation cost whenever a feasible operational solution exists. To illustrate the complete procedure, we will solve a small capacitated 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 the master problem, the auxiliary Farkas problem used to generate feasibility cuts, and the dual transportation subproblem used to generate optimality cuts. Finally, we will combine these components into the complete iterative Benders procedure. 4.1 Loading and preparing the instance We begin by loading the capacitated facility location instance directly from the GitHub repository and converting the JSON data into the sets and dictionaries that will be used throughout the Benders implementation. In addition to facility costs, capacities, customer demands, and transportation costs, we also construct the service neighborhoods describing which facilities are allowed to serve each customer. !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 JSON_URL = ( "https://raw.githubusercontent.com/ceche1212/" "Benders_Tutorials_TDS/refs/heads/main/Data/" "cflp_benders_5x20_instance.json" ) instance = requests.get(JSON_URL, timeout=30).json() facilities = [f["id"] for f in instance["facilities"]] customers = [c["id"] for c in instance["customers"]] fixed_cost = {f["id"]: f["fixed_cost"] for f in instance["facilities"]} capacity = {f["id"]: f["capacity"] for f in instance["facilities"]} demand = {c["id"]: c["demand"] for c in instance["customers"]} facility_coordinates = { f["id"]: (f["x"], f["y"]) for f in instance["facilities"] } customer_coordinates = { c["id"]: (c["x"], c["y"]) for c in instance["customers"] } transportation_cost = { (a["facility"], a["customer"]): a["cost_per_unit"] for a in instance["transportation_costs"] } arcs = list(transportation_cost) facilities_for_customer = { j: [i for i in facilities if (i, j) in transportation_cost] for j in customers } customers_for_facility = { i: [j for j in customers if (i, j) in transportation_cost] for i in facilities } assert all(facilities_for_customer[j] for j in customers) print(f"Facilities: {len(facilities)}") print(f"Customers: {len(customers)}") print(f"Transportation arcs: {len(arcs)}") print(f"Total demand: {sum(demand.values())}") print(f"Total capacity: {sum(capacity.values())}") The resulting data structures closely follow the mathematical notation introduced earlier. The dictionaries fixed_cost, capacity, and demand contain the parameters (fif_i), (KiK_i), and (djd_j), while transportation_cost stores (cijc_{ij}) only for permitted facility-customer connections. The neighborhood dictionaries are particularly important because the transportation network is not complete. facilities_for_customer[j] identifies the facilities capable of serving customer (jj), while customers_for_facility[i] identifies the customers reachable from facility (ii). For this instance, the network contains five candidate facilities, twenty customers, 57 permitted transportation arcs, a total demand of 49 units, and 89 units of potential capacity. Having more total capacity than demand does not guarantee feasibility, since the master may still open the wrong combination of facilities or leave insufficient capacity available to particular groups of customers. 4.2 Building the initial master problem We can now construct the Benders master problem. As in the mathematical formulation, the binary variables (xix_i) determine which facilities are opened, while (θ\theta) represents the master’s current approximation of the transportation cost. We also create two initially empty constraint lists that will store the feasibility and optimality cuts generated throughout the algorithm. Importantly, we deliberately do not impose any capacity or connectivity requirements at this stage. The initial master knows only the facility-opening costs and that (θ\theta) cannot be negative. The purpose of the feasibility cuts will be precisely to teach the master which combinations of facilities can support a feasible transportation plan. master = pyo.ConcreteModel() master.x = pyo.Var(facilities, domain=pyo.Binary) master.theta = pyo.Var(domain=pyo.NonNegativeReals) master.FeasibilityCuts = pyo.ConstraintList() master.OptimalityCuts = pyo.ConstraintList() master.TotalCost = pyo.Objective( expr=sum(fixed_cost[i] * master.x[i] for i in facilities) + master.theta, sense=pyo.minimize ) solver = SolverFactory("appsi_highs") results = solver.solve(master) if results.solver.termination_condition != TerminationCondition.optimal: raise RuntimeError("The master problem was not solved to optimality.") x_solution = { i: int(round(pyo.value(master.x[i]))) for i in facilities } theta_value = pyo.value(master.theta) master_objective = pyo.value(master.TotalCost) print("Initial master solution:", x_solution) print(f"Theta: {theta_value:,.2f}") print(f"Lower bound: {master_objective:,.2f}") Because the initial master contains no Benders cuts, the cheapest solution is to open no facilities and set (θ=0\theta=0). The resulting master objective is therefore zero, providing our first lower bound. Of course, this solution cannot possibly serve any customer. That is intentional. Rather than manually adding constraints that prevent such configurations, we will send the master solution to the feasibility subproblem and let Farkas’ lemma generate a cut explaining why it cannot work. This is the main difference with the implementation in Part I, where the structure of the model guaranteed that every master solution produced a feasible assignment subproblem. 4.3 Building the normalized Farkas subproblem The next component is the auxiliary problem used to detect infeasibility and generate feasibility cuts. For a fixed master solution (x‾\bar{x}), we search for nonnegative Farkas multipliers (rjr_j) and (sis_i) satisfying the conditions derived earlier. Because these multipliers define a cone, any valid certificate can be multiplied by an arbitrary positive constant, so we add a normalization constraint to obtain a bounded linear program. The objective measures the violation associated with the current master solution. If its optimal value is strictly positive, the resulting multipliers provide a certificate that the transportation subproblem is infeasible and can immediately be converted into a Benders feasibility cut. def build_farkas_subproblem(x_bar): m = pyo.ConcreteModel() m.r = pyo.Var(customers, domain=pyo.NonNegativeReals) m.s = pyo.Var(facilities, domain=pyo.NonNegativeReals) m.FarkasCondition = pyo.Constraint( arcs, rule=lambda m, i, j: m.r[j] - m.s[i] 1e-8: print(f"r[{j}] = {value:.6f}") for i, value in s_solution.items(): if value > 1e-8: print(f"s[{i}] = {value:.6f}") master.FeasibilityCuts.add( sum(capacity[i] * s_solution[i] * master.x[i] for i in facilities) >= sum(demand[j] * r_solution[j] for j in customers) ) The constraints (rj−si≤0r_j-s_i\leq0) reproduce the Farkas condition associated with every permitted transportation arc. The objective then compares the weighted customer demand with the capacity made available by the current facility configuration (x‾\bar{x}). For the initial master solution, the optimal objective is positive, confirming that no feasible transportation plan exists. The corresponding values of (rjr_j) and (sis_i) form the Farkas certificate. We then replace the fixed values (x‾i\bar{x}_i) by the master variables (xix_i) and add the resulting inequality directly to master.FeasibilityCuts. The master has therefore received its first piece of operational information. Instead of simply being told that its current solution is infeasible, it now receives a constraint describing a family of facility configurations that must be avoided. 4.4 Building the dual transportation subproblem Once the Farkas problem confirms that the current facility configuration admits a feasible transportation plan, we can evaluate its operating cost. As in Part I, we do this through the dual of the transportation subproblem. The variables (αj\alpha_j) are associated with the customer-demand constraints, while (βi\beta_i) correspond to the transformed facility-capacity constraints. For a fixed master solution (x‾\bar{x}), solving this dual gives the minimum transportation cost through strong duality. Its optimal multipliers will then be used to generate the usual Benders optimality cut. def build_dual_subproblem(x_bar): m = pyo.ConcreteModel() m.ARCS = pyo.Set(dimen=2, initialize=arcs) m.alpha = pyo.Var(customers, domain=pyo.NonNegativeReals) m.beta = pyo.Var(facilities, domain=pyo.NonNegativeReals) m.DualFeasibility = pyo.Constraint( m.ARCS, rule=lambda m, i, j: m.alpha[j] - m.beta[i] tolerance: r = {j: pyo.value(farkas.r[j]) for j in customers} s = {i: pyo.value(farkas.s[i]) for i in facilities} master.FeasibilityCuts.add( sum(capacity[i] * s[i] * master.x[i] for i in facilities) >= sum(demand[j] * r[j] for j in customers) ) transportation = None cut_type = "Feasibility" gap = float("inf") else: # 3. Evaluate feasible solution dual = build_dual_subproblem(x_bar) solver.solve(dual) transportation = pyo.value(dual.TransportationCost) total_cost = fixed + transportation if total_cost = sum(demand[j] * alpha[j] for j in customers) - sum( capacity[i] * beta[i] * master.x[i] for i in facilities ) ) cut_type = "Optimality" iteration_results.append({ "iteration": iteration, "open facilities": ", ".join( i for i in facilities if x_bar[i] ) or "None", "theta": theta, "Farkas objective": certificate, "transportation cost": transportation, "cut added": cut_type, "lower bound": lower_bound, "upper bound": upper_bound }) print( f"{iteration:2d} | " f"open={[i for i in facilities if x_bar[i]]} | " f"cut={cut_type} | " f"LB={lower_bound:,.2f} | " f"UB={upper_bound:,.2f}" ) if certificate = demand[j] ) transport.Capacity = pyo.Constraint( facilities, rule=lambda m, i: sum(m.y[i, j] for j in customers_for_facility[i]) 1e-6 } print(f"Transportation cost: {pyo.value(transport.Cost):,.2f}") print(f"Total cost: {best_total_cost:,.2f}") # Plot the solution fig, ax = plt.subplots(figsize=(12, 8)) max_flow = max(flows.values()) for (i, j), flow in flows.items(): xi, yi = facility_coordinates[i] xj, yj = customer_coordinates[j] ax.plot( [xi, xj], [yi, yj], linewidth=1 + 2 * flow / max_flow, alpha=0.5 ) ax.scatter( *zip(*(customer_coordinates[j] for j in customers)), s=70, label="Customers" ) open_facilities = [i for i in facilities if best_solution[i]] closed_facilities = [i for i in facilities if not best_solution[i]] ax.scatter( *zip(*(facility_coordinates[i] for i in open_facilities)), marker="*", s=350, label="Open facilities" ) ax.scatter( *zip(*(facility_coordinates[i] for i in closed_facilities)), marker="X", s=150, alpha=0.45, label="Closed facilities" ) for j in customers: ax.annotate( f"{j}\nd={demand[j]}", customer_coordinates[j], xytext=(5, 5), textcoords="offset points", fontsize=8 ) for i in facilities: ax.annotate( i, facility_coordinates[i], xytext=(7, 7), textcoords="offset points", fontweight="bold" ) ax.set( title="Optimal Capacitated Facility Location Solution", xlabel="X coordinate", ylabel="Y coordinate" ) ax.set_aspect("equal") ax.grid(alpha=0.25) ax.legend() plt.tight_layout() plt.show() The final transportation model fixes the facility decisions obtained by Benders and optimizes only the continuous flows. We retain the positive (yijy_{ij}) values and use them to draw the active transportation connections, with thicker lines representing larger flows. For this instance, the optimal solution opens facilities F1, F3, F4, and F5, while F2 remains closed. The final objective value is 14,757, combining 9,300 in facility-opening costs with 5,457 in transportation costs. Solving the original monolithic CFLP in the complete notebook produces exactly the same objective value, providing a useful numerical check that the Benders implementation has converged to the correct optimum. 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 second article, we extended the classical Benders decomposition framework to situations in which the master can propose decisions that make the operational subproblem infeasible. Using a capacitated facility location problem, we saw how Farkas’ lemma provides a certificate of infeasibility and how this certificate can be transformed into a Benders feasibility cut. This distinction is important. An optimality cut tells the master that a feasible decision is more expensive than it currently believes. A feasibility cut tells the master that a decision cannot work at all. More importantly, the Farkas certificate allows us to identify some of the structural requirements that feasible solutions must satisfy, so a single cut can potentially eliminate many master solutions that would lead to the same operational impossibility. There is, however, another way of excluding infeasible integer solutions. If we simply want to prevent the master from selecting a particular discrete solution again, we can use a no-good cut. Unlike the Farkas cuts studied here, a no-good cut does not necessarily explain the structural reason for infeasibility. It can simply remove a specific combination of integer decisions from the search space. This may sound less sophisticated, but it becomes particularly useful when the subproblem itself is no longer a linear program. So far, our decomposition has followed a very convenient structure. The difficult integer variables remain in the master, while fixing them leaves a continuous linear subproblem. When that LP is feasible, its dual variables provide the optimality cuts. When it is infeasible, Farkas’ lemma provides the feasibility cuts. In both cases, LP duality gives us the mathematical feedback required by the master. But what happens when fixing the master variables still leaves an integer or mixed-integer optimization problem? In that case, there is generally no LP dual solution that we can use to construct the classical Benders optimality cuts developed in the first two articles. This is where Logic-Based Benders Decomposition enters the picture. Instead of deriving every cut from linear programming duality, we use the structure and logic of the subproblem itself to communicate useful information back to the master. No-good cuts are one of the simplest tools available in this setting, although stronger problem-specific cuts are usually far more valuable. That will be the subject of the next article in this series. We will move from facility location to a parallel-machine scheduling problem with sequence-dependent setup times (a difficult optimization problem, that appears a lot in production and manufacturing)” . The master will decide how jobs are assigned to machines, but determining the best sequence of jobs on each machine will remain a combinatorial optimization problem. This will force us to leave the comfortable world of LP duality and construct Benders cuts using a different kind of reasoning. I sincerely hope you found this article useful and that feasibility cuts and Farkas’ lemma now feel a little less mysterious. 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 III.