{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP2 - Modeling using docplex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. The `docplex` python package\n",
    "\n",
    "`DOcplex` is a python package developped by IBM — It provides easy-to-use API for IBM solvers Cplex and Cpoptimizer.\n",
    "\n",
    "DOcplex documentation for mathematical programming can be found here: http://ibmdecisionoptimization.github.io/docplex-doc/#mathematical-programming-modeling-for-python-using-docplex-mp-docplex-mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Solving TSP using `docplex`\n",
    "\n",
    "### 2.1. TSP model using `docplex`\n",
    "\n",
    "**Exercice:** Using `docplex`, create a model for the travelling salesman problem using the MTZ or Flow formulation and compare them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from docplex.mp.model import Model\n",
    "import tsp.data as data\n",
    "\n",
    "distances = data.grid42\n",
    "N = len(distances)\n",
    "\n",
    "tsp = Model(\"TSP\")\n",
    "tsp.log_output = True\n",
    "\n",
    "...  # TODO\n",
    "\n",
    "solution = tsp.solve()\n",
    "print(\"z* =\", solution.objective_value)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The largest set of distances contains 42 nodes, and should be easily solved by `docplex`.\n",
    "\n",
    "### 2.2. Generating random TSP instances\n",
    "\n",
    "**Question:** What method could you implement to generate a realistic set of distances for $n$ customers?\n",
    "\n",
    "**Exercice:** Implement the method proposed above and test it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.spatial.distance\n",
    "\n",
    "\n",
    "def generate_distances(n: int):\n",
    "    ...  # TODO\n",
    "\n",
    "\n",
    "from docplex.mp.model import Model\n",
    "\n",
    "distances = generate_distances(50)\n",
    "print(distances)\n",
    "\n",
    "N = len(distances)\n",
    "\n",
    "tsp = Model(\"TSP\")\n",
    "tsp.log_output = True\n",
    "\n",
    "...  # TODO: Copy your model from the first question here.\n",
    "\n",
    "solution = tsp.solve()\n",
    "print(\"z* =\", solution.objective_value)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Solving Warehouse Allocation using Benders decomposition with `docplex`\n",
    "\n",
    "### 3.1. The warehouse problem\n",
    "\n",
    "A company needs to supply a set of $n$ clients and needs to open new warehouses (from a\n",
    "set of $m$ possible warehouses).\n",
    "Opening a warehouse $j$ costs $f_j$ and supplying a client $i$ from a warehouse $j$ costs $c_{ij}$ per supply unit.\n",
    "Which warehouses should be opened in order to satisfy all clients while minimizing the total cost?\n",
    "\n",
    "### 3.2. Solving the warehouse problem with a single ILP\n",
    "\n",
    "- $y_{j} = 1$ if and only if warehouse $j$ is opened.\n",
    "- $x_{ij}$ is the fraction supplied from warehouse $j$ to customer $i$.\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "  \\text{min.} \\quad & \\sum_{i=1}^{n} \\sum_{j=1}^{m} c_{ij} x_{ij} + \\sum_{j=1}^{m} f_{j} y_{j} & \\\\\n",
    "  \\text{s.t.} \\quad & \\sum_{j=1}^{m} x_{ij} = 1, & \\forall i \\in\\{1,\\ldots,n\\}\\\\\n",
    "                    & x_{ij} \\leq y_{j}, & \\forall i\\in\\{1,\\ldots,n\\},\\forall j\\in\\{1,\\ldots,m\\}\\\\\n",
    "                    & y_{j} \\in \\left\\{0,~1\\right\\}, & \\forall j \\in \\left\\{1,~\\ldots,~m\\right\\}\\\\\n",
    "                    & x_{ij} \\geq 0, & \\forall i \\in \\left\\{1,~\\ldots,~n\\right\\}, \\forall j \\in \\left\\{1,~\\ldots,~m\\right\\}\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "\n",
    "**Exercice:** Implement the ILP model for the warehouse allocation problem and test it on the given instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from docplex.mp.model import Model\n",
    "\n",
    "# We will start with a small instances with 3 warehouses and 3 clients:\n",
    "N = 3\n",
    "M = 3\n",
    "\n",
    "# Opening and distribution costs:\n",
    "f = [20, 20, 20]\n",
    "c = [[15, 1, 2], [1, 16, 3], [4, 1, 17]]\n",
    "\n",
    "wa = Model(\"Warehouse Allocation\")\n",
    "wa.log_output = True\n",
    "\n",
    "...  # TODO: Model for the warehouse allocation.\n",
    "\n",
    "solution = wa.solve()\n",
    "print(\"z* =\", solution.objective_value)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Benders' decomposition for the Warehouse Allocation problem\n",
    "\n",
    "We are going to use Benders' decomposition to solve the Warehouse Allocation problem. \n",
    "\n",
    "#### Dual subproblem\n",
    "\n",
    "$\n",
    "\\begin{align*}\n",
    "\\text{max.} \\quad & \\sum_{i=1}^{n} v_{i} - \\sum_{i=1}^{n}\\sum_{j=1}^{m} \\bar{y}_{j} u_{ij} & \\\\\n",
    "\\text{s.t.} \\quad & v_{i} - u_{ij} \\leq c_{ij}, & \\forall i\\in\\{1,\\ldots,n\\},\\forall j\\in\\{1,\\ldots,m\\}\\\\\n",
    "                  & v_{i} \\in\\mathbb{R},\\ u_{ij} \\geq 0 & \\forall i \\in\\{1,\\ldots,n\\}, \\forall j\\in\\{1,\\ldots,m\\}\n",
    "\\end{align*}\n",
    "$\n",
    "\n",
    "#### Master problem\n",
    "\n",
    "$\n",
    "\\begin{align*}\n",
    "  \\text{min.} \\quad & \\sum_{j=1}^{m} f_j y_j + z & \\\\\n",
    "  \\text{s.t.} \\quad & z \\geq \\sum_{i=1}^{n}v_i^p - \\sum_{i=1}^{n} \\sum_{j=1}^{m} u_{ij}^p y_j, & \\forall p\\in l_1\\\\\n",
    "                    & 0 \\geq \\sum_{i=1}^{n}v_i^r - \\sum_{i=1}^{n} \\sum_{j=1}^{n} u_{ij}^r y_j, & \\forall r\\in l_2\\\\\n",
    "                    & y_{j} \\in\\{0,1\\}, & \\forall j\\in\\{1,\\ldots,m\\}\n",
    "\\end{align*}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice:** Implement the method `create_master_problem` that creates the initial master problem (without feasibility or optimality constraints) for the warehouse allocation problem.\n",
    "\n",
    "<div class=\"alert alert-info alert-block\">\n",
    "\n",
    "You can use `print(m.export_as_lp_string())` to display a textual representation of a `docplex` model `m`.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from docplex.mp.model import Model\n",
    "from docplex.mp.linear import Var\n",
    "from docplex.mp.constr import AbstractConstraint\n",
    "from typing import List, Sequence, Tuple\n",
    "\n",
    "\n",
    "def create_master_problem(\n",
    "    N: int, M: int, f: Sequence[float], c: Sequence[Sequence[float]]\n",
    ") -> Tuple[Model, Var, Sequence[Var]]:\n",
    "    \"\"\"\n",
    "    Creates the initial Benders master problem for the Warehouse Allocation problem.\n",
    "\n",
    "    Args:\n",
    "        N: Number of clients.\n",
    "        M: Number of warehouses.\n",
    "        f: Array-like containing costs of opening warehouses.\n",
    "        c: 2D-array like containing transport costs from client to warehouses.\n",
    "\n",
    "    Returns:\n",
    "        A 3-tuple containing the docplex problem, the z variable and the y variables.\n",
    "    \"\"\"\n",
    "\n",
    "    wa = Model(\"Warehouse Allocation - Benders master problem\")\n",
    "\n",
    "    ...  # TODO\n",
    "\n",
    "    return wa, z, y\n",
    "\n",
    "\n",
    "# Check your method:\n",
    "wa, z, y = create_master_problem(N, M, f, c)\n",
    "print(wa.export_as_lp_string())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice:** Implement the method `add_optimality_constraints` that add optimality constraints to the Benders master problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_optimality_constraint(\n",
    "    N: int,\n",
    "    M: int,\n",
    "    wa: Model,\n",
    "    z: Var,\n",
    "    y: Sequence[Var],\n",
    "    v: Sequence[float],\n",
    "    u: Sequence[Sequence[float]],\n",
    ") -> List[AbstractConstraint]:\n",
    "    \"\"\"\n",
    "    Adds an optimality constraints to the given Warehouse Allocation model\n",
    "    using the given optimal values from the Benders dual subproblem.\n",
    "\n",
    "    Args:\n",
    "        N: Number of clients.\n",
    "        M: Number of warehouses.\n",
    "        wa: The Benders master problem (docplex.mp.model.Model).\n",
    "        z: The z variable of the master problem.\n",
    "        y: The y variables of the master problem.\n",
    "        v: The optimal values for the v variables of the Benders dual subproblem.\n",
    "        u: The optimal values for the u variables of the Benders dual subproblem.\n",
    "\n",
    "    Return: The optimality constraint added.\n",
    "    \"\"\"\n",
    "    return ...  # TODO\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice:** Implement the method `add_feasibility_constraints` that add feasibility constraints to the Benders master problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_feasibility_constraints(\n",
    "    N: int,\n",
    "    M: int,\n",
    "    wa: Model,\n",
    "    z: Var,\n",
    "    y: Sequence[Var],\n",
    "    v: Sequence[float],\n",
    "    u: Sequence[Sequence[float]],\n",
    ") -> List[AbstractConstraint]:\n",
    "    \"\"\"\n",
    "    Adds an optimality constraints to the given Warehouse Allocation model\n",
    "    using the given optimal values from the Benders dual subproblem.\n",
    "\n",
    "    Args:\n",
    "      - N: Number of clients.\n",
    "      - M: Number of warehouses.\n",
    "      - wa: The Benders master problem (docplex.mp.model.Model).\n",
    "      - z: The z variable of the master problem.\n",
    "      - y: The y variables of the master problem.\n",
    "      - v: The extreme rays for the v variables of the Benders dual subproblem.\n",
    "      - u: The extreme rays for the u variables of the Benders dual subproblem.\n",
    "\n",
    "    Returns:\n",
    "        The feasibility constraint added.\n",
    "    \"\"\"\n",
    "    # TODO:\n",
    "    return\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice:** Implement the method `create_dual_subproblem` that, given a solution `y` of the master problem, create the corresponding Benders dual subproblem.\n",
    "\n",
    "$\n",
    "\\begin{align*}\n",
    "\\text{max.} \\quad & \\sum_{i=1}^{n} v_{i} - \\sum_{i=1}^{n}\\sum_{j=1}^{m} \\bar{y}_{j} u_{ij} & \\\\\n",
    "\\text{s.t.} \\quad & v_{i} - u_{ij} \\leq c_{ij}, & \\forall i\\in\\{1,\\ldots,n\\},\\forall j\\in\\{1,\\ldots,m\\}\\\\\n",
    "                  & v_{i} \\in\\mathbb{R},\\ u_{ij} \\geq 0 & \\forall i \\in\\{1,\\ldots,n\\}, \\forall j\\in\\{1,\\ldots,m\\}\n",
    "\\end{align*}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from docplex.mp.model import Model\n",
    "\n",
    "\n",
    "def create_dual_subproblem(\n",
    "    N: int, M: int, f: Sequence[float], c: Sequence[Sequence[float]], y: Sequence[int]\n",
    ") -> Tuple[Model, Sequence[Var], Sequence[Sequence[Var]]]:\n",
    "    \"\"\"\n",
    "    Creates a Benders dual subproblem for the Warehouse Allocation problem corresponding\n",
    "    to the given master solution.\n",
    "\n",
    "    Args:\n",
    "        N: Number of clients.\n",
    "        M: Number of warehouses.\n",
    "        f: Array-like containing costs of opening warehouses.\n",
    "        c: 2D-array like containing transport costs from client to warehouses.\n",
    "        y: Values of the y variables from the Benders master problem.\n",
    "\n",
    "    Returns:\n",
    "        A 3-tuple containing the docplex problem, the v variable and the u variables.\n",
    "    \"\"\"\n",
    "\n",
    "    dsp = Model(\"Warehouse Allocation - Benders dual subproblem\")\n",
    "\n",
    "    # We disable pre-solve to be able to retrieve a meaningful status in the main\n",
    "    # algorithm:\n",
    "    dsp.parameters.preprocessing.presolve.set(0)\n",
    "\n",
    "    ...  # TODO\n",
    "\n",
    "    return dsp, v, u\n",
    "\n",
    "\n",
    "# Check your method (assuming y = [1 1 1 ... 1]):\n",
    "dsp, v, u = create_dual_subproblem(N, M, f, c, [1] * M)\n",
    "print(dsp.export_as_lp_string())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice:** Using the methods you implemented, write the Benders decomposition algorithm for the warehouse allocation problem.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "The `get_extreme_rays` function can be used to retrieve the extreme rays associated with an unbounded solution of the dual subproblem.\n",
    "    \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    \n",
    "You can use `model.get_solve_status()` to obtain the status of the resolution and compare it to members of `JobSolveStatus`:\n",
    "    \n",
    "```python\n",
    "if model.get_solve_status() == JobSolveStatus.OPTIMAL_SOLUTION:\n",
    "    pass\n",
    "```\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from docplex.mp.model import Model\n",
    "from docplex.util.status import JobSolveStatus\n",
    "\n",
    "\n",
    "def get_extreme_rays(\n",
    "    N: int, M: int, model: Model, v: Sequence[Var], u: Sequence[Sequence[Var]]\n",
    ") -> Tuple[Sequence[float], Sequence[Sequence[float]]]:\n",
    "    \"\"\"\n",
    "    Retrieves the extreme rays associated to the dual subproblem.\n",
    "\n",
    "    Args:\n",
    "        N: Number of clients.\n",
    "        M: Number of warehouses.\n",
    "        model: The Benders dual subproblem model (docplex.mp.model.Model).\n",
    "        v: 1D array containing the v variables of the subproblem.\n",
    "        u: Either a 2D array of a tuple-index dictionary containing the u variables\n",
    "            of the subproblem.\n",
    "\n",
    "    Returns:\n",
    "        A 2-tuple containing the list of extreme rays correspondig to v,\n",
    "        and the 2D-list of extreme rays corresponding to u.\n",
    "    \"\"\"\n",
    "    ray = model.get_engine().get_cplex().solution.advanced.get_ray()\n",
    "\n",
    "    if isinstance(u, dict):\n",
    "\n",
    "        def get_uij(i, j):\n",
    "            return u[i, j]\n",
    "\n",
    "    else:\n",
    "\n",
    "        def get_uij(i, j):\n",
    "            return u[i][j]\n",
    "\n",
    "    return (\n",
    "        [ray[v[i].index] for i in range(N)],\n",
    "        [[ray[get_uij(i, j).index] for j in range(M)] for i in range(N)],\n",
    "    )\n",
    "\n",
    "\n",
    "# We will start with a small instances with 3 warehouses and 3 clients:\n",
    "N = 3\n",
    "M = 3\n",
    "\n",
    "# Opening and distribution costs:\n",
    "f = [20, 20, 20]\n",
    "c = [[15, 1, 2], [1, 16, 3], [4, 1, 17]]\n",
    "\n",
    "# We stop iterating if the new solution is less than epsilon\n",
    "# better than the previous one:\n",
    "epsilon = 1e-6\n",
    "\n",
    "wa, z, y = create_master_problem(N, M, f, c)\n",
    "\n",
    "n = 0\n",
    "while True:\n",
    "\n",
    "    # Print iteration:\n",
    "    n = n + 1\n",
    "    print(\"Iteration {}\".format(n))\n",
    "\n",
    "    ...  # TODO\n",
    "\n",
    "print(\"Done.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4. Generating instances for the Warehouse Allocation problem\n",
    "\n",
    "**Exercice:** Using the TSP instances contained in `tsp.data` or the `generate_distances` method, create instances for the warehouse allocation problem with randomized opening costs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\"></div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}