notebooks/04.04-Maintenance-Planning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9O_tHTY_PL_s"
   },
   "source": [
    "# Maintenance Planning\n",
    "Keywords: scheduling, cbc usage, gdp, disjunctive programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "import shutil\n",
    "import sys\n",
    "import os.path\n",
    "\n",
    "if not shutil.which(\"pyomo\"):\n",
    "    !pip install -q pyomo\n",
    "    assert(shutil.which(\"pyomo\"))\n",
    "\n",
    "if not (shutil.which(\"cbc\") or os.path.isfile(\"cbc\")):\n",
    "    if \"google.colab\" in sys.modules:\n",
    "        !apt-get install -y -qq coinor-cbc\n",
    "    else:\n",
    "        try:\n",
    "            !conda install -c conda-forge coincbc \n",
    "        except:\n",
    "            pass\n",
    "\n",
    "assert(shutil.which(\"cbc\") or os.path.isfile(\"cbc\"))\n",
    "import pyomo.environ as pyo\n",
    "import pyomo.gdp as gdp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Y_yzzUJQDJtr"
   },
   "source": [
    "## Problem statement\n",
    "\n",
    "\n",
    "A process unit is operating over a maintenance planning horizon from $1$ to $T$ days.  On day $t$ the unit makes a profit $c[t]$ which is known in advance. The unit needs to shut down for $P$ maintenance periods during the planning period.  Once started, a maintenance period takes $M$ days to finish.\n",
    "\n",
    "Find a maintenance schedule that allows the maximum profit to be produced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RrSgeD0ewWpE"
   },
   "source": [
    "## Modeling with disjunctive constraints\n",
    "\n",
    "The model is comprised of two sets of the binary variables indexed 1 to $T$. Binary variables $x_t$ correspond to the operating mode of the process unit, with $x_t=1$ indicating the unit is operating on day $t$ and able to earn a profit $c_t$. Binary variable $y_t=1$ indicates the first day of a maintenance period during which the unit is not operating and earning $0$ profit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RrSgeD0ewWpE"
   },
   "source": [
    "### Objective\n",
    "\n",
    "The planning objective is to maximize profit realized during the days the plant is operational. \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mbox{Profit} & = \\max_{x, y} \\sum_{t=1}^T c_t x_t\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "subject to completing $P$ maintenance periods.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RrSgeD0ewWpE"
   },
   "source": [
    "### Constraints\n",
    "\n",
    "**Number of planning periods is equal to P.**\n",
    "\n",
    "Completing $P$ maintenance periods requires a total of $P$ starts.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{t=1}^T y_t &  = P \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "**No more than one maintenance period can start in any consecutive set of M days.**\n",
    "\n",
    "No more than one maintenance period can start in any consecutive set of M days.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{s=0}^{M-1}y_{t+s} & \\leq 1 \\qquad \\forall t = 1, 2, \\ldots, T-M+1\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "This last requirement could be modified if some period of time should occur between maintenance periods.\n",
    "\n",
    "**The unit must shut down for M days following a maintenance start.**\n",
    "\n",
    "The final requirement is a disjunctive constraint that says either $y_t = 0$ or the sum $\\sum_{s}^{M-1}x_{t+s} = 0$, but not both. Mathematically, this forms a set of constraints reading\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\left(y_t = 0\\right) \\lor \\left(\\sum_{s=0}^{M-1}x_{t+s} = 0\\right)\\qquad \\forall t = 1, 2, \\ldots, T-M+1\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $\\lor$ denotes a disjunction.\n",
    "\n",
    "Disjunctive constraints of this nature are frequently encountered in scheduling problems. In this particular case, the disjunctive constraints can be replaced by a set of linear inequalities using the Big-M method.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{s=0}^{M-1}x_{t+s} \\leq M(1-y_t) \\qquad \\forall t = 1, 2, \\ldots, T-M+1\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "In this case, the $M$ appearing on the right-hand side of this constraint happens to be the same as the length of each maintenance period."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pyomo solution using the big-M method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Bae6-dR_lYkm"
   },
   "source": [
    "### Parameter values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "L5TaVPbkEPZ0"
   },
   "outputs": [],
   "source": [
    "# problem parameters\n",
    "T = 90        # planning period from 1..T\n",
    "M = 3         # length of maintenance period\n",
    "P = 4         # number of maintenance periods\n",
    "\n",
    "# daily profits\n",
    "c = {k:np.random.uniform() for k in range(1, T+1)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "e28FWA1GlNyQ"
   },
   "source": [
    "### Pyomo model\n",
    "\n",
    "The disjunctive constraints can be represented directly in Pyomo using the [Generalized Disjunctive Programming](https://pyomo.readthedocs.io/en/latest/modeling_extensions/gdp.html) extension. The GDP extension transforms the disjunctive constraints to an MILP using convex hull and cutting plane methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1014
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1543,
     "status": "ok",
     "timestamp": 1557947668999,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "cUARDFzP9fla",
    "outputId": "ba0f2867-3a14-4314-c09f-bb4a76ebd1a7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# ==========================================================\n",
      "# = Solver Results                                         =\n",
      "# ==========================================================\n",
      "# ----------------------------------------------------------\n",
      "#   Problem Information\n",
      "# ----------------------------------------------------------\n",
      "Problem: \n",
      "- Name: unknown\n",
      "  Lower bound: 41.42942891\n",
      "  Upper bound: 41.42942891\n",
      "  Number of objectives: 1\n",
      "  Number of constraints: 176\n",
      "  Number of variables: 178\n",
      "  Number of binary variables: 180\n",
      "  Number of integer variables: 180\n",
      "  Number of nonzeros: 90\n",
      "  Sense: maximize\n",
      "# ----------------------------------------------------------\n",
      "#   Solver Information\n",
      "# ----------------------------------------------------------\n",
      "Solver: \n",
      "- Status: ok\n",
      "  User time: -1.0\n",
      "  System time: 0.02\n",
      "  Wallclock time: 0.02\n",
      "  Termination condition: optimal\n",
      "  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.\n",
      "  Statistics: \n",
      "    Branch and bound: \n",
      "      Number of bounded subproblems: 0\n",
      "      Number of created subproblems: 0\n",
      "    Black box: \n",
      "      Number of iterations: 0\n",
      "  Error rc: 0\n",
      "  Time: 0.049736976623535156\n",
      "# ----------------------------------------------------------\n",
      "#   Solution Information\n",
      "# ----------------------------------------------------------\n",
      "Solution: \n",
      "- number of solutions: 0\n",
      "  number of solutions displayed: 0\n"
     ]
    }
   ],
   "source": [
    "def maintenance_planning_bigm(c, T, M, P):\n",
    "    m = pyo.ConcreteModel()\n",
    "\n",
    "    m.T = pyo.RangeSet(1, T)\n",
    "    m.Y = pyo.RangeSet(1, T - M + 1)\n",
    "    m.S = pyo.RangeSet(0, M - 1)\n",
    "\n",
    "    m.c = pyo.Param(m.T, initialize = c)\n",
    "    m.x = pyo.Var(m.T, domain=pyo.Binary)\n",
    "    m.y = pyo.Var(m.T, domain=pyo.Binary)\n",
    "\n",
    "    # objective\n",
    "    m.profit = pyo.Objective(expr = sum(m.c[t]*m.x[t] for t in m.T), sense=pyo.maximize)\n",
    "\n",
    "    # required number P of maintenance starts\n",
    "    m.sumy = pyo.Constraint(expr = sum(m.y[t] for t in m.Y) == P)\n",
    "\n",
    "    # no more than one maintenance start in the period of length M\n",
    "    m.sprd = pyo.Constraint(m.Y, rule = lambda m, t: sum(m.y[t+s] for s in m.S) <= 1)\n",
    "\n",
    "    # disjunctive constraints\n",
    "    m.bigm = pyo.Constraint(m.Y, rule = lambda m, t: sum(m.x[t+s] for s in m.S) <= M*(1 - m.y[t]))\n",
    "    \n",
    "    return m\n",
    "   \n",
    "m = maintenance_planning_bigm(c, T, M, P)\n",
    "pyo.SolverFactory('cbc').solve(m).write()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAEYCAYAAAA+gNBwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfVklEQVR4nO3de7RkZX2n8ecrDXLTIDYqDd0cVNKRYAQWEhKJQWGWKFdvARKEmDCsrJGAGRGB5YwxY1aISRxiNGYISkAUUESDODEmoEvjGKSb9sJFlEtDNzTQREFACSK/+aM2WBbnUqf6nFN1dj2ftc7q2vff3uc9Vd/e7967UlVIkiSpHZ427AIkSZI0dwx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFO0tAk+Yck7+ljvuuTHNA1vDbJQfNZWz+SrEyyJsmDSU7urVOShmHJsAuQpJlU1S8Pu4YpnAZ8qar2aobf3z0xyVrghKr614UuTNL48sydJE0iST//+d0FuH6+a5Gk2TDcSVowSfZKcm3TjXkJsGXXtNOT3NJMuyHJa7umTdoNm+TtST7VM+5vkpw9xfbXJjmjWf8PkpyXZMue6e9I8i3g4SRLkrwoyZeS3N90ux7ezHsV8ArgA0keSvKL3XUm+SiwAvhsM/20KWo6JsnqJA80+39An4dTkiZluJO0IJJsAXwG+CiwPfBJ4PVds9wC/AbwC8C7gQuT7DjDai8EDk6yXbONJcBRzTam8jvAq4AXAL8IvLNn+jHAIcB2QIDPAl8AngP8IfCxJCur6pXAV4CTqmrbqvpu90qq6k3AHcBhzfT39haS5G3N9v8r8CzgSGDtDPssSdMy3ElaKPsBmwNnV9VPqupS4JonJlbVJ6vqrqp6vKouAb4H7DvdCqtqA/Bl4I3NqIOB+6pq9TSLfaCq1lXV94E/pRPmur2/mf7jpuZtgbOq6tGqugq4YpJlZi3JDsC7gN+uqmub/f52Va1tpr8iycSmbkfS+DHcSVooy4A7q6q6xt3+xIskxyX5RtP9eT+wB7C0j/WeDxzbvD6W6c/aAazr2f6yaaYvA9ZV1eM9y+zUR10zOQj4dlV9c4rpv0fnzKEkzYrhTtJC2QDslKQ7sKwASLIL8PfAScCzq2o74Dr6CzefAX4lyR7AocDHZph/ec/27+qZ3h0+7wKWJ+l+r1wB3NlHXb3r6rU9cP9kE5rr+g4DzktyXJ/bkiTAcCdp4XwNeAw4ublR4XX8rNt1GzpBaCNAkjfTOXM3o6p6BLgU+Djw9aq6Y4ZF3pJk5yTbA2cCl0wz79XAw8BpSTZvbnY4DLi4n9qAe4DnTzFtDbB/kpekY7ckL2qmXQGsqaoDquqCPrclSYDhTtICqapHgdcBvwv8gM6ND5c1024A/opOALwHeDHw1Vms/vxmmZm6ZKETAr8A3Nr8TPkQ5abmw4FXA/cBfwscV1Xf6bOuPwPe2XQ1n9qz7v/XbPsK4EHg08BWzeQXAjf1uQ1J+jn5+ctfJGnxSbIC+A7wvKr64TTzrWURPFQ4yZHARFWdPeRSJC1CnrmTtKg118P9d+Di6YLdIvNd4ISpntcnSdPx68ckLVpJtqHTjXs7ncegtELTTd3XNYeS1MtuWUmSpBaZsVs2yUeS3JvkuoUoSJIkSYOb8cxdkpcDDwEXVFVf3QRLly6tiYmJTa9OkiRJT7F69er7qmqHyabNeM1dVX15tl+BMzExwapVq2aziCRJkvqU5Papps3Z3bJJTkyyKsmqjRs3ztVqJUmSNAtzdrdsVZ0DnAOwzz77eJeGZm3i9M9NOn7tWYcscCXScEz2N/BE+59umiR181Eo2iR+4EiS2qQNn2tDDXdtOICShsP3D0mjYtR6nmYMd0kuAg4AliZZD7yrqj4834VJktrNgC7Nj37ulj1mIQqRJEnSpvO7ZSVJklrEGypGgF0TkiRprnjmTpIkqUU8cydJWjSmuyvRXhCpw3An9VioD4hRu3VekkaZ4b1/dstKkiS1iGfuJGkGU50x8Oyr5oLtSHPNcCdJktSHxdI13Jpw5/98pOHz71CaW4slTGi0tCbcSRptfkgtDuMW0Mdtf0eBx3z+Ge4kSVoA/gdHvebrel7DnaSh80NPkuaOj0KRJElqEcOdJElSi4xkt6wXW44Wu8ykxWkh/3Z9nxgdfobOrO3HaCTDnTSK2v5mMBt+kGu+zXUbs81qnBjuFimDhqTFwPcqzYYhfG4Y7gT4BizNNT+k1BZ+Piw+iy7c+YYpTc434PHk733hjcox9/NQU/FuWUmSpBYx3EmSJLXIouuWHXWeJtd8s41JkqZjuGshP/wlzQffW6TFwXA3hVG5YFYzm68vXpYkaTEa+3A36v8THfX6ND/8vc8PA780t0b9vWrU65svYx/uND+m+4Ma1z82SYuL71VarAx30jzzA0KS5s+ov8cOoz7DnSRJGvmQpP4Z7iQJP9gktYfhbgB+CEijzb9RSePMcLdARuUuPT/0JGl8+RkwHuY93NmQNIpslzPzGEnaFL6HDI9n7qQxNipnlKW2MNBoFBjutCj4hjmeDJ+SNHtPG3YBkiRJmjtjcebOsz4ad/4NjBZ/H5LmU19n7pIcnOSmJDcnOX2+i5IkSdJgZjxzl2Qz4IPAfwHWA9ckubyqbpjv4iRJ88Ozh1J79dMtuy9wc1XdCpDkYuAIwHAnNfyglCSNilTV9DMkbwAOrqoTmuE3Ab9aVSf1zHcicGIzuBK4ae7L1SK1FLhv2EVopNgm1Ms2oV62ientUlU7TDahnzN3mWTcUxJhVZ0DnDPLwjQGkqyqqn2GXYdGh21CvWwT6mWbGFw/N1SsB5Z3De8M3DU/5UiSJGlT9BPurgF2S7Jrki2Ao4HL57csSZIkDWLGbtmqeizJScA/A5sBH6mq6+e9MrWJ3fXqZZtQL9uEetkmBjTjDRWSJElaPPz6MUmSpBYx3EmSJLWI4U5zJsnyJF9McmOS65Oc0ozfPsm/JPle8++zhl2rFlaSzZKsSXJFM2ybGGNJtktyaZLvNO8Xv2abGG9J/qj53LguyUVJtrRNDM5wp7n0GPC2qnoRsB/wliS7A6cDV1bVbsCVzbDGyynAjV3Dtonx9tfA56vql4CX0GkbtokxlWQn4GRgn6rag87Nm0djmxiY4U5zpqo2VNW1zesH6bxh70Tn6+rOb2Y7HzhyKAVqKJLsDBwCnNs12jYxppI8E3g58GGAqnq0qu7HNjHulgBbJVkCbE3nebq2iQEZ7jQvkkwAewFXA8+tqg3QCYDAc4ZYmhbe2cBpwONd42wT4+v5wEbgvKar/twk22CbGFtVdSfwl8AdwAbggar6AraJgRnuNOeSbAt8CnhrVf1w2PVoeJIcCtxbVauHXYtGxhJgb+BDVbUX8DB2t4215lq6I4BdgWXANkmOHW5Vi5vhTnMqyeZ0gt3HquqyZvQ9SXZspu8I3Dus+rTgXgYcnmQtcDHwyiQXYpsYZ+uB9VV1dTN8KZ2wZ5sYXwcBt1XVxqr6CXAZ8OvYJgZmuNOcSRI619HcWFXv65p0OXB88/p44B8XujYNR1WdUVU7V9UEnQukr6qqY7FNjK2quhtYl2RlM+pA4AZsE+PsDmC/JFs3nyMH0rlm2zYxIL+hQnMmyf7AV4Bv87Prq86kc93dJ4AVdP6I31hV3x9KkRqaJAcAp1bVoUmejW1ibCXZk84NNlsAtwJvpnOywTYxppK8GziKzlMX1gAnANtimxiI4U6SJKlF7JaVJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U7Sk5Jcn+SAYdcxF0ZpX5KsTXLQHK3rH5K8Z6G3K2nxMNxJelJV/XJVfQkWVzCYrNbufZGkcWK4kzTSkiwZdg2StJgY7qSWSVJJXtg1/GQ3XnOG69Qk30ryQJJLkmzZNe/aJAcl+SiwAvhskoeSnDbFtl6U5EtJ7m+6QQ/vWdcZSW5I8oMk5/Vsa1mSTyXZmOS2JCf3LPuOJN8CHk6yJMnpSW5J8mCzztc2805aa/fZvD72e+8ka5p1f7KZPmnXZ1PXnc28NyU5sGva8iSXNfv0H0k+0LXontNsf7pjsVeSa5vtXQJs2bXOaX/fk9Q+5XYmmfe9ST7dNfwXSa5MsvlUy0gaDYY7afz8FnAwsCvwK8Dv9s5QVW8C7gAOq6ptq+q9vfM0H/KfBb4APAf4Q+BjSVZ2zfY7wKuAFwC/CLyzWfZpzbLfBHYCDgTemuRVXcseAxwCbFdVjwG3AL8B/ALwbuDCJDv2U+t0+51kC+DTwD8A2wMXAa+dbAXNvp0EvLSqntHs29pm2mbAFcDtwESzXxf3sf0pj0VT22eAjza1fRJ4/RT7N60+j3m3PwdekWTPJH/Q1P66qvrJINuXtHAMd9L4eX9V3VVV36fzYb/ngOvZD9gWOKuqHq2qq+iEm2O65vlAVa1rtvWnXdNeCuxQVX/SLHsr8PfA0T11rquqHwNU1Sebuh+vqkuA7wH7zqLeqfZ7P2BJM/0nVXUZ8PUp1vFT4OnA7kk2r6q1VXVLM21fYBnw9qp6uKoeqap/62P70x2L/YDNgbOb2i4FrpnFPnfr55g/qar+AzgbuAA4A3hNVT0AkOQVSSYGrEPSPDPcSePn7q7XP6IT0AaxDFhXVY93jbudzlmhJ6zrmbaseb0LsKzpzr0/yf3AmcBzp1iWJMcl+UbX/HsAS2dR71T7vQy4s6pqqm0/oapuBt4K/DFwb5KLkzyxT8uB25uzjLPZ/nTHYrLabp9yD6fXzzHvtQZ4MXBGVXUfk98DMmAdkuaZ4U5qnx8BW3cNP2/A9dQM0+8CljfdfU9YAdzZNby8Z9pdzet1wG1VtV3XzzOq6jWTbT/JLnTOMp0EPLuqtgOu42cBY6Zap7MB2ClJd1hZPtXMVfXxqtqfTlgqOt2XT+zTisz+BpDpjsVkta3oWb7f33c/x/xJSV4MfAg4n06Ye2L84cBhwHlJjpvNjkpaGIY7qX2+Afx2ks2SHAz85oDruQd4/jTTrwYeBk5Lsnk6z5Q7jJ+/zuwtSXZOsj2ds0SXNOO/DvywuTlhq6bWPZK8dIptbUMnSG0ESPJmOmfu+q11Ol+j0916UnPjxhFM0d2bZGWSVyZ5OvAI8ONm2Sf2aQNwVpJtkmyZ5GV9bH+6Y/E14DHg5Ka2101S2zfo7/fd9zFPshOdruM/AP4b8OL87JmBVwBrquqAqrqgj/2TtMAMd1L7nEInZN1P54aGzwy4nj8D3tl04Z3aO7GqHgUOB14N3Af8LXBcVX2na7aP07nh4tbm5z3Nsj9tatwTuK1Z/lw6N0s8RVXdAPwVnbBzD52uwq/2W+t0mv14HfD7dI7ZsXQCzH9OMvvTgbOaeu+mcyPJmT379EI6N3isB47qY/tTHouu2n4X+EGzvst6VtHX77vfY57kmcD/Bd5XVZdX1Y+Av6BzzSTN/t00035JGp78/KUckjQ3kqwFTqiqfx12LbOV5Grg76rqvGHXMmqSHAlMVNXZQy5F0hQ8cydp7CX5zSTPa7o+j6fzqJLPD7uuEfVd4IQkZw+7EEmT88nvkgQrgU/QuYP1FuANVbVhuCWNpqaLfI8ZZ5Q0NHbLSpIktciM3bJJPpLk3iTXLURBkiRJGtyMZ+6SvBx4CLigqvo6Fb906dKamJjY9OokSZL0FKtXr76vqnaYbNqM19xV1Zdn+zUzExMTrFq1ajaLSJIkqU9Jpvy2mjm7oSLJicCJACtW/OwB6hOnf+4p864965Bpp002frppMy0z6hbqGA17faNQw7iub66Nwj5Nt7655t+N6xu1Glzf6NcwX+vrx5w9CqWqzqmqfapqnx12mPQsoSRJkuaZz7mTJElqEcOdJElSi/TzKJSL6Hyf48ok65P8/vyXJUmSpEH0c7fsMQtRiCRJkjad3bKSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSi/QV7pIcnOSmJDcnOX2+i5IkSdJgZgx3STYDPgi8GtgdOCbJ7vNdmCRJkmavnzN3+wI3V9WtVfUocDFwxPyWJUmSpEGkqqafIXkDcHBVndAMvwn41ao6qWe+E4ETm8GVwE1zX64WqaXAfcMuQiPFNqFetgn1sk1Mb5eq2mGyCUv6WDiTjHtKIqyqc4BzZlmYxkCSVVW1z7Dr0OiwTaiXbUK9bBOD66dbdj2wvGt4Z+Cu+SlHkiRJm6KfcHcNsFuSXZNsARwNXD6/ZUmSJGkQM3bLVtVjSU4C/hnYDPhIVV0/75WpTeyuVy/bhHrZJtTLNjGgGW+okCRJ0uLhN1RIkiS1iOFOkiSpRQx3mjNJlif5YpIbk1yf5JRm/PZJ/iXJ95p/nzXsWrWwkmyWZE2SK5ph28QYS7JdkkuTfKd5v/g128R4S/JHzefGdUkuSrKlbWJwhjvNpceAt1XVi4D9gLc0X1V3OnBlVe0GXNkMa7ycAtzYNWybGG9/DXy+qn4JeAmdtmGbGFNJdgJOBvapqj3o3Lx5NLaJgRnuNGeqakNVXdu8fpDOG/ZOdL6u7vxmtvOBI4dSoIYiyc7AIcC5XaNtE2MqyTOBlwMfBqiqR6vqfmwT424JsFWSJcDWdJ6na5sYkOFO8yLJBLAXcDXw3KraAJ0ACDxniKVp4Z0NnAY83jXONjG+ng9sBM5ruurPTbINtomxVVV3An8J3AFsAB6oqi9gmxiY4U5zLsm2wKeAt1bVD4ddj4YnyaHAvVW1eti1aGQsAfYGPlRVewEPY3fbWGuupTsC2BVYBmyT5NjhVrW4Ge40p5JsTifYfayqLmtG35Nkx2b6jsC9w6pPC+5lwOFJ1gIXA69MciG2iXG2HlhfVVc3w5fSCXu2ifF1EHBbVW2sqp8AlwG/jm1iYIY7zZkkoXMdzY1V9b6uSZcDxzevjwf+caFr03BU1RlVtXNVTdC5QPqqqjoW28TYqqq7gXVJVjajDgRuwDYxzu4A9kuydfM5ciCda7ZtEwPyGyo0Z5LsD3wF+DY/u77qTDrX3X0CWEHnj/iNVfX9oRSpoUlyAHBqVR2a5NnYJsZWkj3p3GCzBXAr8GY6JxtsE2MqybuBo+g8dWENcAKwLbaJgRjuJEmSWsRuWUmSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJM27JNcnOWDYdUjSODDcSWMoyW5JHkly4UJsr6p+uaq+1M+8SdYmOWieSxoJc7Gv43S8JPXHcCeNpw8C1wy7CA0uyZJh1yBpNBnupDGT5GjgfuDKGeZbm+TtSb6V5OEkH07y3CT/lOTBJP+a5Fld85+e5JZm2g1JXtuzroO6Xp/arPeBJJck2bKZ9lFgBfDZJA8lOa0ZvyzJp5JsTHJbkpN71j3p+vqsa7pllye5rNnufyT5wEz1THIc35Hkzmb7NyU5cKp9na7WrnrfkeRbwMNJLppkHZNur2c92yb5aZIdu8btkWRDkmdM2SgkLQ5V5Y8//ozJD/BM4LvAcuCPgQunmXct8O/Ac4GdgHuBa4G9gKcDVwHv6pr/jcAyOv9pPAp4GNixa10Hdb3+ejPv9sCNwB/0bPegruGnAauB/wlsATwfuBV4VZ/rm6muSZcFNgO+CfxvYBtgS2D/merpOYYrgXXAsmZ4AnjBNPs6Za1d83+j+f1tNcmxnXZ7PbVdDxzSNXwF8IfDbqP++OPPpv945k4aL/8L+HBVretz/r+pqnuq6k7gK8DVVbWmqv4T+DSdoAdAVX2yqu6qqser6hLge8C+U6z3/c283wc+C+w5TQ0vBXaoqj+pqker6lbg74Gj+1lfH3VNtey+dILW26vq4ap6pKr+rc96nvBTOkF49ySbV9Xaqrplqh3t8xi+v6rWVdWPN3F71wB7AyR5ObA78H+a4VckmZiqTkmjzXAnjYkkewIH0TkT1a97ul7/eJLhbbvWf1ySbyS5P8n9wB7A0inWe3fX6x91r2cSuwDLnlhvs+4z6ZxRnHF9fdQ11bLLgdur6rEB6gGgqm4G3krnLOm9SS5OsmyqHe3zGE4ZzGe5vSfDHfBe4H9U1aPN8O8BmWo7kkab4U4aHwfQ6aa7I8ndwKnA65Ncu6krTrILnbNXJwHPrqrtgOsYLCBUz/A64Laq2q7r5xlV9Zp5rmsdsGKSGxdmVU9Vfbyq9qcTCgv488n2dRa19h6fnxueYXvdrgH2TvJ6YCvgoqaOw4HDgPOSHDfFspJGmOFOGh/nAC+g0+24J/B3wOeAV83BurehEyQ2AiR5M52zToO4h851bE/4OvDD5kaBrZJs1lz8/9J5ruvrwAbgrCTbJNkyyctmU0+SlUlemeTpwCN0znb+dIp9HbTWJ9fRx/a6fRN4HvBXwOlV9Xgz/gpgTVUdUFUX9LF9SSPGcCeNiar6UVXd/cQP8BDwSFVtnIN130AnJHyNTth4MfDVAVf3Z8A7m67JU6vqp3TOJO0J3AbcB5wL/MJ81tW13RcCdwDrgaNmWc/TgbOaee4GnkOnC/cp+wq8ZsBau9dx1Azb696//wS+Daytqn/qmvRC4KY+titpRKWq9wy/JKntkmwB3Az8VlX9e9f4I4GJqjp7SKVJ2kSeuZOk8fQu4Kvdwa7xXeCEJGcvfEmS5oJn7iRpjCTZG/gi8C3gtVV135BLkjTHDHeSJEktMmO3bJKPJLk3yXULUZAkSZIGN+OZu+bJ5Q8BF1RVX48QWLp0aU1MTGx6dZIkSXqK1atX31dVO0w2rffhnE9RVV+e7dfQTExMsGrVqtksIkmSpD4luX2qaTOGu1ls5ETgRIAVK1bM1WrnxMTpn3vKuLVnHTIy65Pm27i12XHbX0nqNmePQqmqc6pqn6raZ4cdJj1LKEmSpHnmc+4kSZJaxHAnSZLUIv08CuUiOt91uDLJ+iS/P/9lSZIkaRD93C17zEIUIkmSpE1nt6wkSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFukr3CU5OMlNSW5Ocvp8FyVJkqTBzBjukmwGfBB4NbA7cEyS3ee7MEmSJM1eP2fu9gVurqpbq+pR4GLgiPktS5IkSYNIVU0/Q/IG4OCqOqEZfhPwq1V1Us98JwInNoMrgZvmvlwtUkuB+4ZdhEaKbUK9bBPqZZuY3i5VtcNkE5b0sXAmGfeURFhV5wDnzLIwjYEkq6pqn2HXodFhm1Av24R62SYG10+37HpgedfwzsBd81OOJEmSNkU/4e4aYLckuybZAjgauHx+y5IkSdIgZuyWrarHkpwE/DOwGfCRqrp+3itTm9hdr162CfWyTaiXbWJAM95QIUmSpMXDb6iQJElqEcOdJElSixjuNGeSLE/yxSQ3Jrk+ySnN+O2T/EuS7zX/PmvYtWphJdksyZokVzTDtokxlmS7JJcm+U7zfvFrtonxluSPms+N65JclGRL28TgDHeaS48Bb6uqFwH7AW9pvqrudODKqtoNuLIZ1ng5Bbixa9g2Md7+Gvh8Vf0S8BI6bcM2MaaS7AScDOxTVXvQuXnzaGwTAzPcac5U1YaqurZ5/SCdN+yd6Hxd3fnNbOcDRw6lQA1Fkp2BQ4Bzu0bbJsZUkmcCLwc+DFBVj1bV/dgmxt0SYKskS4Ct6TxP1zYxIMOd5kWSCWAv4GrguVW1AToBEHjOEEvTwjsbOA14vGucbWJ8PR/YCJzXdNWfm2QbbBNjq6ruBP4SuAPYADxQVV/ANjEww53mXJJtgU8Bb62qHw67Hg1PkkOBe6tq9bBr0chYAuwNfKiq9gIexu62sdZcS3cEsCuwDNgmybHDrWpxM9xpTiXZnE6w+1hVXdaMvifJjs30HYF7h1WfFtzLgMOTrAUuBl6Z5EJsE+NsPbC+qq5uhi+lE/ZsE+PrIOC2qtpYVT8BLgN+HdvEwAx3mjNJQuc6mhur6n1dky4Hjm9eHw/840LXpuGoqjOqaueqmqBzgfRVVXUstomxVVV3A+uSrGxGHQjcgG1inN0B7Jdk6+Zz5EA612zbJgbkN1RoziTZH/gK8G1+dn3VmXSuu/sEsILOH/Ebq+r7QylSQ5PkAODUqjo0ybOxTYytJHvSucFmC+BW4M10TjbYJsZUkncDR9F56sIa4ARgW2wTAzHcSZIktYjdspIkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktcj/B4n+nkmy4GsgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "def plot_schedule(m):\n",
    "    fig,ax = plt.subplots(3,1, figsize=(9,4))\n",
    "    \n",
    "    ax[0].bar(m.T, [m.c[t] for t in m.T])\n",
    "    ax[0].set_title('daily profit $c_t$')\n",
    "    \n",
    "    ax[1].bar(m.T, [m.x[t]() for t in m.T], label='normal operation')\n",
    "    ax[1].set_title('unit operating schedule $x_t$')\n",
    "    \n",
    "    ax[2].bar(m.Y, [m.y[t]() for t in m.Y])\n",
    "    ax[2].set_title(str(P) + ' maintenance starts $y_t$')\n",
    "    for a in ax:\n",
    "        a.set_xlim(0.1, len(m.T)+0.9)\n",
    "        \n",
    "    plt.tight_layout()\n",
    "\n",
    "plot_schedule(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uC2EbGWaDCJ4"
   },
   "source": [
    "## Pyomo solution using the generalized disjunctive constraints extension\n",
    "\n",
    "Disjunctive constraints can be represented directly in Pyomo using the [Generalized Disjunctive Programming](https://pyomo.readthedocs.io/en/latest/modeling_extensions/gdp.html) extension. The advantage of using the extension is that constraints can be transformed to an MILP using alternatives to the big-M, such as convex hull and cutting plane methods.\n",
    "\n",
    "The following cell replaces the Big-M constraints with disjunctions. Disjunctions are represented by lists of mutually exclusive constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1014
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 5519,
     "status": "ok",
     "timestamp": 1557947666034,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "bHrlX2c2lyDc",
    "outputId": "1c0c76bf-2d6f-496c-d259-4340c8684565"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# ==========================================================\n",
      "# = Solver Results                                         =\n",
      "# ==========================================================\n",
      "# ----------------------------------------------------------\n",
      "#   Problem Information\n",
      "# ----------------------------------------------------------\n",
      "Problem: \n",
      "- Name: unknown\n",
      "  Lower bound: 41.42942891\n",
      "  Upper bound: 41.42942891\n",
      "  Number of objectives: 1\n",
      "  Number of constraints: 440\n",
      "  Number of variables: 266\n",
      "  Number of binary variables: 356\n",
      "  Number of integer variables: 356\n",
      "  Number of nonzeros: 90\n",
      "  Sense: maximize\n",
      "# ----------------------------------------------------------\n",
      "#   Solver Information\n",
      "# ----------------------------------------------------------\n",
      "Solver: \n",
      "- Status: ok\n",
      "  User time: -1.0\n",
      "  System time: 0.03\n",
      "  Wallclock time: 0.03\n",
      "  Termination condition: optimal\n",
      "  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.\n",
      "  Statistics: \n",
      "    Branch and bound: \n",
      "      Number of bounded subproblems: 0\n",
      "      Number of created subproblems: 0\n",
      "    Black box: \n",
      "      Number of iterations: 0\n",
      "  Error rc: 0\n",
      "  Time: 0.05097818374633789\n",
      "# ----------------------------------------------------------\n",
      "#   Solution Information\n",
      "# ----------------------------------------------------------\n",
      "Solution: \n",
      "- number of solutions: 0\n",
      "  number of solutions displayed: 0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAEYCAYAAAA+gNBwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfVklEQVR4nO3de7RkZX2n8ecrDXLTIDYqDd0cVNKRYAQWEhKJQWGWKFdvARKEmDCsrJGAGRGB5YwxY1aISRxiNGYISkAUUESDODEmoEvjGKSb9sJFlEtDNzTQREFACSK/+aM2WBbnUqf6nFN1dj2ftc7q2vff3uc9Vd/e7967UlVIkiSpHZ427AIkSZI0dwx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFO0tAk+Yck7+ljvuuTHNA1vDbJQfNZWz+SrEyyJsmDSU7urVOShmHJsAuQpJlU1S8Pu4YpnAZ8qar2aobf3z0xyVrghKr614UuTNL48sydJE0iST//+d0FuH6+a5Gk2TDcSVowSfZKcm3TjXkJsGXXtNOT3NJMuyHJa7umTdoNm+TtST7VM+5vkpw9xfbXJjmjWf8PkpyXZMue6e9I8i3g4SRLkrwoyZeS3N90ux7ezHsV8ArgA0keSvKL3XUm+SiwAvhsM/20KWo6JsnqJA80+39An4dTkiZluJO0IJJsAXwG+CiwPfBJ4PVds9wC/AbwC8C7gQuT7DjDai8EDk6yXbONJcBRzTam8jvAq4AXAL8IvLNn+jHAIcB2QIDPAl8AngP8IfCxJCur6pXAV4CTqmrbqvpu90qq6k3AHcBhzfT39haS5G3N9v8r8CzgSGDtDPssSdMy3ElaKPsBmwNnV9VPqupS4JonJlbVJ6vqrqp6vKouAb4H7DvdCqtqA/Bl4I3NqIOB+6pq9TSLfaCq1lXV94E/pRPmur2/mf7jpuZtgbOq6tGqugq4YpJlZi3JDsC7gN+uqmub/f52Va1tpr8iycSmbkfS+DHcSVooy4A7q6q6xt3+xIskxyX5RtP9eT+wB7C0j/WeDxzbvD6W6c/aAazr2f6yaaYvA9ZV1eM9y+zUR10zOQj4dlV9c4rpv0fnzKEkzYrhTtJC2QDslKQ7sKwASLIL8PfAScCzq2o74Dr6CzefAX4lyR7AocDHZph/ec/27+qZ3h0+7wKWJ+l+r1wB3NlHXb3r6rU9cP9kE5rr+g4DzktyXJ/bkiTAcCdp4XwNeAw4ublR4XX8rNt1GzpBaCNAkjfTOXM3o6p6BLgU+Djw9aq6Y4ZF3pJk5yTbA2cCl0wz79XAw8BpSTZvbnY4DLi4n9qAe4DnTzFtDbB/kpekY7ckL2qmXQGsqaoDquqCPrclSYDhTtICqapHgdcBvwv8gM6ND5c1024A/opOALwHeDHw1Vms/vxmmZm6ZKETAr8A3Nr8TPkQ5abmw4FXA/cBfwscV1Xf6bOuPwPe2XQ1n9qz7v/XbPsK4EHg08BWzeQXAjf1uQ1J+jn5+ctfJGnxSbIC+A7wvKr64TTzrWURPFQ4yZHARFWdPeRSJC1CnrmTtKg118P9d+Di6YLdIvNd4ISpntcnSdPx68ckLVpJtqHTjXs7ncegtELTTd3XNYeS1MtuWUmSpBaZsVs2yUeS3JvkuoUoSJIkSYOb8cxdkpcDDwEXVFVf3QRLly6tiYmJTa9OkiRJT7F69er7qmqHyabNeM1dVX15tl+BMzExwapVq2aziCRJkvqU5Papps3Z3bJJTkyyKsmqjRs3ztVqJUmSNAtzdrdsVZ0DnAOwzz77eJeGZm3i9M9NOn7tWYcscCXScEz2N/BE+59umiR181Eo2iR+4EiS2qQNn2tDDXdtOICShsP3D0mjYtR6nmYMd0kuAg4AliZZD7yrqj4834VJktrNgC7Nj37ulj1mIQqRJEnSpvO7ZSVJklrEGypGgF0TkiRprnjmTpIkqUU8cydJWjSmuyvRXhCpw3An9VioD4hRu3VekkaZ4b1/dstKkiS1iGfuJGkGU50x8Oyr5oLtSHPNcCdJktSHxdI13Jpw5/98pOHz71CaW4slTGi0tCbcSRptfkgtDuMW0Mdtf0eBx3z+Ge4kSVoA/gdHvebrel7DnaSh80NPkuaOj0KRJElqEcOdJElSi4xkt6wXW44Wu8ykxWkh/3Z9nxgdfobOrO3HaCTDnTSK2v5mMBt+kGu+zXUbs81qnBjuFimDhqTFwPcqzYYhfG4Y7gT4BizNNT+k1BZ+Piw+iy7c+YYpTc434PHk733hjcox9/NQU/FuWUmSpBYx3EmSJLXIouuWHXWeJtd8s41JkqZjuGshP/wlzQffW6TFwXA3hVG5YFYzm68vXpYkaTEa+3A36v8THfX6ND/8vc8PA780t0b9vWrU65svYx/uND+m+4Ma1z82SYuL71VarAx30jzzA0KS5s+ov8cOoz7DnSRJGvmQpP4Z7iQJP9gktYfhbgB+CEijzb9RSePMcLdARuUuPT/0JGl8+RkwHuY93NmQNIpslzPzGEnaFL6HDI9n7qQxNipnlKW2MNBoFBjutCj4hjmeDJ+SNHtPG3YBkiRJmjtjcebOsz4ad/4NjBZ/H5LmU19n7pIcnOSmJDcnOX2+i5IkSdJgZjxzl2Qz4IPAfwHWA9ckubyqbpjv4iRJ88Ozh1J79dMtuy9wc1XdCpDkYuAIwHAnNfyglCSNilTV9DMkbwAOrqoTmuE3Ab9aVSf1zHcicGIzuBK4ae7L1SK1FLhv2EVopNgm1Ms2oV62ientUlU7TDahnzN3mWTcUxJhVZ0DnDPLwjQGkqyqqn2GXYdGh21CvWwT6mWbGFw/N1SsB5Z3De8M3DU/5UiSJGlT9BPurgF2S7Jrki2Ao4HL57csSZIkDWLGbtmqeizJScA/A5sBH6mq6+e9MrWJ3fXqZZtQL9uEetkmBjTjDRWSJElaPPz6MUmSpBYx3EmSJLWI4U5zJsnyJF9McmOS65Oc0ozfPsm/JPle8++zhl2rFlaSzZKsSXJFM2ybGGNJtktyaZLvNO8Xv2abGG9J/qj53LguyUVJtrRNDM5wp7n0GPC2qnoRsB/wliS7A6cDV1bVbsCVzbDGyynAjV3Dtonx9tfA56vql4CX0GkbtokxlWQn4GRgn6rag87Nm0djmxiY4U5zpqo2VNW1zesH6bxh70Tn6+rOb2Y7HzhyKAVqKJLsDBwCnNs12jYxppI8E3g58GGAqnq0qu7HNjHulgBbJVkCbE3nebq2iQEZ7jQvkkwAewFXA8+tqg3QCYDAc4ZYmhbe2cBpwONd42wT4+v5wEbgvKar/twk22CbGFtVdSfwl8AdwAbggar6AraJgRnuNOeSbAt8CnhrVf1w2PVoeJIcCtxbVauHXYtGxhJgb+BDVbUX8DB2t4215lq6I4BdgWXANkmOHW5Vi5vhTnMqyeZ0gt3HquqyZvQ9SXZspu8I3Dus+rTgXgYcnmQtcDHwyiQXYpsYZ+uB9VV1dTN8KZ2wZ5sYXwcBt1XVxqr6CXAZ8OvYJgZmuNOcSRI619HcWFXv65p0OXB88/p44B8XujYNR1WdUVU7V9UEnQukr6qqY7FNjK2quhtYl2RlM+pA4AZsE+PsDmC/JFs3nyMH0rlm2zYxIL+hQnMmyf7AV4Bv87Prq86kc93dJ4AVdP6I31hV3x9KkRqaJAcAp1bVoUmejW1ibCXZk84NNlsAtwJvpnOywTYxppK8GziKzlMX1gAnANtimxiI4U6SJKlF7JaVJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U7Sk5Jcn+SAYdcxF0ZpX5KsTXLQHK3rH5K8Z6G3K2nxMNxJelJV/XJVfQkWVzCYrNbufZGkcWK4kzTSkiwZdg2StJgY7qSWSVJJXtg1/GQ3XnOG69Qk30ryQJJLkmzZNe/aJAcl+SiwAvhskoeSnDbFtl6U5EtJ7m+6QQ/vWdcZSW5I8oMk5/Vsa1mSTyXZmOS2JCf3LPuOJN8CHk6yJMnpSW5J8mCzztc2805aa/fZvD72e+8ka5p1f7KZPmnXZ1PXnc28NyU5sGva8iSXNfv0H0k+0LXontNsf7pjsVeSa5vtXQJs2bXOaX/fk9Q+5XYmmfe9ST7dNfwXSa5MsvlUy0gaDYY7afz8FnAwsCvwK8Dv9s5QVW8C7gAOq6ptq+q9vfM0H/KfBb4APAf4Q+BjSVZ2zfY7wKuAFwC/CLyzWfZpzbLfBHYCDgTemuRVXcseAxwCbFdVjwG3AL8B/ALwbuDCJDv2U+t0+51kC+DTwD8A2wMXAa+dbAXNvp0EvLSqntHs29pm2mbAFcDtwESzXxf3sf0pj0VT22eAjza1fRJ4/RT7N60+j3m3PwdekWTPJH/Q1P66qvrJINuXtHAMd9L4eX9V3VVV36fzYb/ngOvZD9gWOKuqHq2qq+iEm2O65vlAVa1rtvWnXdNeCuxQVX/SLHsr8PfA0T11rquqHwNU1Sebuh+vqkuA7wH7zqLeqfZ7P2BJM/0nVXUZ8PUp1vFT4OnA7kk2r6q1VXVLM21fYBnw9qp6uKoeqap/62P70x2L/YDNgbOb2i4FrpnFPnfr55g/qar+AzgbuAA4A3hNVT0AkOQVSSYGrEPSPDPcSePn7q7XP6IT0AaxDFhXVY93jbudzlmhJ6zrmbaseb0LsKzpzr0/yf3AmcBzp1iWJMcl+UbX/HsAS2dR71T7vQy4s6pqqm0/oapuBt4K/DFwb5KLkzyxT8uB25uzjLPZ/nTHYrLabp9yD6fXzzHvtQZ4MXBGVXUfk98DMmAdkuaZ4U5qnx8BW3cNP2/A9dQM0+8CljfdfU9YAdzZNby8Z9pdzet1wG1VtV3XzzOq6jWTbT/JLnTOMp0EPLuqtgOu42cBY6Zap7MB2ClJd1hZPtXMVfXxqtqfTlgqOt2XT+zTisz+BpDpjsVkta3oWb7f33c/x/xJSV4MfAg4n06Ye2L84cBhwHlJjpvNjkpaGIY7qX2+Afx2ks2SHAz85oDruQd4/jTTrwYeBk5Lsnk6z5Q7jJ+/zuwtSXZOsj2ds0SXNOO/DvywuTlhq6bWPZK8dIptbUMnSG0ESPJmOmfu+q11Ol+j0916UnPjxhFM0d2bZGWSVyZ5OvAI8ONm2Sf2aQNwVpJtkmyZ5GV9bH+6Y/E14DHg5Ka2101S2zfo7/fd9zFPshOdruM/AP4b8OL87JmBVwBrquqAqrqgj/2TtMAMd1L7nEInZN1P54aGzwy4nj8D3tl04Z3aO7GqHgUOB14N3Af8LXBcVX2na7aP07nh4tbm5z3Nsj9tatwTuK1Z/lw6N0s8RVXdAPwVnbBzD52uwq/2W+t0mv14HfD7dI7ZsXQCzH9OMvvTgbOaeu+mcyPJmT379EI6N3isB47qY/tTHouu2n4X+EGzvst6VtHX77vfY57kmcD/Bd5XVZdX1Y+Av6BzzSTN/t00035JGp78/KUckjQ3kqwFTqiqfx12LbOV5Grg76rqvGHXMmqSHAlMVNXZQy5F0hQ8cydp7CX5zSTPa7o+j6fzqJLPD7uuEfVd4IQkZw+7EEmT88nvkgQrgU/QuYP1FuANVbVhuCWNpqaLfI8ZZ5Q0NHbLSpIktciM3bJJPpLk3iTXLURBkiRJGtyMZ+6SvBx4CLigqvo6Fb906dKamJjY9OokSZL0FKtXr76vqnaYbNqM19xV1Zdn+zUzExMTrFq1ajaLSJIkqU9Jpvy2mjm7oSLJicCJACtW/OwB6hOnf+4p864965Bpp002frppMy0z6hbqGA17faNQw7iub66Nwj5Nt7655t+N6xu1Glzf6NcwX+vrx5w9CqWqzqmqfapqnx12mPQsoSRJkuaZz7mTJElqEcOdJElSi/TzKJSL6Hyf48ok65P8/vyXJUmSpEH0c7fsMQtRiCRJkjad3bKSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSi/QV7pIcnOSmJDcnOX2+i5IkSdJgZgx3STYDPgi8GtgdOCbJ7vNdmCRJkmavnzN3+wI3V9WtVfUocDFwxPyWJUmSpEGkqqafIXkDcHBVndAMvwn41ao6qWe+E4ETm8GVwE1zX64WqaXAfcMuQiPFNqFetgn1sk1Mb5eq2mGyCUv6WDiTjHtKIqyqc4BzZlmYxkCSVVW1z7Dr0OiwTaiXbUK9bBOD66dbdj2wvGt4Z+Cu+SlHkiRJm6KfcHcNsFuSXZNsARwNXD6/ZUmSJGkQM3bLVtVjSU4C/hnYDPhIVV0/75WpTeyuVy/bhHrZJtTLNjGgGW+okCRJ0uLhN1RIkiS1iOFOkiSpRQx3mjNJlif5YpIbk1yf5JRm/PZJ/iXJ95p/nzXsWrWwkmyWZE2SK5ph28QYS7JdkkuTfKd5v/g128R4S/JHzefGdUkuSrKlbWJwhjvNpceAt1XVi4D9gLc0X1V3OnBlVe0GXNkMa7ycAtzYNWybGG9/DXy+qn4JeAmdtmGbGFNJdgJOBvapqj3o3Lx5NLaJgRnuNGeqakNVXdu8fpDOG/ZOdL6u7vxmtvOBI4dSoIYiyc7AIcC5XaNtE2MqyTOBlwMfBqiqR6vqfmwT424JsFWSJcDWdJ6na5sYkOFO8yLJBLAXcDXw3KraAJ0ACDxniKVp4Z0NnAY83jXONjG+ng9sBM5ruurPTbINtomxVVV3An8J3AFsAB6oqi9gmxiY4U5zLsm2wKeAt1bVD4ddj4YnyaHAvVW1eti1aGQsAfYGPlRVewEPY3fbWGuupTsC2BVYBmyT5NjhVrW4Ge40p5JsTifYfayqLmtG35Nkx2b6jsC9w6pPC+5lwOFJ1gIXA69MciG2iXG2HlhfVVc3w5fSCXu2ifF1EHBbVW2sqp8AlwG/jm1iYIY7zZkkoXMdzY1V9b6uSZcDxzevjwf+caFr03BU1RlVtXNVTdC5QPqqqjoW28TYqqq7gXVJVjajDgRuwDYxzu4A9kuydfM5ciCda7ZtEwPyGyo0Z5LsD3wF+DY/u77qTDrX3X0CWEHnj/iNVfX9oRSpoUlyAHBqVR2a5NnYJsZWkj3p3GCzBXAr8GY6JxtsE2MqybuBo+g8dWENcAKwLbaJgRjuJEmSWsRuWUmSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJM27JNcnOWDYdUjSODDcSWMoyW5JHkly4UJsr6p+uaq+1M+8SdYmOWieSxoJc7Gv43S8JPXHcCeNpw8C1wy7CA0uyZJh1yBpNBnupDGT5GjgfuDKGeZbm+TtSb6V5OEkH07y3CT/lOTBJP+a5Fld85+e5JZm2g1JXtuzroO6Xp/arPeBJJck2bKZ9lFgBfDZJA8lOa0ZvyzJp5JsTHJbkpN71j3p+vqsa7pllye5rNnufyT5wEz1THIc35Hkzmb7NyU5cKp9na7WrnrfkeRbwMNJLppkHZNur2c92yb5aZIdu8btkWRDkmdM2SgkLQ5V5Y8//ozJD/BM4LvAcuCPgQunmXct8O/Ac4GdgHuBa4G9gKcDVwHv6pr/jcAyOv9pPAp4GNixa10Hdb3+ejPv9sCNwB/0bPegruGnAauB/wlsATwfuBV4VZ/rm6muSZcFNgO+CfxvYBtgS2D/merpOYYrgXXAsmZ4AnjBNPs6Za1d83+j+f1tNcmxnXZ7PbVdDxzSNXwF8IfDbqP++OPPpv945k4aL/8L+HBVretz/r+pqnuq6k7gK8DVVbWmqv4T+DSdoAdAVX2yqu6qqser6hLge8C+U6z3/c283wc+C+w5TQ0vBXaoqj+pqker6lbg74Gj+1lfH3VNtey+dILW26vq4ap6pKr+rc96nvBTOkF49ySbV9Xaqrplqh3t8xi+v6rWVdWPN3F71wB7AyR5ObA78H+a4VckmZiqTkmjzXAnjYkkewIH0TkT1a97ul7/eJLhbbvWf1ySbyS5P8n9wB7A0inWe3fX6x91r2cSuwDLnlhvs+4z6ZxRnHF9fdQ11bLLgdur6rEB6gGgqm4G3krnLOm9SS5OsmyqHe3zGE4ZzGe5vSfDHfBe4H9U1aPN8O8BmWo7kkab4U4aHwfQ6aa7I8ndwKnA65Ncu6krTrILnbNXJwHPrqrtgOsYLCBUz/A64Laq2q7r5xlV9Zp5rmsdsGKSGxdmVU9Vfbyq9qcTCgv488n2dRa19h6fnxueYXvdrgH2TvJ6YCvgoqaOw4HDgPOSHDfFspJGmOFOGh/nAC+g0+24J/B3wOeAV83BurehEyQ2AiR5M52zToO4h851bE/4OvDD5kaBrZJs1lz8/9J5ruvrwAbgrCTbJNkyyctmU0+SlUlemeTpwCN0znb+dIp9HbTWJ9fRx/a6fRN4HvBXwOlV9Xgz/gpgTVUdUFUX9LF9SSPGcCeNiar6UVXd/cQP8BDwSFVtnIN130AnJHyNTth4MfDVAVf3Z8A7m67JU6vqp3TOJO0J3AbcB5wL/MJ81tW13RcCdwDrgaNmWc/TgbOaee4GnkOnC/cp+wq8ZsBau9dx1Azb696//wS+Daytqn/qmvRC4KY+titpRKWq9wy/JKntkmwB3Az8VlX9e9f4I4GJqjp7SKVJ2kSeuZOk8fQu4Kvdwa7xXeCEJGcvfEmS5oJn7iRpjCTZG/gi8C3gtVV135BLkjTHDHeSJEktMmO3bJKPJLk3yXULUZAkSZIGN+OZu+bJ5Q8BF1RVX48QWLp0aU1MTGx6dZIkSXqK1atX31dVO0w2rffhnE9RVV+e7dfQTExMsGrVqtksIkmSpD4luX2qaTOGu1ls5ETgRIAVK1bM1WrnxMTpn3vKuLVnHTIy65Pm27i12XHbX0nqNmePQqmqc6pqn6raZ4cdJj1LKEmSpHnmc+4kSZJaxHAnSZLUIv08CuUiOt91uDLJ+iS/P/9lSZIkaRD93C17zEIUIkmSpE1nt6wkSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFukr3CU5OMlNSW5Ocvp8FyVJkqTBzBjukmwGfBB4NbA7cEyS3ee7MEmSJM1eP2fu9gVurqpbq+pR4GLgiPktS5IkSYNIVU0/Q/IG4OCqOqEZfhPwq1V1Us98JwInNoMrgZvmvlwtUkuB+4ZdhEaKbUK9bBPqZZuY3i5VtcNkE5b0sXAmGfeURFhV5wDnzLIwjYEkq6pqn2HXodFhm1Av24R62SYG10+37HpgedfwzsBd81OOJEmSNkU/4e4aYLckuybZAjgauHx+y5IkSdIgZuyWrarHkpwE/DOwGfCRqrp+3itTm9hdr162CfWyTaiXbWJAM95QIUmSpMXDb6iQJElqEcOdJElSixjuNGeSLE/yxSQ3Jrk+ySnN+O2T/EuS7zX/PmvYtWphJdksyZokVzTDtokxlmS7JJcm+U7zfvFrtonxluSPms+N65JclGRL28TgDHeaS48Bb6uqFwH7AW9pvqrudODKqtoNuLIZ1ng5Bbixa9g2Md7+Gvh8Vf0S8BI6bcM2MaaS7AScDOxTVXvQuXnzaGwTAzPcac5U1YaqurZ5/SCdN+yd6Hxd3fnNbOcDRw6lQA1Fkp2BQ4Bzu0bbJsZUkmcCLwc+DFBVj1bV/dgmxt0SYKskS4Ct6TxP1zYxIMOd5kWSCWAv4GrguVW1AToBEHjOEEvTwjsbOA14vGucbWJ8PR/YCJzXdNWfm2QbbBNjq6ruBP4SuAPYADxQVV/ANjEww53mXJJtgU8Bb62qHw67Hg1PkkOBe6tq9bBr0chYAuwNfKiq9gIexu62sdZcS3cEsCuwDNgmybHDrWpxM9xpTiXZnE6w+1hVXdaMvifJjs30HYF7h1WfFtzLgMOTrAUuBl6Z5EJsE+NsPbC+qq5uhi+lE/ZsE+PrIOC2qtpYVT8BLgN+HdvEwAx3mjNJQuc6mhur6n1dky4Hjm9eHw/840LXpuGoqjOqaueqmqBzgfRVVXUstomxVVV3A+uSrGxGHQjcgG1inN0B7Jdk6+Zz5EA612zbJgbkN1RoziTZH/gK8G1+dn3VmXSuu/sEsILOH/Ebq+r7QylSQ5PkAODUqjo0ybOxTYytJHvSucFmC+BW4M10TjbYJsZUkncDR9F56sIa4ARgW2wTAzHcSZIktYjdspIkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktcj/B4n+nkmy4GsgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def maintenance_planning_gdp(c, T, M, P):\n",
    "    m = pyo.ConcreteModel()\n",
    "\n",
    "    m.T = pyo.RangeSet(1, T)\n",
    "    m.Y = pyo.RangeSet(1, T - M + 1)\n",
    "    m.S = pyo.RangeSet(0, M - 1)\n",
    "\n",
    "    m.c = pyo.Param(m.T, initialize = c)\n",
    "    m.x = pyo.Var(m.T, domain=pyo.Binary)\n",
    "    m.y = pyo.Var(m.T, domain=pyo.Binary)\n",
    "\n",
    "    # objective\n",
    "    m.profit = pyo.Objective(expr = sum(m.c[t]*m.x[t] for t in m.T), sense=pyo.maximize)\n",
    "\n",
    "    # required number P of maintenance starts\n",
    "    m.sumy = pyo.Constraint(expr = sum(m.y[t] for t in m.Y) == P)\n",
    "\n",
    "    # no more than one maintenance start in the period of length M\n",
    "    m.sprd = pyo.Constraint(m.Y, rule = lambda m, t: sum(m.y[t+s] for s in m.S) <= 1)\n",
    "\n",
    "    # disjunctive constraints\n",
    "    m.disj = gdp.Disjunction(m.Y, rule = lambda m, t: [m.y[t]==0, sum(m.x[t+s] for s in m.S)==0])\n",
    "\n",
    "    # transformation and soluton\n",
    "    pyo.TransformationFactory('gdp.hull').apply_to(m)\n",
    "    \n",
    "    return m\n",
    "\n",
    "m = maintenance_planning_gdp(c, T, M, P)\n",
    "pyo.SolverFactory('cbc').solve(m).write()\n",
    "plot_schedule(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gcy3yD4J_6OC"
   },
   "source": [
    "## Ramping constraints\n",
    "\n",
    "Prior to maintenance shutdown, a large processing unit may take some time to safely ramp down from full production. And then require more time to safely ramp back up to full production following maintenance. To provide for ramp-down and ramp-up periods, we modify the problem formation in the following ways.\n",
    "\n",
    "* The variable denoting unit operation, $x_t$ is changed from a binary variable to a continuous variable $0 \\leq x_t \\leq 1$ denoting the fraction of total capacity at which the unit is operating on day $t$.\n",
    "\n",
    "* Two new variable sequences, $0 \\leq u_t^+ \\leq u_t^{+,\\max}$ and $0\\leq u_t^- \\leq u_t^{-,\\max}$, are introduced which denote the fraction increase or decrease in unit capacity to completed on day $t$.\n",
    "\n",
    "* An additional sequence of equality constraints is introduced relating $x_t$ to $u_t^+$ and $u_t^-$.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "x_{t} & = x_{t-1} + u^+_t - u^-_t\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "We begin the Pyomo model by specifying the constraints, then modifying the Big-M formulation to add the features described above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "9fRtZc5nCInh"
   },
   "outputs": [],
   "source": [
    "upos_max = 0.3334\n",
    "uneg_max = 0.5000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3307,
     "status": "ok",
     "timestamp": 1557947674808,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "ESHZWcEaA9gP",
    "outputId": "030e96a8-1ae1-47f5-f4c6-d83378818838"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAEYCAYAAAA+gNBwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfZklEQVR4nO3de7RkZX2n8ecrDXLTIDYqDd0cVNKRYAQWEhKJQWGWKFdvARKEmDCsrJGAGRGB5YwxY1aISRxiNGYISkAUUESDODEmoEvjGKQbvHAR5dLQDQ00URBQgshv/qjdWBbnUqf6nFN1dj2ftc7q2vff3uc9Vd/e7967UlVIkiSpHZ427AIkSZI0dwx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFO0tAk+Yck7+ljvhuSHNA1vCbJQfNZWz+SrExyXZKHkpzcW6ckDcOSYRcgSTOpql8edg1TOA34UlXt1Qy/v3tikjXACVX1rwtdmKTx5Zk7SZpEkn7+87sLcMN81yJJs2G4k7RgkuyV5NqmG/MSYMuuaacnubWZdmOS13ZNm7QbNsnbk3yqZ9zfJDl7iu2vSXJGs/4fJDkvyZY909+R5FvAI0mWJHlRki8leaDpdj28mfcq4BXAB5I8nOQXu+tM8lFgBfDZZvppU9R0TJLVSR5s9v+APg+nJE3KcCdpQSTZAvgM8FFge+CTwOu7ZrkV+A3gF4B3Axcm2XGG1V4IHJxku2YbS4Cjmm1M5XeAVwEvAH4ReGfP9GOAQ4DtgACfBb4APAf4Q+BjSVZW1SuBrwAnVdW2VfXd7pVU1ZuAO4HDmunv7S0kydua7f9X4FnAkcCaGfZZkqZluJO0UPYDNgfOrqqfVNWlwDUbJ1bVJ6vq7qp6oqouAb4H7DvdCqtqPfBl4I3NqIOB+6tq9TSLfaCq1lbV94E/pRPmur2/mf7jpuZtgbOq6rGqugq4YpJlZi3JDsC7gN+uqmub/f52Va1ppr8iycSmbkfS+DHcSVooy4C7qqq6xt2x8UWS45J8o+n+fADYA1jax3rPB45tXh/L9GftANb2bH/ZNNOXAWur6omeZXbqo66ZHAR8u6q+OcX036Nz5lCSZsVwJ2mhrAd2StIdWFYAJNkF+HvgJODZVbUdcD39hZvPAL+SZA/gUOBjM8y/vGf7d/dM7w6fdwPLk3S/V64A7uqjrt519doeeGCyCc11fYcB5yU5rs9tSRJguJO0cL4GPA6c3Nyo8Dp+1u26DZ0gtAEgyZvpnLmbUVU9ClwKfBz4elXdOcMib0myc5LtgTOBS6aZ92rgEeC0JJs3NzscBlzcT23AvcDzp5h2HbB/kpekY7ckL2qmXQFcV1UHVNUFfW5LkgDDnaQFUlWPAa8Dfhf4AZ0bHy5rpt0I/BWdAHgv8GLgq7NY/fnNMjN1yUInBH4BuK35mfIhyk3NhwOvBu4H/hY4rqq+02ddfwa8s+lqPrVn3f+v2fYVwEPAp4GtmskvBG7ucxuS9HPy85e/SNLik2QF8B3geVX1w2nmW8MieKhwkiOBiao6e8ilSFqEPHMnaVFrrof778DF0wW7Rea7wAlTPa9Pkqbj149JWrSSbEOnG/cOOo9BaYWmm7qvaw4lqZfdspIkSS0yY7dsko8kuS/J9QtRkCRJkgY345m7JC8HHgYuqKq+ugmWLl1aExMTm16dJEmSnmL16tX3V9UOk02b8Zq7qvrybL8CZ2JiglWrVs1mEUmSJPUpyR1TTZuzu2WTnJhkVZJVGzZsmKvVSpIkaRbm7G7ZqjoHOAdgn3328S4NzdrE6Z+bdPyasw5Z4Eqk4Zjsb2Bj+59umiR181Eo2iR+4EiS2qQNn2tDDXdtOICShsP3D0mjYtR6nmYMd0kuAg4AliZZB7yrqj4834VJktrNgC7Nj37ulj1mIQqRJEnSpvO7ZSVJklrEGypGgF0TkiRprnjmTpIkqUU8cydJWjSmuyvRXhCpw3An9VioD4hRu3VekkaZ4b1/dstKkiS1iGfuJGkGU50x8Oyr5oLtSHPNcCdJktSHxdI13Jpw5/98pOHz71CaW4slTGi0tCbcSRptfkgtDuMW0Mdtf0eBx3z+Ge4kSVoA/gdHvebrel7DnaSh80NPkuaOj0KRJElqEcOdJElSi4xkt6wXW44Wu8ykxWkh/3Z9nxgdfobOrO3HaCTDnTSK2v5mMBt+kGu+zXUbs81qnBjuFimDhqTFwPcqzYYhfG4Y7gT4BizNNT+k1BZ+Piw+iy7c+YYpTc434PHk733hjcox9/NQU/FuWUmSpBYx3EmSJLXIouuWHXWeJtd8s41JkqZjuGshP/wlzQffW6TFwXA3hVG5YFYzm68vXpYkaTEa+3A36v8THfX6ND/8vc8PA780t0b9vWrU65svYx/uND+m+4Ma1z82SYuL71VarAx30jzzA0KS5s+ov8cOoz7DnSRJGvmQpP4Z7iQJP9gktYfhbgB+CEijzb9RSePMcLdARuUuPT/0JGl8+RkwHuY93NmQNIpslzPzGEnaFL6HDI9n7qQxNipnlKW2MNBoFBjutCj4hjmeDJ+SNHtPG3YBkiRJmjtjcebOsz4ad/4NjBZ/H5LmU19n7pIcnOTmJLckOX2+i5IkSdJgZjxzl2Qz4IPAfwHWAdckubyqbpzv4iRJ88Ozh1J79dMtuy9wS1XdBpDkYuAIwHAnNfyglCSNilTV9DMkbwAOrqoTmuE3Ab9aVSf1zHcicGIzuBK4ee7L1SK1FLh/2EVopNgm1Ms2oV62ientUlU7TDahnzN3mWTcUxJhVZ0DnDPLwjQGkqyqqn2GXYdGh21CvWwT6mWbGFw/N1SsA5Z3De8M3D0/5UiSJGlT9BPurgF2S7Jrki2Ao4HL57csSZIkDWLGbtmqejzJScA/A5sBH6mqG+a9MrWJ3fXqZZtQL9uEetkmBjTjDRWSJElaPPz6MUmSpBYx3EmSJLWI4U5zJsnyJF9MclOSG5Kc0ozfPsm/JPle8++zhl2rFlaSzZJcl+SKZtg2McaSbJfk0iTfad4vfs02Md6S/FHzuXF9kouSbGmbGJzhTnPpceBtVfUiYD/gLUl2B04Hrqyq3YArm2GNl1OAm7qGbRPj7a+Bz1fVLwEvodM2bBNjKslOwMnAPlW1B52bN4/GNjEww53mTFWtr6prm9cP0XnD3onO19Wd38x2PnDkUArUUCTZGTgEOLdrtG1iTCV5JvBy4MMAVfVYVT2AbWLcLQG2SrIE2JrO83RtEwMy3GleJJkA9gKuBp5bVeuhEwCB5wyxNC28s4HTgCe6xtkmxtfzgQ3AeU1X/blJtsE2Mbaq6i7gL4E7gfXAg1X1BWwTAzPcac4l2Rb4FPDWqvrhsOvR8CQ5FLivqlYPuxaNjCXA3sCHqmov4BHsbhtrzbV0RwC7AsuAbZIcO9yqFjfDneZUks3pBLuPVdVlzeh7k+zYTN8RuG9Y9WnBvQw4PMka4GLglUkuxDYxztYB66rq6mb4UjphzzYxvg4Cbq+qDVX1E+Ay4NexTQzMcKc5kyR0rqO5qare1zXpcuD45vXxwD8udG0ajqo6o6p2rqoJOhdIX1VVx2KbGFtVdQ+wNsnKZtSBwI3YJsbZncB+SbZuPkcOpHPNtm1iQH5DheZMkv2BrwDf5mfXV51J57q7TwAr6PwRv7Gqvj+UIjU0SQ4ATq2qQ5M8G9vE2EqyJ50bbLYAbgPeTOdkg21iTCV5N3AUnacuXAecAGyLbWIghjtJkqQWsVtWkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJT0pyQ5IDhl3HXBilfUmyJslBc7Suf0jynoXerqTFw3An6UlV9ctV9SVYXMFgslq790WSxonhTtJIS7Jk2DVI0mJiuJNaJkkleWHX8JPdeM0ZrlOTfCvJg0kuSbJl17xrkhyU5KPACuCzSR5OctoU23pRki8leaDpBj28Z11nJLkxyQ+SnNezrWVJPpVkQ5Lbk5zcs+w7knwLeCTJkiSnJ7k1yUPNOl/bzDtprd1n8/rY772TXNes+5PN9Em7Ppu67mrmvTnJgV3Tlie5rNmn/0jyga5F95xm+9Mdi72SXNts7xJgy651Tvv7nqT2KbczybzvTfLpruG/SHJlks2nWkbSaDDcSePnt4CDgV2BXwF+t3eGqnoTcCdwWFVtW1Xv7Z2n+ZD/LPAF4DnAHwIfS7Kya7bfAV4FvAD4ReCdzbJPa5b9JrATcCDw1iSv6lr2GOAQYLuqehy4FfgN4BeAdwMXJtmxn1qn2+8kWwCfBv4B2B64CHjtZCto9u0k4KVV9Yxm39Y00zYDrgDuACaa/bq4j+1PeSya2j4DfLSp7ZPA66fYv2n1ecy7/TnwiiR7JvmDpvbXVdVPBtm+pIVjuJPGz/ur6u6q+j6dD/s9B1zPfsC2wFlV9VhVXUUn3BzTNc8Hqmpts60/7Zr2UmCHqvqTZtnbgL8Hju6pc21V/Rigqj7Z1P1EVV0CfA/Ydxb1TrXf+wFLmuk/qarLgK9PsY6fAk8Hdk+yeVWtqapbm2n7AsuAt1fVI1X1aFX9Wx/bn+5Y7AdsDpzd1HYpcM0s9rlbP8f8SVX1H8DZwAXAGcBrqupBgCSvSDIxYB2S5pnhTho/93S9/hGdgDaIZcDaqnqia9wddM4KbbS2Z9qy5vUuwLKmO/eBJA8AZwLPnWJZkhyX5Btd8+8BLJ1FvVPt9zLgrqqqqba9UVXdArwV+GPgviQXJ9m4T8uBO5qzjLPZ/nTHYrLa7phyD6fXzzHvdR3wYuCMquo+Jr8HZMA6JM0zw53UPj8Ctu4aft6A66kZpt8NLG+6+zZaAdzVNby8Z9rdzeu1wO1VtV3XzzOq6jWTbT/JLnTOMp0EPLuqtgOu52cBY6Zap7Me2ClJd1hZPtXMVfXxqtqfTlgqOt2XG/dpRWZ/A8h0x2Ky2lb0LN/v77ufY/6kJC8GPgScTyfMbRx/OHAYcF6S42azo5IWhuFOap9vAL+dZLMkBwO/OeB67gWeP830q4FHgNOSbJ7OM+UO4+evM3tLkp2TbE/nLNElzfivAz9sbk7Yqql1jyQvnWJb29AJUhsAkryZzpm7fmudztfodLee1Ny4cQRTdPcmWZnklUmeDjwK/LhZduM+rQfOSrJNki2TvKyP7U93LL4GPA6c3NT2uklq+wb9/b77PuZJdqLTdfwHwH8DXpyfPTPwCuC6qjqgqi7oY/8kLTDDndQ+p9AJWQ/QuaHhMwOu58+AdzZdeKf2Tqyqx4DDgVcD9wN/CxxXVd/pmu3jdG64uK35eU+z7E+bGvcEbm+WP5fOzRJPUVU3An9FJ+zcS6er8Kv91jqdZj9eB/w+nWN2LJ0A85+TzP504Kym3nvo3EhyZs8+vZDODR7rgKP62P6Ux6Krtt8FftCs77KeVfT1++73mCd5JvB/gfdV1eVV9SPgL+hcM0mzfzfPtF+Shic/fymHJM2NJGuAE6rqX4ddy2wluRr4u6o6b9i1jJokRwITVXX2kEuRNAXP3Ekae0l+M8nzmq7P4+k8quTzw65rRH0XOCHJ2cMuRNLkfPK7JMFK4BN07mC9FXhDVa0fbkmjqeki32PGGSUNjd2ykiRJLTJjt2ySjyS5L8n1C1GQJEmSBjfjmbskLwceBi6oqr5OxS9durQmJiY2vTpJkiQ9xerVq++vqh0mmzbjNXdV9eXZfs3MxMQEq1atms0ikiRJ6lOSKb+tZs5uqEhyInAiwIoVP3uA+sTpn3vKvGvOOmTaaZONn27aoMvMdQ2ub/5q0GhpU7sc9fWNQg2ub/RqcH2jX8N8ra8fc/YolKo6p6r2qap9dthh0rOEkiRJmmc+506SJKlFDHeSJEkt0s+jUC6i832OK5OsS/L781+WJEmSBtHP3bLHLEQhkiRJ2nR2y0qSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLXIjN9QIWk8TZz+uaeMW3PWIUOoRJI0G565kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iI9C0aLgYzkkSeqPZ+4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRn3MnjbHJnh8IPkNQkhYzz9xJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRfoKd0kOTnJzkluSnD7fRUmSJGkwM4a7JJsBHwReDewOHJNk9/kuTJIkSbPXz5m7fYFbquq2qnoMuBg4Yn7LkiRJ0iBSVdPPkLwBOLiqTmiG3wT8alWd1DPficCJzeBK4Oa5L1eL1FLg/mEXoZFim1Av24R62Samt0tV7TDZhCV9LJxJxj0lEVbVOcA5syxMYyDJqqraZ9h1aHTYJtTLNqFetonB9dMtuw5Y3jW8M3D3/JQjSZKkTdFPuLsG2C3Jrkm2AI4GLp/fsiRJkjSIGbtlq+rxJCcB/wxsBnykqm6Y98rUJnbXq5dtQr1sE+plmxjQjDdUSJIkafHwGyokSZJaxHAnSZLUIoY7zZkky5N8MclNSW5Ickozfvsk/5Lke82/zxp2rVpYSTZLcl2SK5ph28QYS7JdkkuTfKd5v/g128R4S/JHzefG9UkuSrKlbWJwhjvNpceBt1XVi4D9gLc0X1V3OnBlVe0GXNkMa7ycAtzUNWybGG9/DXy+qn4JeAmdtmGbGFNJdgJOBvapqj3o3Lx5NLaJgRnuNGeqan1VXdu8fojOG/ZOdL6u7vxmtvOBI4dSoIYiyc7AIcC5XaNtE2MqyTOBlwMfBqiqx6rqAWwT424JsFWSJcDWdJ6na5sYkOFO8yLJBLAXcDXw3KpaD50ACDxniKVp4Z0NnAY80TXONjG+ng9sAM5ruurPTbINtomxVVV3AX8J3AmsBx6sqi9gmxiY4U5zLsm2wKeAt1bVD4ddj4YnyaHAfVW1eti1aGQsAfYGPlRVewGPYHfbWGuupTsC2BVYBmyT5NjhVrW4Ge40p5JsTifYfayqLmtG35tkx2b6jsB9w6pPC+5lwOFJ1gAXA69MciG2iXG2DlhXVVc3w5fSCXu2ifF1EHB7VW2oqp8AlwG/jm1iYIY7zZkkoXMdzU1V9b6uSZcDxzevjwf+caFr03BU1RlVtXNVTdC5QPqqqjoW28TYqqp7gLVJVjajDgRuxDYxzu4E9kuydfM5ciCda7ZtEwPyGyo0Z5LsD3wF+DY/u77qTDrX3X0CWEHnj/iNVfX9oRSpoUlyAHBqVR2a5NnYJsZWkj3p3GCzBXAb8GY6JxtsE2MqybuBo+g8deE64ARgW2wTAzHcSZIktYjdspIkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZp3SW5IcsCw65CkcWC4k8ZQkt2SPJrkwoXYXlX9clV9qZ95k6xJctA8lzQS5mJfx+l4SeqP4U4aTx8Erhl2ERpckiXDrkHSaDLcSWMmydHAA8CVM8y3Jsnbk3wrySNJPpzkuUn+KclDSf41ybO65j89ya3NtBuTvLZnXQd1vT61We+DSS5JsmUz7aPACuCzSR5OclozflmSTyXZkOT2JCf3rHvS9fVZ13TLLk9yWbPd/0jygZnqmeQ4viPJXc32b05y4FT7Ol2tXfW+I8m3gEeSXDTJOibdXs96tk3y0yQ7do3bI8n6JM+YslFIWhyqyh9//BmTH+CZwHeB5cAfAxdOM+8a4N+B5wI7AfcB1wJ7AU8HrgLe1TX/G4FldP7TeBTwCLBj17oO6nr99Wbe7YGbgD/o2e5BXcNPA1YD/xPYAng+cBvwqj7XN1Ndky4LbAZ8E/jfwDbAlsD+M9XTcwxXAmuBZc3wBPCCafZ1ylq75v9G8/vbapJjO+32emq7ATika/gK4A+H3Ub98cefTf/xzJ00Xv4X8OGqWtvn/H9TVfdW1V3AV4Crq+q6qvpP4NN0gh4AVfXJqrq7qp6oqkuA7wH7TrHe9zfzfh/4LLDnNDW8FNihqv6kqh6rqtuAvweO7md9fdQ11bL70glab6+qR6rq0ar6tz7r2eindILw7kk2r6o1VXXrVDva5zF8f1Wtraofb+L2rgH2BkjycmB34P80w69IMjFVnZJGm+FOGhNJ9gQOonMmql/3dr3+8STD23at/7gk30jyQJIHgD2ApVOs956u1z/qXs8kdgGWbVxvs+4z6ZxRnHF9fdQ11bLLgTuq6vEB6gGgqm4B3krnLOl9SS5OsmyqHe3zGE4ZzGe5vSfDHfBe4H9U1WPN8O8BmWo7kkab4U4aHwfQ6aa7M8k9wKnA65Ncu6krTrILnbNXJwHPrqrtgOsZLCBUz/Ba4Paq2q7r5xlV9Zp5rmstsGKSGxdmVU9Vfbyq9qcTCgv488n2dRa19h6fnxueYXvdrgH2TvJ6YCvgoqaOw4HDgPOSHDfFspJGmOFOGh/nAC+g0+24J/B3wOeAV83BurehEyQ2ACR5M52zToO4l851bBt9Hfhhc6PAVkk2ay7+f+k81/V1YD1wVpJtkmyZ5GWzqSfJyiSvTPJ04FE6Zzt/OsW+Dlrrk+voY3vdvgk8D/gr4PSqeqIZfwVwXVUdUFUX9LF9SSPGcCeNiar6UVXds/EHeBh4tKo2zMG6b6QTEr5GJ2y8GPjqgKv7M+CdTdfkqVX1UzpnkvYEbgfuB84FfmE+6+ra7guBO4F1wFGzrOfpwFnNPPcAz6HThfuUfQVeM2Ct3es4aobtde/ffwLfBtZU1T91TXohcHMf25U0olLVe4ZfktR2SbYAbgF+q6r+vWv8kcBEVZ09pNIkbSLP3EnSeHoX8NXuYNf4LnBCkrMXviRJc8Ezd5I0RpLsDXwR+Bbw2qq6f8glSZpjhjtJkqQWmbFbNslHktyX5PqFKEiSJEmDm/HMXfPk8oeBC6qqr0cILF26tCYmJja9OkmSJD3F6tWr76+qHSab1vtwzqeoqi/P9mtoJiYmWLVq1WwWkSRJUp+S3DHVtBnD3Sw2ciJwIsCKFSvmarXSjCZO/9xTxq056xDXJ0kaS3P2KJSqOqeq9qmqfXbYYdKzhJIkSZpnPudOkiSpRQx3kiRJLdLPo1AuovNdhyuTrEvy+/NfliRJkgbRz92yxyxEIZIkSdp0dstKkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJapG+wl2Sg5PcnOSWJKfPd1GSJEkazIzhLslmwAeBVwO7A8ck2X2+C5MkSdLs9XPmbl/glqq6raoeAy4GjpjfsiRJkjSIVNX0MyRvAA6uqhOa4TcBv1pVJ/XMdyJwYjO4Erh57svVIrUUuH/YRWik2CbUyzahXraJ6e1SVTtMNmFJHwtnknFPSYRVdQ5wziwL0xhIsqqq9hl2HRodtgn1sk2ol21icP10y64DlncN7wzcPT/lSJIkaVP0E+6uAXZLsmuSLYCjgcvntyxJkiQNYsZu2ap6PMlJwD8DmwEfqaob5r0ytYnd9eplm1Av24R62SYGNOMNFZIkSVo8/IYKSZKkFjHcSZIktYjhTnMmyfIkX0xyU5IbkpzSjN8+yb8k+V7z77OGXasWVpLNklyX5Ipm2DYxxpJsl+TSJN9p3i9+zTYx3pL8UfO5cX2Si5JsaZsYnOFOc+lx4G1V9SJgP+AtzVfVnQ5cWVW7AVc2wxovpwA3dQ3bJsbbXwOfr6pfAl5Cp23YJsZUkp2Ak4F9qmoPOjdvHo1tYmCGO82ZqlpfVdc2rx+i84a9E52vqzu/me184MihFKihSLIzcAhwbtdo28SYSvJM4OXAhwGq6rGqegDbxLhbAmyVZAmwNZ3n6domBmS407xIMgHsBVwNPLeq1kMnAALPGWJpWnhnA6cBT3SNs02Mr+cDG4Dzmq76c5Nsg21ibFXVXcBfAncC64EHq+oL2CYGZrjTnEuyLfAp4K1V9cNh16PhSXIocF9VrR52LRoZS4C9gQ9V1V7AI9jdNtaaa+mOAHYFlgHbJDl2uFUtboY7zakkm9MJdh+rqsua0fcm2bGZviNw37Dq04J7GXB4kjXAxcArk1yIbWKcrQPWVdXVzfCldMKebWJ8HQTcXlUbquonwGXAr2ObGJjhTnMmSehcR3NTVb2va9LlwPHN6+OBf1zo2jQcVXVGVe1cVRN0LpC+qqqOxTYxtqrqHmBtkpXNqAOBG7FNjLM7gf2SbN18jhxI55pt28SA/IYKzZkk+wNfAb7Nz66vOpPOdXefAFbQ+SN+Y1V9fyhFamiSHACcWlWHJnk2tomxlWRPOjfYbAHcBryZzskG28SYSvJu4Cg6T124DjgB2BbbxEAMd5IkSS1it6wkSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS3y/wEZsJ5PbI15gQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def maintenance_planning_ramp(c, T, M, P):\n",
    "    m = pyo.ConcreteModel()\n",
    "\n",
    "    m.T = pyo.RangeSet(1, T)\n",
    "    m.Y = pyo.RangeSet(1, T - M + 1)\n",
    "    m.S = pyo.RangeSet(0, M - 1)\n",
    "\n",
    "    m.c = pyo.Param(m.T, initialize = c)\n",
    "    m.x = pyo.Var(m.T, bounds=(0, 1))\n",
    "    m.y = pyo.Var(m.T, domain=pyo.Binary)\n",
    "    m.upos = pyo.Var(m.T, bounds=(0, upos_max))\n",
    "    m.uneg = pyo.Var(m.T, bounds=(0, uneg_max))\n",
    "\n",
    "    # objective\n",
    "    m.profit = pyo.Objective(expr = sum(m.c[t]*m.x[t] for t in m.T), sense=pyo.maximize)\n",
    "    \n",
    "    # ramp constraint\n",
    "    m.ramp = pyo.Constraint(m.T, rule = lambda m, t: \n",
    "         m.x[t] == m.x[t-1] + m.upos[t] - m.uneg[t] if t > 1 else pyo.Constraint.Skip)\n",
    "      \n",
    "    # required number P of maintenance starts\n",
    "    m.sumy = pyo.Constraint(expr = sum(m.y[t] for t in m.Y) == P)\n",
    "\n",
    "    # no more than one maintenance start in the period of length M\n",
    "    m.sprd = pyo.Constraint(m.Y, rule = lambda m, t: sum(m.y[t+s] for s in m.S) <= 1)\n",
    "    \n",
    "    # disjunctive constraints\n",
    "    m.disj = gdp.Disjunction(m.Y, rule = lambda m, t: [m.y[t]==0, sum(m.x[t+s] for s in m.S)==0])\n",
    "    \n",
    "    # transformation and soluton\n",
    "    pyo.TransformationFactory('gdp.hull').apply_to(m)\n",
    "    \n",
    "    return m\n",
    "  \n",
    "m = maintenance_planning_ramp(c, T, M, P)\n",
    "pyo.SolverFactory('cbc').solve(m)\n",
    "plot_schedule(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HM_k7HkbThAe"
   },
   "source": [
    "## Specifying the minimum number of operational days between maintenance periods\n",
    "\n",
    "Up to this point we have imposed no constraints on the frequency of maintenance periods. Without such constraints, particularly when ramping constraints are imposed, is that maintenance periods will be scheduled back-to-back, which is clearly not a useful result for most situations.\n",
    "\n",
    "The next revision of the model is to incorporate a requirement that $N$ operational days be scheduled between any maintenance periods. This does allow for maintenance to be postponed until the very end of the planning period. The disjunctive constraints read\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\left(y_t = 0\\right) \\lor \\left(\\sum_{s=0}^{(M + N -1) \\wedge (t + s \\leq T)}x_{t+s} = 0\\right)\\qquad \\forall t = 1, 2, \\ldots, T-M+1\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where the upper bound on the summation is needed to handle the terminal condition. \n",
    "\n",
    "Paradoxically, this is an example where the Big-M method provides a much faster solution.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{s=0}^{(M + N -1) \\wedge (t + s \\leq T)}x_{t+s} \\leq (M+N)(1-y_t) \\qquad \\forall t = 1, 2, \\ldots, T-M+1\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The following cell implements both sets of constraints. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "fEOMGHPvUJzJ"
   },
   "outputs": [],
   "source": [
    "N = 10  # minimum number of operational days between maintenance periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3249,
     "status": "ok",
     "timestamp": 1557947681109,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "m93Na6TkUqDb",
    "outputId": "cae87fae-7eba-46c3-ea3b-7e2427c4f529"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAEYCAYAAAA+gNBwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfnElEQVR4nO3de7RkZX2n8ecrDXLTIDYqDd0cVNKRYAQWEhKJQWGWKFdvARKEmDCsrJGAGRGB5YwxY1aISRxiNGYISkAUUESDODEmoEvjGKQbvHAR5dLQDQ00URBQgshv/qjdWBbnUqf6nFN1dj2ftc7q2vff3uc9Vd/e7967UlVIkiSpHZ427AIkSZI0dwx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFO0tAk+Yck7+ljvhuSHNA1vCbJQfNZWz+SrExyXZKHkpzcW6ckDcOSYRcgSTOpql8edg1TOA34UlXt1Qy/v3tikjXACVX1rwtdmKTx5Zk7SZpEkn7+87sLcMN81yJJs2G4k7RgkuyV5NqmG/MSYMuuaacnubWZdmOS13ZNm7QbNsnbk3yqZ9zfJDl7iu2vSXJGs/4fJDkvyZY909+R5FvAI0mWJHlRki8leaDpdj28mfcq4BXAB5I8nOQXu+tM8lFgBfDZZvppU9R0TJLVSR5s9v+APg+nJE3KcCdpQSTZAvgM8FFge+CTwOu7ZrkV+A3gF4B3Axcm2XGG1V4IHJxku2YbS4Cjmm1M5XeAVwEvAH4ReGfP9GOAQ4DtgACfBb4APAf4Q+BjSVZW1SuBrwAnVdW2VfXd7pVU1ZuAO4HDmunv7S0kydua7f9X4FnAkcCaGfZZkqZluJO0UPYDNgfOrqqfVNWlwDUbJ1bVJ6vq7qp6oqouAb4H7DvdCqtqPfBl4I3NqIOB+6tq9TSLfaCq1lbV94E/pRPmur2/mf7jpuZtgbOq6rGqugq4YpJlZi3JDsC7gN+uqmub/f52Va1ppr8iycSmbkfS+DHcSVooy4C7qqq6xt2x8UWS45J8o+n+fADYA1jax3rPB45tXh/L9GftANb2bH/ZNNOXAWur6omeZXbqo66ZHAR8u6q+OcX036Nz5lCSZsVwJ2mhrAd2StIdWFYAJNkF+HvgJODZVbUdcD39hZvPAL+SZA/gUOBjM8y/vGf7d/dM7w6fdwPLk3S/V64A7uqjrt519doeeGCyCc11fYcB5yU5rs9tSRJguJO0cL4GPA6c3Nyo8Dp+1u26DZ0gtAEgyZvpnLmbUVU9ClwKfBz4elXdOcMib0myc5LtgTOBS6aZ92rgEeC0JJs3NzscBlzcT23AvcDzp5h2HbB/kpekY7ckL2qmXQFcV1UHVNUFfW5LkgDDnaQFUlWPAa8Dfhf4AZ0bHy5rpt0I/BWdAHgv8GLgq7NY/fnNMjN1yUInBH4BuK35mfIhyk3NhwOvBu4H/hY4rqq+02ddfwa8s+lqPrVn3f+v2fYVwEPAp4GtmskvBG7ucxuS9HPy85e/SNLik2QF8B3geVX1w2nmW8MieKhwkiOBiao6e8ilSFqEPHMnaVFrrof778DF0wW7Rea7wAlTPa9Pkqbj149JWrSSbEOnG/cOOo9BaYWmm7qvaw4lqZfdspIkSS0yY7dsko8kuS/J9QtRkCRJkgY345m7JC8HHgYuqKq+ugmWLl1aExMTm16dJEmSnmL16tX3V9UOk02b8Zq7qvrybL8CZ2JiglWrVs1mEUmSJPUpyR1TTZuzu2WTnJhkVZJVGzZsmKvVSpIkaRbm7G7ZqjoHOAdgn3328S4NzdrE6Z+bdPyasw5Z4Eqk4Zjsb2Bj+59umiR181Eo2iR+4EiS2qQNn2tDDXdtOICShsP3D0mjYtR6nmYMd0kuAg4AliZZB7yrqj4834VJktrNgC7Nj37ulj1mIQqRJEnSpvO7ZSVJklrEGypGgF0TkiRprnjmTpIkqUU8cydJWjSmuyvRXhCpw3An9VioD4hRu3VekkaZ4b1/dstKkiS1iGfuJGkGU50x8Oyr5oLtSHPNcCdJktSHxdI13Jpw5/98pOHz71CaW4slTGi0tCbcSRptfkgtDuMW0Mdtf0eBx3z+Ge4kSVoA/gdHvebrel7DnaSh80NPkuaOj0KRJElqEcOdJElSi4xkt6wXW44Wu8ykxWkh/3Z9nxgdfobOrO3HaCTDnTSK2v5mMBt+kGu+zXUbs81qnBjuFimDhqTFwPcqzYYhfG4Y7gT4BizNNT+k1BZ+Piw+iy7c+YYpTc434PHk733hjcox9/NQU/FuWUmSpBYx3EmSJLXIouuWHXWeJtd8s41JkqZjuGshP/wlzQffW6TFwXA3hVG5YFYzm68vXpYkaTEa+3A36v8THfX6ND/8vc8PA780t0b9vWrU65svYx/uND+m+4Ma1z82SYuL71VarAx30jzzA0KS5s+ov8cOoz7DnSRJGvmQpP4Z7iQJP9gktYfhbgB+CEijzb9RSePMcLdARuUuPT/0JGl8+RkwHuY93NmQNIpslzPzGEnaFL6HDI9n7qQxNipnlKW2MNBoFBjutCj4hjmeDJ+SNHtPG3YBkiRJmjtjcebOsz4ad/4NjBZ/H5LmU19n7pIcnOTmJLckOX2+i5IkSdJgZjxzl2Qz4IPAfwHWAdckubyqbpzv4iRJ88Ozh1J79dMtuy9wS1XdBpDkYuAIwHAnNfyglCSNilTV9DMkbwAOrqoTmuE3Ab9aVSf1zHcicGIzuBK4ee7L1SK1FLh/2EVopNgm1Ms2oV62ientUlU7TDahnzN3mWTcUxJhVZ0DnDPLwjQGkqyqqn2GXYdGh21CvWwT6mWbGFw/N1SsA5Z3De8M3D0/5UiSJGlT9BPurgF2S7Jrki2Ao4HL57csSZIkDWLGbtmqejzJScA/A5sBH6mqG+a9MrWJ3fXqZZtQL9uEetkmBjTjDRWSJElaPPz6MUmSpBYx3EmSJLWI4U5zJsnyJF9MclOSG5Kc0ozfPsm/JPle8++zhl2rFlaSzZJcl+SKZtg2McaSbJfk0iTfad4vfs02Md6S/FHzuXF9kouSbGmbGJzhTnPpceBtVfUiYD/gLUl2B04Hrqyq3YArm2GNl1OAm7qGbRPj7a+Bz1fVLwEvodM2bBNjKslOwMnAPlW1B52bN4/GNjEww53mTFWtr6prm9cP0XnD3onO19Wd38x2PnDkUArUUCTZGTgEOLdrtG1iTCV5JvBy4MMAVfVYVT2AbWLcLQG2SrIE2JrO83RtEwMy3GleJJkA9gKuBp5bVeuhEwCB5wyxNC28s4HTgCe6xtkmxtfzgQ3AeU1X/blJtsE2Mbaq6i7gL4E7gfXAg1X1BWwTAzPcac4l2Rb4FPDWqvrhsOvR8CQ5FLivqlYPuxaNjCXA3sCHqmov4BHsbhtrzbV0RwC7AsuAbZIcO9yqFjfDneZUks3pBLuPVdVlzeh7k+zYTN8RuG9Y9WnBvQw4PMka4GLglUkuxDYxztYB66rq6mb4UjphzzYxvg4Cbq+qDVX1E+Ay4NexTQzMcKc5kyR0rqO5qare1zXpcuD45vXxwD8udG0ajqo6o6p2rqoJOhdIX1VVx2KbGFtVdQ+wNsnKZtSBwI3YJsbZncB+SbZuPkcOpHPNtm1iQH5DheZMkv2BrwDf5mfXV51J57q7TwAr6PwRv7Gqvj+UIjU0SQ4ATq2qQ5M8G9vE2EqyJ50bbLYAbgPeTOdkg21iTCV5N3AUnacuXAecAGyLbWIghjtJkqQWsVtWkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJT0pyQ5IDhl3HXBilfUmyJslBc7Suf0jynoXerqTFw3An6UlV9ctV9SVYXMFgslq790WSxonhTtJIS7Jk2DVI0mJiuJNaJkkleWHX8JPdeM0ZrlOTfCvJg0kuSbJl17xrkhyU5KPACuCzSR5OctoU23pRki8leaDpBj28Z11nJLkxyQ+SnNezrWVJPpVkQ5Lbk5zcs+w7knwLeCTJkiSnJ7k1yUPNOl/bzDtprd1n8/rY772TXNes+5PN9Em7Ppu67mrmvTnJgV3Tlie5rNmn/0jyga5F95xm+9Mdi72SXNts7xJgy651Tvv7nqT2KbczybzvTfLpruG/SHJlks2nWkbSaDDcSePnt4CDgV2BXwF+t3eGqnoTcCdwWFVtW1Xv7Z2n+ZD/LPAF4DnAHwIfS7Kya7bfAV4FvAD4ReCdzbJPa5b9JrATcCDw1iSv6lr2GOAQYLuqehy4FfgN4BeAdwMXJtmxn1qn2+8kWwCfBv4B2B64CHjtZCto9u0k4KVV9Yxm39Y00zYDrgDuACaa/bq4j+1PeSya2j4DfLSp7ZPA66fYv2n1ecy7/TnwiiR7JvmDpvbXVdVPBtm+pIVjuJPGz/ur6u6q+j6dD/s9B1zPfsC2wFlV9VhVXUUn3BzTNc8Hqmpts60/7Zr2UmCHqvqTZtnbgL8Hju6pc21V/Rigqj7Z1P1EVV0CfA/Ydxb1TrXf+wFLmuk/qarLgK9PsY6fAk8Hdk+yeVWtqapbm2n7AsuAt1fVI1X1aFX9Wx/bn+5Y7AdsDpzd1HYpcM0s9rlbP8f8SVX1H8DZwAXAGcBrqupBgCSvSDIxYB2S5pnhTho/93S9/hGdgDaIZcDaqnqia9wddM4KbbS2Z9qy5vUuwLKmO/eBJA8AZwLPnWJZkhyX5Btd8+8BLJ1FvVPt9zLgrqqqqba9UVXdArwV+GPgviQXJ9m4T8uBO5qzjLPZ/nTHYrLa7phyD6fXzzHvdR3wYuCMquo+Jr8HZMA6JM0zw53UPj8Ctu4aft6A66kZpt8NLG+6+zZaAdzVNby8Z9rdzeu1wO1VtV3XzzOq6jWTbT/JLnTOMp0EPLuqtgOu52cBY6Zap7Me2ClJd1hZPtXMVfXxqtqfTlgqOt2XG/dpRWZ/A8h0x2Ky2lb0LN/v77ufY/6kJC8GPgScTyfMbRx/OHAYcF6S42azo5IWhuFOap9vAL+dZLMkBwO/OeB67gWeP830q4FHgNOSbJ7OM+UO4+evM3tLkp2TbE/nLNElzfivAz9sbk7Yqql1jyQvnWJb29AJUhsAkryZzpm7fmudztfodLee1Ny4cQRTdPcmWZnklUmeDjwK/LhZduM+rQfOSrJNki2TvKyP7U93LL4GPA6c3NT2uklq+wb9/b77PuZJdqLTdfwHwH8DXpyfPTPwCuC6qjqgqi7oY/8kLTDDndQ+p9AJWQ/QuaHhMwOu58+AdzZdeKf2Tqyqx4DDgVcD9wN/CxxXVd/pmu3jdG64uK35eU+z7E+bGvcEbm+WP5fOzRJPUVU3An9FJ+zcS6er8Kv91jqdZj9eB/w+nWN2LJ0A85+TzP504Kym3nvo3EhyZs8+vZDODR7rgKP62P6Ux6Krtt8FftCs77KeVfT1++73mCd5JvB/gfdV1eVV9SPgL+hcM0mzfzfPtF+Shic/fymHJM2NJGuAE6rqX4ddy2wluRr4u6o6b9i1jJokRwITVXX2kEuRNAXP3Ekae0l+M8nzmq7P4+k8quTzw65rRH0XOCHJ2cMuRNLkfPK7JMFK4BN07mC9FXhDVa0fbkmjqeki32PGGSUNjd2ykiRJLTJjt2ySjyS5L8n1C1GQJEmSBjfjmbskLwceBi6oqr5OxS9durQmJiY2vTpJkiQ9xerVq++vqh0mmzbjNXdV9eXZfs3MxMQEq1atms0ikiRJ6lOSKb+tZs5uqEhyInAiwIoVP3uA+sTpn3vKvGvOOmTaaZONn27aTMuMuoU6RsNe3yjU0Ib1jbJROUZtWt+m1DDqRvWYz7S+UahhuvWNulH/uxn19fVjzh6FUlXnVNU+VbXPDjtMepZQkiRJ88zn3EmSJLWI4U6SJKlF+nkUykV0vs9xZZJ1SX5//suSJEnSIPq5W/aYhShEkiRJm85uWUmSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBaZ8RsqNDcmTv/cpOPXnHXIAlcibbrJ2rNtWZPxvW98+T4xPJ65kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1yFg8CsXbsdUmtmf1y7Yynnz8jDxzJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSi4zFc+4Wks+VUlvYljUbtpfx5O99NHnmTpIkqUUMd5IkSS1iuJMkSWoRw50kSVKLGO4kSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS3SV7hLcnCSm5PckuT0+S5KkiRJg5kx3CXZDPgg8Gpgd+CYJLvPd2GSJEmavX7O3O0L3FJVt1XVY8DFwBHzW5YkSZIGkaqafobkDcDBVXVCM/wm4Fer6qSe+U4ETmwGVwI3z325WqSWAvcPuwiNFNuEetkm1Ms2Mb1dqmqHySYs6WPhTDLuKYmwqs4BzpllYRoDSVZV1T7DrkOjwzahXrYJ9bJNDK6fbtl1wPKu4Z2Bu+enHEmSJG2KfsLdNcBuSXZNsgVwNHD5/JYlSZKkQczYLVtVjyc5CfhnYDPgI1V1w7xXpjaxu169bBPqZZtQL9vEgGa8oUKSJEmLh99QIUmS1CKGO0mSpBYx3GnOJFme5ItJbkpyQ5JTmvHbJ/mXJN9r/n3WsGvVwkqyWZLrklzRDNsmxliS7ZJcmuQ7zfvFr9kmxluSP2o+N65PclGSLW0TgzPcaS49Drytql4E7Ae8pfmqutOBK6tqN+DKZljj5RTgpq5h28R4+2vg81X1S8BL6LQN28SYSrITcDKwT1XtQefmzaOxTQzMcKc5U1Xrq+ra5vVDdN6wd6LzdXXnN7OdDxw5lAI1FEl2Bg4Bzu0abZsYU0meCbwc+DBAVT1WVQ9gmxh3S4CtkiwBtqbzPF3bxIAMd5oXSSaAvYCrgedW1XroBEDgOUMsTQvvbOA04ImucbaJ8fV8YANwXtNVf26SbbBNjK2qugv4S+BOYD3wYFV9AdvEwAx3mnNJtgU+Bby1qn447Ho0PEkOBe6rqtXDrkUjYwmwN/ChqtoLeAS728Zacy3dEcCuwDJgmyTHDreqxc1wpzmVZHM6we5jVXVZM/reJDs203cE7htWfVpwLwMOT7IGuBh4ZZILsU2Ms3XAuqq6uhm+lE7Ys02Mr4OA26tqQ1X9BLgM+HVsEwMz3GnOJAmd62huqqr3dU26HDi+eX088I8LXZuGo6rOqKqdq2qCzgXSV1XVsdgmxlZV3QOsTbKyGXUgcCO2iXF2J7Bfkq2bz5ED6VyzbZsYkN9QoTmTZH/gK8C3+dn1VWfSue7uE8AKOn/Eb6yq7w+lSA1NkgOAU6vq0CTPxjYxtpLsSecGmy2A24A30znZYJsYU0neDRxF56kL1wEnANtimxiI4U6SJKlF7JaVJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U6SJKlFDHeSJEktYriTJElqEcOdJElSixjuJEmSWsRwJ0mS1CKGO0mSpBYx3EmSJLWI4U7SvEtyQ5IDhl2HJI0Dw500hpLsluTRJBcuxPaq6per6kv9zJtkTZKD5rmkkTAX+zpOx0tSfwx30nj6IHDNsIvQ4JIsGXYNkkaT4U4aM0mOBh4ArpxhvjVJ3p7kW0keSfLhJM9N8k9JHkryr0me1TX/6UlubabdmOS1Pes6qOv1qc16H0xySZItm2kfBVYAn03ycJLTmvHLknwqyYYktyc5uWfdk66vz7qmW3Z5ksua7f5Hkg/MVM8kx/EdSe5qtn9zkgOn2tfpau2q9x1JvgU8kuSiSdYx6fZ61rNtkp8m2bFr3B5J1id5xpSNQtLiUFX++OPPmPwAzwS+CywH/hi4cJp51wD/DjwX2Am4D7gW2At4OnAV8K6u+d8ILKPzn8ajgEeAHbvWdVDX6683824P3AT8Qc92D+oafhqwGvifwBbA84HbgFf1ub6Z6pp0WWAz4JvA/wa2AbYE9p+pnp5juBJYCyxrhieAF0yzr1PW2jX/N5rf31aTHNtpt9dT2w3AIV3DVwB/OOw26o8//mz6j2fupPHyv4APV9XaPuf/m6q6t6ruAr4CXF1V11XVfwKfphP0AKiqT1bV3VX1RFVdAnwP2HeK9b6/mff7wGeBPaep4aXADlX1J1X1WFXdBvw9cHQ/6+ujrqmW3ZdO0Hp7VT1SVY9W1b/1Wc9GP6UThHdPsnlVramqW6fa0T6P4furam1V/XgTt3cNsDdAkpcDuwP/pxl+RZKJqeqUNNoMd9KYSLIncBCdM1H9urfr9Y8nGd62a/3HJflGkgeSPADsASydYr33dL3+Ufd6JrELsGzjept1n0nnjOKM6+ujrqmWXQ7cUVWPD1APAFV1C/BWOmdJ70tycZJlU+1on8dwymA+y+09Ge6A9wL/o6oea4Z/D8hU25E02gx30vg4gE433Z1J7gFOBV6f5NpNXXGSXeicvToJeHZVbQdcz2ABoXqG1wK3V9V2XT/PqKrXzHNda4EVk9y4MKt6qurjVbU/nVBYwJ9Ptq+zqLX3+Pzc8Azb63YNsHeS1wNbARc1dRwOHAacl+S4KZaVNMIMd9L4OAd4AZ1uxz2BvwM+B7xqDta9DZ0gsQEgyZvpnHUaxL10rmPb6OvAD5sbBbZKsllz8f9L57murwPrgbOSbJNkyyQvm009SVYmeWWSpwOP0jnb+dMp9nXQWp9cRx/b6/ZN4HnAXwGnV9UTzfgrgOuq6oCquqCP7UsaMYY7aUxU1Y+q6p6NP8DDwKNVtWEO1n0jnZDwNTph48XAVwdc3Z8B72y6Jk+tqp/SOZO0J3A7cD9wLvAL81lX13ZfCNwJrAOOmmU9TwfOaua5B3gOnS7cp+wr8JoBa+1ex1EzbK97//4T+Dawpqr+qWvSC4Gb+9iupBGVqt4z/JKktkuyBXAL8FtV9e9d448EJqrq7CGVJmkTeeZOksbTu4Cvdge7xneBE5KcvfAlSZoLnrmTpDGSZG/gi8C3gNdW1f1DLknSHDPcSZIktciM3bJJPpLkviTXL0RBkiRJGtyMZ+6aJ5c/DFxQVX09QmDp0qU1MTGx6dVJkiTpKVavXn1/Ve0w2bTeh3M+RVV9ebZfQzMxMcGqVatms4gkSZL6lOSOqabNGO5msZETgRMBVqxYMVernXcTp3/uKePWnHXIECqRNs1kbRlsz5qc733jyd/7eJizR6FU1TlVtU9V7bPDDpOeJZQkSdI88zl3kiRJLWK4kyRJapF+HoVyEZ3vOlyZZF2S35//siRJkjSIfu6WPWYhCpEkSdKms1tWkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUosY7iRJklrEcCdJktQihjtJkqQWMdxJkiS1iOFOkiSpRQx3kiRJLWK4kyRJahHDnSRJUov0Fe6SHJzk5iS3JDl9vouSJEnSYGYMd0k2Az4IvBrYHTgmye7zXZgkSZJmr58zd/sCt1TVbVX1GHAxcMT8liVJkqRBpKqmnyF5A3BwVZ3QDL8J+NWqOqlnvhOBE5vBlcDNc1+uFqmlwP3DLkIjxTahXrYJ9bJNTG+XqtphsglL+lg4k4x7SiKsqnOAc2ZZmMZAklVVtc+w69DosE2ol21CvWwTg+unW3YdsLxreGfg7vkpR5IkSZuin3B3DbBbkl2TbAEcDVw+v2VJkiRpEDN2y1bV40lOAv4Z2Az4SFXdMO+VqU3srlcv24R62SbUyzYxoBlvqJAkSdLi4TdUSJIktYjhTpIkqUUMd5ozSZYn+WKSm5LckOSUZvz2Sf4lyfeaf5817Fq1sJJsluS6JFc0w7aJMZZkuySXJvlO837xa7aJ8Zbkj5rPjeuTXJRkS9vE4Ax3mkuPA2+rqhcB+wFvab6q7nTgyqraDbiyGdZ4OQW4qWvYNjHe/hr4fFX9EvASOm3DNjGmkuwEnAzsU1V70Ll582hsEwMz3GnOVNX6qrq2ef0QnTfsneh8Xd35zWznA0cOpUANRZKdgUOAc7tG2ybGVJJnAi8HPgxQVY9V1QPYJsbdEmCrJEuArek8T9c2MSDDneZFkglgL+Bq4LlVtR46ARB4zhBL08I7GzgNeKJrnG1ifD0f2ACc13TVn5tkG2wTY6uq7gL+ErgTWA88WFVfwDYxMMOd5lySbYFPAW+tqh8Oux4NT5JDgfuqavWwa9HIWALsDXyoqvYCHsHutrHWXEt3BLArsAzYJsmxw61qcTPcaU4l2ZxOsPtYVV3WjL43yY7N9B2B+4ZVnxbcy4DDk6wBLgZemeRCbBPjbB2wrqquboYvpRP2bBPj6yDg9qraUFU/AS4Dfh3bxMAMd5ozSULnOpqbqup9XZMuB45vXh8P/ONC16bhqKozqmrnqpqgc4H0VVV1LLaJsVVV9wBrk6xsRh0I3IhtYpzdCeyXZOvmc+RAOtds2yYG5DdUaM4k2R/4CvBtfnZ91Zl0rrv7BLCCzh/xG6vq+0MpUkOT5ADg1Ko6NMmzsU2MrSR70rnBZgvgNuDNdE422CbGVJJ3A0fReerCdcAJwLbYJgZiuJMkSWoRu2UlSZJaxHAnSZLUIoY7SZKkFjHcSZIktYjhTpIkqUUMd5IkSS1iuJMkSWqR/w8QKp5POY4ebgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def maintenance_planning_ramp_operational(c, T, M, P, N):\n",
    "    m = pyo.ConcreteModel()\n",
    "\n",
    "    m.T = pyo.RangeSet(1, T)\n",
    "    m.Y = pyo.RangeSet(1, T - M + 1)\n",
    "    m.S = pyo.RangeSet(0, M - 1)\n",
    "    m.W = pyo.RangeSet(0, M + N - 1)\n",
    "\n",
    "    m.c = pyo.Param(m.T, initialize = c)\n",
    "    m.x = pyo.Var(m.T, bounds=(0, 1))\n",
    "    m.y = pyo.Var(m.T, domain=pyo.Binary)\n",
    "    m.upos = pyo.Var(m.T, bounds=(0, upos_max))\n",
    "    m.uneg = pyo.Var(m.T, bounds=(0, uneg_max))\n",
    "\n",
    "    # objective\n",
    "    m.profit = pyo.Objective(expr = sum(m.c[t]*m.x[t] for t in m.T), sense=pyo.maximize)\n",
    "    \n",
    "    # ramp constraint\n",
    "    m.ramp = pyo.Constraint(m.T, rule = lambda m, t: \n",
    "         m.x[t] == m.x[t-1] + m.upos[t] - m.uneg[t] if t > 1 else pyo.Constraint.Skip)\n",
    "      \n",
    "    # required number P of maintenance starts\n",
    "    m.sumy = pyo.Constraint(expr = sum(m.y[t] for t in m.Y) == P)\n",
    "\n",
    "    # no more than one maintenance start in the period of length M\n",
    "    m.sprd = pyo.Constraint(m.Y, rule = lambda m, t: sum(m.y[t+s] for s in m.W if t + s <= T) <= 1)   \n",
    "    \n",
    "    # Choose one or the other the following methods. Comment out the method not used.\n",
    "    \n",
    "    # disjunctive constraints, big-M method.\n",
    "    m.bigm = pyo.Constraint(m.Y, rule = lambda m, t: sum(m.x[t+s] for s in m.S) <= (M+N)*(1 - m.y[t]))\n",
    "    \n",
    "    # disjunctive constraints, GDP programming method\n",
    "    #m.disj = gdp.Disjunction(m.Y, rule = lambda m, t: [m.y[t]==0, sum(m.x[t+s] for s in m.W if t + s <= T)==0])\n",
    "    #pyo.TransformationFactory('gdp.hull').apply_to(m)\n",
    "    \n",
    "    return m\n",
    "\n",
    "m = maintenance_planning_ramp_operational(c, T, M, P, N)\n",
    "pyo.SolverFactory('cbc').solve(m)\n",
    "plot_schedule(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3j0fLxBOwSch"
   },
   "source": [
    "## Exercises\n",
    "\n",
    "1. Rather than specify how many maintenance periods must be accomodated, modify the model so that the process unit can operate no more than $N$ days without a maintenance shutdown.  (Hint. You may to introduce an additional set of binary variables, $z_t$ to denote the start of an operational period.)\n",
    "\n",
    "2. Do a systematic comparison of the Big-M, Convex Hull, and Cutting Plane techniques for implementing the disjunctive constraints. Your comparison should include a measure of complexity (such as the number of decision variables and constraints in the resulting transformed problems), computational effort, and the effect of solver (such as glpk vs cbc)."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "04.03-Maintenance-Planning.ipynb",
   "provenance": [
    {
     "file_id": "https://github.com/jckantor/ND-Pyomo-Cookbook/blob/master/notebooks/04.03-Maintenance-Planning.ipynb",
     "timestamp": 1556576712896
    }
   ],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}