AMPL Development Tutorial 3/6 – Benders Decomposition via AMPL scripting
Description: In this third installment of our six-part series, we continue our exploration by addressing the complexities introduced by the stochastic programming formulation presented in part two.
While this approach offers superior results compared to the deterministic model, it significantly increases the problem’s complexity.
In this section, we delve into the Benders decomposition technique, applying it to our two-stage stochastic problem using AMPL’s advanced features.
This method aims to provide a ‘robust’ solution that is adaptable to any future customer demand scenario, yet avoids the complexity of solving a large, monolithic problem by breaking it down into more manageable subproblems.
As we move to the fourth notebook, our focus will shift from leveraging AMPL’s language capabilities to utilizing AMPL’s Python API.
We’ll explore three additional methods for implementing Benders decomposition, with the amplpy module managing the algorithm’s control flow.
This transition will further allow us to showcase the parallel processing of subproblems.
Lastly, we will introduce the AMPLS library, which is designed to boost the efficiency of the Benders decomposition by enabling the export of AMPL model instances as persistent solver representations.
Tags: amplpy, ampl, mip, stochastic, facility location, benders, decomposition
Notebook author: 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)
Introduction
This notebook begins by formulating a general two-stage stochastic programming problem using matrix and vector notation, setting the foundation for our exploration.
We then delve into the theory behind Benders Decomposition, providing a comprehensive understanding of its principles.
Following the theoretical framework, we describe the Benders Decomposition algorithm, tailored to a general two-stage stochastic programming problem, again utilizing matrix and vector notation for consistency.
This sets the stage for applying the decomposition approach to the extensive form of the stochastic facility location problem, as introduced in the second part of this tutorial series.
Once the theoretical groundwork is laid, we transition to practical application by developing the AMPL model for our decomposed problem.
We introduce an AMPL script designed to implement the Benders Decomposition algorithm, leveraging AMPL’s robust language features.
These include defining multiple problems within a single model, dynamically modifying parameter values and their indices, and utilizing AMPL’s loop constructs to iterate through different problems until an optimal solution is discovered.
Although this approach necessitates starting the problem-solving process from scratch in each iteration, it provides a valuable rapid prototyping avenue for stochastic programming problems that are well-suited to decomposition strategies.
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 & & & \
& p^i \in P,; w^j \in P & & & & \
\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}
$
AMPL implementation
Converting the mathematical formulation of our decomposition scheme into an AMPL model is straightforward.
Similar to the extensive form, AMPL enables the modeling of problems in a way that closely aligns with their conceptual structure.
In AMPL, sets and parameters are globally accessible, allowing for seamless integration across different problem components.
AMPL facilitates the assignment of variables, objective functions, and constraints to specific named problems using the syntax problem <problemName>: <list of model entities>.
In the model file presented below, we’ve organized the content to clearly indicate which elements pertain to the master problem and which to the subproblems.
This organization is primarily for clarity and ease of understanding; the AMPL interpreter treats the entire file as a single problem until it encounters explicit problem statements.
Beyond the entities familiar from the extensive form, we introduce additional parameters like sub_facility_open, which convey the master problem’s solutions to the subproblems.
We also incorporate parameters for the dual values from the subproblems, namely customer_price and facility_price.
New constraints, often referred to as cuts, are added to the master problem based on the subproblems’ feasibility and optimality relative to the master problem’s objectives.
Furthermore, a new set of variables, sub_variable_cost, is added to the master problem. These variables account for the variable costs identified in each scenario’s subproblem, contributing to the overall total cost calculation.
Data file
For this implementation, we can utilize the same data file that was employed for the extensive form.
This is feasible because the values needed for our newly introduced parameters will be dynamically derived from the solutions of the subproblems.
The Benders algorithm script
Our AMPL script, named floc_bend.run, initiates the process by loading the necessary model and data files.
To ensure the retrieval of dual rays from the solver, we disable both AMPL’s and the solver’s presolve features.
We configure the solver to employ the dual simplex method, which is particularly suited for our requirements.
We introduce the .dunbdd suffix to access the dual rays.
It’s worth noting that the .dual suffix for constraints is predefined in AMPL and doesn’t require explicit definition.
Following this setup, we leverage the problem statement syntax to delineate the Master and Sub problems.
Before entering the main loop, we establish several parameters to manage the loop’s logic and set the initial values of the master problem’s variable sub_variable_cost to zero.
The looping logic is straightforward.
It involves a counter nITER for tracking the number of iterations, a parameter for identifying the current scenario or sub-problem, and a no_violation parameter to monitor any optimality or feasibility violations during each iteration.
The loop continues until the solutions for all sub-problems are free of violations, indicating that the problem has been solved to optimality.
Key AMPL constructs utilized in this script include:
The problem command, which enables the definition of the master and sub-problems.
The let command for modifying data and variable values.
Looping constructs such as for and repeat.
Conditional statements like if <condition> then <statement> to handle various scenarios within the loop.
The break command to exit from the loop.
Running the script
In this Jupyter Notebook environment, we’ll leverage AMPL’s Python API to execute the script outlined above.
For those not using a notebook, the script can be run directly from a terminal, provided the three files floc_bend.mod, floc_ef.dat, and floc_bend.run are located in the same working directory.
To do this, simply enter the command ampl floc_bend.run in the terminal.
ITERATION 1
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 11621793.46
4 simplex iterations
0 barrier iterations
1: Optimality cut added for scenario Low
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 14779784.86
5 simplex iterations
0 barrier iterations
1: Optimality cut added for scenario Medium
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 21050711.2
6 simplex iterations
0 barrier iterations
1: Optimality cut added for scenario High
SOLVING MASTER PROBLEM
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 16688018.6
0 simplex iterations
1 branching nodes
ITERATION 2
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 11621793.46
0 simplex iterations
0 barrier iterations
2: No cut needed for scenario Low
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 14966984.86
2 simplex iterations
0 barrier iterations
2: Optimality cut added for scenario Medium
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 21570711.2
0 simplex iterations
0 barrier iterations
2: No cut needed for scenario High
SOLVING MASTER PROBLEM
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 16758018.6
1 simplex iterations
1 branching nodes
ITERATION 3
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 11621793.46
0 simplex iterations
0 barrier iterations
3: No cut needed for scenario Low
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 14779784.86
1 simplex iterations
0 barrier iterations
3: No cut needed for scenario Medium
HiGHS 1.6.0: pre:solve = off
alg:simplex = 1
HiGHS 1.6.0: optimal solution; objective 21050711.2
1 simplex iterations
0 barrier iterations
3: No cut needed for scenario High
OPTIMAL SOLUTION FOUND
Expected cost = 16758018.596250
facility_open [*] :=
Baton_Rouge_LA 1
Baytown_TX 1
Beaumont_TX 1
;
Conclusion
Implementing the Benders decomposition algorithm was straightforward, largely because we built upon the extensive form of our stochastic programming problem.
However, there are a few important caveats to consider in our approach:
API Integration: In production environments, AMPL is often used through an API. It would be beneficial to shift the algorithm’s control flow to the API level for better integration and flexibility.
Parallelization of Subproblems: Our analysis highlighted that the subproblems are independent, allowing for the possibility of parallel execution. For complex subproblems or when dealing with a large number of scenarios, parallelizing these computations could markedly improve efficiency.
Solver Instances: Each time we issue a solve command for the master or subproblems, a new solver instance is created. While AMPL provides initial values, the solver might not utilize them, essentially solving each problem from scratch. This approach can overlook valuable information that could be leveraged if the problems were persistent.
In the following section, we will explore how to utilize AMPL’s Python API to move the algorithm’s control flow to the API level.