AMPL Development Tutorial 4/6 – Benders Decomposition via PYTHON scripting
Description: In this fourth installment of our six-part series, we advance our exploration by demonstrating how to adapt our AMPL script for use with AMPL’s Python API.
This adaptation facilitates easier integration into production environments by transitioning the control flow of the Benders algorithm to the API level.
We will explore how to import data into AMPL from Pandas data frames and detail the process of extracting data from AMPL for use within Python.
This approach not only streamlines the interaction between AMPL and Python but also sets the stage for implementing parallel processing of subproblems in subsequent sections.
In our final notebook, we will introduce the AMPLS library.
This powerful tool enhances the efficiency of Benders decomposition by allowing for the export of AMPL model instances as persistent solver representations, further optimizing the solution process.
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
To ensure a cohesive understanding and for ease of reference, we will revisit the foundational concepts of stochastic programming and the theory behind Benders decomposition in this notebook.
This repetition will help reinforce these critical concepts and provide a seamless transition into the PYTHON script implementing Benders decomposition of our stochastic facility location problem.
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}
$
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.
Overwriting floc_bend.mod
Data as Python objects
In this notebook, we continue to work with the same dataset as in the previous installments.
However, the format in which we store our data differs; we now utilize Python lists and Pandas data frames. Initially, we create these Python objects and subsequently employ amplpy methods to seamlessly transfer the data into AMPL.
The Benders algorithm script in Python
The Python snippet below begins the process by importing the required model file and configuring the necessary settings.
Additionally, it loads data from the previously defined Python objects.
For lists, we can directly assign data to AMPL sets using the syntax ampl.set['<set name>'] = <Python list holding data>.
When dealing with parameters represented in Pandas dataframes, we employ the set_data() method.
This method matches the dataframe’s column names with corresponding entities in the AMPL model, ensuring proper assignment of values.
It’s important to ensure that the dataframe’s indexing aligns with the AMPL model’s indexing scheme.
To define suffixes and problem structures, we utilize the eval() method.
Unlike in the AMPL-only approach, we no longer need to define most loop logic parameters in AMPL; they will be represented as Python objects instead.
The control flow of the algorithm is managed through Python’s control structures.
The notable exception is the nITER parameter, which remains defined in AMPL to allow us to index constraints and price parameters adequately.
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 THE 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
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
SOLVING THE 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
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
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
OPTIMAL SOLUTION FOUND
Expected cost = 16758018.596250001
facility_open [*] :=
Baton_Rouge_LA 1
Baytown_TX 1
Beaumont_TX 1
;
Conclusion
Transitioning control flow of the Benders algorithm to Python has not altered the core solution of our problem.
However, by shifting the control of the solve executions to Python code, we’ve significantly enhanced the flexibility of our algorithm’s control mechanism.
This newfound flexibility opens the door to sophisticated enhancements, such as parallelizing the scenario iterations within the main while loop.
With the Benders algorithm’s control flow now in Python, implementing such parallelization becomes far more feasible and straightforward.