Benders Decomposition in CPLEX

Author's note: This post is getting a bit long in the tooth. For users with CPLEX 12.7 or later, you might want to read this more recent post (after reading this one).

I finally got around to constructing a reasonably small example of Benders decomposition for solving a mixed-integer linear program (MIP) in CPLEX. What I am about to describe works in CPLEX 12.4 (possibly 12.2 or 12.3) and hopefully will work in the future (since I have no ambition to revise the post). The significance of the CPLEX version is twofold: IBM refactored the API for cuts, emphasizing the distinction between "user cuts" (cuts that tighten bounds without eliminating integer-feasible solutions) and "lazy constraints" (cuts that eliminate feasible solutions); and they decided that lazy constraint callbacks would be called whenever CPLEX found an incumbent, regardless of when and how it found the incumbent (node LP solution, application of heuristics, lucky guess ...), which eliminated the need for incumbent and branch callbacks. If you are still using an earlier version of CPLEX, I strongly recommend you upgrade. Among various reasons, you will need a bunch more code to pull off Benders correctly with an old version.

My test problem is simple: a transportation problem with fixed charges for the use of supply sources. (In the code, I refer to sources as "warehouses" and sinks as "customers"). The fixed-charge transportation problem is well know and rather easy to solve -- too easy to warrant Benders decomposition, in fact. I chose it because it is sufficiently straightforward that it should not interfere with understanding how the code works. (Subsequent to writing the first draft of this post, I discovered a nice document by Erwin Kalvelagen that demonstrates Benders decomposition in GAMS. He too used a fixed-charge transportation problem, but in his version there is a fixed charge for each arc used, rather than for each source.)

The mathematical model is as follows: \begin{gather*} \textrm{minimize }c'x+f'y\\ \textrm{s.t. }\sum_{i}x_{ij}\ge d_{j}\quad\forall j\\ \sum_{j}x_{ij}\le s_{i}y_{i}\quad\forall i\\ x\ge0\\ y_{i}\in\{0,1\}\quad\forall i. \end{gather*} Continuous variable $x_{ij}$ represents the flow from source $i$ to sink $j$; binary variable $y_{i}$ signals whether source $i$ is used (1) or not (0). Parameters $c\ge0$ and $f\ge0$ are unit flow costs and facility use charges respectively, while $s$ and $d$ are supplies (source capacities) and demands.

My Benders decomposition puts the facility decisions in a master problem (MIP) and the flow decisions in a subproblem (linear program). The subproblem is just \begin{gather*} \textrm{minimize }z=c'x\\ \textrm{s.t. }\sum_{i}x_{ij}\ge d_{j}\quad\forall j\\ \sum_{j}x_{ij}\le s_{i}y_{i}\quad\forall i\\ x\ge0 \end{gather*} in which the $y_{i}$ are treated as constants. The master problem is \begin{gather*} \textrm{minimize }f'y+\hat{z}\\ \textrm{s.t. }a^{(k)\prime}y-\hat{z}\le b^{(k)}\quad\forall k\in\mathcal{O}\\ a^{(k)\prime}y\le b^{(k)}\quad\forall k\in\mathcal{F}\\ y_{i}\in\{0,1\}\quad\forall i\\ \hat{z}\ge0 \end{gather*} where $\hat{z}$ is a surrogate for the flow cost component ($z$) of the original objective function (i.e., the value of the subproblem objective), $a^{(k)}$ and $b^{(k)}$ are coefficient vector and constant term for a cut generated during the solution of the problem, $\mathcal{O}$ is the index set of "optimality" cuts (cuts that attempt to correct the master underestimating the true flow costs for a given set of decisions $y$, i.e., $\hat{z}<z$) and $\mathcal{F}$ is the index set of "feasibility" cuts (cuts that eliminate use decisions $y$ that provide too little capacity to satisfy demand). Initially $\mathcal{O}=\emptyset=\mathcal{F}$, and the first master problem solution is both trivial ($y=0$, $\hat{z}=0$) and infeasible in the original model.

To generate cuts, we will need to work indirectly with the dual of the subproblem: \begin{gather*} \textrm{maximize }z=d'\lambda+\sum_{i}(s_{i}\mu_{i})y_{i}\\ \textrm{s.t. }\lambda_{j}+\mu_{i}\le c_{ij}\quad\forall i,j\\ \lambda\ge0\\ \mu\le0 \end{gather*} where $\lambda$ is the vector of shadow prices of the demand constraints and $\mu$ the vector of shadow prices of the supply constraints. Given the assumption that all flow costs are nonnegative ($c\ge0$), the dual is obviously feasible ($\lambda=0=\mu$ works). If, for a given choice of $y$, the LP subproblem has an optimal solution, then so will the dual (with equal objective value $z$). Should the master problem generate a solution $(y,\hat{z})$ for which the subproblem is feasible with optimal value $z=\hat{z}$, we accept that solution as the new incumbent. If, instead, the subproblem is feasible but $\hat{z}<z$, then we add a new optimality cut $\hat{z}\ge d'\lambda+\sum_{i}(s_{i}\mu_{i})y_{i}$ where $(\lambda,\mu)$ is the dual solution to the subproblem (so that $a_{i}^{(k)}=s_{i}\mu_{i},b^{(k)}=-d'\lambda$).

Should the master problem generate a solution $(y,\hat{z})$ for which the subproblem is infeasible, then the dual will be unbounded, and there will be a dual solution $(\lambda,\mu)$ such that $t(\lambda,\mu)$ is feasible in the dual for all $t>0$ and $d'\lambda+\sum_{i}(s_{i}\mu_{i})y_{i}>0$. Note that the master problem variables ($y$) have no role in defining the dual problem's feasible region; in the dual, they occur only in the objective. So the dual ray we just obtained will remain a recession direction for the dual regardless of $y$, and to get a feasible subproblem we need to modify $y$ so that $d'\lambda+\sum_{i}(s_{i}\mu_{i})y_{i}\le0$ (i.e., so that the ray now points in an undesirable direction). This gives us our feasibility cut ($a_{i}^{(k)}=s_{i}\mu_{i},b^{(k)}=-d'\lambda$). Observe that optimality and feasibility cuts have the same constant term and coefficient vector; the only difference is the inclusion of $\hat{z}$ in optimality cuts and its exclusion in feasibility cuts.

