dp01bis_tbc_intro_with_howard_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Special lecture 2: dynamic programming</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass on optimal transport and economic applications</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018-2022 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed [here](http://www.math-econ-code.org/team).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass on optimal transport and economic applications, January 2022. https://github.com/math-econ-code/mec_optim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "* Basics of (finite-horizon, discrete) dynamic programming: Bellman's equation; forward induction, backward induction\n",
    "\n",
    "* Markov decision processes\n",
    "\n",
    "* Dynamic programming as linear programming: interpretation of duality\n",
    "\n",
    "* Vectorization, Kronecker products, multidimensional arrays\n",
    "\n",
    "### References\n",
    "\n",
    "* Ford Jr, L. R., & Fulkerson, D. R. (1958). Constructing maximal dynamic flows from static flows. *Operations research*, 6(3), 419-433.\n",
    "\n",
    "* Howard, R. (1960). *Dynamic programming and Markov processes*. Wiley.\n",
    "\n",
    "* Schrijver, A. (2003). *Combinatorial optimization: polyhedra and efficiency* Vol. A. Springer. Section 12.5.c.\n",
    "\n",
    "* Bertsekas, D. (2011), *Dynamic Programming and Optimal Control*, Vols. I and II. 3rd ed. Athena.\n",
    "\n",
    "* Ljungqvist, L., & Sargent, T. (2012), *Recursive Macroeconomic Theory* 3rd ed. MIT.\n",
    "\n",
    "* Rust (1987), Optimal replacement of GMC bus engines: an empirical model of Harold Zurcher. *Econometrica*.\n",
    "\n",
    "* Nazareth, J., & Kulkarni, R. (1986). Linear programming formulations of Markov decision processes. *Operations Research Letters*.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the libraries we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as spr\n",
    "from tabulate import tabulate\n",
    "# !python -m pip install -i https://pypi.gurobi.com gurobipy ## only if Gurobi not here\n",
    "import gurobipy as grb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Movitation\n",
    "\n",
    "Howard (1960, chapter 5) describes the following problem.\n",
    "\n",
    "``The examples of the policy-iteration method presented up to this point have been somewhat far removed from the realm of practical problems. It would be extremely interesting to see the method applied to a problem that is of major importance to industry. As an example of such a practical application, the replacement problem was chosen. This is the problem of when to replace a piece of capital equipment that deteriorates with time. The question to be answered is this: If we now own a machine of a certain age, should we keep it or should we trade it in; further, if we trade it in, how old a machine should we buy?\n",
    "\n",
    "In order to fix ideas, let us consider the problem of automobile replacement over a time interval of ten years. We agree to review our current situation every three months and to make a decision on keeping our present car or trading it in at that time. The state of the system, x, is described by the age of the car in three-month periods; $x$ may run from 0 to 39. In order to keep the number of states finite, a car of age 39 remains a car of age 39 forever (it is considered to be essentially worn out). The alternatives available in each state are these: The first alternative, y = 0, is to keep the present car for another quarter. The other alternatives, y > 0, are to buy a car of age y - 1, where y - 1 may be as large as 39. We have then 40 states with 41 alternatives in each state, with the result that there are 41^40 possible policies.``\n",
    "\n",
    "(we have replaced the notations by $x$ and $y$ the notations $i$ and $k$ from the original text for consistency with the rest of these notes, and have changed the ranges to start from 0 as is the convention in Python)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{aligned}\n",
    "C_x &: \\text{the cost of buying a car of age } x \\\\\n",
    "T_x &: \\text{the trade-in value of a car of age } x \\\\\n",
    "E_x &: \\text{the expected cost of operating a car of age } x \\text{ until it reaches age } x + 1 \\\\\n",
    "p_x &: \\text{the probability that a car of age } x \\text{ will survive to be age } x + 1 \\text{ without incurring a prohibitively expensive repair.}\n",
    "\\end{aligned}\n",
    "\n",
    "The probability defined here is necessary to limit the number of states. A car of any age that has a hopeless breakdown is immediately sent to state 39. Naturally, $ p_{39} = 0 $.\n",
    "\n",
    "The basic equations governing the system when it is in state $ x $ are the following:\n",
    "\n",
    "If $y = 0$ (keep present car),\n",
    "$$\n",
    "g + v_x = - E_x + p_{x}v_{x+1} + (1 - p_{x})v_{39}\n",
    "$$\n",
    "\n",
    "If $y > 0$ (trade for car of age $y - 1$),\n",
    "$$\n",
    "g + v_x = T_x - C_{y-1} - E_{y-1} + p_{y-1}v_{y-1} + (1 - p_{y-1})v_{39}\n",
    "$$\n",
    "\n",
    "Phrased in terms of earlier notation, the equations are:\n",
    "\n",
    "$q_x^y = -E_x$ for $y = 0$, and  $q_x^y = T_x - C_{y-1} - E_{y-1}$\n",
    "for $y > 0$\n",
    "\n",
    "For the probabilities, for $y = 0$:\n",
    "$$\n",
    "P_{x^\\prime|xy} = \n",
    "\\begin{cases} \n",
    "p_x & x^\\prime = x + 1 \\\\\n",
    "1 - p_x & x^\\prime = 39 \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "and for $y > 0$, $y<39$:\n",
    "$$\n",
    "P_{x^\\prime|xy} = \n",
    "\\begin{cases} \n",
    "p_{y-1} & x^\\prime = y \\\\\n",
    "1 - p_{y-1} & x^\\prime = 39 \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "and for $y \\geq 39$:\n",
    "$$\n",
    "P_{x^\\prime|xy} = \n",
    "\\begin{cases} \n",
    "1  & x^\\prime = 39 \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def howard_ex_data():\n",
    "    C_x = np.array([2000, 1840, 1680, 1560, 1300, 1220, 1150, 1080, 900, 840,\n",
    "                   780, 730, 600, 560, 520, 480, 440, 420, 400, 380, 360, 345,\n",
    "                   330, 315, 300, 290, 280, 265, 250, 240, 230, 220, 210, 200,\n",
    "                   190, 180, 170, 160, 150, 140, 130])\n",
    "    \n",
    "    T_x = np.array([1600, 1460, 1340, 1230, 1050, 980, 910, 840, 710, 650, 600,\n",
    "                   550, 480, 430, 390, 360, 330, 310, 290, 270, 255, 240, 225,\n",
    "                   210, 200, 190, 180, 170, 160, 150, 145, 140, 135, 130, 120,\n",
    "                   115, 110, 105, 95, 87, 80])\n",
    "\n",
    "\n",
    "    E_x = np.array([50, 53, 56, 59, 62, 65, 68, 71, 75, 78, 81, 84, 87, 90, 93,\n",
    "                   96, 100, 103, 106, 109, 112, 115, 118, 121, 125, 129, 133,\n",
    "                   137, 141, 145, 150, 155, 160, 167, 175, 182, 190, 205, 220,\n",
    "                   235, 250])\n",
    "    \n",
    "    p_x = np.array([1, 0.999, 0.998, 0.997, 0.996, 0.994, 0.991, 0.988, 0.985, 0.983,\n",
    "                   0.98, 0.975, 0.97, 0.965, 0.96, 0.955, 0.95, 0.945, 0.94, 0.935,\n",
    "                   0.93, 0.925, 0.919, 0.91, 0.9, 0.89, 0.88, 0.865, 0.85, 0.82,\n",
    "                   0.79, 0.76, 0.73, 0.66, 0.59, 0.51, 0.43, 0.3, 0.2, 0.1, 0])\n",
    "\n",
    "\n",
    "                    \n",
    "    return (C_x[1:],T_x[1:],E_x[1:],p_x[1:])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(40, 40, 40, 40)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_x,T_x,E_x,p_x = howard_ex_data()\n",
    "len(p_x),len(E_x),len(T_x),len(C_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(730, 550, 84, 0.975)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = 11\n",
    "C_x[x],T_x[x],E_x[x],p_x[x],"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "def howard_ex_objects():\n",
    "    C_x,T_x,E_x,p_x = howard_ex_data()\n",
    "    nbx,nby = 40,40,41\n",
    "    P_xp_x_y = np.zeros((nbx,nbx,nby))\n",
    "    P_xp_x_y[-1,:,0] = 1 - p_x\n",
    "    P_xp_x_y[:,:,0] += np.diag(p_x[:-1],k=-1)\n",
    "    for y in range(1,nby-2):\n",
    "        P_xp_x_y[y,:,y] = p_x[y-1]\n",
    "        P_xp_x_y[-1,:,y] = 1 - p_x[y-1]\n",
    "    \n",
    "    for y in range(nby-2,nby):\n",
    "        P_xp_x_y[-1,:,y] = 1\n",
    "    \n",
    "    \n",
    "    u_x_y = np.zeros((nbx,nby))\n",
    "    u_x_y[:,0] = -E_x\n",
    "    u_x_y[:,1:]= T_x[:,None] - C_x[None,:-1]-E_x[None,:-1]\n",
    "    return P_xp_x_y\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 0.],\n",
       "       [1., 1., 1., ..., 1., 1., 1.]])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tst =howard_ex_objects()\n",
    "tst[:,:,40]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.   , 0.   , 0.   , ..., 0.   , 0.   , 0.   ],\n",
       "       [0.999, 0.   , 0.   , ..., 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.998, 0.   , ..., 0.   , 0.   , 0.   ],\n",
       "       ...,\n",
       "       [0.   , 0.   , 0.   , ..., 0.   , 0.   , 0.   ],\n",
       "       [0.   , 0.   , 0.   , ..., 0.2  , 0.   , 0.   ],\n",
       "       [0.001, 0.002, 0.003, ..., 0.8  , 1.   , 1.   ]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tst[:,:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Mdp:\n",
    "    def __init__(self,)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear dynamic programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## States and actions\n",
    "\n",
    "\n",
    "Consider a finite set of states $x\\in\\mathcal{X}$; and a set of possible actions $y\\in\\mathcal{Y}$.\n",
    "\n",
    "It is assumed that at initial time, the mass of agents in state $x$ is $q_{x}$ (exogenous). \n",
    "\n",
    "Let $\\pi_{xy}^{t}$ be the (endogenous) number of individuals who are in state $x$ and choose $y$ (\"policy variable\").\n",
    "\n",
    "Because $\\pi_{xy}^{t}$ is our decision variable, we shall need to worry about how the computer represents it in the memory: by stacking columns or by stacking rows? this brings us to the important question of *vectorization*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bus example**. In our bus maintenance example, one faces a maintenance decision, which is captured by $y\\in\\mathcal{Y}=\\left\\{0,1\\right\\}$, where $y = 0$ is to keep going, and $y=1$ is to perform overhaul. The state $x\\in\\mathcal{X}=\\left\\{x_{0},...,x_{I}\\right\\}$ represents the mileage level of a bus. The transition between states is as follows:\n",
    "* When no overhaul is performed, these states undergo random transitions (depending on how much the bus is used): $x\\rightarrow x^{\\prime}$ with some probability $P_{x^{\\prime}|x}$, where state $x^{\\prime}$ marks more mileage than $x$.\n",
    "* When overhaul is performed on a bus, the state is restored to the zero state $x_{0}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vectorization and Kronecker products\n",
    "\n",
    "We will need to represent matrices (such as $U_{x}^{t}$) and 3-dimensional arrays (such as $u_{xy}^{t}$). Under the row-major order (a.k.a. 'C order') used in `C` and by default in `numpy`, we will represent a matrix $M_{ij}$ by varying the last index first, i.e. a $2\\times2$ matrix will be represented as $vec_C\\left(M\\right) = M_{11}, M_{12}, M_{21}, M_{22}.$ Likewise, a 2x2x2 3-dimensional array $A$ will be represented by varying the first index first, then the second, i.e.\n",
    "\n",
    "$vec_C\\left(A\\right) = A_{111}, A_{112}, A_{121}, A_{122}, A_{211}, A_{212}, A_{221}, A_{222}$.\n",
    "\n",
    "In `numpy`, this is implemented by `reshape(...)`.\n",
    "\n",
    "A very important identity is\n",
    "\\begin{align*}\n",
    "vec_C\\left(AXB\\right) = \\left(  A\\otimes B^\\top\\right)  vec_C\\left(X\\right),\n",
    "\\end{align*}\n",
    "where $vec_C$ is the vectorization under the C (row-major) order, and where the Kronecker product $\\otimes$ is defined as follows for 2x2 matrices (with obvious generalization):\n",
    "\n",
    "\\begin{align*}\n",
    "A\\otimes B=\n",
    "\\begin{pmatrix}\n",
    "a_{11}B & a_{12}B\\\\\n",
    "a_{21}B & a_{22}B\n",
    "\\end{pmatrix}.\n",
    "\\end{align*}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Representation of the decision variable\n",
    "\n",
    " We shall adopt the convention to store $\\pi^t_{xy}$ by taking the indexing ordering to be $(t,x,y)$ under the row-major order. That way we shall store $(\\pi^1_{xy})_{xy}$ on top of  $(\\pi^2_{xy})_{xy}$, etc, and whithin the $(\\pi^1_{xy})_{xy}$ block, $(\\pi^1_{x_1y})_{y}$ on top of $(\\pi^1_{x_2y})_{y}$, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choice-counting equation\n",
    "\n",
    "Define $n_{x}^{t}$ be the (endogenous) number of individuals in state $x$ at time $t$. The choice-counting equation expresses that the sum over $y$ of individuals of type $x$ who made choice $y$ is equal to $n_x$, that is:\n",
    "\n",
    "$\n",
    "n_{x}^{t} = \\sum_{y\\in\\mathcal{Y}}\\pi_{xy}^{t}.\n",
    "$\n",
    "\n",
    "Clearly, this can we rewritten for each $t$\n",
    "\n",
    "$\n",
    "n^t = \\begin{pmatrix}\n",
    "1_{Y}^{\\top } & 0 & \\cdots  & 0 \\\\ \n",
    "0 & 1_{Y}^{\\top } & \\ddots  & \\vdots  \\\\ \n",
    "\\vdots  & \\ddots  & \\ddots  & 0 \\\\ \n",
    "0 & \\cdots  & 0 & 1_{Y}^{\\top }\n",
    "\\end{pmatrix} vec_C(\\pi^{t}) =\\left( I_{X}\\otimes 1_{Y}^{\\top } \\right) vec_C(\\pi^{t}),\n",
    "$\n",
    "\n",
    "and stacking this over $t$ yields\n",
    "\n",
    "$\n",
    "n = \\begin{pmatrix}\n",
    "I_{X}\\otimes 1_{Y}^{\\top } & 0 & \\cdots  & 0 \\\\ \n",
    "0 & I_{X}\\otimes 1_{Y}^{\\top } & \\ddots  & \\vdots  \\\\ \n",
    "\\vdots  & \\ddots  & \\ddots  & 0 \\\\ \n",
    "0 & \\cdots  & 0 & I_{X}\\otimes 1_{Y}^{\\top }%\n",
    "\\end{pmatrix}%\n",
    "vec_C(\\pi) = \\left( I_{T}\\otimes I_{X}\\otimes 1_{Y}^{\\top } \\right) vec_C(\\pi).  \n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore the choice-counting equation reads:\n",
    "\n",
    "\\begin{equation}\n",
    "\\left( I_{T}\\otimes I_{X}\\otimes 1_{Y}^{\\top } \\right) vec_C(\\pi)  = n.\n",
    "\\end{equation}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov transitions\n",
    "\n",
    "Now assume that if an agent is in state $x$ and takes decision $y$ at time $t-1$, then the agent's state will transition to state $x'$ at time $t$ with probability $P_{x^{\\prime}|xy}$, which expresses that among the individual in state $x$ who choose $y$ at time $t-1$, a fraction $P_{x^{\\prime}|xy}$ shall transit to state $x^{\\prime}$ at time $t$.\n",
    "\n",
    "$P$ is represented by a matrix whose rows are indexed by $x^{\\prime}$ and whose columns are indexed by $xy$. Under the row-major order, the index ordering will be $(x^{\\prime},x,y)$.\n",
    "\n",
    "For $t$ such that $1\\leq t\\leq T-1$, the Markov transitions are given by\n",
    "\n",
    "$\n",
    "n_{x^{\\prime}}^{t} = \\sum_{x\\in\\mathcal{X},~y\\in\\mathcal{Y}}P_{x^{\\prime}|xy}\\pi_{xy}^{t-1}\n",
    "$, that is\n",
    "\n",
    "$\n",
    "n^t = P \\pi^{t-1}\n",
    "$, and  $n^{1}=q$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bus example**. In the bus example, assume that the Markov transitions laws are given by:\n",
    "\n",
    "* If no overhaul is performed, each state but the last one has a probability $25\\%$ of transiting to the next one, and probability $75\\%$ of remaining the same. The last state transits to $x=0$ with probability $25\\%$ (overhaul is performed when beyond last state).\n",
    "\n",
    "* If overhaul is performed, the state transits to $0$ for sure.\n",
    "\n",
    "\n",
    "Define the matrix $L^{\\mathcal{X}}_{x^{\\prime}|x}$ associated with a transition to the next state, and a reset at the last state, and $R^{\\mathcal{X}}_{x^{\\prime}|x}$ associated with a reset to the initial state.\n",
    "\n",
    "for the transitions $x\\to x^{\\prime}$ by:\n",
    "\\begin{align*}\n",
    "L^{\\mathcal{X}}=%\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & \\cdots & 0 & 1\\\\\n",
    "1 & 0 & \\ddots & \\ddots & 0\\\\\n",
    "0 & 1 & \\ddots & \\ddots & 0\\\\\n",
    "\\vdots & \\ddots & \\ddots & \\ddots & 0 \\\\\n",
    "0 & 0 & 0 & 1 & 0\n",
    "\\end{pmatrix}\n",
    "\\text{ and }\n",
    "R^{\\mathcal{X}}=%\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & \\cdots & 1\\\\\n",
    "0 & \\vdots & \\ddots  & \\vdots\\\\\n",
    "0 & \\vdots &  \\ddots & \\vdots\\\\\n",
    "0 & 0 & \\cdots & 0\n",
    "\\end{pmatrix}\n",
    "\\end{align*}\n",
    "Now:\n",
    "* If $y=0$, the transition matrix $x\\to x^{\\prime}$ will be given by $0.75I_{\\mathcal{X}}+0.25L_{\\mathcal{X}}$, while\n",
    "* If $y=1$ it will be given by $R^{\\mathcal{X}}$\n",
    "\n",
    "so that $P$ is given by\n",
    "\\begin{align*}\n",
    "P=\\left(  0.75I_{\\mathcal{X}}+0.25L^{\\mathcal{X}}\\right) \\otimes (1,0)\n",
    "+ R^{\\mathcal{X}} \\otimes (0,1)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To fix ideas, we shall assume that $\\mathcal{X}=\\left\\{0,1,2\\right\\}$, that $\\mathcal{Y}=\\left\\{0,1\\right\\}$ as argued before, and  $\\mathcal{T}=\\left\\{0,...,3\\right\\}$,, and we implement this as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class dynamic_problem:\n",
    "    def __init__(self,nbT,nbX):\n",
    "        self.nbT = nbT\n",
    "        self.nbX = nbX\n",
    "        self.nbY = 2\n",
    "\n",
    "\n",
    "small_example = dynamic_problem(4,3)\n",
    "        \n",
    "def P_x_xy(self):\n",
    "    LX = spr.csr_matrix((np.ones(self.nbX), (list(range(1,self.nbX))+[0], range(self.nbX))), shape = (self.nbX,self.nbX))\n",
    "    RX = spr.csr_matrix((np.ones(self.nbX), ([0 for i in range(self.nbX)], range(self.nbX))), shape = (self.nbX,self.nbX))\n",
    "    return spr.kron((0.75 * spr.identity(self.nbX) + 0.25 * LX), np.array([[1,0]])) + spr.kron(RX, np.array([[0,1]])) \n",
    "\n",
    "dynamic_problem.P_x_xy = P_x_xy\n",
    "\n",
    "small_example.P_x_xy().toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov forward equation\n",
    "\n",
    "Denoting $N_{T}$ the $T \\times T$ matrix defined by:\n",
    "\n",
    "$\n",
    "N_{T}=\n",
    "\\begin{pmatrix}\n",
    "0  & 0 & \\cdots & \\cdots & 0\\\\\n",
    "1  & \\ddots & \\cdots & \\cdots & 0\\\\\n",
    "0 & \\ddots & \\ddots &  & \\vdots\\\\\n",
    "\\vdots  & \\ddots & \\ddots & \\ddots & \\vdots\\\\\n",
    "0  & \\cdots & 0 & 1 & 0\n",
    "\\end{pmatrix}\n",
    "$\n",
    "\n",
    "this rewrites\n",
    "\n",
    "$\n",
    "n = (N_{T} \\otimes P) \\pi + b \n",
    "$\n",
    "\n",
    "where $b^t_x$ is the vector obtained by stacking column vector $q_x$ on top of $(T-1)X$ zeros.\n",
    "\n",
    "Therefore the Markov forward equation reads:\n",
    "\n",
    "\\begin{equation}\n",
    "\\left(N_{T} \\otimes P\\right) vec_C(\\pi) + b = n.\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def N_t_t(self):\n",
    "    return(spr.csr_matrix((np.ones(self.nbT-1), (list(range(1,self.nbT)), range(self.nbT-1))), shape = (self.nbT,self.nbT)))\n",
    "\n",
    "dynamic_problem.N_t_t = N_t_t\n",
    "\n",
    "small_example.N_t_t().toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Population dynamics\n",
    "\n",
    "Combining the choice-counting equation with the Markov forward equation to substitute out $n$ yields\n",
    "\n",
    "$$\n",
    "A ~ vec_C(\\pi)  = b.\n",
    "$$\n",
    "\n",
    "where we have defined \n",
    "\n",
    "$$\n",
    "A :=  I_{T}\\otimes I_{X}\\otimes 1_{Y}^{\\top }   - N_{T} \\otimes P . \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implementation is straightforward:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def A_tx_txy(self):\n",
    "    return(spr.kron(spr.identity(self.nbT),spr.kron(spr.identity(self.nbX), np.ones(self.nbY))) - spr.kron(N_t_t(self), P_x_xy(self)) )\n",
    "\n",
    "dynamic_problem.A_tx_txy = A_tx_txy\n",
    "\n",
    "def b_tx(self,q=np.array([])): \n",
    "    if len(q)==0:\n",
    "        q = np.ones(self.nbX)\n",
    "    return( np.concatenate((q,np.zeros(self.nbX * (self.nbT-1 )))))\n",
    "\n",
    "dynamic_problem.b_tx = b_tx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"b=\",small_example.b_tx())\n",
    "small_example.A_tx_txy()[0:5,0:8].toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Payoffs\n",
    "\n",
    "Let $u_{xy}^{t}$ be the (exogenous) payoff associated with choice $y\\in\\mathcal{Y}$ at time $t\\in\\mathcal{T}=\\left\\{  1,...,T\\right\\}  $ in state $x\\in\\mathcal{X}$, discounted in order to be expressed in zeroth-period equivalent (for example: $u_{xy}^{t}=\\beta^{t}u_{xy}$ where $\\beta$ is a constant discount factor).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example**. Back to our bus example, we assume there is a fixed cost $C$ associated with overhaul (independent of mileage), while operations costs $c\\left(  x\\right)$ increase with mileage (maintenance, fuel, insurance and costs of unexpected breakdowns). Specifically, assume the following:\n",
    "* The cost of replacing an engine is $C=\\$8,000$ (in $1985$ dollars).\n",
    "\n",
    "* The operations cost is $c\\left(  x\\right)  =  \\frac {C x} {\\left| \\mathcal{X} \\right|} .$\n",
    "\n",
    "* The discount factor is $\\beta=0.9$.\n",
    "\n",
    "Next, we build $u_{xyt}$ \n",
    "* First the $u_{xy}$'s so that $u_{x0}=-x\\times5.10^{2}$ for $x<\\bar{x}$, and $u_{\\bar{x}0}=-C$, while $u_{x1}=-C$ for all $x$.\n",
    "\n",
    "* Next the $u_{xyt}$ so that $u_{xyt}=u_{xy}\\beta^{t}=vec_F\\left(\\left(\\beta^{1},...,\\beta^{T}\\right)  \\otimes u_{xy}\\right)$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u_txy(self, overhaulCost = 8000 , beta = 0.9):\n",
    "    maintCost = lambda x: x * overhaulCost / self.nbX\n",
    "    discount = lambda x: beta**x\n",
    "    u_xy = np.array(list(zip([-maintCost(x) for x in range(self.nbX)],[-overhaulCost for x in range(self.nbX)] ))).reshape(1,-1)\n",
    "    return(np.kron(np.array([beta ** i for i in range(self.nbT)]),u_xy) )\n",
    "\n",
    "dynamic_problem.u_txy = u_txy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_example.u_txy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intertemporal optimization problem (finite horizon)\n",
    "\n",
    "## Duality\n",
    "\n",
    "We now have everything ready to write down the intertemporal (primal) optimization problem, which expresses\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{\\pi\\geq0}  &  \\, u^{\\top}\\pi\\\\\n",
    "s.t.~  &  A\\pi=b~\\left[U\\right]\n",
    "\\end{align*}\n",
    "\n",
    "while the dual problem is given by\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{U}  & \\, b^{\\top}U\\\\\n",
    "s.t.~  &  A^{\\top} U\\geq u~\\left[\\pi \\geq 0\\right]  .\n",
    "\\end{align*}\n",
    "\n",
    "We shall write explicitely both problems to interpret the dual, but first, let's implement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_lp(self):\n",
    "    m = grb.Model()\n",
    "    pi_txy = m.addMVar(shape=self.nbT*self.nbX*self.nbY )\n",
    "    m.setObjective(self.u_txy() @ pi_txy, grb.GRB.MAXIMIZE)\n",
    "    m.addConstr( self.A_tx_txy() @ pi_txy == self.b_tx() )\n",
    "    m.optimize()\n",
    "    piopt_t_x_y = np.array(m.getAttr('X')).reshape((self.nbT,self.nbX,self.nbY))\n",
    "    Uopt_t_x = np.array(m.getAttr('pi')).reshape((self.nbT,self.nbX))\n",
    "    return({'pi_t_x_y':piopt_t_x_y,'U_t_x':Uopt_t_x})\n",
    "\n",
    "dynamic_problem.solve_lp = solve_lp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_lp = small_example.solve_lp()\n",
    "sol_lp['U_t_x'][0,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Primal problem: central planner's problem\n",
    "\n",
    "The central planner's problem is:\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{\\pi_{xy}^{t}\\geq0}  &  \\sum_{x\\in\\mathcal{X},~y\\in\\mathcal{Y},~t\\in\\mathcal{T}}\\pi_{xy}^{t}u_{xy}^{t} \\\\\n",
    "s.t.  &  \n",
    "\\sum_{y^{\\prime}\\in\\mathcal{Y}}\\pi_{xy^{\\prime}}^{0}=q_{x}~\\left[U_{x}^{1}\\right] \\\\\n",
    "& \\sum_{y^{\\prime}\\in\\mathcal{Y}}\\pi_{x^{\\prime}y^{\\prime}}^{t}=\\sum_{x\\in\\mathcal{X},~y\\in\\mathcal{Y}}P_{x^{\\prime}|xy}\\pi_{xy}^{t-1}~\\forall t\\in\\mathcal{T}\\backslash\\left\\{  0\\right\\}  ~\\left[\n",
    "U_{x^{\\prime}}^{t}\\right]\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dual problem: dynamic programming\n",
    "\n",
    "We have introduced $U_{x}^{t}$ the Lagrange multiplier associated with the constraints at time $t$. The dual problem is:\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{U_{x}^{t},~t\\in\\mathcal{T},~x\\in\\mathcal{X}}  &  \\sum_{x\\in\\mathcal{X}}q_{x}U_{x}^{0} \\\\\n",
    "s.t.~  &  U_{x}^{t}\\geq u_{xy}^{t}+\\sum_{x^{\\prime}}U_{x^{\\prime}}^{t+1}P_{x^{\\prime}|xy}~\\forall x\\in\\mathcal{X},~y\\in\\mathcal{Y},~t\\in\\mathcal{T}\\backslash\\left\\{  T - 1\\right\\} \\\\\n",
    "&  U_{x}^{T-1}\\geq u_{xy}^{T-1}~\\forall x\\in\\mathcal{X},y\\in\\mathcal{Y}\n",
    "\\end{align*}\n",
    "\n",
    "We shall see that $U_{x}^t$ represents the *intertemporal payoff* of being in state $x$ at time $t$, while the constraints is a *Bellman equation*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complementary slackness and Bellman's equation\n",
    "By complementary slackness, we have\n",
    "\n",
    "\\begin{align*}\n",
    "\\pi_{xy}^{t}>0\\Longrightarrow U_{x}^{t}=u_{xy}^{t}+\\sum_{x^{\\prime}}U_{x^{\\prime}}^{t+1}P_{x^{\\prime}|xy}\n",
    "\\end{align*}\n",
    "\n",
    "whose interpretation is immediate: if $y$ is the optimal choice in state $x$ at time $t$, then the intertemporal payoff of $x$ at $t$ is the sum of her myopic payoff $u_{xy}^{t}$ and her expected payoff at the next step.\n",
    "\n",
    "As a result, the dual variable is called *intertemporal payoff* in the vocable of dynamic programming. The relation yields *Bellman's equation*, verified as soon as $n^t_x>0$:\n",
    "\n",
    "<a name=\"bellman\"></a>\n",
    "\\begin{align*}\n",
    "U_{x}^{t}=\\max_{y\\in\\mathcal{Y}}\\left\\{  u_{xy}^{t}+\\sum_{x^{\\prime}}U_{x^{\\prime}}^{t+1}P_{x^{\\prime}|xy}\\right\\},\n",
    "\\end{align*}\n",
    "\n",
    "It is easy to see that if $U$ satisfies Bellman's equation for all $(x,t)$, then it solves the dual program.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Backward-forward induction\n",
    "\n",
    "But there is in fact a much faster way to compute the primal and dual solutions without having to use the full power of a linear programming solver. Along with the fact that $U^{T}=0$, [Bellman's equation](#bellman) implies that there is a particularly simple method to obtain the dual variables $U^{t}$, by solving recursively backward in time, from $t=T-1$ to $t=0$. This method is called *backward induction*:\n",
    "\n",
    "---\n",
    "**Algorithm [Backward induction]**\n",
    "1. Set $U^{T}=0$\n",
    "\n",
    "2. For $t=T-1$ down to $0$, set $U_{x}^{t}:=\\max_{y\\in\\mathcal{Y}}\\left\\{u_{xy}^{t}+\\sum_{x^{\\prime}}U_{x^{\\prime}}^{t+1}P_{x^{\\prime}|xy}\\right\\}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_induction(self):\n",
    "    U_t_x = np.zeros((self.nbT,self.nbX))\n",
    "    contVals_x = self.u_txy().reshape(self.nbT,self.nbX,self.nbY)[self.nbT-1,:,:].max(axis=1)\n",
    "    U_t_x[ self.nbT-1,:] = contVals_x\n",
    "    for t in range(self.nbT-2, -1, -1):\n",
    "        myopic_x_y = self.u_txy().reshape(self.nbT,self.nbX,self.nbY)[t,:,:]\n",
    "        EcontVals_x_y = (contVals_x.transpose() @ self.P_x_xy()).reshape((self.nbX,self.nbY))\n",
    "        contVals_x = (myopic_x_y + EcontVals_x_y).max(axis = 1)\n",
    "        U_t_x[ t,:] = contVals_x\n",
    "    return(U_t_x)\n",
    "\n",
    "dynamic_problem.backward_induction = backward_induction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_bi = small_example.backward_induction()\n",
    "\n",
    "print(tabulate(zip(sol_bi.flatten(),sol_lp['U_t_x'].flatten()),headers=[\"bi\",\"lp\" ]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The primal variables $\\pi^{t}$ are then deduced also by recursion, but this time forward in time from $t=1$ to $t=T-1$, by the so-called *forward induction* method:\n",
    "\n",
    "---\n",
    "**Algorithm [Forward induction]**\n",
    "1. Set $b^{0}=q$ and compute $\\left(  U^{t}\\right)$ by backward induction.\n",
    "\n",
    "2. For $t=0$ up to $T-1$, pick $\\pi^{t}$ such that $\\pi_{xy}^{t}/n_{x}^{t}$ is a probability measure supported in the set\n",
    "\n",
    "\\begin{align*}\n",
    "\\left\\{  y:U_{x}^{t}=u_{xy}^{t}+\\sum_{x^{\\prime}}U_{x^{\\prime}}^{t+1}P_{x^{\\prime}|xy}\\right\\}  .\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "3. Set $n_{x^{\\prime}}^{t+1}:=\\sum_{x\\in\\mathcal{X},~y\\in\\mathcal{Y}}P_{x^{\\prime}|xy}\\pi_{xy}^{t-1}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_forward_induction(self, q=np.array([]) ):\n",
    "    U_t_x = self.backward_induction()\n",
    "    n_t_x = self.b_tx(q).reshape((self.nbT,self.nbX))\n",
    "    pi_t_x_y = np.zeros((self.nbT,self.nbX,self.nbY))\n",
    "    self.u_txy().reshape(self.nbT,self.nbX,self.nbY)[self.nbT-1,:,:].max(axis=1)\n",
    "    for t in range(self.nbT-1):\n",
    "        myopic_x_y = self.u_txy().reshape(self.nbT,self.nbX,self.nbY)[t,:,:]\n",
    "        EcontVals_x_y = (U_t_x[t+1,:].transpose() @ self.P_x_xy()).reshape((self.nbX,self.nbY))\n",
    "        Y_x = (myopic_x_y + EcontVals_x_y).argmax(axis = 1)\n",
    "        for x in range(self.nbX):\n",
    "            pi_t_x_y[t,x,Y_x[x]] = n_t_x[t,x] \n",
    "        n_t_x[t+1,:] = (self.P_x_xy() @ pi_t_x_y[t,:,:].flatten())\n",
    "    myopic_x_y = self.u_txy().reshape(self.nbT,self.nbX,self.nbY)[self.nbT-1,:,:]\n",
    "    Y_x = myopic_x_y.argmax(axis = 1)\n",
    "    for x in range(self.nbX):\n",
    "        pi_t_x_y[self.nbT-1,x,Y_x[x]] = n_t_x[self.nbT-1,x] \n",
    "    return({'pi_t_x_y':pi_t_x_y,'U_t_x':U_t_x})\n",
    "\n",
    "dynamic_problem.backward_forward_induction = backward_forward_induction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_bf = small_example.backward_forward_induction()\n",
    "#print(Value )\n",
    "print('Value sol bf=',(small_example.u_txy() * (sol_bf['pi_t_x_y'].flatten())).sum(), ' ; value sol lp=',(small_example.u_txy() * (sol_lp['pi_t_x_y'].flatten())).sum())\n",
    "\n",
    "print(tabulate(zip(sol_bf['pi_t_x_y'].flatten(),sol_lp['pi_t_x_y'].flatten()),headers=[\"pi (bi)\",\"pi (lp)\" ]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks \n",
    "\n",
    "1. The dual variable is $U$ not necessarily unique (if $(x,t)$ is not visited, $U^t_x$ can take typically several values); the primal variable is not either, as there may be ties between several states.\n",
    "\n",
    "2. The computation by the combination of the backward and forward algorithms is much faster than the computation by a black-box linear programming solver.\n",
    "\n",
    "3. However, as soon as we introduce capacity constraints, the computation by backward induction no longer works, and the linear programming formulation is necessary, as we shall now see."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Population constraints\n",
    "\n",
    "Let us now assume that for each category $x$ at most $k_y$ individuals can take choice $y$ at each time, where $\\sum_y k_y \\geq \\sum_x q_x$. An additional constraint is therefore that for all $t\\in\\mathcal{T}$, $x \\in \\mathcal{X}$ and $y \\in \\mathcal{Y}$, one should have:\n",
    "\n",
    "$$\n",
    "\\pi^t_{xy} \\leq k_y\n",
    "$$\n",
    "\n",
    "---\n",
    "**Exercise**. \n",
    "1. Write down (in matrix form) how the primal problem is modified.\n",
    "2. In the example above, assume $k_0=+\\infty$ and $k_1 = 1$ and compute the problem using LP.\n",
    "3. Can we use the backward-forward algorithm?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution to the exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a constraint of the form $B\\pi\\leq m$, the primal problem then writes\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{\\pi\\geq0}  &  u^{\\top}\\pi\\\\\n",
    "s.t.~  &  A \\pi=b~\\left[  U\\right] \\\\\n",
    "&  B \\pi\\leq m~\\left[  \\Lambda\\right]\n",
    "\\end{align*}\n",
    "\n",
    "whose dual is\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{U,\\Lambda\\geq0}  &  b^{\\top}U+m^{\\top}\\Lambda\\\\\n",
    "s.t.~  &  A^\\top U+B^\\top\\Lambda\\geq u~\\left[  \\pi\\right]\n",
    "\\end{align*}\n",
    "\n",
    "The dual becomes\n",
    "\\begin{align*}\n",
    "\\min_{U_{x}^{t},\\lambda_{y}^{t}\\geq0}  &  \\sum_{x\\in\\mathcal{X}}n_{x}U_{x}%\n",
    "^{1}+\\sum_{x\\in\\mathcal{X}}\\sum_{t\\in\\mathcal{T}}m_{y}\\lambda_{y}^{t}\\\\\n",
    "s.t.~  &  U_{x}^{t}\\geq u_{xy}^{t}-\\lambda_{y}^{t}+\\sum_{x^{\\prime}%\n",
    "}U_{x^{\\prime}}^{t+1}P_{x^{\\prime}|xy}~\\forall x\\in\\mathcal{X},~y\\in\n",
    "\\mathcal{Y},~t\\in\\mathcal{T}\\backslash\\left\\{  T\\right\\} \\nonumber\\\\\n",
    "&  U_{x}^{T}\\geq u_{xy}^{T}~\\forall x\\in\\mathcal{X},y\\in\\mathcal{Y}\\nonumber\n",
    "\\end{align*}\n",
    "\n",
    "and $\\lambda_{y}^{t}$ interprets as the shadow price of alternative $y$ at\n",
    "time $t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution to the exercise**. The constraints $\\pi^t_{xy} \\leq k_y$  express as\n",
    "$$vec_C(\\pi) \\leq 1_{XT} \\otimes k$$\n",
    "However, when the upper bound is $+\\infty$, we'd rather drop the corresponding constraints. In which case, the matrix form becomes $$(I_{XT} \\otimes (0,1)) vec_C(\\pi) \\leq 1_{XT}$$ which we code as: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Abis_tx_txy = spr.kron(spr.identity(small_example.nbT * small_example.nbX),np.array([0,1]))\n",
    "bbis_tx = np.ones(small_example.nbX * small_example.nbT )\n",
    "\n",
    "m = grb.Model()\n",
    "pi_txy = m.addMVar(shape=small_example.nbT*small_example.nbX*small_example.nbY )\n",
    "m.setObjective(small_example.u_txy() @ pi_txy, grb.GRB.MAXIMIZE)\n",
    "m.addConstr( small_example.A_tx_txy() @ pi_txy == small_example.b_tx())\n",
    "m.addConstr( Abis_tx_txy @ pi_txy <= bbis_tx) \n",
    "m.optimize()\n",
    "pioptbis_t_x_y = np.array(m.getAttr('X')).reshape((small_example.nbT,small_example.nbX,small_example.nbY))\n",
    "Uoptbis_t_x = np.array(m.getAttr('pi'))[:(small_example.nbX * small_example.nbT)].reshape((small_example.nbT,small_example.nbX))\n",
    "                        \n",
    "Uoptbis_t_x[0,]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another plausible constraint is to assume that the maximum of buses that can go on maintenance all together is $k$, say e.g. $k=1.5$. This constraint now writes\n",
    "$$\\left( I_{T}\\otimes 1_{X}^{\\top }\\otimes (0,1)\\right) vec\\left( \\pi \\right)\n",
    "\\leq k 1_{T},$$\n",
    "and is coded as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Abis_tx_txy = spr.kron(spr.kron(spr.identity(small_example.nbT ), np.ones(small_example.nbX)),np.array([0,1]))\n",
    "bbis_tx = 1.5 * np.ones( small_example.nbT )\n",
    "\n",
    "m = grb.Model()\n",
    "pi_txy = m.addMVar(shape=small_example.nbT*small_example.nbX*small_example.nbY )\n",
    "m.setObjective(small_example.u_txy() @ pi_txy, grb.GRB.MAXIMIZE)\n",
    "m.addConstr( small_example.A_tx_txy() @ pi_txy == small_example.b_tx())\n",
    "m.addConstr( Abis_tx_txy @ pi_txy <= bbis_tx) \n",
    "m.optimize()\n",
    "pioptbis_t_x_y = np.array(m.getAttr('X')).reshape((small_example.nbT,small_example.nbX,small_example.nbY))\n",
    "Uoptbis_t_x = np.array(m.getAttr('pi'))[:(small_example.nbX * small_example.nbT)].reshape((small_example.nbT,small_example.nbX))\n",
    "                        \n",
    "Uoptbis_t_x[0,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In these cases, the linear programming solver is our only way to compute the solution. The backward induction algorithm would not apply because of the capacity constraint."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}