electrical_power_generation/electrical_power_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Electrical Power Generation 2\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "This example (which is an extension of the earlier ‘Electrical Power Generation 1 ‘ example) will teach you how to choose an optimal set of power stations to turn on in order to satisfy anticipated power demand over a 24-hour time horizon – but gives you the option of using hydroelectric power plants to satisfy that demand.\n",
    "\n",
    "This model is example 16 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 271-272 and 326-327.\n",
    "\n",
    "This example is at the intermediate level, where we assume that you know Python and the Gurobi Python API and that you have some knowledge of building mathematical optimization models.\n",
    "\n",
    "**Download the Repository** <br />\n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). \n",
    "\n",
    "---\n",
    "## Problem Description\n",
    "\n",
    "In this problem, thermal power generation units are grouped into three distinct types, with different characteristics for each type (power output, cost per megawatt hour, startup cost, etc.).  A unit can be on or off, with a startup cost associated with transitioning from off to on, and power output that can fall anywhere between a specified minimum and maximum value when the unit is on.  There are also two hydroelectric plants available, also with different cost and power generation characteristics.  A 24-hour time horizon is divided into 5 discrete time periods, each with an expected total power demand.  The model decides which units to turn on, and when, in order to satisfy demand for each time period.  The model also captures a reserve requirement, where the selected power plants must be capable of increasing their output, while still respecting their maximum output, in order to cope with the situation where actual demand exceeds predicted demand.\n",
    "\n",
    "A set of generators is available to satisfy power demand for the following day.  Anticipated demand is as follows:\n",
    "\n",
    "| Time Period | Demand (megawatts) |\n",
    "| --- | --- |\n",
    "| 12 pm to 6 am | 15000 |\n",
    "| 6 am to 9 am | 30000 |\n",
    "| 9 am to 3 pm | 25000 |\n",
    "| 3 pm to 6 pm | 40000 |\n",
    "| 6 pm to 12 pm | 27000 |\n",
    "\n",
    "Thermal generators are grouped into three types, with the following minimum and maximum output for each type (when they are on):\n",
    "\n",
    "| Type | Number available | Minimum output (MW) | Maximum output (MW) |\n",
    "| --- | --- | --- | --- |\n",
    "| 0 | 12 |  850 | 2000 |\n",
    "| 1 | 10 | 1250 | 1750 |\n",
    "| 2 | 5 | 1500 | 4000 |\n",
    "\n",
    "There are costs associated with using a thermal generator: a cost per hour when the generator is on (and generating its minimum output), a cost per megawatt hour above its minimum, and a startup cost for turning a generator on:\n",
    "\n",
    "| Type | Cost per hour (when on) | Cost per MWh above minimum | Startup cost |\n",
    "| --- | --- | --- | --- |\n",
    "| 0 | $\\$1000$ | $\\$2.00$ | $\\$2000$ |\n",
    "| 1 | $\\$2600$ | $\\$1.30$ | $\\$1000$ |\n",
    "| 2 | $\\$3000$ | $\\$3.00$ | $\\$500$ |\n",
    "\n",
    "Two hydroelectric generators are also available, each with a fixed power output (when on):\n",
    "\n",
    "| Hydro plant | Output (MW) |\n",
    "| --- | --- |\n",
    "| A | 900 |\n",
    "| B | 1400 |\n",
    "\n",
    "The costs associated with using a hydro plant are slightly different.  There's an hourly cost, but it is much smaller than the hourly cost of a thermal generator.  The real cost for a hydroelectric plant comes from depletion of the water in the reservoir, which happens at different rates for the two units.  The reservoir must be replenished before the end of the time horizon by pumping water into it, which consumes electricity.  A hydroelectric plant also has a startup cost.\n",
    "\n",
    "| Hydro plant | Cost per hour (when on) | Startup cost | Reservoir depth reduction (m/hr) |\n",
    "| --- | --- | --- | --- |\n",
    "| A | $\\$90$ | $\\$1500$ | 0.31 |\n",
    "| B | $\\$150$ | $\\$1200$ | 0.47 |\n",
    "\n",
    "Pumping water into the reservoir consumes electricity at a rate of 3000 MWh of electricity per meter of height.  The height of the reservoir at the end of the time horizon must be equal to the height at the beginning.\n",
    "\n",
    "Generators must meet predicted demand, but they must also have sufficient reserve capacity to be able to cope with the situation where actual demand exceeds predicted demand.  For this model, the set of selected thermal generators plus the set of hydro generators must be able to produce as much as 115% of predicted demand.\n",
    "\n",
    "Which generators should be committed to meet anticipated demand in order to minimize total cost?\n",
    "\n",
    "---\n",
    "## Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "\n",
    "$t \\in \\text{Types}=\\{0,1,2\\}$: Types of thermal generators.\n",
    "\n",
    "$h \\in \\text{HydroUnits}=\\{0,1\\}$: Two hydro generators.\n",
    "\n",
    "$p \\in \\text{Periods}=\\{0,1,2,3,4\\}$: Time periods.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$\\text{period_hours}_{p} \\in \\mathbb{N}^+$: Number of hours in each time period.\n",
    "\n",
    "$\\text{demand}_p \\in \\mathbb{R}^+$: Total power demand for time period $p$.\n",
    "\n",
    "$\\text{generators}_t \\in \\mathbb{N}^+$: Number of thermal generators of type $t$.\n",
    "\n",
    "$\\text{start0} \\in \\mathbb{N}^+$: Number of thermal generators that are on at the beginning of the time horizon (and available in time period 0 without paying a startup cost).\n",
    "\n",
    "$\\text{min_output}_t \\in \\mathbb{R}^+$: Minimum output for thermal generator type $t$ (when on).\n",
    "\n",
    "$\\text{max_output}_t \\in \\mathbb{R}^+$: Maximum output for thermal generator type $t$.\n",
    "\n",
    "$\\text{base_cost}_t \\in \\mathbb{R}^+$: Minimum operating cost (per hour) for a thermal generator of type $t$.\n",
    "\n",
    "$\\text{per_mw_cost}_t \\in \\mathbb{R}^+$: Cost to generate one additional MW (per hour) for a thermal generator of type $t$.\n",
    "\n",
    "$\\text{startup_cost}_t \\in \\mathbb{R}^+$: Startup cost for thermal generator of type $t$.\n",
    "\n",
    "$\\text{hydro_load}_h \\in \\mathbb{R}^+$: Output for hydro generator $h$.\n",
    "\n",
    "$\\text{hydro_cost}_h \\in \\mathbb{R}^+$: Cost for operating hydro generator $h$.\n",
    "\n",
    "$\\text{hydro_startup_cost}_h \\in \\mathbb{R}^+$: Startup cost for hydro generator $h$.\n",
    "\n",
    "$\\text{hydro_height_reduction}_h \\in \\mathbb{R}^+$: Hourly reduction in reservoir height from operating hydro generator $h$.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$\\text{ngen}_{t,p} \\in \\mathbb{N}^+$: Number of thermal generators of type $t$ that are on in time period $p$.\n",
    "\n",
    "$\\text{output}_{t,p} \\in \\mathbb{R}^+$: Total power output from thermal generators of type $t$ in time period $p$.\n",
    "\n",
    "$\\text{nstart}_{t,p} \\in \\mathbb{N}^+$: Number of thermal generators of type $t$ to start in time period $p$.\n",
    "\n",
    "$\\text{hydro}_{h,p} \\in [0,1]$: Indicates whether hydro generators $h$ is on in time period $p$.\n",
    "\n",
    "$\\text{hydro_start}_{h,p} \\in [0,1]$: Indicates whether hydro generator $h$ starts in time period $p$.\n",
    "\n",
    "$\\text{height}_{p} \\in \\mathbb{R}^+$: Height of reservoir in time period $p$.\n",
    "\n",
    "$\\text{pumping}_{p} \\in \\mathbb{R}^+$: Power used to replenish reservoir in time period $p$.\n",
    "\n",
    "\n",
    "### Objective Function\n",
    "\n",
    "- **Cost**: Minimize the cost (in USD) to satisfy the predicted electricity demand.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Minimize} \\quad Z_{on} + Z_{extra} + Z_{startup} + Z_{hydro} + Z_{hydro\\_startup}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "Z_{on} = \\sum_{(t,p) \\in \\text{Types} \\times \\text{Periods}}{\\text{base_cost}_t*\\text{ngen}_{t,p}}\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "Z_{extra} = \\sum_{(t,p) \\in \\text{Types} \\times \\text{Periods}}{\\text{per_mw_cost}_t*(\\text{output}_{t,p} - \\text{min_load}_t})\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "Z_{startup} = \\sum_{(t,p) \\in \\text{Types} \\times \\text{Periods}}{\\text{startup_cost}_t*\\text{nstart}_{t,p}}\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "Z_{hydro} = \\sum_{(h,p) \\in \\text{HydroUnits} \\times \\text{Periods}}{\\text{hydro_cost}_h*\\text{hydro}_{h,p}}\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "Z_{hydro\\_startup} = \\sum_{(h,p) \\in \\text{HydroUnits} \\times \\text{Periods}}{\\text{hydro_startup_cost}_h*\\text{hydro_start}_{h,p}}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "### Constraints\n",
    "\n",
    "- **Available generators**: Number of generators used must be less than or equal to the number available.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{ngen}_{t,p} \\leq \\text{generators}_t  \\quad \\forall (t,p) \\in \\text{Types} \\times \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- **Demand**: Total power generated across all generator types must meet anticipated demand plus pumping for each time period $p$.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{t \\in \\text{Types}}{\\text{output}_{t,p}} +\n",
    "\\sum_{h \\in \\text{HydroUnits}}{\\text{hydro_load}_h*\\text{hydro}_{h,p}} \\geq\n",
    "\\text{demand}_p + \\text{pumping}_p \\quad \\forall p \\in \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- **Min/max generation**: Power generation must respect thermal generator min/max values.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{output}_{t,p} \\geq \\text{min_output}_t*\\text{ngen}_{t,p} \\quad \\forall (t,p) \\in \\text{Types} \\times \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{output}_{t,p} \\leq \\text{max_output}_t*\\text{ngen}_{t,p} \\quad \\forall (t,p) \\in \\text{Types} \\times \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- **Reserve**: Selected thermal generators plus hydro units must be able to satisfy demand that is as much as 15% above the prediction.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{t \\in \\text{Types}}{\\text{max_output}_t*\\text{ngen}_{t,p}} +\n",
    "\\sum_{h \\in \\text{HydroUnits}}{\\text{hydro_load}_h} \\geq 1.15 * \\text{demand}_p \\quad \\forall p \\in \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- **Thermal startup**: Establish relationship between number of active thermal generators and number of startups (use $start0$ for period before the time horizon starts)\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{ngen}_{t,p} \\leq \\text{ngen}_{t,p-1} + \\text{nstart}_{t,p} \\quad \\forall (t,p) \\in \\text{Types} \\times \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- **Hydro startup**: Establish relationship between hydro generator state and number of hydro startups (assume hydro plants are off at the beginning of the horizon)\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{hydro}_{h,p} \\leq \\text{hydro}_{h,p-1} + \\text{hydro_start}_{h,p} \\quad \\forall (h,p) \\in \\text{HydroUnits} \\times \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- **Reservoir height**: Track reservoir height.   Note that the height at the end of the final time period must equal the height at the beginning of the first.\n",
    "\n",
    "- Reservoir level constraints: Height increases due to pumping activity and decreases due to hydroelectric generation.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{height}_{p} = \\text{height}_{p-1}  + \\text{period_hours}_{p}*\\text{pumping}_{p}/3000 -\n",
    "\\sum_{h \\in \\text{HydroUnits}}{\\text{period_hours}_{p}*\\text{hydro_height_reduction}_{h}*\\text{hydro}_{h,p}} \\quad \\forall p \\in \\text{Periods}\n",
    "\\end{equation}\n",
    "\n",
    "- Cyclic constraint: Height at the first period must be equal to height at the last period.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{height}_{pfirst} = \\text{height}_{plast}  + \\text{period_hours}_{pfirst}*\\text{pumping}_{pfirst}/3000 -\n",
    "\\sum_{h \\in \\text{HydroUnits}}{\\text{period_hours}_{pfirst}*\\text{hydro_height_reduction}_{h}*\\text{hydro}_{h,pfirst}}\n",
    "\\end{equation}\n",
    "\n",
    "---\n",
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.11 & Gurobi 11.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input Data\n",
    "We define all the input data of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "\n",
    "ntypes = 3\n",
    "nperiods = 5\n",
    "maxstart0 = 5\n",
    "hydrounits = 2\n",
    "\n",
    "generators = [12, 10, 5]\n",
    "period_hours = [6, 3, 6, 3, 6]\n",
    "demand = [15000, 30000, 25000, 40000, 27000]\n",
    "min_load = [850, 1250, 1500]\n",
    "max_load = [2000, 1750, 4000]\n",
    "base_cost = [1000, 2600, 3000]\n",
    "per_mw_cost = [2, 1.3, 3]\n",
    "startup_cost = [2000, 1000, 500]\n",
    "\n",
    "hydro_load = [900, 1400]\n",
    "hydro_cost = [90, 150]\n",
    "hydro_height_reduction = [0.31, 0.47]\n",
    "hydro_startup_cost = [1500, 1200]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "We create a model and the variables. For each time period, we have: an integer decision variable to capture the number of active generators of each type (ngen), an integer variable to capture the number of generators of each type we must start (nstart), a continuous decision variable to capture the total power output for each generator type (output), a binary decision variable that indicates whether a hydro unit is active (hydro), a binary decision variable that indicates whether a hydro unit must be started (hydrstart), a continuous decision variable that captures the amount of enery used to replenish the reservoir (pumping), and a continuous decision variable that captures the height of the reservoir (height)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n"
     ]
    }
   ],
   "source": [
    "model = gp.Model('PowerGeneration2')\n",
    "\n",
    "ngen = model.addVars(ntypes, nperiods, vtype=GRB.INTEGER, name=\"ngen\")\n",
    "nstart = model.addVars(ntypes, nperiods, vtype=GRB.INTEGER, name=\"nstart\")\n",
    "output = model.addVars(ntypes, nperiods, vtype=GRB.CONTINUOUS, name=\"genoutput\")\n",
    "\n",
    "hydro = model.addVars(hydrounits, nperiods, vtype=GRB.BINARY, name=\"hydro\")\n",
    "hydrostart = model.addVars(hydrounits, nperiods, vtype=GRB.BINARY, name=\"hydrostart\")\n",
    "pumping = model.addVars(nperiods, vtype=GRB.CONTINUOUS, name=\"pumping\")\n",
    "height = model.addVars(nperiods, vtype=GRB.CONTINUOUS, name=\"height\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we insert the constraints:\n",
    "\n",
    "The number of active generators can't exceed the number of generators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generator count\n",
    "\n",
    "numgen = model.addConstrs(ngen[type, period] <= generators[type]\n",
    "                         for type in range(ntypes) for period in range(nperiods))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total power output for a thermal generator type depends on the number of generators of that type that are active."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Respect minimum and maximum output per generator type\n",
    "\n",
    "min_output = model.addConstrs((output[type, period] >= min_load[type] * ngen[type, period])\n",
    "                              for type in range(ntypes) for period in range(nperiods))\n",
    "\n",
    "max_output = model.addConstrs((output[type, period] <= max_load[type] * ngen[type, period])\n",
    "                              for type in range(ntypes) for period in range(nperiods))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total generator output (thermal plus hydro) for each time period must meet predicted demand plus pumping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Meet demand\n",
    "\n",
    "meet_demand = model.addConstrs(gp.quicksum(output[type, period] for type in range(ntypes)) +\n",
    "                               gp.quicksum(hydro_load[unit]*hydro[unit,period] for unit in range(hydrounits))\n",
    "                               >= demand[period] + pumping[period]\n",
    "                               for period in range(nperiods))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maintain appropriate reservoir levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reservoir levels\n",
    "\n",
    "reservoir = model.addConstrs(height[period] == height[period-1] + period_hours[period]*pumping[period]/3000 -\n",
    "                             gp.quicksum(hydro_height_reduction[unit]*period_hours[period]*hydro[unit,period] for unit in range(hydrounits))\n",
    "                             for period in range(1,nperiods))\n",
    "\n",
    "# cyclic - height at end must equal height at beginning\n",
    "reservoir0 = model.addConstr(height[0] == height[nperiods-1] + period_hours[0]*pumping[0]/3000 -\n",
    "                             gp.quicksum(hydro_height_reduction[unit]*period_hours[0]*hydro[unit,0]\n",
    "                             for unit in range(hydrounits)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Selected generators must be able to cope with an excess of demand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Provide sufficient reserve capacity\n",
    "\n",
    "reserve = model.addConstrs(gp.quicksum(max_load[type]*ngen[type, period] for type in range(ntypes)) >=\n",
    "                           1.15*demand[period] - sum(hydro_load)\n",
    "                           for period in range(nperiods))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Connect the decision variables that capture active thermal generators with the decision variables that count startups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Startup constraint\n",
    "\n",
    "startup0 = model.addConstrs((ngen[type,0] <= maxstart0 + nstart[type,0])\n",
    "                            for type in range(ntypes))\n",
    "\n",
    "startup = model.addConstrs((ngen[type,period] <= ngen[type,period-1] + nstart[type,period])\n",
    "                           for type in range(ntypes) for period in range(1,nperiods))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Connect hydro decision variables with hydro startup decision variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hydro startup constraint\n",
    "\n",
    "hydro_startup0 = model.addConstrs((hydro[unit,0] <= hydrostart[unit,0])\n",
    "                                    for unit in range(hydrounits))\n",
    "\n",
    "hydro_startup = model.addConstrs((hydro[unit,period] <= hydro[unit,period-1] + hydrostart[unit,period])\n",
    "                           for unit in range(hydrounits) for period in range(1,nperiods))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Objective: minimize total cost.  Cost consists of five components: the cost for running active thermal generation units, the cost to generate power beyond the minimum for each unit, the cost to start up thermal generation units, the cost to operate hydro units, and the cost to start up hydro units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective: minimize total cost\n",
    "\n",
    "active = gp.quicksum(base_cost[type]*period_hours[period]*ngen[type,period]\n",
    "                    for type in range(ntypes) for period in range(nperiods))\n",
    "\n",
    "per_mw = gp.quicksum(per_mw_cost[type]*period_hours[period]*(output[type,period] - min_load[type]*ngen[type,period])\n",
    "                       for type in range(ntypes) for period in range(nperiods))\n",
    "\n",
    "startup_obj = gp.quicksum(startup_cost[type]*nstart[type,period]\n",
    "                         for type in range(ntypes) for period in range(nperiods))\n",
    "\n",
    "hydro_active = gp.quicksum(hydro_cost[unit]*period_hours[period]*hydro[unit,period]\n",
    "                           for unit in range(hydrounits) for period in range(nperiods))\n",
    "\n",
    "hydro_startup = gp.quicksum(hydro_startup_cost[unit]*hydrostart[unit,period]\n",
    "                            for unit in range(hydrounits) for period in range(nperiods))\n",
    "\n",
    "model.setObjective(active + per_mw + startup_obj + hydro_active + hydro_startup)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we start the optimization and Gurobi finds the optimal solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 85 rows, 75 columns and 215 nonzeros\n",
      "Model fingerprint: 0x1b27f759\n",
      "Variable types: 25 continuous, 50 integer (20 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e-03, 4e+03]\n",
      "  Objective range  [4e+00, 9e+03]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [5e+00, 4e+04]\n",
      "Presolve removed 18 rows and 3 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 67 rows, 72 columns, 194 nonzeros\n",
      "Variable types: 25 continuous, 47 integer (18 binary)\n",
      "\n",
      "Root relaxation: objective 9.990143e+05, 28 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 999014.286    0    7          - 999014.286      -     -    0s\n",
      "H    0     0                    1106190.0000 999014.286  9.69%     -    0s\n",
      "H    0     0                    1030900.0000 999014.286  3.09%     -    0s\n",
      "H    0     0                    1011330.0000 999014.286  1.22%     -    0s\n",
      "H    0     0                    1008935.0000 999014.286  0.98%     -    0s\n",
      "H    0     0                    1006220.0000 999014.286  0.72%     -    0s\n",
      "H    0     0                    1005470.0000 999014.286  0.64%     -    0s\n",
      "H    0     0                    1003274.0000 999014.286  0.42%     -    0s\n",
      "H    0     0                    1000630.0000 999955.477  0.07%     -    0s\n",
      "     0     0 999955.477    0   10 1000630.00 999955.477  0.07%     -    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 1\n",
      "  MIR: 4\n",
      "\n",
      "Explored 1 nodes (39 simplex iterations) in 0.04 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 7: 1.00063e+06 1.00327e+06 1.00547e+06 ... 1.10619e+06\n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 1.000630000000e+06, best bound 1.000611428571e+06, gap 0.0019%\n"
     ]
    }
   ],
   "source": [
    "model.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Analysis\n",
    "\n",
    "The anticipated demand for electricity over the 24-hour time window can be met for a total cost of $\\$1,000,630$, as compared with the $\\$1,002,540$ that was required to meet the same demand without the hydro plants.  The detailed plan for each time period is as follows.\n",
    "\n",
    "### Unit Commitments\n",
    "\n",
    "The following table shows the number of thermal generators of each type that are active in each time period in the optimal solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Thermal0</th>\n",
       "      <td>12.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>12.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Thermal1</th>\n",
       "      <td>3.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Thermal2</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             0     1     2     3     4\n",
       "Thermal0  12.0  12.0  12.0  12.0  12.0\n",
       "Thermal1   3.0   9.0   9.0   9.0   9.0\n",
       "Thermal2   0.0   0.0   0.0   1.0   0.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = [\"Thermal\" + str(t) for t in range(ntypes)]\n",
    "units = pd.DataFrame(columns=range(nperiods), index=rows, data=0.0)\n",
    "\n",
    "for t in range(ntypes):\n",
    "    for p in range(nperiods):\n",
    "        units.loc[\"Thermal\"+str(t), p] = ngen[t,p].x\n",
    "units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following shows the number of generators of each type that must be started in each time period to achieve this plan (recall that the model assumes that 5 are available at the beginning of the time horizon):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>HydroA</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>HydroB</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0    1    2    3    4\n",
       "HydroA  0.0  0.0  0.0  0.0  0.0\n",
       "HydroB  0.0  0.0  0.0  1.0  1.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = [\"HydroA\", \"HydroB\"]\n",
    "hydrotable = pd.DataFrame(columns=range(nperiods), index=rows, data=0.0)\n",
    "\n",
    "for p in range(nperiods):\n",
    "    hydrotable.loc[\"HydroA\", p] = int(hydro[0,p].x)\n",
    "    hydrotable.loc[\"HydroB\", p] = int(hydro[1,p].x)\n",
    "hydrotable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The hydro plants are lightly used in the schedule - we only use plant B for the last two time periods.\n",
    "\n",
    "The following shows the pumping that must be done in order to support this level of hydro activity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Pumping</th>\n",
       "      <td>815.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>950.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>350.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             0    1      2    3      4\n",
       "Pumping  815.0  0.0  950.0  0.0  350.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pumptable = pd.DataFrame(columns=range(nperiods), index=[\"Pumping\"], data=0.0)\n",
    "\n",
    "for p in range(nperiods):\n",
    "        pumptable.loc[\"Pumping\", p] = pumping[p].x\n",
    "pumptable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is interesting to note that the plan simultaneously operates HydroB and pumps in time period 4.\n",
    "While it might appear that costs could be reduced by operating the hydro plant at a lower load, in this model the hydro plants have fixed output.  If it makes economic sense to turn a hydro plant on, then we have to replenish the water drained, even if it means using some of the generated electricity to pump."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}