AMPL Development Tutorial 6/6 – Implementing Benders Decomposition with ampls
Description: This concluding notebook in our six-part series delves into enhancing the efficiency of our decomposition algorithm by utilizing AMPL Solver Libraries (ampls).
We showcase the use of ampls to implement our Benders decomposition approach for solving our stochastic facility location problem.
AMPL’s intuitive syntax allows for straightforward allocation of variables and constraints to master and subproblems.
However, each iteration typically involves solving the problems from scratch.
To enhance efficiency, we introduce a method using ampls to export the master problem to its solver representation.
This allows us to add cuts directly to the solver’s matrix, enabling solvers to efficiently utilize most of the previously computed solution.
This technique significantly boosts the overall process efficiency allowing to readily embed algorithmic development into production applications.
Tags: amplpy, ampl, ampls, mip, stochastic, facility location, benders
Notebook author: Christian Valente <ccv@ampl.com>, Gyorgy Matyasfalvi <gyorgy@ampl.com>
References:
AMPL a Modeling Language for Mathematical Programming – Robert Fourer et al.
Intro to Linear Optimization – Dimitris Bertsimas et al.
SCIP Optimization Suite example projects – Stephen J. Maher (ZIB)
Two-stage stochastic programming problem formulation
Consider a decision-making process occurring over two sequential stages.
Initially, a decision vector $x$ is selected.
Following this, new information becomes available, prompting the selection of a second-stage decision vector $y$.
We introduce $S$ possible scenarios for the newly acquired information, with the actual scenario being revealed only after the selection of $x$.
Each scenario is indexed by $s$, and its occurrence probability is denoted by $\alpha^s$, which we assume to be positive. Given the sequential nature of decision-making, the choice of $y$ can be scenario-dependent, denoted as $y^s$ to reflect this.
The decision vectors $x$ and $y^s$ incur costs represented by vectors $f$ and $q$, respectively.
The first-stage decisions are constrained by:
\begin{equation}
\begin{array}{crlllllll}
& A x & = & {b} \
& x & \geq & 0.
\end{array}
\tag{1}
\end{equation}
Furthermore, the combination of first and second-stage decisions must satisfy the scenario-specific constraints:
\begin{equation}
\begin{array}{crcllllll}
& {B}^s x + {D}^s y^s & = & {d}^s \
& y^s & \geq & 0,
\end{array}
\tag{2}
\end{equation}
for every scenario $s$.
The uncertainty in the model is encapsulated by $\xi := (q, B, D, d)$, where a constraint is deemed part of the second stage if it involves any nonzero coefficients of variables ${y}$.
The matrices ${B}^s$ are referred to as technology matrices, which determine the influence of the decision at the first stage on the second stage.
The objective is to determine the choices of $x$ and $y^s$ for all $s=1, \ldots, S$ that minimize the expected total cost:
\begin{equation}
\begin{array}{cllllllll}
& f^T x & + & \alpha_1 {(q^1)^T} y^1 & + & \alpha^2 {(q^2)^T} y^2 & \ldots & + & \alpha^S {(q^S)^T} y^S
\end{array}
\tag{3}
\end{equation}
This formulation leads to the following two-stage stochastic optimization problem:
\begin{equation}
\begin{array}{crcrcrccrcl}
\min & f^T x & + & \alpha^1 {(q^1)^T} y_1 & + & \alpha^2 {(q^1)^T} y^2 & \ldots & + & \alpha^s {(q^s)^T} y^S & & \
\textrm{s.t.} & A x & & & & & & & & &= {b} \
& {B^1} x & + & {D^1} y^1 & & & & & & &= {d^1} \
& {B^2} x & + & & & {D^2} y^2 & & & & &= {d^2} \
& \vdots & & & & & \ddots & & & &\vdots \
& {B^S} x & + & & & & & & {D^S} y^S & &= {d^S} \
& x \geq 0, & & y^1 \geq 0, & \ldots & & & & y^S \geq 0 & &
\end{array}
\tag{4}
\end{equation}
we will refer to the above problem as the master problem.
Notice that even if the number of possible scenarios $S$ is moderate, this formulation can be a large linear programming problem.
However, a decomposition method can help.
Benders decomposition for two-stage problems
Given a vector $\hat{x}$, using the dual simplex method solve the subproblem (for all $s \in S$):
\begin{equation}
\begin{array}{crl}
\min & (q^s)^T y^s & \
\textrm{s.t.} & D^s y^s & = d^s - B^s \hat{x} \
& y^s & \geq 0
\end{array}
\tag{5}
\end{equation}
The dual of (5) takes the following form:
\begin{equation}
\begin{array}{crl}
\max & (p^s)^T (d^s - B^s \hat{x}) & \
\textrm{s.t.} & (p^s)^T D^s & <= (q^s)^T
\end{array}
\tag{6}
\end{equation}
Let
\begin{equation}
\begin{array}{crrl}
& P ;; = & { p & | ;; p^T D \leq q^T } & \
\end{array}
\tag{7}
\end{equation}
We assume that $P$ is nonempty and has at least one extreme point (so we can apply duality theory). Let $p_i$, $i = 1, \ldots, I$, be the extreme points, and let $w^j$, $j = 1, \ldots, J$, be a complete set of extreme rays of $P$.
Define $z^s(\hat{x})^*$ as the optimal cost obtained from problem (5) described above. Returning to the optimization of $x$, we encounter the following problem:
\begin{equation}
\begin{array}{crcrcrccrcl}
\min & f^T x & + & \sum_{s=1}^S \alpha^s z^s(x) & \
\textrm{s.t.} & A x & & & = {b} \
& x & \geq 0 & &
\end{array}
\tag{8}
\end{equation}
Under our assumption that the set $P$ is nonempty, either the dual problem (6) has an optimal solution and $z^s(\hat{x})$ is finite, or the optimal dual cost is infinite (indicating that the
primal (5) is infeasible).
In particular, $z^s(x) \leq \infty \leftrightarrow (w^j)^T (d^s - B^s x) \leq 0, ; \forall j$ and whenever $z^s(x)$ is finite the optimum of (6) must be attained at an extreme point of the set $P$ meaning $z^s(x)$ is the smallest number $z^s$ such that $(p^i)^T (d^s - B^s x) \leq z^s, ; \forall i$, leading us to the following master-problem formulation:
\begin{equation}
\begin{array}{crcrcrccrcl}
\min & f^T x & + & \sum_{s=1}^S \alpha^s z^s & & \
\textrm{s.t.} & A x & & & = b & \
& (p^i)^T (d^s - B^s x) & & & \leq z^s & \forall ; i, s \
& (w^j)^T (d^s - B^s x) & & & \leq 0 & \forall ; j, s \
& x & \geq 0 & & &
\end{array}
\tag{9}
\end{equation}
In our algorithm we will only generate constraints that we find to be violated by the current solution.
The algorithm
Initialization step: initialize $(\hat{x}, \hat{z}^1, \ldots, \hat{z}^S)$ to zero.
Sub-problem step: using the dual-simplex method solve the sub-problems individually for every $s \in S$:
\begin{equation}
\begin{array}{crl}
\min & (q^s)^T y^s & \
\textrm{s.t.} & D^s y^s & = d^s - B^s \hat{x} \
& y^s & \geq 0
\end{array}
\tag{5}
\end{equation}
Check optimality of master-problem: If for every $s$, the corresponding sub-problem is feasible and $(q^s)^T \hat{y}^s \leq \hat{z}^s$ then all constraints are satisfied and we have an optimal solution to the master problem.
Add optimality cut: If the sub-problem corresponding to some $s$ has an optimal solution such that $(q^s)^T \hat{y}^s > \hat{z}^s$ the following cut is added to the master problem:
\begin{equation}
(p^s)^T (d^s - B^s x) \leq z^s,
\tag{10}
\end{equation}
where $p^s$ is an optimal basic feasible solution to the dual of the sub-problem.
Add feasibility cut: If the sub-problem corresponding to some $s$ is infeasible, its dual has infinite cost and the following cut is added to the master problem:
\begin{equation}
(w^s)^T (d^s - B^s x) \leq 0,
\tag{11}
\end{equation}
where $w^s$ is a positive cost extreme ray.
Stochastic facility location a concrete example
Given:
Task:
Variables
$x_i \in {0, 1} \quad \forall i \in I$
$y_{ij}^s \geq 0 \quad \forall i \in I, \forall j \in J, \forall s \in S$
Parameters:
$\alpha^s$: the probability of scenario $s$.
$f_i$: the fixed cost for opening facility $i$,
$q_{ij}$: the cost of servicing customer $j$ from facility $i$,
$\lambda_j^s$: the demand of customer $j$ in scenario $s$,
$k_i:$ the capacity of facility $i$.
The sub-problem
Given the above formulation for the extensive form (12) we can express our scenario specific sub-problems as follows
$
\begin{equation}
\begin{array}{lcrllll}
z_s(\hat{x}) & = & \min \quad & \sum_{i \in I} \sum_{j \in J} q_{ij} y_{ij}^s & & & \
& & & & & & \
& & \textrm{s.t.} \quad & \sum_{i \in I} y_{ij}^s & \geq & \lambda_j^s & \forall j \in J \
& & & \sum_{j \in J} y_{ij}^s & \leq & k_i \hat{x_i} & \forall i \in I \
& & & & & & \
& & & y_{ij}^s \geq 0 & & & \forall ; i \in I, \forall j \in J
\end{array} \tag{13}
\end{equation}
$
The master-problem
The master-problem takes the following form
$
\begin{equation}
\begin{array}{rllll}
\min \quad & \sum_{i \in I} f_i x_i + \sum_{s \in S} \alpha^s z^s & & & \
& & & & \
\textrm{s.t.} \quad & \sum_{i \in I} k_i x_i & \geq & \max_{s \in S} \sum_{j \in J} \lambda_j^s & \
& \sum_{j \in J} p_j^s \lambda_j^s + \sum_{i \in I} p_i^s k_i x_i & \leq & z^s & \forall s \in S, \forall p \in P^s \
& \sum_{j \in J} w_j^s \lambda_j^s + \sum_{i \in I} w_i^s k_i x_i & \leq & 0 & \forall s \in S, \forall w \in P^s \
& & & & \
& x_i \in {0, 1} & & & \forall i \in I \
\end{array} \tag{12}
\end{equation}
$
Create amplpy instances of master and sub-problems
create_common(): this function initializes a base AMPL object. It incorporates the model’s sets and parameters that are common to both the master and the sub-problem.
create_sub_problem(): this function is dedicated to setting up the AMPL instance for the sub-problem. It begins by invoking create_common(). Following this, it employs the eval() method to dynamically load the additional model entities and parameters specifically required for the sub-problem.
create_master_problem(): similar to the sub-problem setup, this function also starts with create_common() to instantiate an AMPL object with the basic shared model elements. It then utilizes the eval() method to incorporate additional model entities and parameters that are exclusive to the master problem.
Data
We can use the same data file for both the master and the sub-problem.
Utility functions
We can declare a few python functions that will make the subsequent code more succint.
Firstly, to retain the ability of assigning data from AMPL data files but also allow assigning it directlry from python, we declare get_data_from_ampl that extracts the data from AMPL and creates the python data structures that will be needed when executing the algorithm. Note that in real life usage, the data will mostly come from Python itself.
When exporting the model from AMPL to ampls we will lose the semantics of the model, because ampls stores its representation at solver level. Intuitively, instead of having the AMPL structured model, made of sets, parameters, variables and constraints, we have the Jacobian matrix, a cost vector and a vector of right hand sides.
To be able to add cuts (done via ampls) from values we get by solving the subproblems (via AMPL), we need some helper functions and some maps.
Finally, we write a function that adds the Benders cuts using ampls; the optimality and feasibility cuts are very similar to each other:
$
\begin{equation}
\begin{array}{rllll}
& & & & \
\sum_{j \in J} p_j^s \lambda_j^s + \sum_{i \in I} p_i^s k_i x_i & \leq & z^s & \forall s \in S, \forall p \in P^s \
\sum_{j \in J} w_j^s \lambda_j^s + \sum_{i \in I} w_i^s k_i x_i & \leq & 0 & \forall s \in S, \forall w \in P^s \
& & & & \
\end{array}
\end{equation}
$
which, applied to this problem, correspond to the following (structured) AMPL:
s.t. optimality_cut {k in 1..nITER, s in SCENARIOS}:
sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s] +
sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i] <= sub_variable_cost[s];
s.t. feasibility_cut {k in 1..nITER, s in SCENARIOS}:
sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s] +
sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i] <= 0;
In ampls we lose the structure, so we need to pass the indices (in the Jacobian matrix) and the coefficients of each non-zero element of the cut.
Also, the (numeric) right hand side will have to be calculated explicitly.
Let’s start by reordering the formulation above: the RHS is actually sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s], as this is the only numerical data.
Let’s then move the $z^s$ (which in AMPL is sub_variable_cost[s]) to the left. What we will pass to ampls is akin to:
sub_variable_cost[s] - sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i]
>= sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s];
We therefore need the indices of the variables sub_variable_cost[s] in the master problem, along with the indices of all the facility_open[i]. We can get them from the maps we pass to the function: index_sub_variable_cost and index_facility_open respectively. Also note the negative coefficients for the product facility_price*facility_capacity*facility_open.
Benders decomposition
The main function is, at this point, easy to implement.
We create two instances of AMPL for the master and the sub problem respectively, clearly using the same data and set some options. At this stage, we can also get the data out from one of the AMPL instances to have it readily available in Python for the algorithm implementation. We also initialize some entities that we’ll use in the algorithm.
We use the function AMPL.to_ampls to export an instance of the master problem to the ampls library. Note that at this point, solver_master is a solver-level representation of the master problem, and we lose direct access to AMPL’s structured entities. This is why we also get the maps linking entities names to solver index.
Now the master loop, which implements Benders.
At every algorithm iteration, we:
Set the scenario index in the subproblem and solve it. We check the solution:
Infeasible -> We extract the dual rays and use them to construct a feasibility cut
Feasible but violating the optimality condition -> We extract the dual values and use them to construct an optimality cut
Feasible and not violating -> We mark this scenario as non-violating
If no scenarios had violations, we reached optimality and quit
Otherwise, we resolve the master problem (with the cuts added in step 1)
We use the solution to assign:
the parameter sub_facility_open in the subproblems
the dictionary sub_variable_cost used at the next iteration to check for violations
****** Iteration 1 *******
Scenario Low Objective=11621793.455
Added optimality cut: +1.000000*X[0]>= 11621793.455000;
Scenario Medium Objective=14779784.865
Added optimality cut: +308961.500000*X[3]+1.000000*X[1]>= 15088746.365000;
Scenario High Objective=21050711.200000003
Added optimality cut: +1548961.500000*X[3]+520000.000000*X[4]+1.000000*X[2]>= 23119672.700000;
SOLVING MASTER PROBLEM
****** Iteration 2 *******
Version identifier: 22.1.1.0 | 2022-11-28 | 9160aff4d
Scenario Low Objective=11621793.455
CPXPARAM_Simplex_Display 0
CPXPARAM_MIP_Display 0
CPXPARAM_Barrier_Display 0
Scenario Medium Objective=14966984.865
Added optimality cut: +1548961.500000*X[3]+520000.000000*X[4]+1.000000*X[1]>= 16515946.365000;
Scenario High Objective=21570711.200000003
SOLVING MASTER PROBLEM
Version identifier: 22.1.1.0 | 2022-11-28 | 9160aff4d
****** Iteration 3 *******
CPXPARAM_Simplex_Display 0
CPXPARAM_MIP_Display 0
CPXPARAM_Barrier_Display 0
Scenario Low Objective=11621793.455
Scenario Medium Objective=14779784.865
Scenario High Objective=21050711.200000003
Import solution back to AMPL
In the following part of code, we have the correct solution in the object solver_master. We could settle at this, but often it is convenient to have the solution back in AMPL.
We therefore use the method AMPL.import_solution to import the solution values back to the AMPL object.
For demonstrational purposes, we also use the line:
master.eval(solver_master.get_recorded_entities())
to show how to import back into AMPL all the entities that have been added via ampls (in our case the cuts) and we print them in AMPL.
Optimal solution found, cost: 16758018.596249998
CPLEX 22.1.1: optimal solution; objective 16758018.59625
0 simplex iterations
CPLEX 22.1.1: optimal solution; objective 16758018.59625
0 simplex iterations
subject to sufficient_production_capacity:
1550*facility_open['Baytown_TX'] + 650*facility_open['Beaumont_TX'] +
1750*facility_open['Baton_Rouge_LA'] >= 3223;
subject to ampls_con_1:
sub_variable_cost['Low'] >= 11621800;
subject to ampls_con_2:
sub_variable_cost['Medium'] + 308962*facility_open['Baytown_TX'] >=
15088700;
subject to ampls_con_3:
sub_variable_cost['High'] + 1548960*facility_open['Baytown_TX'] +
520000*facility_open['Beaumont_TX'] >= 23119700;
subject to ampls_con_4:
sub_variable_cost['Medium'] + 1548960*facility_open['Baytown_TX'] +
520000*facility_open['Beaumont_TX'] >= 16515900;
Check solution
The following is just to double check the result: we solve the extensive form of the problem and check that its result coincides to the one we get from our implementation of Benders.
CPLEX 22.1.1.0: optimal integer solution; objective 16758018.6
18 MIP simplex iterations
0 branch-and-bound nodes
Optimal solution EF form: 16758018.596250001
Optimal solution Benders: 16758018.596249994
EF and Benders objective values match!