authors/gomfy/dev/4_benders_in_python_stoch_floc.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# AMPL Development Tutorial 4/6 -- Benders Decomposition via PYTHON scripting\n",
        "[![4_benders_in_python_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/4_benders_in_python_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/4_benders_in_python_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/4_benders_in_python_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/4_benders_in_python_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/4_benders_in_python_stoch_floc.ipynb) [![Hits](https://h.ampl.com/https://github.com/ampl/colab.ampl.com/blob/master/authors/gomfy/dev/4_benders_in_python_stoch_floc.ipynb)](https://colab.ampl.com)\n",
        "\n",
        "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. \n",
        "This adaptation facilitates easier integration into production environments by transitioning the control flow of the Benders algorithm to the API level.\n",
        "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. \n",
        "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.\n",
        "In our final notebook, we will introduce the AMPLS library. \n",
        "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.\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": 7,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Install dependencies\n",
        "%pip install -q amplpy pandas"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "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": 9,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Import all necessary libraries\n",
        "import pandas as pd"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Introduction\n",
        "\n",
        "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. \n",
        "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."
      ]
    },
    {
      "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",
        "    \\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",
        "### 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": 10,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Overwriting 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 as Python objects\n",
        "\n",
        "In this notebook, we continue to work with the same dataset as in the previous installments. \n",
        "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."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Define the data\n",
        "\n",
        "# Define sets\n",
        "FACILITIES = [\"Baytown_TX\", \"Beaumont_TX\", \"Baton_Rouge_LA\"]\n",
        "CUSTOMERS = [\"San_Antonio_TX\", \"Dallas_TX\", \"Jackson_MS\", \"Birmingham_AL\"]\n",
        "SCENARIOS = [\"Low\", \"Medium\", \"High\"]\n",
        "\n",
        "# Define parameters indexed over FACILITIES\n",
        "fixed_cost_df = pd.DataFrame(\n",
        "    [400000, 200000, 600000], index=FACILITIES, columns=[\"fixed_cost\"]\n",
        ")\n",
        "\n",
        "facility_capacity_df = pd.DataFrame(\n",
        "    [1550, 650, 1750], index=FACILITIES, columns=[\"facility_capacity\"]\n",
        ")\n",
        "\n",
        "# Define parameters indexed over SCENARIOS\n",
        "prob_df = pd.DataFrame([0.25, 0.5, 0.25], index=SCENARIOS, columns=[\"prob\"])\n",
        "\n",
        "# Define parameters indexed over CUSTOMERS and SCENARIOS\n",
        "customer_demand_df = pd.DataFrame(\n",
        "    [450, 650, 887, 910, 1134, 1456, 379, 416, 673, 91, 113, 207],\n",
        "    index=pd.MultiIndex.from_product(\n",
        "        [CUSTOMERS, SCENARIOS], names=[\"CUSTOMERS\", \"SCENARIOS\"]\n",
        "    ),\n",
        "    columns=[\"customer_demand\"],\n",
        ")\n",
        "\n",
        "# Define parameters indexed over FACILITIES and CUSTOMERS\n",
        "variable_cost_df = pd.DataFrame(\n",
        "    [\n",
        "        5739.725,\n",
        "        6539.725,\n",
        "        8650.40,\n",
        "        22372.1125,\n",
        "        6055.05,\n",
        "        6739.055,\n",
        "        8050.40,\n",
        "        21014.225,\n",
        "        8650.40,\n",
        "        7539.055,\n",
        "        4539.72,\n",
        "        15024.325,\n",
        "    ],\n",
        "    index=pd.MultiIndex.from_product(\n",
        "        [FACILITIES, CUSTOMERS], names=[\"FACILITIES\", \"CUSTOMERS\"]\n",
        "    ),\n",
        "    columns=[\"variable_cost\"],\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### The Benders algorithm script in Python\n",
        "\n",
        "The Python snippet below begins the process by importing the required model file and configuring the necessary settings. \n",
        "Additionally, it loads data from the previously defined Python objects.\n",
        "\n",
        "For lists, we can directly assign data to AMPL sets using the syntax `ampl.set['<set name>'] = <Python list holding data>`. \n",
        "When dealing with parameters represented in Pandas dataframes, we employ the `set_data()` method. \n",
        "This method matches the dataframe's column names with corresponding entities in the AMPL model, ensuring proper assignment of values. \n",
        "It's important to ensure that the dataframe's indexing aligns with the AMPL model's indexing scheme.\n",
        "\n",
        "To define suffixes and problem structures, we utilize the `eval()` method.\n",
        "\n",
        "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. \n",
        "The control flow of the algorithm is managed through Python's control structures. \n",
        "The notable exception is the `nITER` parameter, which remains defined in AMPL to allow us to index constraints and price parameters adequately."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "\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 THE 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",
            "\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",
            "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",
            "\n",
            "SOLVING THE 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",
            "\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",
            "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",
            "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",
            "\n",
            "\n",
            "OPTIMAL SOLUTION FOUND\n",
            "Expected cost = 16758018.596250001\n",
            "\n",
            "facility_open [*] :=\n",
            "Baton_Rouge_LA  1\n",
            "    Baytown_TX  1\n",
            "   Beaumont_TX  1\n",
            ";\n",
            "\n"
          ]
        }
      ],
      "source": [
        "ampl = AMPL()  # Instantiate an AMPL object\n",
        "file_path = \"floc_bend.mod\"  # Define the file path\n",
        "ampl.read(file_path)  # Read the model from file\n",
        "\n",
        "# Set options\n",
        "ampl.option[\"solver\"] = \"highs\"  # Select the solver\n",
        "# Ensure that the solver uses the dual simplex method and returns dual rays that's why presolve is turned off:\n",
        "ampl.option[\"highs_options\"] = \"pre:solve off alg:simplex 1\"\n",
        "# Also turn off presolve in AMPL so we get info from the solver:\n",
        "ampl.option[\"presolve\"] = 0\n",
        "ampl.option[\"display_eps\"] = 0.000001  # Set the display precision\n",
        "ampl.option[\"omit_zero_rows\"] = 1  # Omit zero rows in the output\n",
        "\n",
        "# Load data\n",
        "# Load sets\n",
        "ampl.set[\"FACILITIES\"] = FACILITIES\n",
        "ampl.set[\"CUSTOMERS\"] = CUSTOMERS\n",
        "ampl.set[\"SCENARIOS\"] = SCENARIOS\n",
        "# Load parameters\n",
        "ampl.set_data(fixed_cost_df)\n",
        "ampl.set_data(variable_cost_df)\n",
        "ampl.set_data(customer_demand_df)\n",
        "ampl.set_data(facility_capacity_df)\n",
        "ampl.set_data(prob_df)\n",
        "\n",
        "# Set up AMPL and Python for Benders decomposition\n",
        "# Define suffix for storing dual rays\n",
        "ampl.eval(r\"suffix dunbdd;\")\n",
        "\n",
        "# Define Master and Sub problems\n",
        "ampl.eval(\n",
        "    r\"\"\"\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",
        ")\n",
        "\n",
        "# Loop logic variables\n",
        "epsilon = 0.00001\n",
        "n_iter = 0\n",
        "ampl.param[\"nITER\"] = n_iter\n",
        "violation_ctr = 0\n",
        "\n",
        "# Start Bender's loop\n",
        "while True:\n",
        "    violation_ctr = 0\n",
        "    n_iter += 1\n",
        "    ampl.param[\"nITER\"] = n_iter\n",
        "    print(f\"\\n\\nITERATION {n_iter}\\n\")\n",
        "\n",
        "    # Loop through the scenarios and solve each subproblem and update dual prices for each scenario\n",
        "    for s in SCENARIOS:\n",
        "        ampl.param[\"sub_scenario\"] = s\n",
        "        ampl.solve(problem=\"Sub\")\n",
        "        # Check for feasibility\n",
        "        if ampl.solve_result == \"infeasible\":\n",
        "            ampl.param[\"cut_type\"][n_iter, s] = \"feas\"\n",
        "            # Update cut_type, customer_price, and facility_price inside AMPL\n",
        "            ampl.eval(\n",
        "                r\"\"\"\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",
        "                \"\"\"\n",
        "            )\n",
        "            print(f\"{n_iter}: Feasibility cut added for scenario {s}\")\n",
        "        # Check for optimality\n",
        "        elif (\n",
        "            ampl.obj[\"operating_cost\"].value()\n",
        "            > ampl.var[\"sub_variable_cost\"][s].value() + epsilon\n",
        "        ):\n",
        "            ampl.param[\"cut_type\"][n_iter, s] = \"opt\"\n",
        "            # Update cut_type, customer_price, and facility_price inside AMPL\n",
        "            ampl.eval(\n",
        "                r\"\"\" \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",
        "                \"\"\"\n",
        "            )\n",
        "            print(f\"{n_iter}: Optimality cut added for scenario {s}\")\n",
        "        else:\n",
        "            violation_ctr += 1\n",
        "            ampl.param[\"cut_type\"][n_iter, s] = \"none\"\n",
        "    # If no cuts are added we can  break the loop as we have found the optimal solution\n",
        "    if violation_ctr == len(SCENARIOS):\n",
        "        break\n",
        "    # Solve the master problem\n",
        "    print(f\"\\nSOLVING THE MASTER PROBLEM\\n\")\n",
        "    ampl.solve(problem=\"Master\")\n",
        "    ampl.eval(r\"let {i in FACILITIES} sub_facility_open[i] := facility_open[i];\")\n",
        "\n",
        "\n",
        "print(f\"\\n\\nOPTIMAL SOLUTION FOUND\\nExpected cost = {ampl.obj['total_cost'].value()}\\n\")\n",
        "ampl.eval(r\"display facility_open;\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Conclusion\n",
        "\n",
        "Transitioning control flow of the Benders algorithm to Python has not altered the core solution of our problem. \n",
        "However, by shifting the control of the solve executions to Python code, we've significantly enhanced the flexibility of our algorithm's control mechanism. \n",
        "This newfound flexibility opens the door to sophisticated enhancements, such as parallelizing the scenario iterations within the main `while` loop. \n",
        "With the Benders algorithm's control flow now in Python, implementing such parallelization becomes far more feasible and straightforward."
      ]
    }
  ],
  "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
}