tutorial/westeros/westeros_baseline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Welcome to the MESSAGEix framework\n",
    "\n",
    "### *Integrated Assessment Modeling for the 21st Century*\n",
    "\n",
    "Please refer to the [user guidelines](https://github.com/iiasa/message_ix/blob/master/NOTICE.rst)\n",
    "for additional information on using *MESSAGEix*, including the recommended citation and how to name new models.\n",
    "\n",
    "This tutorial is based on a presentation by Matthew Gidden ([@gidden](https://github.com/gidden))\n",
    "for a summer school at the the **Centre National de la Recherche Scientifique (CNRS)**\n",
    "on *Integrated Assessment Modeling* in June 2018."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Scope of this tutorial: Building a Simple Energy Model\n",
    "\n",
    "The goal of this tutorial is to build a simple energy model using *MESSAGEix* with minimal features that can be expanded in future tutorials.  \n",
    "\n",
    "We will build the model component by component, focusing on both the **how** (code implementation) and **why** (mathematical formulation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Online documentation\n",
    "\n",
    "The full framework documentation is available at [https://messageix.iiasa.ac.at](https://messageix.iiasa.ac.at)\n",
    "    \n",
    "<img src='_static/doc_page.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Simple installation\n",
    "\n",
    "When working with Python, you can easily install `MESSAGEix` yourself and get all the tutorials:\n",
    "\n",
    "```shell\n",
    "$ conda install -c conda-forge message-ix # install the MESSAGEix package and all dependencies\n",
    "\n",
    "$ messageix-dl # download all tutorials to your current directory\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A stylized reference energy system model for Westeros\n",
    "\n",
    "This tutorial is based on the country of Westeros from the TV show \"Game of Thrones\".\n",
    "\n",
    "<table align='center'><tr><td><img src='_static/westeros.jpg' width='150'></td><td><img src='_static/base_res.png'></td></tr></table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MESSAGEix: the mathematical paradigm\n",
    "\n",
    "At its core, *MESSAGEix* is an optimization problem:\n",
    "\n",
    "> $\\min \\quad ~c^T \\cdot x$  \n",
    "> $~s.t. \\quad A \\cdot x \\leq b$\n",
    "\n",
    "More explicitly, the model...\n",
    "- optimizes an **objective function**, nominally minimizing total **system costs**\n",
    "- under a system of **constraints** equations (inequalities or equality conditions)\n",
    "\n",
    "The mathematical implementation includes a number of features that make it particularly geared towards the modelling of *energy-water-land systems* in the context of *climate change mitigation and sustainable development*.\n",
    "\n",
    "Throughout this document, the mathematical formulation follows the convention that\n",
    "- decision **VARIABLES** ($x$) are capitalized\n",
    "- **parameters** ($A$, $b$) are lower case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MESSAGEix: connected to the *ix modeling platform*\n",
    "\n",
    "The *modeling platform for integrated and cross-cutting analysis* (ixmp) provides a powerful framework for working with scenarios, including a database infrastucture for data version control and interfaces to scientific programming languages.\n",
    "\n",
    "<img src='_static/message_ixmp.png' width='700'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ready, steady, go!\n",
    "\n",
    "First, we import all the packages we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import ixmp\n",
    "import message_ix\n",
    "\n",
    "from message_ix.utils import make_df\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The `Platform` is your connection to a database for storing model input data and scenario results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mp = ixmp.Platform(dbtype='HSQLDB')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Once connected, we create a new `Scenario` to build our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario = message_ix.Scenario(mp, model='Westeros Electrified', \n",
    "                               scenario='baseline', version='new')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model Structure\n",
    "\n",
    "We start by defining basic characteristics of the model, including time, space, and the energy system structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The model horizon will span 3 decades. Let's assume that we're far in the future after the events of A Song of Ice and Fire (which occur ~300 years after Aegon the conqueror).\n",
    "\n",
    "| Math Notation | Model Meaning                |\n",
    "|---------------|------------------------------|\n",
    "| $y \\in Y^H$   | time periods in history      |\n",
    "| $y \\in Y^M$   | time periods in model horizon|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "history = [690]\n",
    "model_horizon = [700, 710, 720]\n",
    "scenario.add_horizon({'year': history + model_horizon,  \n",
    "                      'firstmodelyear': model_horizon[0]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our model will have a single `node`, i.e., its spatial dimension.\n",
    "\n",
    "\n",
    "| Math Notation | Model Meaning|\n",
    "|---------------|--------------|\n",
    "| $n \\in N$     | node         |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "country = 'Westeros'\n",
    "scenario.add_spatial_sets({'country': country})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And we fill in the energy system's `commodities`, `levels`, `technologies`, and `mode` (defining how certain technologies operate). \n",
    "\n",
    "\n",
    "| Math Notation | Model Meaning|\n",
    "|---------------|--------------|\n",
    "| $c \\in C$     | commodity    |\n",
    "| $l \\in L$     | level        |\n",
    "| $t \\in T$     | technology   |\n",
    "| $m \\in M$     | mode         |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.add_set(\"commodity\", [\"electricity\", \"light\"])\n",
    "scenario.add_set(\"level\", [\"secondary\", \"final\", \"useful\"])\n",
    "scenario.add_set(\"technology\", ['coal_ppl', 'wind_ppl', 'grid', 'bulb'])\n",
    "scenario.add_set(\"mode\", \"standard\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supply and Demand (or Balancing Commodities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The fundamental premise of the model is to satisfy demand for energy (services).\n",
    "To first order, demand for services like electricity track with economic productivity (GDP).\n",
    "We define a GDP profile similar to first-world GDP growth from [1900-1930](https://en.wikipedia.org/wiki/List_of_regions_by_past_GDP):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gdp_profile = pd.Series([1., 1.5, 1.9],\n",
    "                        index=pd.Index(model_horizon, name='Time'))\n",
    "gdp_profile.plot(title='Demand')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The `COMMODITY_BALANCE` equation ensures that `demand` for each `commodity` is met at each `level` in the energy system.\n",
    "\n",
    "$\\sum_{\\substack{n^L,t,m \\\\ y^V \\leq y}} \\text{output}_{n^L,t,y^V,y,m,n,c,l} \\cdot \\text{ACT}_{n^L,t,y^V,y,m}$\n",
    "$- \\sum_{\\substack{n^L,t,m, \\\\ y^V \\leq y}} \\text{input}_{n^L,t,y^V,y,m,n,c,l} \\cdot \\text{ACT}_{n^L,t,m,y}$  \n",
    "$\\geq \\text{demand}_{n,c,l,y} \\quad \\forall \\ l \\in L$\n",
    "\n",
    "While `demand` must be met, it can also be *exceeded* allowing the model to plan for future periods of demand for storable commodities.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "First we establish demand. Let's assume\n",
    "\n",
    "- 40 million people in [300 AC](https://atlasoficeandfireblog.wordpress.com/2016/03/06/the-population-of-the-seven-kingdoms/)\n",
    "- similar population growth to Earth in the same time frame [(~factor of 12)](https://en.wikipedia.org/wiki/World_population_estimates)\n",
    "- a per capita demand for electricity of 1000 kWh\n",
    "- and 8760 hours in a year (of course!)\n",
    "\n",
    "Then we can add the demand parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Note present day: [~72000 GWh in Austria](https://www.iea.org/statistics/?country=AUSTRIA&year=2016&category=Energy%20consumption&indicator=undefined&mode=chart&dataTable=INDICATORS) with population [~8.7M](http://www.austria.org/population/) which is ~8300 kWh per capita"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "demand_per_year = 40 * 12 * 1000 / 8760\n",
    "light_demand = pd.DataFrame({\n",
    "        'node': country,\n",
    "        'commodity': 'light',\n",
    "        'level': 'useful',\n",
    "        'year': model_horizon,\n",
    "        'time': 'year',\n",
    "        'value': (100 * gdp_profile).round(),\n",
    "        'unit': 'GWa',\n",
    "    })"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "`light_demand` illustrates the data format for MESSAGEix parameters. It is a `pandas.DataFrame` containing three types of information in a specific format:\n",
    "\n",
    "- A \"value\" column containing the numerical values for this parameter.\n",
    "- A \"unit\" column.\n",
    "- Other columns (\"node\", \"commodity\", \"level\", \"time\") that indicate the key to which each value applies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "light_demand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.add_par(\"demand\", light_demand)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In order to add values for the reference energy system, we define some common keys.\n",
    "\n",
    "- **Input** quantities require `_origin` keys that specify where the inputs are *received from*.\n",
    "- **Output** quantities require `_dest` keys that specify where the outputs are *transferred to*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "year_df = scenario.vintage_and_active_years()\n",
    "vintage_years, act_years = year_df['year_vtg'], year_df['year_act']\n",
    "\n",
    "base = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': vintage_years,\n",
    "    'year_act': act_years,\n",
    "    'mode': 'standard',\n",
    "    'time': 'year',\n",
    "    'unit': '-',\n",
    "}\n",
    "\n",
    "base_input = make_df(base, node_origin=country, time_origin='year')\n",
    "base_output = make_df(base, node_dest=country, time_dest='year')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Working backwards along the Reference Energy System, we can add connections for the `bulb`. A light bulb…\n",
    "\n",
    "- receives *input* in the form of the \"electricity\" *commodity* at the \"final [energy]\" *level*, and\n",
    "- *outputs* the commodity \"light\" at the \"useful [energy]\" level.\n",
    "\n",
    "Because there are no other commodities or forms of useful energy in Westeros, we specify that 100% of the input and output are from/to these commodities/levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "bulb_out = make_df(base_output, technology='bulb', commodity='light', \n",
    "                   level='useful', value=1.0)\n",
    "scenario.add_par('output', bulb_out)\n",
    "\n",
    "bulb_in = make_df(base_input, technology='bulb', commodity='electricity',  \n",
    "                  level='final', value=1.0)\n",
    "scenario.add_par('input', bulb_in)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Next, the electrical `grid`, which…\n",
    "\n",
    "- receives electricity at the \"secondary\" energy level.\n",
    "- also outputs electricity, but at the \"final\" energy level (to be used by the light bulb).\n",
    "\n",
    "Because the grid has transmission losses, only 90% of the input electricity is available as output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "grid_efficiency = 0.9\n",
    "grid_out = make_df(base_output, technology='grid', commodity='electricity', \n",
    "                   level='final', value=grid_efficiency)\n",
    "scenario.add_par('output', grid_out)\n",
    "\n",
    "grid_in = make_df(base_input, technology='grid', commodity='electricity',\n",
    "                  level='secondary', value=1.0)\n",
    "scenario.add_par('input', grid_in)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And finally, our power plants. The model does not include the fossil resources used as `input` for coal plants; however, costs of coal extraction are included in the parameter $variable\\_cost$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "coal_out = make_df(base_output, technology='coal_ppl', commodity='electricity', \n",
    "                   level='secondary', value=1.)\n",
    "scenario.add_par('output', coal_out)\n",
    "\n",
    "wind_out = make_df(base_output, technology='wind_ppl', commodity='electricity', \n",
    "                   level='secondary', value=1.)\n",
    "scenario.add_par('output', wind_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Operational Constraints and Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The model has a number of \"reality\" constraints, which relate built *capacity* (`CAP`) to available power, or the *activity* (`ACT`) of that technology.\n",
    "\n",
    "The **capacity constraint** limits the total time a each technology can operate.\n",
    "\n",
    "$$\\sum_{m} \\text{ACT}_{n,t,y^V,y,m,h}\n",
    "    \\leq \\text{duration_time}_{h} \\cdot \\text{capacity_factor}_{n,t,y^V,y,h} \\cdot \\text{CAP}_{n,t,y^V,y}\n",
    "    \\quad t \\ \\in \\ T^{INV}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This requires us to provide the `capacity_factor` for each technology. Here, we call `make_df()` and `add_par()` in a loop to execute similar code for three technologies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_capacity_factor = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': vintage_years,\n",
    "    'year_act': act_years,\n",
    "    'time': 'year',\n",
    "    'unit': '-',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "capacity_factor = {\n",
    "    'coal_ppl': 1,\n",
    "    'wind_ppl': 1,\n",
    "    'bulb': 1, \n",
    "}\n",
    "\n",
    "for tec, val in capacity_factor.items():\n",
    "    df = make_df(base_capacity_factor, technology=tec, value=val)\n",
    "    scenario.add_par('capacity_factor', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The model can further be provided `technical_lifetime`s in order to properly manage deployed capacity and related costs via the **capacity maintenance** constraint:\n",
    "\n",
    "$\\text{CAP}_{n,t,y^V,y} \\leq \\text{remaining_capacity}_{n,t,y^V,y} \\cdot \\text{value} \\quad \\forall \\quad t \\in T^{INV}$\n",
    "\n",
    "where `value` can take different forms depending on what time period is considered:\n",
    "\n",
    "| Value                               | Condition                                           |\n",
    "|-------------------------------------|-----------------------------------------------------|\n",
    "| $\\Delta_y \\text{historical_new_capacity}_{n,t,y^V}$ | $y$ is first model period           |\n",
    "| $\\Delta_y \\text{CAP_NEW}_{n,t,y^V}$ | $y = y^V$                                           |\n",
    "| $\\text{CAP}_{n,t,y^V,y-1}$          | $0 < y - y^V < \\text{technical_lifetime}_{n,t,y^V}$ |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_technical_lifetime = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': model_horizon,\n",
    "    'unit': 'y',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lifetime = {\n",
    "    'coal_ppl': 20,\n",
    "    'wind_ppl': 20,\n",
    "    'bulb': 1,\n",
    "}\n",
    "\n",
    "for tec, val in lifetime.items():\n",
    "    df = make_df(base_technical_lifetime, technology=tec, value=val)\n",
    "    scenario.add_par('technical_lifetime', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Technological Diffusion and Contraction\n",
    "\n",
    "We know from historical precedent that energy systems can not be transformed instantaneously. Therefore, we use a family of constraints on activity and capacity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$\\sum_{y^V \\leq y,m} \\text{ACT}_{n,t,y^V,y,m,h} \\leq$  \n",
    "$\\text{initial_activity_up}_{n,t,y,h}\n",
    "        \\cdot \\frac{ \\Big( 1 + growth\\_activity\\_up_{n,t,y,h} \\Big)^{|y|} - 1 }\n",
    "                   { growth\\_activity\\_up_{n,t,y,h} }+ \\Big( 1 + growth\\_activity\\_up_{n,t,y,h} \\Big)^{|y|} \\cdot \\Big( \\sum_{y^V \\leq y-1,m} ACT_{n,t,y^V,y-1,m,h} + \\sum_{m} historical\\_activity_{n,t,y-1,m,h}\\Big)$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This example limits the ability for technologies to **grow**. To do so, we  need to provide `growth_activity_up` values for each technology that we want to model as being diffusion constrained. Here, we set this constraint at 10% per year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_growth = {\n",
    "    'node_loc': country,\n",
    "    'year_act': model_horizon,\n",
    "    'time': 'year',\n",
    "    'unit': '-',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "growth_technologies = [\n",
    "    \"coal_ppl\", \n",
    "    \"wind_ppl\", \n",
    "]\n",
    "\n",
    "for tec in growth_technologies:\n",
    "    df = make_df(base_growth, technology=tec, value=0.1) \n",
    "    scenario.add_par('growth_activity_up', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Defining an Energy Mix\n",
    "\n",
    "To model the transition of an energy system, one must start with the existing system which are defined by the parameters `historical_activity` and `historical_capacity`. These parameters define the energy mix before the model horizon.  \n",
    "\n",
    "We begin by defining a few key values:\n",
    "\n",
    "- how much useful energy was needed\n",
    "- how much final energy was generated\n",
    "- and the mix for different technologies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "historic_demand = 0.85 * demand_per_year\n",
    "historic_generation = historic_demand / grid_efficiency\n",
    "coal_fraction = 0.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_capacity = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': history,\n",
    "    'unit': 'GWa',\n",
    "}\n",
    "\n",
    "base_activity = {\n",
    "    'node_loc': country,\n",
    "    'year_act': history,\n",
    "    'mode': 'standard',\n",
    "    'time': 'year',\n",
    "    'unit': 'GWa',\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Then, we can define the **activity** and **capacity** in the historic period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "old_activity = {\n",
    "    'coal_ppl': coal_fraction * historic_generation,\n",
    "    'wind_ppl': (1 - coal_fraction) * historic_generation,\n",
    "}\n",
    "\n",
    "for tec, val in old_activity.items():\n",
    "    df = make_df(base_activity, technology=tec, value=val)\n",
    "    scenario.add_par('historical_activity', df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "act_to_cap = {\n",
    "    'coal_ppl': 1 / 10 / capacity_factor['coal_ppl'] / 2, # 20 year lifetime\n",
    "    'wind_ppl': 1 / 10 / capacity_factor['wind_ppl'] / 2,\n",
    "}\n",
    "\n",
    "for tec in act_to_cap:\n",
    "    value = old_activity[tec] * act_to_cap[tec]\n",
    "    df = make_df(base_capacity, technology=tec, value=value)\n",
    "    scenario.add_par('historical_new_capacity', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Objective Function\n",
    "\n",
    "The objective function drives the purpose of the optimization. Do we wish to seek maximum utility of the social planner, minimize carbon emissions, or something else? Classical IAMs seek to minimize total discounted system cost over space and time. \n",
    "\n",
    "$$\\min \\sum_{n,y \\in Y^{M}} \\text{interestrate}_{y} \\cdot \\text{COST_NODAL}_{n,y}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "First, let's add the interest rate parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "rate = [0.05] * len(model_horizon)\n",
    "unit = ['-'] * len(model_horizon)\n",
    "scenario.add_par(\"interestrate\", key=model_horizon, val=rate, unit=unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "`COST_NODAL` is comprised of a variety of costs related to the use of different technologies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Investment Costs\n",
    "\n",
    "Capital, or investment, costs are invoked whenever a new plant or unit is built\n",
    "\n",
    "$$\\text{inv_cost}_{n,t,y} \\cdot \\text{construction_time_factor}_{n,t,y} \\cdot \\text{CAP_NEW}_{n,t,y}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_inv_cost = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': model_horizon,\n",
    "    'unit': 'USD/GWa',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# in $ / kW\n",
    "costs = {\n",
    "    'coal_ppl': 500,\n",
    "    'wind_ppl': 1500,\n",
    "    'bulb': 5,\n",
    "}\n",
    "\n",
    "for tec, val in costs.items():\n",
    "    df = make_df(base_inv_cost, technology=tec, value=val)\n",
    "    scenario.add_par('inv_cost', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Fixed O&M Costs\n",
    "\n",
    "Fixed cost are only relevant as long as the capacity is active. This formulation allows to include the potential cost savings from early retirement of installed capacity.\n",
    "\n",
    "$$\\sum_{y^V \\leq y} \\text{fix_cost}_{n,t,y^V,y} \\cdot \\text{CAP}_{n,t,y^V,y}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_fix_cost = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': vintage_years,\n",
    "    'year_act': act_years,\n",
    "    'unit': 'USD/GWa',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# in $ / kW\n",
    "costs = {\n",
    "    'coal_ppl': 30,\n",
    "    'wind_ppl': 10,\n",
    "}\n",
    "\n",
    "for tec, val in costs.items():\n",
    "    df = make_df(base_fix_cost, technology=tec, value=val)\n",
    "    scenario.add_par('fix_cost', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Variable O&M Costs\n",
    "\n",
    "Variable Operation and Maintence costs are associated with the costs of actively running the plant. Thus, they are not applied if a plant is on standby (i.e., constructed, but not currently in use).\n",
    "\n",
    "$$\\sum_{\\substack{y^V \\leq y \\\\ m,h}} \\text{var_cost}_{n,t,y^V,y,m,h} \\cdot \\text{ACT}_{n,t,y^V,y,m,h} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "base_var_cost = {\n",
    "    'node_loc': country,\n",
    "    'year_vtg': vintage_years,\n",
    "    'year_act': act_years,\n",
    "    'mode': 'standard',\n",
    "    'time': 'year',\n",
    "    'unit': 'USD/GWa',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# in $ / MWh\n",
    "costs = {\n",
    "    'coal_ppl': 30,\n",
    "    'grid': 50,\n",
    "}\n",
    "\n",
    "for tec, val in costs.items():\n",
    "    df = make_df(base_var_cost, technology=tec, value=val)\n",
    "    scenario.add_par('var_cost', df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A full model will also have costs associated with\n",
    "\n",
    "- costs associated with technologies (investment, fixed, variable costs)\n",
    "- resource extraction: $\\sum_{c,g} \\ resource\\_cost_{n,c,g,y} \\cdot EXT_{n,c,g,y} $\n",
    "- emissions\n",
    "- land use (emulator): $\\sum_{s} land\\_cost_{n,s,y} \\cdot LAND_{n,s,y}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Time to Solve the Model\n",
    "\n",
    "First, we *commit* the model structure and input data (sets and parameters).\n",
    "In the `ixmp` backend, this creates a new model version in the database, which is assigned a number automatically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ixmp.logger().info('version number prior to commit: {}'.format(scenario.version))\n",
    "\n",
    "scenario.commit(comment='basic model of Westerosi electrification')\n",
    "\n",
    "ixmp.logger().info('version number prior committing to the database: {}'.format(scenario.version))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "An `ixmp` database can contain many scenarios, and possibly multiple versions of the same model and scenario name.\n",
    "These are distinguished by unique version numbers.\n",
    "\n",
    "To make it easier to retrieve the \"correct\" version (e.g., the latest one), you can set a specific scenario as the default version to use if the \"Westeros Electrified\" model is loaded from the `ixmp` database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.set_as_default()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "name": "solve-objective-value"
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "scenario.var('OBJ')['lvl']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plotting Results\n",
    "\n",
    "We make use of some custom code; see `tools.py` in the tutorial directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from tools import Plots\n",
    "p = Plots(scenario, country, firstyear=model_horizon[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Activity\n",
    "\n",
    "How much energy is generated in each time period from the different potential sources?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "p.plot_activity(baseyear=True, subset=['coal_ppl', 'wind_ppl'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Capacity\n",
    "\n",
    "Given how many new plants are built, how many are actually used?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "p.plot_capacity(baseyear=True, subset=['coal_ppl', 'wind_ppl'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Electricity Price\n",
    "\n",
    "And how much does the electricity cost? These prices are taken from the **dual variables** of the solution and are given the name **shadow prices**. They reflect the marginal price of electricity, taken from the most expensive producer.  \n",
    "\n",
    "Note that the price drop when the most expensive technology is no longer in the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "p.plot_prices(subset=['light'], baseyear=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Close the connection to the database\n",
    "\n",
    "When working with local HSQLDB database instances, you cannot connect to one database from multipe Jupyter notebooks (or processes) at the same time.\n",
    "\n",
    "If you want to easily switch between notebooks with connections to the same `ixmp` database, you need to close the connection in one notebook before initializing the platform using `ixmp.Platform(dbtype='HSQLDB')` in another notebook.\n",
    "\n",
    "After having closed the database connection, you can reopen it using\n",
    "```\n",
    "mp.open_db()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "mp.close_db()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Congratulations! \n",
    "\n",
    "You have built and run your very first MESSAGEix model. Welcome to the community!\n",
    "\n",
    "The next tutorials will introduce you to other features of the framework, including energy system constraints, emissions taxes, and other policy options.\n",
    "\n",
    "Check us out on Github https://github.com/iiasa/message_ix  \n",
    "and get in touch with us online https://groups.google.com/forum/message-ix ..."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}