On to the code. Let me start with my usual disclaimer: I am not a professional coder, and I make no claim that my code is particularly efficient, nor that I took the fastest/most straightforward route to the solution. As usual, the code is in Java. For the most part, conversion to other APIs that expose the callback functionality should be straightforward, but there are a few wrinkles to note.

• If you code the comparisons $\hat{z}<z$ and $z=\hat{z}$ literally, the gods of computing will wait until your back is turned and then bite you in your derrière. There will be rounding errors in both $z$ and $\hat{z}$, so asking for exact equality is inviting disaster. Insert some rounding tolerance into your comparisons. In my code, look for a constant named FUZZ.
• When the LP subproblem is feasible, we can get the dual values directly with a function named getDuals. When the LP is infeasible, it is tempting to use the getRay function to get the dual ray, but to do that you would have to formulate the dual itself, solve it directly, and then invoke getRay on the dual. It is easier to ask for a Farkas certificate, which provides the ray without further ado.
• I found it necessary to turn off presolving of the subproblem. If the presolver detects that the subproblem is infeasible, it reports that fact with no option to get a Farkas certificate. So I sacrifice whatever speed advantage the presolver confers in order to have access to the Farkas certificate when the subproblem is infeasible.
• Update 11/30/16: I recently tripped over another thing that I should have realized but somehow overlooked. In order to access the Farkas certificate, the LP subproblem must be solved by one of the simplex solvers (primal, dual or network). I'd been leaving the choice up to CPLEX, and on one instance it used the barrier method, decided the LP was infeasible, and then threw an exception when the Farkas certificate was requested. There's also a glitch in the CPLEX documentation for the dualFarkas method. It turns out that dualFarkas only works if the LP was solved by the dual simplex algorithm. So the key is to set the RootAlg parameter for the LP to 2 (dual simplex).
• Forming a cut (either optimality or feasibility) involves multiplying the right-hand side of a subproblem constraint, expressed as a linear function of the master variables, by the corresponding dual value. In this particular example, the right-hand side always has a single term, either a constant ($d_{j}$) in a demand constraint or a multiple of a single master variable ($s_{i}y_{i}$) in a supply constraint.
• More generally, though, each right-hand side could contain multiple terms. I deliberately over-engineered my code a bit so that it will generalize to other problems. What I do in the code is form and store an instance of IloNumExpr for each right-hand side, at the time that I am creating the subproblem. Inside the lazy constraint callback, I multiply those expressions by the corresponding dual values and sum to get the necessary expression ($a^{(k)\prime}y-b^{(k)}$) for the new cut. Somewhat counterintuitively, and for rather arcane reasons, this cannot be done with IloLinearNumExpr, even though the right-hand sides are in fact linear.
• Since any master variable can occur in multiple subproblem right-hand sides, the method I just described for assembling the cut expression will typically result in multiple terms involving the same variable ($\dots+\alpha_{3}y_{7}+\dots+\alpha_{12}y_{7}+\dots+\alpha_{23}y_{7}+\dots$). Not to worry: when you actually add the cut, CPLEX will combine like terms ($\dots+(\alpha_{3}+\alpha_{12}+\alpha_{23})y_{7}+\dots$) for you.
• Interpreting the Farkas certificate is a bit tricky, since CPLEX returns the dual values in an arbitrary order (see my earlier post "Infeasible LPs and Farkas Certificates"). To deal with this, I store the subproblem constraints in two arrays (one for supply constraints, one for demand constraints), and then use a HashMap to map each constraint to the corresponding right-hand side expression.
• The C++ Standard Template Library has some analogous container class you can use (don't ask me which; I'm blissfully C++-free and mean to stay that way). With Python, you would create a dict. With other APIs, you are on your own.
• The reason I use two arrays, rather than storing all constraints in a single array, is that it facilitates updating right-hand sides in the subproblem each time the callback is invoked. I need to recompute the RHS of each supply constraint ($s_{i}y_{i}$, using the new value of $y_{i}$), but not the RHS of each demand constraint ($d_{j}$). Having the supply and demand constraints separate makes this a bit easier, but in a more general application (where most or all constraints require an update), you may want to use more than two arrays, or just a single array.
• Once the master problem has been solved, you have the optimal set of warehouse use decisions ($y$). There's a fairly good chance you will also want to know the optimal flows ($x$). The subproblem may well have been solved more recently than when the callback accepted what proved to be the optimal $y$; after that node, CPLEX may have generated some incumbents that your callback tested and rejected. So the most recent subproblem solution may not correspond to the optimal $y$. One resolution is to solve the subproblem one last time, using the optimal $y$. The other (which I use in the code) is to store the subproblem solution any time an incumbent is accepted, overwriting any previously stored solution. When the master is solved, the stored subproblem solution will be the correct flows.

That's it. To avoid inducing comas among readers, I will not go through the code line by line here. Instead, I will just provide this link for anyone who wants to download it. [UPDATE: A newer version of CPLEX broke the code. I have created a repository for the updated code; please see this more recent post for details.]