authors/gomfy/dev/3_benders_stoch_floc.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# AMPL Development Tutorial 3/6 -- Benders Decomposition via AMPL scripting\n",
        "[![3_benders_stoch_floc.ipynb](https://img.shields.io/badge/github-%23121011.svg?logo=github)](https://github.com/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/3_benders_stoch_floc.ipynb) [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/3_benders_stoch_floc.ipynb) [![Kaggle](https://kaggle.com/static/images/open-in-kaggle.svg)](https://kaggle.com/kernels/welcome?src=https://github.com/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/3_benders_stoch_floc.ipynb) [![Gradient](https://assets.paperspace.io/img/gradient-badge.svg)](https://console.paperspace.com/github/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/3_benders_stoch_floc.ipynb) [![Open In SageMaker Studio Lab](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/3_benders_stoch_floc.ipynb) [![Hits](https://h.ampl.com/https://github.com/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/3_benders_stoch_floc.ipynb)](https://colab.ampl.com)\n",
        "\n",
        "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. \n",
        "While this approach offers superior results compared to the deterministic model, it significantly increases the problem's complexity.\n",
        "In this section, we delve into the Benders decomposition technique, applying it to our two-stage stochastic problem using AMPL's advanced features. \n",
        "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.\n",
        "As we move to the fourth notebook, our focus will shift from leveraging AMPL's language capabilities to utilizing AMPL's Python API. \n",
        "We'll explore three additional methods for implementing Benders decomposition, with the `amplpy` module managing the algorithm's control flow. \n",
        "This transition will further allow us to showcase the parallel processing of subproblems.\n",
        "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.\n",
        "\n",
        "Tags: amplpy, ampl, mip, stochastic, facility location, benders, decomposition\n",
        "\n",
        "Notebook author: Gyorgy Matyasfalvi <<gyorgy@ampl.com>>\n",
        "\n",
        "References:\n",
        "\n",
        "* AMPL a Modeling Language for Mathematical Programming -- Robert Fourer et al.\n",
        "* Intro to Linear Optimization -- Dimitris Bertsimas et al.\n",
        "* SCIP Optimization Suite example projects -- Stephen J. Maher (ZIB)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Install dependencies\n",
        "%pip install -q amplpy pandas"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Google Colab & Kaggle integration\n",
        "from amplpy import AMPL, ampl_notebook\n",
        "\n",
        "ampl = ampl_notebook(\n",
        "    modules=[\"open\"],  # modules to install\n",
        "    license_uuid=\"default\",  # license to use\n",
        ")  # instantiate AMPL object and register magics"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Import all necessary libraries\n",
        "import pandas as pd"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Introduction\n",
        "\n",
        "This notebook begins by formulating a general two-stage stochastic programming problem using matrix and vector notation, setting the foundation for our exploration. \n",
        "We then delve into the theory behind Benders Decomposition, providing a comprehensive understanding of its principles.\n",
        "\n",
        "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. \n",
        "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.\n",
        "\n",
        "Once the theoretical groundwork is laid, we transition to practical application by developing the AMPL model for our decomposed problem. \n",
        "We introduce an AMPL script designed to implement the Benders Decomposition algorithm, leveraging AMPL's robust language features. \n",
        "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.\n",
        "\n",
        "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."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Two-stage stochastic programming problem formulation\n",
        "\n",
        "Consider a decision-making process occurring over two sequential stages. \n",
        "Initially, a decision vector $x$ is selected. \n",
        "Following this, new information becomes available, prompting the selection of a second-stage decision vector $y$. \n",
        "We introduce $S$ possible scenarios for the newly acquired information, with the actual scenario being revealed only after the selection of $x$. \n",
        "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.\n",
        "\n",
        "The decision vectors $x$ and $y^s$ incur costs represented by vectors $f$ and $q$, respectively. \n",
        "The first-stage decisions are constrained by:\n",
        "\n",
        "\\begin{equation}\n",
        "\\begin{array}{crlllllll}\n",
        "\t& A x \t\t&  =  \t& {b} \\\\\n",
        "\t& x\t\t\t& \\geq\t& 0.\n",
        "\\end{array}\n",
        "\\tag{1}\n",
        "\\end{equation}\n",
        "\n",
        "Furthermore, the combination of first and second-stage decisions must satisfy the scenario-specific constraints:\n",
        "\n",
        "\\begin{equation}\n",
        "\\begin{array}{crcllllll}\n",
        "\t& {B}^s x + {D}^s y^s \t\t& =\t\t& {d}^s \\\\\n",
        "\t& y^s\t\t\t\t\t\t\t\t\t\t& \\geq \t& \t0,\n",
        "\\end{array}\n",
        "\\tag{2}\n",
        "\\end{equation}\n",
        "\n",
        "for every scenario $s$. \n",
        "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}$. \n",
        "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. \n",
        "The objective is to determine the choices of $x$ and $y^s$ for all $s=1, \\ldots, S$ that minimize the expected total cost:\n",
        "\n",
        "\\begin{equation}\n",
        "\\begin{array}{cllllllll}\n",
        "\t& 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 \n",
        "\\end{array}\n",
        "\\tag{3}\n",
        "\\end{equation}\n",
        "\n",
        "This formulation leads to the following two-stage stochastic optimization problem:\n",
        "\n",
        "\n",
        "\\begin{equation}\n",
        "\\begin{array}{crcrcrccrcl}\n",
        "\\min\t\t\t& f^T x \t\t\t& + & \\alpha^1 {(q^1)^T} y_1 & + \t& \\alpha^2 {(q^1)^T} y^2 \t& \\ldots  \t& + & \\alpha^s {(q^s)^T} y^S\t&   & \\\\\n",
        "\\textrm{s.t.}\t& A x \t\t\t\t&   &\t\t \t\t\t   &\t& \t\t\t\t \t\t& \t \t\t& \t&\t\t\t\t\t\t&   &= {b} \\\\\n",
        "\t\t\t\t& {B^1} x \t\t\t& + & {D^1} y^1 \t \t   &    &\t\t \t\t\t \t&\t \t\t&\t&\t\t\t\t\t\t& \t&= {d^1} \\\\\n",
        "\t\t\t\t& {B^2} x \t\t\t& + &\t\t\t\t\t\t&   & {D^2} y^2 \t \t\t&\t\t\t&\t&\t\t\t\t\t\t& \t&= {d^2} \\\\\n",
        "\t\t\t\t& \\vdots \t\t\t&   &  \t\t\t  \t\t  \t&   & \t\t\t\t\t\t & \\ddots \t&\t&\t\t\t\t\t\t& \t&\\vdots \\\\\n",
        "\t\t\t\t& {B^S} x \t\t\t& + & \t  \t\t\t\t  \t&   & \t\t\t\t\t\t &\t \t  \t& \t& {D^S} y^S\t\t\t\t&\t&= {d^S} \\\\\n",
        "\t\t\t\t& x \\geq 0, \t\t& \t& y^1 \\geq 0,\t\t  & \\ldots\t&\t\t\t \t\t&\t\t\t& \t& y^S \\geq 0\t\t\t&   &\n",
        "\\end{array}\n",
        "\\tag{4}\n",
        "\\end{equation}\n",
        "\n",
        "we will refer to the above problem as the **master problem**. \n",
        "Notice that even if the number of possible scenarios $S$ is moderate, this formulation can be a large linear programming problem. \n",
        "However, a decomposition method can help."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Benders decomposition for two-stage problems\n",
        "\n",
        "1. Given a vector $\\hat{x}$, using the dual simplex method solve the subproblem (for all $s \\in S$):\n",
        "\n",
        "    \\begin{equation}\n",
        "    \\begin{array}{crl}\n",
        "        \\min\t    &  (q^s)^T y^s\t& \\\\\n",
        "    \\textrm{s.t.}   &  D^s y^s \t\t& = d^s - B^s \\hat{x} \\\\\n",
        "    \t            &  y^s       \t& \\geq 0 \n",
        "    \\end{array}\n",
        "    \\tag{5}\n",
        "    \\end{equation}\n",
        "\n",
        "2. The dual of (5) takes the following form:\n",
        "\n",
        "    \\begin{equation}\n",
        "    \\begin{array}{crl}\n",
        "        \\max\t    & (p^s)^T (d^s - B^s \\hat{x}) & \\\\\n",
        "    \\textrm{s.t.}   & (p^s)^T D^s           \t\t& <=  (q^s)^T\n",
        "    \\end{array}\n",
        "    \\tag{6}\n",
        "    \\end{equation}\n",
        "\n",
        "    Let  \n",
        "    \n",
        "    \\begin{equation}\n",
        "    \\begin{array}{crrl}\n",
        "        \t    & P \\;\\; = & \\{ p & | \\;\\; p^T D \\leq q^T \\}  & \\\\\n",
        "    \\end{array}\n",
        "    \\tag{7}\n",
        "    \\end{equation}\n",
        "\n",
        "    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$.\n",
        "\n",
        "3. 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:\n",
        "\n",
        "    \\begin{equation}\n",
        "    \\begin{array}{crcrcrccrcl}\n",
        "    \\min\t\t\t& f^T x \t    \t& + & \\sum_{s=1}^S \\alpha^s z^s(x) &    \\\\\n",
        "    \\textrm{s.t.}\t& A x \t\t\t\t&   &\t\t \t\t\t\t\t& = {b} \\\\\n",
        "    \t\t\t\t& x          \t    &\t\\geq 0 &              \t\t&\n",
        "    \\end{array}\n",
        "    \\tag{8}\n",
        "    \\end{equation}\n",
        "\n",
        "    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 \n",
        "    primal (5) is infeasible).\n",
        "\n",
        "4. 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:\n",
        "\n",
        "\n",
        "    \\begin{equation}\n",
        "    \\begin{array}{crcrcrccrcl}\n",
        "    \\min\t\t\t& f^T x \t    \t    & + & \\sum_{s=1}^S \\alpha^s z^s &           & \\\\\n",
        "    \\textrm{s.t.}\t& A x \t\t\t\t    &   &\t\t \t\t\t\t\t& = b       & \\\\\n",
        "                \t& (p^i)^T (d^s - B^s x)\t&   &\t\t \t\t\t\t\t& \\leq z^s  & \\forall \\; i, s \\\\\n",
        "                \t& (w^j)^T (d^s - B^s x)\t&   &\t\t \t\t\t\t\t& \\leq 0    & \\forall \\; j, s \\\\\n",
        "    \t\t\t\t& x          \t        &\t\\geq 0 &              \t\t&           & \\\\\n",
        "                    & p^i \\in P,\\; w^j \\in P &   &                           &           & \\\\\n",
        "\n",
        "    \\end{array}\n",
        "    \\tag{9}\n",
        "    \\end{equation}\n",
        "\n",
        "5. **In our algorithm we will only generate constraints that we find to be violated by the current solution.**\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The algorithm\n",
        "\n",
        "0. **Initialization step**: initialize $(\\hat{x}, \\hat{z}^1, \\ldots, \\hat{z}^S)$ to zero.\n",
        "\n",
        "1. **Sub-problem step**: using the dual-simplex method solve the sub-problems individually for every $s \\in S$:\n",
        "\n",
        "\\begin{equation}\n",
        "\\begin{array}{crl}\n",
        "    \\min\t    &  (q^s)^T y^s\t& \\\\\n",
        "\\textrm{s.t.}   &  D^s y^s \t\t& = d^s - B^s \\hat{x} \\\\\n",
        "\t            &  y^s       \t& \\geq 0 \n",
        "\\end{array}\n",
        "\\tag{5}\n",
        "\\end{equation}\n",
        "\n",
        "2. **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.\n",
        "\n",
        "3. **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:  \n",
        "\\begin{equation}\n",
        "(p^s)^T (d^s - B^s x) \\leq z^s,\n",
        "\\tag{10}\n",
        "\\end{equation}\n",
        "where $p^s$ is an optimal basic feasible solution to the dual of the sub-problem. \n",
        "\n",
        "4. **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:  \n",
        "\\begin{equation}\n",
        "(w^s)^T (d^s - B^s x) \\leq 0, \n",
        "\\tag{11}\n",
        "\\end{equation}\n",
        "where $w^s$ is a positive cost extreme ray. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Stochastic facility location a concrete example\n",
        "\n",
        "**Given:** \n",
        "* A set of facilities: $I$.\n",
        "* A set of customers: $J$.\n",
        "* Set of scenarios: $S$ (representing different customer demands).\n",
        "\n",
        "**Task:** \n",
        "* Find the minimum cost facilities to open such that the customer demand can be satisfied in all scenarios.\n",
        "\n",
        "### Variables\n",
        "* $x_i \\in \\{0, 1\\} \\quad \\forall i \\in I$\n",
        "    * $x_i = 1$ if facility $i$ is opened.\n",
        "* $y_{ij}^s \\geq 0 \\quad \\forall i \\in I, \\forall j \\in J, \\forall s \\in S$\n",
        "    * $y_{ij}^s$ is the level of demand for customer $j$ satisfied by facility $i$ in scenario $s$.\n",
        "\n",
        "### Parameters:\n",
        "* $\\alpha^s$: the probability of scenario $s$.\n",
        "* $f_i$: the fixed cost for opening facility $i$,\n",
        "* $q_{ij}$: the cost of servicing customer $j$ from facility $i$,\n",
        "* $\\lambda_j^s$: the demand of customer $j$ in scenario $s$,\n",
        "* $k_i:$ the capacity of facility $i$.\n",
        "\n",
        "\n",
        "### The extensive form\n",
        "The extensive form of our stochastic program can be formulated as follows:\n",
        "\n",
        "$\n",
        "\\begin{equation}\n",
        "\\begin{array}{rrlrrll}\n",
        " \\min \\quad & \\sum_{i \\in I} f_i x_i & + & \\sum_{s \\in S} \\sum_{i \\in I} \\sum_{j \\in J} \\alpha^s q_{ij} y_{ij}^s & & & \\\\\n",
        " & & & & & & \\\\\n",
        " \\textrm{s.t.} \\quad & & & \\sum_{i \\in I} y_{ij}^s & \\geq & \\lambda_j^s & \\forall j \\in J, \\forall s \\in S \\\\\n",
        " & & & \\sum_{j \\in J} y_{ij}^s & \\leq & k_i x_i & \\forall i \\in I, \\forall s \\in S \\\\\n",
        " & \\sum_{i \\in I} k_i x_i & & & \\geq & \\max_{s \\in S} \\sum_{j \\in J} \\lambda_j^s & \\\\\n",
        " & & & & & & \\\\\n",
        " & x_i  & \\in \\{0, 1\\} & & & & \\forall i \\in I \\\\\n",
        " & y_{ij}^s & \\geq 0 & & & & \\forall i \\in I, \\forall j \\in J, \\forall s \\in S  \n",
        "\\end{array} \\tag{12}\n",
        "\\end{equation}\n",
        "$\n",
        "\n",
        "Rearranged to fit the general form introduced above:\n",
        "\n",
        "$\n",
        "\\begin{equation}\n",
        "\\begin{array}{rrlrrll}\n",
        " \\min \\quad & \\sum_{i \\in I} f_i x_i & + & \\sum_{s \\in S} \\sum_{i \\in I} \\sum_{j \\in J} \\alpha^s q_{ij} y_{ij}^s & & & \\\\\n",
        " & & & & & & \\\\\n",
        " \\textrm{s.t.} \\quad & \\sum_{i \\in I} k_i x_i & & & \\geq & \\max_{s \\in S} \\sum_{j \\in J} \\lambda_j^s & \\\\\n",
        " & - k_i x_i & & \\sum_{j \\in J} y_{ij}^s & \\leq & 0 & \\forall i \\in I, \\forall s \\in S \\\\\n",
        " & & & \\sum_{i \\in I} y_{ij}^s & \\geq & \\lambda_j^s & \\forall j \\in J, \\forall s \\in S \\\\\n",
        " & & & & & & \\\\\n",
        " & x_i  & \\in \\{0, 1\\} & & & & \\forall i \\in I \\\\\n",
        " & y_{ij}^s & \\geq 0 & & & & \\forall i \\in I, \\forall j \\in J, \\forall s \\in S  \n",
        "\\end{array} \\tag{12}\n",
        "\\end{equation}\n",
        "$\n",
        "\n",
        "### The sub-problem\n",
        "Given the above formulation for the extensive form (12) we can express our scenario specific sub-problems as follows\n",
        "\n",
        "$\n",
        "\\begin{equation}\n",
        "\\begin{array}{lcrllll}\n",
        " z_s(\\hat{x}) & = & \\min \\quad &  \\sum_{i \\in I} \\sum_{j \\in J} q_{ij} y_{ij}^s & & &  \\\\\n",
        " & & & & & & \\\\\n",
        " & & \\textrm{s.t.} \\quad & \\sum_{i \\in I} y_{ij}^s & \\geq  & \\lambda_j^s & \\forall j \\in J \\\\\n",
        " & & & \\sum_{j \\in J} y_{ij}^s & \\leq & k_i \\hat{x_i} & \\forall i \\in I \\\\\n",
        " & & & & & &  \\\\\n",
        " & & & y_{ij}^s \\geq 0 & & & \\forall \\; i \\in I, \\forall j \\in J\n",
        "\\end{array} \\tag{13}\n",
        "\\end{equation}\n",
        "$\n",
        "\n",
        "### The master-problem\n",
        "The master-problem takes the following form\n",
        "\n",
        "$\n",
        "\\begin{equation}\n",
        "\\begin{array}{rllll}\n",
        " \\min \\quad          & \\sum_{i \\in I} f_i x_i + \\sum_{s \\in S} \\alpha^s z^s & & & \\\\\n",
        " & & & & \\\\\n",
        " \\textrm{s.t.} \\quad & \\sum_{i \\in I} k_i x_i & \\geq & \\max_{s \\in S} \\sum_{j \\in J} \\lambda_j^s & \\\\\n",
        "                     & \\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 \\\\\n",
        "                     & \\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 \\\\\n",
        " & & & & \\\\\n",
        "                     & x_i \\in \\{0, 1\\} & & & \\forall i \\in I \\\\\n",
        "\\end{array} \\tag{12}\n",
        "\\end{equation}\n",
        "$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## AMPL implementation\n",
        "\n",
        "Converting the mathematical formulation of our decomposition scheme into an AMPL model is straightforward. \n",
        "Similar to the extensive form, AMPL enables the modeling of problems in a way that closely aligns with their conceptual structure. \n",
        "In AMPL, sets and parameters are globally accessible, allowing for seamless integration across different problem components.\n",
        "\n",
        "AMPL facilitates the assignment of variables, objective functions, and constraints to specific named problems using the syntax `problem <problemName>: <list of model entities>`. \n",
        "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. \n",
        "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.\n",
        "\n",
        "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. \n",
        "We also incorporate parameters for the dual values from the subproblems, namely `customer_price` and `facility_price`. \n",
        "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. \n",
        "\n",
        "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."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Writing floc_bend.mod\n"
          ]
        }
      ],
      "source": [
        "%%writefile floc_bend.mod\n",
        "\n",
        "# Global sets\n",
        "set FACILITIES; # set of facilities\n",
        "set CUSTOMERS;  # set of customers\n",
        "set SCENARIOS;  # set of scenarios\n",
        "\n",
        "# Global parameters\n",
        "param customer_demand{CUSTOMERS, SCENARIOS} >= 0;   # customer_demand of customer j in scenario s\n",
        "param facility_capacity{FACILITIES} >= 0;           # facility_capacity of facility_open i\n",
        "\n",
        "\n",
        "### SUBPROBLEM ###\n",
        "\n",
        "# Subproblem parameters\n",
        "param variable_cost{FACILITIES, CUSTOMERS} >= 0;    # variable cost of satisfying customer_demand of customer j from facility_open i\n",
        "\n",
        "# Bender's parameters\n",
        "param sub_facility_open{FACILITIES} binary, default 1; # 1 if facility i is open, 0 otherwise\n",
        "param sub_scenario symbolic in SCENARIOS;\n",
        "\n",
        "# Subproblem variables\n",
        "var production{FACILITIES, CUSTOMERS, SCENARIOS} >= 0;  # production from facility i to satisfy customer demand j in scenario s\n",
        "\n",
        "# Subproblem objective\n",
        "minimize operating_cost: \n",
        "    sum{i in FACILITIES, j in CUSTOMERS} variable_cost[i,j]*production[i,j,sub_scenario]; \n",
        "\n",
        "# Subproblem constraints\n",
        "s.t. satisfying_customer_demand{j in CUSTOMERS}:\n",
        "    sum{i in FACILITIES} production[i,j,sub_scenario] >= customer_demand[j,sub_scenario];\n",
        "\n",
        "s.t. facility_capacity_limits{i in FACILITIES}:\n",
        "    sum{j in CUSTOMERS} production[i,j,sub_scenario] <= facility_capacity[i] * sub_facility_open[i];\n",
        "\n",
        "\n",
        "\n",
        "### MASTER PROBLEM ### \n",
        "\n",
        "# Master parameters\n",
        "param fixed_cost{FACILITIES} >= 0;               # fixed cost of opening facility_open i\n",
        "param prob{SCENARIOS} default 1/card(SCENARIOS); # probability of scenario s\n",
        "\n",
        "# Bender's parameters\n",
        "param nITER >= 0 integer;\n",
        "param cut_type {1..nITER, SCENARIOS} symbolic within {\"opt\", \"feas\", \"none\"};\n",
        "param customer_price{CUSTOMERS, SCENARIOS, 1..nITER} >= 0;  # dual variable for the demand constraint\n",
        "param facility_price{FACILITIES, SCENARIOS, 1..nITER} <= 0; # dual variable for the capacity constraint\n",
        "\n",
        "# Master variables\n",
        "var sub_variable_cost{SCENARIOS} >= 0; # subproblem objective bound\n",
        "var facility_open{FACILITIES} binary;  # 1 if facility i is open, 0 otherwise\n",
        "\n",
        "# Master objective\n",
        "minimize total_cost: \n",
        "    sum{i in FACILITIES} fixed_cost[i]*facility_open[i] + sum{s in SCENARIOS} prob[s]*sub_variable_cost[s];\n",
        "\n",
        "# Master constraints\n",
        "s.t. sufficient_production_capacity:\n",
        "    sum{i in FACILITIES} facility_capacity[i]*facility_open[i] >= max{s in SCENARIOS} sum{j in CUSTOMERS} customer_demand[j,s];\n",
        "\n",
        "s.t. optimality_cut {k in 1..nITER, s in SCENARIOS: cut_type[k,s] == \"opt\"}:\n",
        "        sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s] + \n",
        "        sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i] <= sub_variable_cost[s]; \n",
        "\n",
        "s.t. feasibility_cut {k in 1..nITER, s in SCENARIOS: cut_type[k,s] == \"feas\"}: \n",
        "        sum{j in CUSTOMERS} customer_price[j,s,k]*customer_demand[j,s] + \n",
        "        sum{i in FACILITIES} facility_price[i,s,k]*facility_capacity[i]*facility_open[i] <= 0; "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Data file\n",
        "\n",
        "For this implementation, we can utilize the same data file that was employed for the extensive form. \n",
        "This is feasible because the values needed for our newly introduced parameters will be dynamically derived from the solutions of the subproblems."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Overwriting floc_ef.dat\n"
          ]
        }
      ],
      "source": [
        "%%writefile floc_ef.dat\n",
        "set FACILITIES  := Baytown_TX Beaumont_TX Baton_Rouge_LA;\n",
        "set CUSTOMERS   := San_Antonio_TX Dallas_TX Jackson_MS Birmingham_AL;\n",
        "set SCENARIOS   := Low Medium High;\n",
        "\n",
        "param prob := Low 0.25 Medium 0.5 High 0.25;\n",
        "\n",
        "param fixed_cost := Baytown_TX 400000 Beaumont_TX 200000 Baton_Rouge_LA 600000;\n",
        "\n",
        "param facility_capacity := Baytown_TX 1550 Beaumont_TX 650 Baton_Rouge_LA 1750;\n",
        "\n",
        "param variable_cost:        San_Antonio_TX  Dallas_TX    Jackson_MS   Birmingham_AL :=\n",
        "             Baytown_TX     5739.725        6539.725     8650.40      22372.1125\n",
        "             Beaumont_TX    6055.05         6739.055     8050.40      21014.225 \n",
        "             Baton_Rouge_LA 8650.40         7539.055     4539.72      15024.325;\n",
        "\n",
        "param customer_demand:          Low    Medium   High :=\n",
        "               San_Antonio_TX   450    650      887 \n",
        "               Dallas_TX        910    1134     1456      \n",
        "               Jackson_MS       379    416      673\n",
        "               Birmingham_AL    91     113      207;"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### The Benders algorithm script\n",
        "\n",
        "Our AMPL script, named `floc_bend.run`, initiates the process by loading the necessary model and data files. \n",
        "To ensure the retrieval of dual rays from the solver, we disable both AMPL's and the solver's presolve features. \n",
        "We configure the solver to employ the dual simplex method, which is particularly suited for our requirements.\n",
        "\n",
        "We introduce the `.dunbdd` suffix to access the dual rays. \n",
        "It's worth noting that the `.dual` suffix for constraints is predefined in AMPL and doesn't require explicit definition.\n",
        "\n",
        "Following this setup, we leverage the problem statement syntax to delineate the `Master` and `Sub` problems. \n",
        "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.\n",
        "\n",
        "The looping logic is straightforward. \n",
        "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. \n",
        "The loop continues until the solutions for all sub-problems are free of violations, indicating that the problem has been solved to optimality.\n",
        "\n",
        "Key AMPL constructs utilized in this script include:\n",
        "- The `problem` command, which enables the definition of the master and sub-problems.\n",
        "- The `let` command for modifying data and variable values.\n",
        "- Looping constructs such as `for` and `repeat`.\n",
        "- Conditional statements like `if <condition> then <statement>` to handle various scenarios within the loop.\n",
        "- The `break` command to exit from the loop."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Writing floc_bend.run\n"
          ]
        }
      ],
      "source": [
        "%%writefile floc_bend.run\n",
        "\n",
        "# Load model and data\n",
        "model floc_bend.mod;\n",
        "data floc_ef.dat;\n",
        "\n",
        "# Set options\n",
        "option presolve 0;\n",
        "# Solver options\n",
        "option solver highs;\n",
        "option highs_options 'pre:solve off alg:simplex 1';\n",
        "# Display optoins\n",
        "option omit_zero_rows 1; option display_eps .000001;\n",
        "\n",
        "# Define suffix for storing dual rays\n",
        "suffix dunbdd;\n",
        "\n",
        "# Define Master and Sub problems\n",
        "problem Master: facility_open, sub_variable_cost, total_cost, optimality_cut, feasibility_cut, sufficient_production_capacity; \n",
        "problem Sub: production, operating_cost, satisfying_customer_demand, facility_capacity_limits;\n",
        "\n",
        "# Loop logic parameters\n",
        "param epsilon := 0.00001;\n",
        "param no_violation{SCENARIOS} binary, default 0;\n",
        "\n",
        "# Initialzie variables and parameters\n",
        "let {s in SCENARIOS} sub_variable_cost[s] := 0;\n",
        "let nITER := 0;\n",
        "\n",
        "# Start Bender's loop\n",
        "repeat { \n",
        "   let nITER := nITER + 1;\n",
        "   printf \"\\nITERATION %d\\n\\n\", nITER;\n",
        "\n",
        "   # Loop through the scenarios and solve each subproblem and update dual prices for each scenario\n",
        "   for {s in SCENARIOS} {\n",
        "      let sub_scenario := s;\n",
        "      solve Sub;\n",
        "      if Sub.result = \"infeasible\" then {\n",
        "         let cut_type[nITER, s] := \"feas\";\n",
        "         let {j in CUSTOMERS} customer_price[j,sub_scenario,nITER] := satisfying_customer_demand[j].dunbdd;\n",
        "         let {i in FACILITIES} facility_price[i,sub_scenario,nITER] := facility_capacity_limits[i].dunbdd;  \n",
        "         printf \"%d: Feasibility cut added for scenario %s\\n\", nITER, s;\n",
        "      }\n",
        "      else if operating_cost > sub_variable_cost[s] + epsilon then {\n",
        "         let cut_type[nITER, s] := \"opt\";\n",
        "         let {j in CUSTOMERS} customer_price[j,sub_scenario,nITER] := satisfying_customer_demand[j].dual;\n",
        "         let {i in FACILITIES} facility_price[i,sub_scenario,nITER] := facility_capacity_limits[i].dual;\n",
        "         printf \"%d: Optimality cut added for scenario %s\\n\", nITER, s;\n",
        "      }\n",
        "      else {\n",
        "         let cut_type[nITER, s] := \"none\";\n",
        "         let no_violation[s] := 1;\n",
        "         printf \"%d: No cut needed for scenario %s\\n\", nITER, s;\n",
        "      }\n",
        "   }\n",
        "\n",
        "   # Check if we have optimal solution\n",
        "   if sum{s in SCENARIOS} no_violation[s] == card(SCENARIOS) then break;\n",
        "   # Reset no violation\n",
        "   let {s in SCENARIOS} no_violation[s] := 0;\n",
        "\n",
        "   # If not optimal, resolve master\n",
        "   printf \"\\nSOLVING MASTER PROBLEM\\n\\n\";\n",
        "   solve Master;\n",
        "\n",
        "   # Update the subproblem right-hand side with master solution\n",
        "   let {i in FACILITIES} sub_facility_open[i] := facility_open[i];\n",
        "};\n",
        "\n",
        "printf \"\\nOPTIMAL SOLUTION FOUND\\nExpected cost = %f\\n\\n\", total_cost;\n",
        "display facility_open;\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Running the script\n",
        "\n",
        "In this Jupyter Notebook environment, we'll leverage AMPL's Python API to execute the script outlined above. \n",
        "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. \n",
        "To do this, simply enter the command `ampl floc_bend.run` in the terminal."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "ITERATION 1\n",
            "\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 11621793.46\n",
            "4 simplex iterations\n",
            "0 barrier iterations\n",
            "1: Optimality cut added for scenario Low\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 14779784.86\n",
            "5 simplex iterations\n",
            "0 barrier iterations\n",
            "1: Optimality cut added for scenario Medium\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 21050711.2\n",
            "6 simplex iterations\n",
            "0 barrier iterations\n",
            "1: Optimality cut added for scenario High\n",
            "\n",
            "SOLVING MASTER PROBLEM\n",
            "\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 16688018.6\n",
            "0 simplex iterations\n",
            "1 branching nodes\n",
            "\n",
            "ITERATION 2\n",
            "\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 11621793.46\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "2: No cut needed for scenario Low\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 14966984.86\n",
            "2 simplex iterations\n",
            "0 barrier iterations\n",
            "2: Optimality cut added for scenario Medium\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 21570711.2\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "2: No cut needed for scenario High\n",
            "\n",
            "SOLVING MASTER PROBLEM\n",
            "\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 16758018.6\n",
            "1 simplex iterations\n",
            "1 branching nodes\n",
            "\n",
            "ITERATION 3\n",
            "\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 11621793.46\n",
            "0 simplex iterations\n",
            "0 barrier iterations\n",
            "3: No cut needed for scenario Low\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 14779784.86\n",
            "1 simplex iterations\n",
            "0 barrier iterations\n",
            "3: No cut needed for scenario Medium\n",
            "HiGHS 1.6.0:   pre:solve = off\n",
            "  alg:simplex = 1\n",
            "\bHiGHS 1.6.0: optimal solution; objective 21050711.2\n",
            "1 simplex iterations\n",
            "0 barrier iterations\n",
            "3: No cut needed for scenario High\n",
            "\n",
            "OPTIMAL SOLUTION FOUND\n",
            "Expected cost = 16758018.596250\n",
            "\n",
            "facility_open [*] :=\n",
            "Baton_Rouge_LA  1\n",
            "    Baytown_TX  1\n",
            "   Beaumont_TX  1\n",
            ";\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ampl = AMPL()  # create an AMPL object\n",
        "ampl.read(\"floc_bend.run\")  # Execute the .run file"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Conclusion\n",
        "\n",
        "Implementing the Benders decomposition algorithm was straightforward, largely because we built upon the extensive form of our stochastic programming problem. \n",
        "However, there are a few important caveats to consider in our approach:\n",
        "\n",
        "1. **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.\n",
        "\n",
        "2. **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.\n",
        "\n",
        "3. **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.\n",
        "\n",
        "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."
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.7"
    },
    "vscode": {
      "interpreter": {
        "hash": "ee4198aac13dae9b6fccd443041680d99c2643c4956ef7c80ff6dcc0057cb523"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}