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 ($x_1$) and ($x_2$) indicate whether the two resources are activated. The continuous variables ($y_1$) and ($y_2$) represent the quantities supplied by each resource. The complete optimization problem is:

$min quad 5x_1+x_2+2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$y_1leq 6x_1,$

$y_2leq 4x_2,$

$x_1+x_2geq 1,$

$x_1,x_2in{0,1},$

$y_1,y_2geq 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 ($x_1$) and ($x_2$) in the master problem. The operational variables ($y_1$) and ($y_2$) are temporarily removed and replaced by the variable ($theta$), which represents an estimate of the operational cost. The initial master problem is:

$min quad 5x_1+x_2+theta$

subject to

$x_1+x_2geq 1,$

$thetageq 0,$

$x_1,x_2in{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:

$bar{x}_1=0, qquad bar{x}_2=1, qquad theta=0.$

The corresponding objective value is, $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 ($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:

$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 ($y_1$) and ($y_2$). The subproblem objective is:

$Q(bar{x})=min quad 2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$y_1leq 6bar{x}_1,$

$y_2leq 4bar{x}_2,$

$y_1,y_2geq 0.$

Substituting the current master solution into the capacity constraints gives:

$y_1leq 6(0)=0,$

$y_2leq 4(1)=4.$

The complete subproblem is therefore:

$Q(0,1)=min quad 2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$y_1leq 0,$

$y_2leq 4,$

$y_1,y_2geq 0.$

The nonnegativity constraint and the upper bound on ($y_1$) force:

$y_1=0.$

The largest possible value of ($y_2$) is:

$y_2=4.$

Consequently, the maximum quantity that can be supplied is:

$y_1+y_2leq 4.$

However, the demand constraint requires:

$y_1+y_2geq 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 ($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:

$Bygeq h,$

$ygeq 0.$

The vector ($y$) contains the variables whose feasibility we want to determine. The matrix ($B$) contains their coefficients, while ($h$) contains the right-hand sides of the constraints.

There are two possible outcomes. The first is that we can find a nonnegative vector ($y$) 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 ($r$) satisfying:

$rgeq 0,$

$B^top rleq 0,$

$h^top r>0.$

The vector ($r$) 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):

$exists, ygeq 0 quad text{such that} quad Bygeq h.$

The second alternative is (the is not feasible):

$exists, rgeq 0 quad text{such that} quad B^top rleq 0$

and

$h^top r>0.$

If a feasible vector ($y$) exists, no Farkas certificate can exist. Conversely, if a Farkas certificate exists, no feasible vector ($y$) can exist.

2.4.2 Why does the certificate prove infeasibility?

The conditions imposed on ($r$) may initially appear abstract. Their meaning becomes clearer if we assume, temporarily, that both a feasible solution ($y$) and a certificate ($r$) exist.

Because we assume ($y$) is feasible, we have:

$Bygeq h.$

Because every component of ($r$) is nonnegative, we can multiply each constraint by its corresponding component of ($r$) and add the resulting inequalities. This gives:

$r^top Bygeq r^top h.$

The left-hand side can also be written as:

$r^top By=y^top B^top r.$

The certificate requires:

$B^top rleq 0.$

Because ($ygeq0$), multiplying the nonpositive vector ($B^top r$) by the nonnegative vector ($y$) gives:

$y^top B^top rleq 0.$

Therefore:

$r^top Byleq 0.$

However, the final condition of the Farkas certificate requires:

$r^top h>0.$

Combining these results gives:

$0geq r^top Bygeq r^top h>0.$

This is impossible. A quantity cannot be both nonpositive and strictly positive.

The contradiction proves that a feasible vector ($y$) and a Farkas certificate ($r$) cannot exist simultaneously. Therefore, finding a vector ($r$) 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 ($bar{x}$):

$Q(bar{x})=min_y quad d^top y$

subject to

$Bygeq b-Abar{x},$

$ygeq 0.$

For the moment, the objective function is not important. Before minimizing the operational cost, we must determine whether there is any vector ($y$) satisfying the constraints.

The feasibility system is:

$Bygeq b-Abar{x},$

$ygeq 0.$

In the notation of Farkas’ lemma, the right-hand side is:

$h(bar{x})=b-Abar{x}.$

If the subproblem is infeasible, Farkas’ lemma guarantees the existence of a vector ($r$) satisfying:

$rgeq 0,$

$B^top rleq 0,$

$left(b-Abar{x}right)^top r>0.$

The first two conditions depend only on the structure of the subproblem. They ensure that ($r$) is a valid combination of its constraints.

The third condition depends on the current master solution ($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 ($r^k$) provides a certificate for the current master solution ($bar{x}$). We therefore have:

$left(b-Abar{x}right)^top r^k>0.$

This inequality proves that ($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 ($x$) that admits a feasible subproblem, the Farkas contradiction must not occur. Consequently, every feasible master decision must satisfy:

$left(b-Axright)^top r^kleq 0.$

This is a Benders feasibility cut. The cut can also be expanded as:

$b^top r^k-x^top A^top r^kleq 0.$

Equivalently:

$x^top A^top r^kgeq b^top r^k.$

The current solution ($bar{x}$) violates this inequality because:

$left(b-Abar{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 ($bar{x}$), consider:

$max_r quad left(b-Abar{x}right)^top r$

subject to

$B^top rleq 0,$

$sum_i r_ileq 1,$

$rgeq 0.$

The normalization constraint $sum_i r_ileq 1$, is necessary because the certificate conditions define a cone. If ($r$) is a certificate, then any positive multiple of ($r$) 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 ($r$) satisfies:

$left(b-Abar{x}right)^top r>0,$

and therefore certifies that the subproblem is infeasible. The corresponding feasibility cut is then:

$left(b-Axright)^top rleq 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 ($y$).

The complete procedure can therefore be summarized as follows. First, fix the current master solution ($bar{x}$). Second, determine whether the operational constraints admit a feasible ($y$). If they do not, find a Farkas certificate ($r$). Finally, replace ($bar{x}$) with the master variables ($x$) in the certificate condition and reverse the strict violation into a valid feasibility requirement:

$left(b-Abar{x}right)^top r>0$

becomes

$left(b-Axright)^top rleq 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 ($bar{x}=(0,1)$) produced the following infeasible subproblem:

$Q(0,1)=min quad 2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$y_1leq 0,$

$y_2leq 4,$

$y_1,y_2geq 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:

$Bygeq h.$

The demand constraint already uses the required greater-than-or-equal-to direction:

$y_1+y_2geq 6.$

The capacity constraints must therefore be multiplied by (-1). Before substituting the current master decision, they can be written as:

$-y_1geq -6x_1,$

$-y_2geq -4x_2.$

The feasibility system for a generic master decision (x) is consequently:

$y_1+y_2geq 6,$

$-y_1geq -6x_1,$

$-y_2geq -4x_2,$

$y_1,y_2geq 0.$

Using matrix notation, the vector of operational variables and The coefficient matrix is:

$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)= begin{bmatrix} 6\ -6x_1\ -4x_2 end{bmatrix}.$

The operational feasibility system can therefore be written as:

$Bygeq h(x),$

$ygeq 0.$

For the current master solution ($bar{x}=(0,1)$), the right-hand side becomes:

$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= begin{bmatrix} r_1\ r_2\ r_3 end{bmatrix}.$

The multiplier ($r_1$) is associated with the demand constraint, ($r_2$) with the capacity of the first resource, and ($r_3$) with the capacity of the second resource.

A valid Farkas certificate must first satisfy $r_1,r_2,r_3geq 0.$ It must also satisfy, $B^top rleq 0.$ Using the matrix ($B$), we obtain:

$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^top r= begin{bmatrix} r_1-r_2\ r_1-r_3 end{bmatrix}.$

The condition ($B^top rleq0$) becomes:

$r_1-r_2leq 0, qquad r_1-r_3leq 0.$

Equivalently:

$r_2geq r_1, qquad r_3geq r_1.$

We now need the certificate to prove that the current right-hand side is inconsistent. This requires:

$h(bar{x})^top r>0.$

Substituting the current master solution gives:

$h(bar{x})^top r = 6r_1-4r_3.$

We must therefore find nonnegative multipliers satisfying:

$r_2geq r_1, qquad r_3geq r_1, qquad 6r_1-4r_3>0.$

The simplest choice is:

$r_1=1, quad r_2=1, quad r_3=1.$

Thus:

$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 ($r$), was obtained for the current master solution ($bar{x}=(0,1)$). To create a Benders cut, we return to the generic right-hand side:

$h(x)= begin{bmatrix} 6\ -6x_1\ -4x_2 end{bmatrix}.$

Using the certificate ($r=(1,1,1)$), we obtain:

$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-6x_1-4x_2leq 0.$

Rearranging the inequality gives:

$6x_1+4x_2geq 6.$

This is the Benders feasibility cut generated by the infeasible subproblem. The current master ($bar{x}=(0,1)$) solution violates the cut because:

$6(0)+4(1)=4<6.$

The cut therefore removes ($bar{x}=(0,1)$) from the master problem. Its meaning is also intuitive: the activated resources must provide at least six units of total capacity. The important point is that this constraint was not added manually based on our interpretation of the problem. It was derived systematically from the Farkas certificate returned by the infeasible subproblem. The master has now learned why its first decision failed and what condition future decisions must satisfy.

2.6 Completing the toy example

The first master solution was infeasible, but its Farkas certificate generated the following feasibility cut, $6x_1+4x_2geq 6.$ The updated master problem is then:

$min quad 5x_1+x_2+theta$

subject to

$x_1+x_2geq 1,$

$6x_1+4x_2geq 6,$

$thetageq 0,$

$x_1,x_2in{0,1}.$

The previous solution ($bar{x}=(0,1)$) is no longer feasible because it provides only four units of capacity. The least expensive remaining master solution is:

$bar{x}_1=1, qquad bar{x}_2=0, qquad theta=0.$

The corresponding master objective value is, $LB=5.$

This value becomes the new lower bound. The feasibility cut guarantees that the decision has enough total capacity, but the master still assumes that the operational cost ($theta$) can be zero. We must therefore send ($bar{x}=(1,0)$) to the subproblem.

2.6.1 Evaluating the new master decision

For the decision ($bar{x}=(1,0)$), the operational subproblem is:

$Q(1,0)=min quad 2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$y_1leq 6,$

$y_2leq 0,$

$y_1,y_2geq 0.$

The second resource is not activated, so $y_2=0.$ The entire demand must therefore be supplied by the first resource $y_1=6.$ The operational cost is:

$Q(1,0)=2(6)+0=12.$

Adding the activation cost gives a feasible solution to the original problem with total cost:

$5(1)+1(0)+12=17.$

The first feasible master decision ($bar{x}=(1,0)$) therefore provides $UB=17.$ We now have:

$LB=5, qquad UB=17.$

The subproblem is feasible, but the master estimated an operational cost ($theta$) of zero instead of twelve. We must therefore generate an optimality cut.

2.6.2 Constructing the dual subproblem

Before substituting the current master solution, the operational subproblem can be written as:

$Q(x)=min quad 2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$-y_1geq -6x_1,$

$-y_2geq -4x_2,$

$y_1,y_2geq 0.$

Let ($u_1$) be the dual variable associated with the demand constraint, while ($u_2$) and ($u_3$) are associated with the capacity constraints of the first and second resources. The dual subproblem is:

$max quad 6u_1-6x_1u_2-4x_2u_3$

subject to

$u_1-u_2leq 2,$

$u_1-u_3leq 1,$

$u_1,u_2,u_3geq 0.$

For the current master solution ($bar{x}=(1,0)$), the dual objective becomes:

$max quad 6u_1-6u_2$

subject to

$u_1-u_2leq 2,$

$u_1-u_3leq 1,$

$u_1,u_2,u_3geq 0.$

An optimal dual solution is:

$u_1=2, qquad u_2=0, qquad u_3=1.$

The corresponding dual objective value is:

$6(2)-6(0)=12.$

This matches the primal subproblem value ($Q(1,0)=12$), as guaranteed by strong duality.

2.6.3 Generating the optimality cut

We substitute the optimal dual solution into the general dual objective:

$thetageq 6u_1-6x_1u_2-4x_2u_3.$

Using ($u_1=2$), ($u_2=0$), and ($u_3=1$), we obtain:

$thetageq 6(2)-6x_1(0)-4x_2(1).$

Therefore, the optimality cut is:

$thetageq 12-4x_2.$

At the current solution ($bar{x}=(1,0)$), this cut requires, $thetageq 12.$

The master can no longer activate only the first resource while pretending that the operational cost ($theta$) is zero. The cut also tells the master that activating the second resource can reduce the operational cost because that resource supplies units more cheaply.

2.6.4 The final master iteration

The master problem now contains both the feasibility cut and the optimality cut:

$min quad 5x_1+x_2+theta$

subject to

$x_1+x_2geq 1,$

$6x_1+4x_2geq 6,$

$thetageq 12-4x_2,$

$thetageq 0,$

$x_1,x_2in{0,1}.$

If only the first resource is activated, the master objective is, $5(1)+1(0)+12=17.$ If both resources are activated, the optimality cut requires, $thetageq 12-4(1)=8.$ The corresponding master objective is:

$5(1)+1(1)+8=14.$

The new master solution is therefore:

$bar{x}_1=1, qquad bar{x}_2=1, qquad theta=8.$

The new lower bound is, $LB=14.$

2.6.5 Evaluating the final solution

For ($bar{x}=(1,1)$), the operational subproblem becomes:

$Q(1,1)=min quad 2y_1+y_2$

subject to

$y_1+y_2geq 6,$

$y_1leq 6,$

$y_2leq 4,$

$y_1,y_2geq 0.$

Because the second resource has the lowest operating cost, the subproblem uses its entire capacity, $y_2=4.$ The remaining two units are supplied by the first resource, $y_1=2.$ The operational cost is:

$Q(1,1)=2(2)+4=8.$

Adding the activation costs ($5x_1 + x_2$) gives:

$5(1)+1(1)+8=14.$

This feasible solution improves the upper bound, $UB=14.$ We now have:

$LB=14, qquad UB=14.$

Because the lower and upper bounds coincide, the algorithm has converged. The optimal solution is:

$x_1^*=1, qquad x_2^*=1, qquad y_1^*=2, qquad y_2^*=4,$

with optimal objective value $z^*=14.$

The toy example illustrates the two distinct forms of feedback used by Benders decomposition. The first master decision produced no feasible operational solution, so Farkas’ lemma generated a feasibility cut that restricted where the master could search. The next decision was feasible, but its operational cost was underestimated, so the dual subproblem generated an optimality cut. Together, these cuts guided the master toward a solution that was both feasible and optimal.

3. Applying feasibility cuts to capacitated facility location

The toy example showed how insufficient resource capacity can make a Benders subproblem infeasible. The capacitated facility location problem follows the same logic on a larger scale. The master decides which facilities to open, while the subproblem determines how customer demand should be supplied by the selected facilities.

Unlike the uncapacitated problem studied in Part I, opening at least one facility does not guarantee feasibility. The opened facilities may provide insufficient total capacity, or their locations may prevent them from serving particular customers. The algorithm must therefore generate feasibility cuts whenever the selected facilities cannot support a valid transportation plan.

3.1 The capacitated facility location problem

Let ($I$) denote the set of candidate facilities and ($J$) the set of customers. We use ($Asubseteq Itimes J$) to represent the permitted connections between facilities and customers.

The parameters are:

• ($f_i$): fixed cost of opening facility ($i$);
• ($K_i$): capacity of facility ($i$);
• ($d_j$): demand of customer ($j$);
• ($c_{ij}$): cost of supplying one unit from facility ($i$) to customer ($j$).

The binary variable ($x_i$) indicates whether facility ($i$) is opened. The continuous variable ($y_{ij}$) represents the quantity supplied from facility ($i$) to customer ($j$). The complete capacitated facility location model is:

$min quad sum_{iin I} f_i x_i + sum_{(i,j)in A} c_{ij}y_{ij}$

subject to

$sum_{i:(i,j)in A} y_{ij}geq d_j qquad forall jin J,$

$sum_{j:(i,j)in A} y_{ij}leq K_i x_i qquad forall iin I,$

$x_iin{0,1} qquad forall iin I,$

$y_{ij}geq 0 qquad forall (i,j)in A.$

The first family of constraints requires every customer’s demand to be satisfied. The second limits the total quantity supplied by each facility to its available capacity and prevents closed facilities from serving customers.

3.2 The Benders master problem

The master problem contains the facility-opening variables and the variable ($theta$), which estimates the transportation cost. The initial master problem is:

$min quad sum_{iin I} f_i x_i+theta$

subject to

$x_iin{0,1} qquad forall iin I,$

$thetageq 0,$

together with the feasibility and optimality cuts generated during the algorithm. At the beginning, the master knows very little about the operational consequences of opening a particular combination of facilities. It may therefore select an inexpensive set of facilities that cannot satisfy all customer demand.

3.3 The operational subproblem

For a fixed master solution ($bar{x}$), the transportation subproblem is:

$Q(bar{x}) = min_y quad sum_{(i,j)in A} c_{ij}y_{ij}$

subject to

$sum_{i:(i,j)in A} y_{ij}geq d_j qquad forall jin J,$

$sum_{j:(i,j)in A} y_{ij}leq K_ibar{x}_i qquad forall iin I,$

$y_{ij}geq 0 qquad forall (i,j)in A.$

Two outcomes are possible. If the opened facilities can satisfy every customer, the subproblem is feasible and returns the minimum transportation cost. If the available capacity or permitted connections are insufficient, the subproblem is infeasible.

3.4 Optimality cuts

To derive the optimality cuts, rewrite the capacity constraints as:

$-sum_{j:(i,j)in A} y_{ij} geq -K_ibar{x}_i qquad forall iin I.$

Let ($alpha_jgeq0$) be the dual variable associated with the demand constraint of customer ($j$), and let ($beta_igeq0$) be associated with the transformed capacity constraint of facility ($i$).

The dual subproblem is:

$max_{alpha,beta} quad sum_{jin J} d_jalpha_j – sum_{iin I} K_ibar{x}_ibeta_i$

subject to

$alpha_j-beta_ileq c_{ij} qquad forall (i,j)in A,$

$alpha_jgeq0 qquad forall jin J,$

$beta_igeq0 qquad forall iin I.$

When the subproblem is feasible, an optimal dual solution ($alpha^k,beta^k$) generates the following optimality cut:

$theta geq sum_{jin J} d_jalpha_j^k – sum_{iin I} K_ibeta_i^k x_i.$

This cut improves the master’s approximation of the transportation cost associated with different facility-opening decisions.

3.5 Feasibility cuts

If the transportation subproblem is infeasible, we cannot obtain an operational cost or update the upper bound. Instead, Farkas’ lemma provides multipliers that certify why the selected facilities cannot satisfy the operational requirements.

Let ($r_jgeq0$) be the Farkas multiplier associated with the demand constraint of customer ($j$), and let ($s_igeq0$) be associated with the transformed capacity constraint of facility ($i$).

A valid Farkas certificate satisfies:

$r_j-s_ileq0 qquad forall (i,j)in A,$

$r_jgeq0 qquad forall jin J,$

$s_igeq0 qquad forall iin I,$

and, for the current master solution ($bar{x}$),

$sum_{jin J} d_jr_j – sum_{iin I} K_ibar{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:

$sum_{jin J} d_jr_j – sum_{iin I} K_i x_i s_i leq0.$

Equivalently, the Benders feasibility cut is:

$sum_{iin I} K_i s_i x_i geq sum_{jin 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:

$sum_{iin I} K_i x_i geq sum_{jin 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 ($f_i$), ($K_i$), and ($d_j$), while transportation_cost stores ($c_{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 ($j$), while customers_for_facility[i] identifies the customers reachable from facility ($i$). 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 ($x_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 ($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 ($bar{x}$), we search for nonnegative Farkas multipliers ($r_j$) and ($s_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] <= 0
    )

    m.Normalization = pyo.Constraint(
        expr=sum(m.r[j] for j in customers)
           + sum(m.s[i] for i in facilities) <= 1
    )

    m.CertificateViolation = pyo.Objective(
        expr=sum(demand[j] * m.r[j] for j in customers)
           - sum(capacity[i] * x_bar[i] * m.s[i] for i in facilities),
        sense=pyo.maximize
    )

    return m


farkas = build_farkas_subproblem(x_solution)
results = solver.solve(farkas)

if results.solver.termination_condition != TerminationCondition.optimal:
    raise RuntimeError("Farkas problem was not solved to optimality.")

certificate_value = pyo.value(farkas.CertificateViolation)

r_solution = {j: pyo.value(farkas.r[j]) for j in customers}
s_solution = {i: pyo.value(farkas.s[i]) for i in facilities}

print(f"Farkas objective: {certificate_value:.6f}")

for j, value in r_solution.items():
    if value > 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 ($r_j-s_ileq0$) 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 ($bar{x}$).

For the initial master solution, the optimal objective is positive, confirming that no feasible transportation plan exists. The corresponding values of ($r_j$) and ($s_i$) form the Farkas certificate. We then replace the fixed values ($bar{x}_i$) by the master variables ($x_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 ($alpha_j$) are associated with the customer-demand constraints, while ($beta_i$) correspond to the transformed facility-capacity constraints.

For a fixed master solution ($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] <= transportation_cost[i, j]
    )

    m.TransportationCost = pyo.Objective(
        expr=sum(demand[j] * m.alpha[j] for j in customers)
           - sum(capacity[i] * x_bar[i] * m.beta[i] for i in facilities),
        sense=pyo.maximize
    )

    return m

The dual feasibility constraints reproduce the relationship ($alpha_j-beta_ileq c_{ij}$) for every permitted transportation arc. The objective combines the demand contribution with the capacity made available by the current master solution ($bar{x}$).

Unlike the Farkas problem, this model is used only after feasibility has been established. Its optimal objective value gives the transportation cost ($Q(bar{x})$), while the values of ($alpha_j$) and ($beta_i$) provide the coefficients required to construct an optimality cut and improve the master’s approximation represented by ($theta$).

4.5 Putting everything together: the complete Benders loop

We can now combine the master problem, the normalized Farkas problem, and the transportation dual into the complete Benders algorithm. Since the first feasibility cut has already been added in the previous subsection, the loop starts by resolving the updated master and then repeatedly determines which type of feedback is required.

At each iteration, the master provides a facility configuration and a lower bound. We first solve the Farkas problem. If its objective is positive, the transportation problem is infeasible and we add a feasibility cut. Otherwise, we solve the transportation dual, obtain a feasible total cost and update the upper bound. If the bounds have not yet converged, the corresponding dual multipliers generate a new optimality cut.

maximum_iterations = 100
tolerance = 1e-6

lower_bound = -float("inf")
upper_bound = float("inf")
best_solution = None
iteration_results = []

for iteration in range(1, maximum_iterations + 1):

    # 1. Solve master
    solver.solve(master)
    x_bar = {
        i: int(round(pyo.value(master.x[i])))
        for i in facilities
    }

    theta = pyo.value(master.theta)
    fixed = sum(fixed_cost[i] * x_bar[i] for i in facilities)
    lower_bound = pyo.value(master.TotalCost)

    # 2. Check feasibility
    farkas = build_farkas_subproblem(x_bar)
    solver.solve(farkas)
    certificate = pyo.value(farkas.CertificateViolation)

    if certificate > 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 < upper_bound:
            upper_bound = total_cost
            best_solution = dict(x_bar)
            best_fixed_cost = fixed
            best_transportation_cost = transportation
            best_total_cost = total_cost

        gap = upper_bound - lower_bound

        # 4. Add optimality cut if needed
        if gap <= tolerance:
            cut_type = "None (converged)"
        else:
            alpha = {j: pyo.value(dual.alpha[j]) for j in customers}
            beta = {i: pyo.value(dual.beta[i]) for i in facilities}

            master.OptimalityCuts.add(
                master.theta >=
                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 <= tolerance and gap <= tolerance:
        break

print("nFinal Benders solution")
print(f"Open facilities: {[i for i in facilities if best_solution[i]]}")
print(f"Fixed cost: {best_fixed_cost:,.2f}")
print(f"Transportation cost: {best_transportation_cost:,.2f}")
print(f"Total cost: {best_total_cost:,.2f}")

Once the code finishes running, we can see it converges on iteration 10, obtaining the minimum cost.

----------------------------------------------------------------------
Iteration: 1
Open facilities: ['F2', 'F4']
Fixed cost: 4,300.00
Theta: 0.00
Farkas objective: 1.500000
Cut added: Feasibility
Lower bound: 4,300.00
Actual transportation cost: infeasible
Upper bound: not available yet
Relative gap: not available yet
----------------------------------------------------------------------
Iteration: 2
Open facilities: ['F4', 'F5']
Fixed cost: 4,400.00
Theta: 0.00
Farkas objective: 1.250000
Cut added: Feasibility
Lower bound: 4,400.00
Actual transportation cost: infeasible
Upper bound: not available yet
Relative gap: not available yet
----------------------------------------------------------------------
Iteration: 3
Open facilities: ['F2', 'F5']
Fixed cost: 4,500.00
Theta: 0.00
Farkas objective: 2.000000
Cut added: Feasibility
Lower bound: 4,500.00
Actual transportation cost: infeasible
Upper bound: not available yet
Relative gap: not available yet
----------------------------------------------------------------------
Iteration: 4
Open facilities: ['F2', 'F4', 'F5']
Fixed cost: 6,600.00
Theta: 0.00
Farkas objective: -0.000000
Cut added: Optimality
Lower bound: 6,600.00
Actual transportation cost: 8,946.00
Upper bound: 15,546.00
Relative gap: 57.545349%
----------------------------------------------------------------------
Iteration: 5
Open facilities: ['F1', 'F3', 'F4', 'F5']
Fixed cost: 9,300.00
Theta: 0.00
Farkas objective: -0.000000
Cut added: Optimality
Lower bound: 9,300.00
Actual transportation cost: 5,457.00
Upper bound: 14,757.00
Relative gap: 36.979061%
----------------------------------------------------------------------
Iteration: 6
Open facilities: ['F2', 'F3', 'F4', 'F5']
Fixed cost: 9,100.00
Theta: 2,194.00
Farkas objective: -0.000000
Cut added: Optimality
Lower bound: 11,294.00
Actual transportation cost: 5,853.00
Upper bound: 14,757.00
Relative gap: 23.466829%
----------------------------------------------------------------------
Iteration: 7
Open facilities: ['F1', 'F2', 'F4', 'F5']
Fixed cost: 9,000.00
Theta: 2,826.00
Farkas objective: -0.000000
Cut added: Optimality
Lower bound: 11,826.00
Actual transportation cost: 6,497.00
Upper bound: 14,757.00
Relative gap: 19.861761%
----------------------------------------------------------------------
Iteration: 8
Open facilities: ['F3', 'F4', 'F5']
Fixed cost: 6,900.00
Theta: 5,853.00
Farkas objective: 0.666667
Cut added: Feasibility
Lower bound: 12,753.00
Actual transportation cost: infeasible
Upper bound: 14,757.00
Relative gap: inf%
----------------------------------------------------------------------
Iteration: 9
Open facilities: ['F1', 'F2', 'F3', 'F4', 'F5']
Fixed cost: 11,500.00
Theta: 2,055.00
Farkas objective: -0.000000
Cut added: Optimality
Lower bound: 13,555.00
Actual transportation cost: 4,208.00
Upper bound: 14,757.00
Relative gap: 8.145287%
----------------------------------------------------------------------
Iteration: 10
Open facilities: ['F1', 'F3', 'F4', 'F5']
Fixed cost: 9,300.00
Theta: 5,457.00
Farkas objective: -0.000000
Cut added: None (converged)
Lower bound: 14,757.00
Actual transportation cost: 5,457.00
Upper bound: 14,757.00
Relative gap: 0.000000%

Benders decomposition converged.

======================================================================
Final Benders solution
======================================================================
Facility F1: Open (capacity = 18)
Facility F2: Closed (capacity = 20)
Facility F3: Open (capacity = 16)
Facility F4: Open (capacity = 18)
Facility F5: Open (capacity = 17)

Fixed opening cost: 9,300.00
Transportation cost: 5,457.00
Total cost: 14,757.00
Iterations after the first feasibility cut: 10

The central difference with Part I is the conditional step immediately after solving the master. A positive Farkas objective means that the current facility configuration cannot support a feasible transportation plan, so no upper bound can be computed and the algorithm adds a feasibility cut instead. Once the Farkas objective falls to zero, the transportation problem is feasible and the algorithm switches to the familiar optimality-cut mechanism.

The master objective provides the lower bound throughout the procedure, while every feasible facility configuration provides a candidate upper bound. The loop terminates only when a feasible subproblem is obtained and these two bounds coincide within the prescribed tolerance. We also store the main information from every iteration so that we can inspect how the feasibility and optimality cuts progressively guide the master toward the optimum.

4.6 Recovering and visualizing the final solution

The Benders algorithm gives us the optimal facility-opening decisions and transportation cost, but the dual subproblem does not directly provide the actual customer flows ($y_{ij}$). We therefore solve the original transportation subproblem one final time with the optimal facility configuration fixed. This recovers the flows used in the final solution, which we can then visualize together with the open and closed facilities.

As an additional validation, the complete notebook also solves the original monolithic capacitated facility location model and confirms that its objective value matches the Benders solution.

# Recover the primal transportation solution
transport = pyo.ConcreteModel()
transport.y = pyo.Var(arcs, domain=pyo.NonNegativeReals)

transport.Demand = pyo.Constraint(
    customers,
    rule=lambda m, j:
        sum(m.y[i, j] for i in facilities_for_customer[j]) >= demand[j]
)

transport.Capacity = pyo.Constraint(
    facilities,
    rule=lambda m, i:
        sum(m.y[i, j] for j in customers_for_facility[i])
        <= capacity[i] * best_solution[i]
)

transport.Cost = pyo.Objective(
    expr=sum(
        transportation_cost[i, j] * transport.y[i, j]
        for i, j in arcs
    ),
    sense=pyo.minimize
)

solver.solve(transport)

flows = {
    (i, j): pyo.value(transport.y[i, j])
    for i, j in arcs
    if pyo.value(transport.y[i, j]) > 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 ($y_{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.

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.