economic_planning/economic_planning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Economic Planning\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "In this example, you’ll discover how mathematical optimization can be used to address a  macroeconomic planning problem that a country may face. We’ll show you how to model and solve an input-output problem encompassing the interrelationships between the different sectors of a country’s economy.\n",
    "\n",
    "This model is example 9 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 263-264 and 316-317.\n",
    "\n",
    "This modeling example is at the intermediate level, where we assume that you know Python and are familiar with the Gurobi Python API. In addition, you should have some knowledge about 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). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "In this problem, we assume that we have an economy with three types of industries:\n",
    "\n",
    "* Coal\n",
    "* Steel\n",
    "* Transport\n",
    "\n",
    "Part of the outputs from these industries are needed as inputs to others. For example, coal is needed to fire the blast furnaces that produce steel, steel is needed in the machinery for extracting coal, etc. We measure all units of production in monetary terms.  That is, each dollar of production by one of the industries requires inputs (in dollars) from possibly its own industry as well as other industries. The required inputs as well as the labor requirements (also measured in dollars) are shown in the following table. \n",
    "\n",
    "| input (t) / <br /> output (t+1) | Coal | Steel | Transport |\n",
    "| --- | --- | --- | --- | \n",
    "| Coal | 0.1 | 0.5 | 0.4 | \n",
    "| Steel | 0.1 | 0.1 | 0.2 |\n",
    "| Transport | 0.2 | 0.1 | 0.2 |\n",
    "| Labor | 0.6 | 0.3 | 0.2 | \n",
    "|**Total**  | **1** | **1** | **1**|\n",
    "\n",
    "There is a time lag in the economy so that the output in year t +1 requires an input in year t.  For example, to produce one-dollar worth of coal requires 0.1 dollars of coal (to provide the necessary power), 0.1 dollars of steel (the steel ‘used up’ in the ‘wear and tear’ on the machinery) and 0.2 dollars of transport (for moving the coal from the mine). In addition, 0.6 dollars of manpower is required. Similarly, the other columns of the table give the inputs required (in dollars) for each dollar of steel and each dollar of transport (trucks, cars, trains, etc.). Notice that the value of each unit of output is exactly matched by the sum of the values of its inputs.\n",
    "\n",
    "Output from an industry may also be used to build productive capacity for itself or other industries in future years. The inputs required to give unit increases (capacity in dollars’ worth of extra production) in productive capacity are given in the following table. Input from an industry in year t results in a (permanent) increase in productive capacity in year t +2.\n",
    "\n",
    "| input (t) / <br /> output (t+2) | Coal | Steel | Transport |\n",
    "| --- | --- | --- | --- | \n",
    "| Coal | 0.1 | 0.7 | 0.9 | \n",
    "| Steel | 0.1 | 0.1 | 0.2 |\n",
    "| Transport | 0.2 | 0.1 | 0.2 |\n",
    "| Labor | 0.4 | 0.2 | 0.1 | \n",
    "\n",
    "Stocks of goods may be held from year to year. At present (the beginning of year 1), the stocks and productive capacities (per year) are given in the following table (in millions of dollars).\n",
    "\n",
    "| Present | Stocks | Productive Capacity |\n",
    "| --- | --- | --- | \n",
    "| Coal | 150 | 300 | \n",
    "| Steel | 80 | 350 | \n",
    "| Transport | 100 | 280 |  \n",
    "\n",
    "We want to maximize the total manpower utilization i.e. employment, over the planning horizon while meeting an exogenous consumption requirement in every year:\n",
    "\n",
    "* 60 million dollars of coal. \n",
    "* 60 million dollars of steel.\n",
    "* 30 million dollars  of transport.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Formulation\n",
    "A widely used type of national economic model is the input–output model representing the interrelationships between the different sectors of a country’s economy. Such models are often referred to as dynamic Leontief models after their originator, who built such a model of the American economy. A similar model has been considered by Wagner (1957).\n",
    "\n",
    "### Sets and indices\n",
    "\n",
    "$i,j \\in \\text{Industries}=\\{\\text{coal}, \\text{steel}, \\text{transport}\\}$\n",
    "\n",
    "$ t \\in \\text{Horizon} = \\{\\text{year1}, \\text{year2}, \\text{year3}, \\text{year4}, \\text{year5}, \\text{year6} \\}$ \n",
    "\n",
    "$ t \\in \\text{fiveYears} = \\{\\text{year1}, \\text{year2}, \\text{year3}, \\text{year4}, \\text{year5} \\}$  \n",
    "\n",
    "\n",
    "$t \\in H_{2,4} = \\{year2, year3, year4 \\} $\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$\\text{demand}_{j} \\in \\mathbb{R}^+$: Exogenous demand of goods in industry $j$.\n",
    "\n",
    "$\\text{initial_stock}_{j} \\in \\mathbb{R}^+$: Stocks of goods in industry $j$ available at the beginning of year 1.\n",
    "\n",
    "$\\text{in_out_prod}_{i,j} \\in \\mathbb{R}^+$: Input of good $i$ required to produce one unit of good $j$ in the following year.\n",
    "\n",
    "$\\text{in_out_cap}_{i,j} \\in \\mathbb{R}^+$: Input of good $i$ in current year results in a permanent increase in productive capacity of good $j$  two years later.\n",
    "\n",
    "$\\text{industry_cap}_{j} \\in \\mathbb{R}^+$: Industry $j$  productive capacity at the beginning of year 1.\n",
    "\n",
    "$\\text{labor_prod}_{j} \\in \\mathbb{R}^+$: Labor requirements for the production of goods in industry  $j$.\n",
    "\n",
    "$\\text{labor_extra_cap}_{j} \\in \\mathbb{R}^+$: Labor requirements to permanently increase capacity for goods in industry $j$.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$\\text{production}_{j,t} \\in \\mathbb{R}^+$: Production of goods in industry $j$ available at year $t$, in millions of dollars.\n",
    "Note: goods of industry j, in millions of dollars, available in year t but produced in previous year.\n",
    "\n",
    "$\\text{stock}_{j,t} \\in \\mathbb{R}^+$: Stock level of industry $j$ at the end of year $t$, in millions of dollars.\n",
    "\n",
    "$\\text{extra_cap}_{j,t} \\in \\mathbb{R}^+$: Extra capacity for industry   $j$ becoming effective at the beginning of year $t$, in millions of dollars.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to build a realistic model, we need to think beyond year 5. Therefore, we make the following assumptions:\n",
    "\n",
    "* External demand is constant up to year 5 and beyond.\n",
    "* Stock level remains constant at year 5 and beyond.\n",
    "* Capacity cannot be increased after year 5.\n",
    "\n",
    "Consequently, for year 6 and beyond, we can assume a static model of the economy. For this static model, the production of goods in each industry $i \\in \\text{Industries}$ can be computed as follows.\n",
    "\n",
    "\\begin{equation}\n",
    "x_{i} = \\text{demand}_{i} + \\sum_{j \\in \\text{Industries} } \\text{in_out_prod}_{i,j} * x_{j}\n",
    "\\end{equation}\n",
    "\n",
    "Where $x_{i}$ is the production of goods in industry $i$.\n",
    "\n",
    "Solving this system of equations gives lower bounds for the production of goods in years beyond year 5. Then, the production of goods of  industry $j \\in \\text{Industries}$ in year 6 and beyond can be defined by the following constraints. \n",
    "\n",
    "$$\n",
    "\\text{production}_{j,t} \\geq x_{j} \\quad \\forall j \\in \\text{Industries}, \\; t=6. \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{extra_cap}_{j,t} = 0 \\quad \\forall j \\in \\text{Industries}, \\; t=6.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constraints\n",
    "\n",
    "**Balance equation year 1**: The initial stocks in industry $i$ should be equal to the total demand (internal and external), plus the extra capacity to build, plus the stock level at the end of year 1 of the goods in that industry.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{initial_stock}_{i} = \\sum_{j \\in \\text{Industries} } \\text{in_out_prod}_{i,j}*\\text{production}_{j,2}\n",
    "+ \\text{demand}_{i} + \\sum_{j \\in \\text{Industries} } \\text{in_out_cap}_{i,j}*\\text{extra_cap}_{j,3} +\n",
    " \\text{stock}_{i,1}\n",
    "\\end{equation}\n",
    "\n",
    "**Balance equation for year 2, year 3, year 4**: The production of goods in industry $i$ at year t plus the stocks at the end of year t-1 should be equal to the total demand (internal and external), plus the extra capacity to build, plus the stock level at the end of year t of the goods in that industry. \n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{production}_{j,t} + \\text{stock}_{i,t-1} =\n",
    "\\sum_{j \\in \\text{Industries} } \\text{in_out_prod}_{i,j}*\\text{production}_{j,t+1} + \\text{demand}_{i} +\n",
    "\\sum_{j \\in \\text{Industries} } \\text{in_out_cap}_{i,j}*\\text{extra_cap}_{j,t+2} \n",
    "+ \\text{stock}_{i,t} \\quad \\forall t \\in H_{2,4}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "**Balance equation year 5**: The production of goods in industry $i$ plus the stocks at the end of year 4 should be equal to the total demand (internal and external), plus the stock level at the end of year 5 of the goods in that industry.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{production}_{j,5} + \\text{stock}_{i,4} =\n",
    "\\sum_{j \\in \\text{Industries} } \\text{in_out_prod}_{i,j}*\\text{production}_{j,6}  +\n",
    "\\text{demand}_{i} + \\text{stock}_{i,5}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "**End of horizon constraints**: \n",
    "\n",
    "$$\n",
    "\\text{production}_{j,t} \\geq x_{j} \\quad \\forall j \\in \\text{Industries}, \\; t=6 \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{extra_cap}_{j,t} = 0 \\quad \\forall j \\in \\text{Industries}, \\; t=6\n",
    "$$\n",
    "\n",
    "Where $x_{j}$ is the solution of the static model.\n",
    "\n",
    "\n",
    "**Productive capacity constraints**: These constraints ensure that the production of goods for each industry $j$ during the planning horizon do not exceed the total production capacity at that year. \n",
    "\n",
    "\\begin{equation}\n",
    "\\text{production}_{j,t} \\leq \\text{base_cap}_{j} + \\sum_{\\tau \\leq t} \\text{extra_cap}_{j,\\tau} \\quad \\forall t \\in \\text{Horizon}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective Function\n",
    "\n",
    "**Labor utilization**: Maximize labor employment.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Maximize} \\quad Z =\n",
    "\\sum_{t \\in \\text{fiveYears} } \\sum_{j \\in \\text{Industries} } \\text{labor_prod}_{j}*\\text{production}_{j,t} +\n",
    "\\sum_{t \\in \\text{fiveYears} } \\sum_{j \\in \\text{Industries} } \\text{labor_extra_cap}_{j}*\\text{extra_cap}_{j,t}\n",
    "\\end{equation}\n",
    "\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",
    "from itertools import product\n",
    "\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.7.0 & Gurobi 9.1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input Data\n",
    "We define all the input data for the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input-output matrix for the production of goods of each industry\n",
    "\n",
    "arcs, inout_prod = gp.multidict({\n",
    "    ('coal', 'coal'): 0.1,\n",
    "    ('coal', 'steel'): 0.5,\n",
    "    ('coal', 'transport'): 0.4,\n",
    "    ('steel', 'coal'): 0.1,\n",
    "    ('steel', 'steel'): 0.1,\n",
    "    ('steel', 'transport'): 0.2,\n",
    "    ('transport', 'coal'): 0.2,\n",
    "    ('transport', 'steel'): 0.1,\n",
    "    ('transport', 'transport'): 0.2\n",
    "})\n",
    "\n",
    "# Labor requirements for the production of goods of each industry\n",
    "labor_prod = dict({'coal': 0.6,\n",
    "                   'steel': 0.3,\n",
    "                   'transport': 0.2})\n",
    "\n",
    "# input-output matrix to create extra capacity for each industry\n",
    "\n",
    "arcs, inout_cap = gp.multidict({\n",
    "    ('coal', 'coal'): 0.1,\n",
    "    ('coal', 'steel'): 0.7,\n",
    "    ('coal', 'transport'): 0.9,\n",
    "    ('steel', 'coal'): 0.1,\n",
    "    ('steel', 'steel'): 0.1,\n",
    "    ('steel', 'transport'): 0.2,\n",
    "    ('transport', 'coal'): 0.2,\n",
    "    ('transport', 'steel'): 0.1,\n",
    "    ('transport', 'transport'): 0.2\n",
    "})\n",
    "\n",
    "# Labor requirements to increase capacity for each industry\n",
    "labor_extra_cap = dict({'coal': 0.4,\n",
    "                  'steel': 0.2,\n",
    "                  'transport': 0.1})\n",
    "\n",
    "# Initial stock, initial capacity, and demand of goods for each industry\n",
    "\n",
    "industries, stock0, capacity0, demand = gp.multidict({\n",
    "    'coal': [250,300,60],\n",
    "    'steel': [180,350,60],\n",
    "    'transport': [200,280,30]\n",
    "})\n",
    "\n",
    "# Time Horizons\n",
    "horizon = [1,2,3,4,5,6]\n",
    "fiveYears = [1,2,3,4,5]\n",
    "years2_4 = [2,3,4]\n",
    "\n",
    "# Computed parameters\n",
    "i2h = set(product(industries, horizon))\n",
    "i2f = set(product(industries, fiveYears))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing\n",
    "\n",
    "We assume a static model of the economy, for year 6 and beyond. This static model is defined by a system of equations. The decision variables are the production of goods of each industry for the static model of the economy. We formulate and solve this model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n",
      "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 3 rows, 3 columns and 9 nonzeros\n",
      "Model fingerprint: 0x5b989667\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e-01, 9e-01]\n",
      "  Objective range  [0e+00, 0e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [3e+01, 6e+01]\n",
      "Presolve removed 3 rows and 3 columns\n",
      "Presolve time: 0.05s\n",
      "Presolve: All rows and columns removed\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    0.0000000e+00   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 0 iterations and 0.05 seconds\n",
      "Optimal objective  0.000000000e+00\n",
      "\n",
      "\n",
      "_________________________________________________________________________________\n",
      "The production of goods by industry for the static model of the economy is:\n",
      "_________________________________________________________________________________\n",
      "Generate $166.40 million dollars of coal \n",
      "Generate $105.67 million dollars of steel \n",
      "Generate $92.31 million dollars of transport \n"
     ]
    }
   ],
   "source": [
    "static = gp.Model('StaticModel')\n",
    "\n",
    "static_prod = static.addVars(industries, name=\"static_prod\")\n",
    "\n",
    "# Static model balance equations\n",
    "\n",
    "static_balance = static.addConstrs(\n",
    "    (static_prod[i] - gp.quicksum(inout_prod[i,j]*static_prod[j] for j in industries) \n",
    "     == demand[i] for i in industries  ), name='static_balance' )\n",
    "\n",
    "# we define a constant objective function to solve the system of balance equations for the static model.\n",
    "\n",
    "static.setObjective(0)\n",
    "\n",
    "# Verify model formulation\n",
    "\n",
    "static.write(\"StaticModel.lp\")\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "static.optimize()\n",
    "\n",
    "# Print solution of static model\n",
    "\n",
    "print(\"\\n\\n_________________________________________________________________________________\")\n",
    "print(f\"The production of goods by industry for the static model of the economy is:\")\n",
    "print(\"_________________________________________________________________________________\")\n",
    "for i in industries:\n",
    "    if (static_prod[i].x > 1e-6):\n",
    "        dollars_static_prod = '${:,.2f}'.format(static_prod[i].x)\n",
    "        print(f\"Generate {dollars_static_prod} million dollars of {i} \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "We create a model and the variables. For each year in the planning horizon and industry, the variables capture the production of goods, the stock level of goods, and the increase in production capacity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = gp.Model('EconomicPlanning')\n",
    "\n",
    "# Decision variables\n",
    "production = model.addVars(i2h, ub=capacity0, name=\"production\")\n",
    "\n",
    "# There is no production at year 1\n",
    "model.setAttr('ub', production.select('*',1), 0)    \n",
    "\n",
    "stock = model.addVars(i2f, name=\"stock\")\n",
    "extra_cap = model.addVars(i2h, name=\"extra_cap\")\n",
    "\n",
    "\n",
    "# No extra capacity can be generated at years 1 and 2\n",
    "for j,t in i2h:\n",
    "    if t < 3:\n",
    "        extra_cap[j,t].ub = 0  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first constraint is a balance equation for year 1 of the planning horizon. The initial stocks in industry  $i$  should be equal to the total demand (internal and external), plus the extra capacity to build, plus the stock level at the end of year 1 of the goods in that industry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Year 1 balance equations \n",
    "\n",
    "balance1 = model.addConstrs( ( stock0[i] == gp.quicksum(inout_prod[i,j]*production[j,2] for j in industries)  \n",
    "                              + gp.quicksum(inout_cap[i,j]*extra_cap[j,3] for j in industries ) \n",
    "                              + demand[i] + stock[i,1] for i in industries ), name='balance1' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following constraints are balance equations for years 2,3, and 4. The production of goods in industry  i  at year t plus the stocks at the end of year t-1 should be equal to the total demand (internal and external), plus the extra capacity to build, plus the stock level at the end of year t of the goods in that industry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Balance equations for years 2,3, and 4\n",
    "\n",
    "balance_t = model.addConstrs(( production[i, year] + stock[i,year-1]  == \n",
    "                              gp.quicksum(inout_prod[i,j]*production[j, year + 1] for j in industries) \n",
    "                              + gp.quicksum(inout_cap[i,j]*extra_cap[j, year + 2] for j in industries ) \n",
    "                              + demand[i] + stock[i, year] for i in industries for year in years2_4 ), name='balance_t' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The constraint for year 5 ensures that the production of goods in industry  i  plus the stocks at the end of year 4 should be equal to the total demand (internal and external), plus the stock level at the end of year 5 of the goods in that industry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Balance equations for year 5\n",
    "\n",
    "balance5 = model.addConstrs( (production[i, 5] + stock[i,4]  == \n",
    "                              gp.quicksum(inout_prod[i,j]*production[j,6] for j in industries) \n",
    "                              + demand[i] + stock[i,5] for i in industries ), name='balance5' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to enforce the constraints of a static model of the economy for year 6 and beyond.\n",
    "* Production of goods at each industry should be larger or equal to the industry production of the static model of the economy.\n",
    "\n",
    "* There is no increased capacity for year 6 and beyond."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Steady state production for year 6 and beyond\n",
    "\n",
    "steadyProduction = model.addConstrs((production[j,6] - static_prod[j].x >= 0 for j in industries ), name='steadyProduction')\n",
    "\n",
    "# Zero increased capacity for year 6 and beyond\n",
    "for j,t in i2h:\n",
    "    if t == 6:\n",
    "        extra_cap[j,t].ub = 0   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The productive capacity constraints ensure that the production of goods for each industry during the planning horizon do not exceed the total production capacity at that year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Productive capacity constraints\n",
    "\n",
    "capacityConstr = model.addConstrs(\n",
    "    (production[industry, year] - gp.quicksum(extra_cap[industry,t] for t in fiveYears if t <= year) \n",
    "     <= capacity0[industry] for industry,year in i2f ), name='capacityConstr' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective is to maximize labor employment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Maximize employment\n",
    "model.setObjective(\n",
    "    (gp.quicksum(labor_prod[j]*production[j,t] for j in industries for t in fiveYears) \n",
    "     + gp.quicksum(labor_extra_cap[j]*extra_cap[j,t] for j in industries for t in fiveYears) ), GRB.MAXIMIZE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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 33 rows, 51 columns and 183 nonzeros\n",
      "Model fingerprint: 0x194e5768\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e-01, 1e+00]\n",
      "  Objective range  [1e-01, 6e-01]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [3e+01, 4e+02]\n",
      "Presolve removed 9 rows and 18 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 24 rows, 33 columns, 126 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    2.5258351e+04   7.671320e+03   0.000000e+00      0s\n",
      "      25    1.9022169e+03   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 25 iterations and 0.01 seconds\n",
      "Optimal objective  1.902216910e+03\n"
     ]
    }
   ],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "model.write('DynamicModel.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Analysis\n",
    "\n",
    "The objective of maximizing employment results in a total manpower utilization of $\\$1,902.22 $ millions. With this objective, the coal industry is boosted in view of its heavy labor requirements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_______________________________________________________________________________________________\n",
      "The production of goods by industry and year for the dynamic Leontief model of the economy is:\n",
      "_______________________________________________________________________________________________\n"
     ]
    },
    {
     "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>coal</th>\n",
       "      <th>steel</th>\n",
       "      <th>transport</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Year1</th>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year2</th>\n",
       "      <td>$300.00</td>\n",
       "      <td>$136.78</td>\n",
       "      <td>$211.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year3</th>\n",
       "      <td>$370.87</td>\n",
       "      <td>$231.91</td>\n",
       "      <td>$217.39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year4</th>\n",
       "      <td>$368.03</td>\n",
       "      <td>$209.57</td>\n",
       "      <td>$275.67</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year5</th>\n",
       "      <td>$961.06</td>\n",
       "      <td>$350.00</td>\n",
       "      <td>$92.31</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          coal    steel transport\n",
       "Year1    $0.00    $0.00     $0.00\n",
       "Year2  $300.00  $136.78   $211.30\n",
       "Year3  $370.87  $231.91   $217.39\n",
       "Year4  $368.03  $209.57   $275.67\n",
       "Year5  $961.06  $350.00    $92.31"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Output reports\n",
    "\n",
    "print(\"_______________________________________________________________________________________________\")\n",
    "print(f\"The production of goods by industry and year for the dynamic Leontief model of the economy is:\")\n",
    "print(\"_______________________________________________________________________________________________\")\n",
    "\n",
    "goods = {}\n",
    " \n",
    "for i in industries:\n",
    "    my_list = []\n",
    "    for t in fiveYears:\n",
    "        my_list.append('${:,.2f}'.format( production[i,t ].x ) )\n",
    "    goods[i] = my_list\n",
    "\n",
    "goods_production = pd.DataFrame(goods, index=[\"Year1\", \"Year2\", \"Year3\", \"Year4\", \"Year5\"])\n",
    "goods_production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_______________________________________________________________________________________________\n",
      "The productive capacity by industry and year for the dynamic Leontief model of the economy is:\n",
      "_______________________________________________________________________________________________\n"
     ]
    },
    {
     "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>coal</th>\n",
       "      <th>steel</th>\n",
       "      <th>transport</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Year1</th>\n",
       "      <td>300</td>\n",
       "      <td>350</td>\n",
       "      <td>280</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year2</th>\n",
       "      <td>300</td>\n",
       "      <td>350</td>\n",
       "      <td>280</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year3</th>\n",
       "      <td>371</td>\n",
       "      <td>350</td>\n",
       "      <td>280</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year4</th>\n",
       "      <td>371</td>\n",
       "      <td>350</td>\n",
       "      <td>280</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year5</th>\n",
       "      <td>961</td>\n",
       "      <td>350</td>\n",
       "      <td>280</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       coal  steel  transport\n",
       "Year1   300    350        280\n",
       "Year2   300    350        280\n",
       "Year3   371    350        280\n",
       "Year4   371    350        280\n",
       "Year5   961    350        280"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"_______________________________________________________________________________________________\")\n",
    "print(f\"The productive capacity by industry and year for the dynamic Leontief model of the economy is:\")\n",
    "print(\"_______________________________________________________________________________________________\")\n",
    "\n",
    "# Compute cummulative capacity\n",
    "totalCap = {}\n",
    " \n",
    "for i in industries:\n",
    "    amount = capacity0[i]\n",
    "    my_list = []\n",
    "    for t in fiveYears:\n",
    "        amount += extra_cap[i,t].x\n",
    "        my_list.append(round(amount))\n",
    "    totalCap[i] = my_list\n",
    "    \n",
    "extra_capacity = pd.DataFrame(totalCap, index=[\"Year1\", \"Year2\", \"Year3\", \"Year4\", \"Year5\"])\n",
    "extra_capacity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "____________________________________________________________________________________\n",
      "Stock level by industry at the end of year for the dynamic Leontif model of the economy is:\n",
      "____________________________________________________________________________________\n"
     ]
    },
    {
     "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>coal</th>\n",
       "      <th>steel</th>\n",
       "      <th>transport</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Year1</th>\n",
       "      <td>$0.00</td>\n",
       "      <td>$26.97</td>\n",
       "      <td>$39.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year2</th>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$80.34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year3</th>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year4</th>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Year5</th>\n",
       "      <td>$794.66</td>\n",
       "      <td>$244.33</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          coal    steel transport\n",
       "Year1    $0.00   $26.97    $39.89\n",
       "Year2    $0.00    $0.00    $80.34\n",
       "Year3    $0.00    $0.00     $0.00\n",
       "Year4    $0.00    $0.00     $0.00\n",
       "Year5  $794.66  $244.33     $0.00"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"____________________________________________________________________________________\")\n",
    "print(f\"Stock level by industry at the end of year for the dynamic Leontif model of the economy is:\")\n",
    "print(\"____________________________________________________________________________________\")\n",
    "\n",
    "inv = {}\n",
    " \n",
    "for i in industries:\n",
    "    my_list = []\n",
    "    for t in fiveYears:\n",
    "        my_list.append('${:,.2f}'.format( stock[i,t ].x ) )\n",
    "    inv[i] = my_list\n",
    "\n",
    "stock_level = pd.DataFrame(inv, index=[\"Year1\", \"Year2\", \"Year3\", \"Year4\", \"Year5\"])\n",
    "stock_level"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Wagner, H.M. (1957) A linear programming solution to dynamic Leontief type models. Management Science, 3, 234–254.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  }
 ],
 "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
}