teaching/01_Electricity_supply_of_an_island.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reminder (Before we start)\n",
    "Check whether the kernel at the top right is set to \"urbs\" in order to be able to run this script."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1: Electricity supply of an island"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "* Translate a mathematical optimization problem into a pyomo ConcreteModel/AbstractModel\n",
    "* Recognize the basic structure of an optimization model\n",
    "* Report the results of pyomo in different formats\n",
    "* Run and edit scripts using pyomo, pandas and matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mathematical formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with a simple example. Let's assume we have a gas power plant ($P_{gas}$ = 100 MW) and a biomass power plant ($P_{bm}$ = 30 MW) supplying an island. The cost of supplying 1 MWh of electricity using the gas power plant is EUR 50, whereas the cost of using biomass is 25 EUR/MWh. These costs include operation and maintenance costs and fuel costs. The efficiency of the power plants is already taken into account.\n",
    "\n",
    "We would like to minimize the cost of operating the system for a given demand of electricity $d(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\min \\quad 50s_{gas}(t) + 25s_{bm}(t)$$\n",
    "$$s.t. \\quad s_{gas}(t) + s_{bm}(t) \\geq d(t)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The supply from the power plant is non-negative:\n",
    "$$s_{gas}(t), s_{bm}(t) \\geq 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It cannot exceed the capacity of the power plants:\n",
    "$$s_{gas}(t) \\leq 100$$\n",
    "$$s_{bm}(t) \\leq 30$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further, we define the demand as follows:\n",
    "$$d(t) = [60, 100, 120, 80, 30]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### <span style=\"color:blue\">Task</span>\n",
    "Try to solve this problem with pen and paper!\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formulation as a pyomo ConcreteModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could solve this problem using a pyomo ConcreteModel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the object \"environ\" from the library \"pyomo\" which is already installed in our urbs environment.\n",
    "# Whenever we will use it, we will call it using its alias \"pyo\"\n",
    "import pyomo.environ as pyo\n",
    "\n",
    "# Let's create a ConcreteModel object and fill it with life!\n",
    "model = pyo.ConcreteModel()\n",
    "model.name = \"Example1\"\n",
    "\n",
    "## Variables\n",
    "# Our variable \"s\" (supply) has two dimensions: time and technology. It is always positive.\n",
    "model.s = pyo.Var([\"t1\", \"t2\", \"t3\", \"t4\", \"t5\"], [\"Gas\", \"Biomass\"], domain=pyo.NonNegativeReals)\n",
    "\n",
    "## Objective function\n",
    "# The objective is also a variable, albeit a special one which we will optimize.\n",
    "model.OBJ = pyo.Objective(expr=50*model.s[\"t1\", \"Gas\"] + 25*model.s[\"t1\", \"Biomass\"] +\\\n",
    "                               50*model.s[\"t2\", \"Gas\"] + 25*model.s[\"t2\", \"Biomass\"] +\\\n",
    "                               50*model.s[\"t3\", \"Gas\"] + 25*model.s[\"t3\", \"Biomass\"] +\\\n",
    "                               50*model.s[\"t4\", \"Gas\"] + 25*model.s[\"t4\", \"Biomass\"] +\\\n",
    "                               50*model.s[\"t5\", \"Gas\"] + 25*model.s[\"t5\", \"Biomass\"])\n",
    "\n",
    "## Constraints\n",
    "# The supply from the Gas power plant cannot exceed its capacity of 100 MW\n",
    "model.ConstraintGasCap1 = pyo.Constraint(expr = model.s[\"t1\", \"Gas\"] <= 100)\n",
    "model.ConstraintGasCap2 = pyo.Constraint(expr = model.s[\"t2\", \"Gas\"] <= 100)\n",
    "model.ConstraintGasCap3 = pyo.Constraint(expr = model.s[\"t3\", \"Gas\"] <= 100)\n",
    "model.ConstraintGasCap4 = pyo.Constraint(expr = model.s[\"t4\", \"Gas\"] <= 100)\n",
    "model.ConstraintGasCap5 = pyo.Constraint(expr = model.s[\"t5\", \"Gas\"] <= 100)\n",
    "\n",
    "# The supply from the Biomass power plant cannot exceed its capacity of 30 MW\n",
    "model.ConstraintBiomassCap1 = pyo.Constraint(expr = model.s[\"t1\", \"Biomass\"] <= 30)\n",
    "model.ConstraintBiomassCap2 = pyo.Constraint(expr = model.s[\"t2\", \"Biomass\"] <= 30)\n",
    "model.ConstraintBiomassCap3 = pyo.Constraint(expr = model.s[\"t3\", \"Biomass\"] <= 30)\n",
    "model.ConstraintBiomassCap4 = pyo.Constraint(expr = model.s[\"t4\", \"Biomass\"] <= 30)\n",
    "model.ConstraintBiomassCap5 = pyo.Constraint(expr = model.s[\"t5\", \"Biomass\"] <= 30)\n",
    "\n",
    "# The supply should at least be equal to the demand\n",
    "model.ConstraintDem1 = pyo.Constraint(expr = model.s[\"t1\", \"Gas\"] + model.s[\"t1\", \"Biomass\"] >= 60)\n",
    "model.ConstraintDem2 = pyo.Constraint(expr = model.s[\"t2\", \"Gas\"] + model.s[\"t2\", \"Biomass\"] >= 100)\n",
    "model.ConstraintDem3 = pyo.Constraint(expr = model.s[\"t3\", \"Gas\"] + model.s[\"t3\", \"Biomass\"] >= 120)\n",
    "model.ConstraintDem4 = pyo.Constraint(expr = model.s[\"t4\", \"Gas\"] + model.s[\"t4\", \"Biomass\"] >= 80)\n",
    "model.ConstraintDem5 = pyo.Constraint(expr = model.s[\"t5\", \"Gas\"] + model.s[\"t5\", \"Biomass\"] >= 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the LP mathematical problem that is solved to a file (optional)\n",
    "# Here, we are reporting the model itself, not its solution\n",
    "model.write(\"01_concrete_a.lp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### <span style=\"color:blue\">Task</span>\n",
    "Open the file \"01_concrete_a.lp\" with a text editor. Can you recognize the variables and constraints?\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try this now\n",
    "model.write(\"01_concrete_b.lp\", io_options={'symbolic_solver_labels': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's solve the model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We first load the solver\n",
    "opt = pyo.SolverFactory('glpk') # glpk: GNU Linear Programming Kit\n",
    "results = opt.solve(model)\n",
    "# First way of reporting the solution\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.Solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> For more on solver status and termination conditions:\n",
    "http://www.pyomo.org/blog/2015/1/8/accessing-solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### <span style=\"color:blue\">Task</span>\n",
    "1. Try to comment one or multiple constraints. What happens?\n",
    "2. Try to maximize instead of minimizing the costs. (Tip: add the option 'sense=pyo.maximize' into the objective function)\n",
    "3. How easy is it to add another power plant? Another time step?\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formulation as a pyomo AbstractModel\n",
    "\n",
    "One way to add flexibility is to write the problem abstractly. For example, the following equations represent a linear program (LP) to find optimal values for the vector $x$ (in our case, the hourly supply from the power plants) with parameters $c_j$ (costs), $a_{i,j}$ and $b_i$ (constraints):\n",
    "\n",
    "$$ \\begin{array}{lll} \\min & \\sum_{j=1}^n c_j x_{j,t} & \\\\\n",
    "s.t. & \\sum_{j=1}^n a_{i,j} x_{j,t} \\geq b_{i,t} & \\forall i = 1 \\ldots m\\\\ & x_{j,t} \\geq 0 & \\forall j = 1 \\ldots n\n",
    "\\end{array} $$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For that, there is the pyomo class AbstractModel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we are loading all the objects within \"pyomo.environ\" and we can use them directly.\n",
    "# The script will recognize their names. This is practical, but not a good coding style.\n",
    "# There could be issues if you use one of the \"reserved names\" that are already defined in the library.\n",
    "from pyomo.environ import *\n",
    "\n",
    "model = AbstractModel()\n",
    "\n",
    "# Sets\n",
    "model.I = Set() # we could define the dimensions, or let pyomo determine them from the data\n",
    "model.J = Set()\n",
    "model.T = Set()\n",
    "\n",
    "# Parameters\n",
    "model.a = Param(model.I, model.J)\n",
    "model.b = Param(model.I, model.T)\n",
    "model.c = Param(model.J)\n",
    "\n",
    "# Variables\n",
    "model.x = Var(model.J, model.T, domain=NonNegativeReals) # the variable is indexed by the set J and the set T\n",
    "\n",
    "# Objective function\n",
    "def obj_expression(model):\n",
    "    sigma = 0\n",
    "    for t in model.T:\n",
    "        for j in model.J:\n",
    "            sigma = sigma + model.c[j] * model.x[(j, t)]\n",
    "    return sigma\n",
    "\n",
    "model.OBJ = Objective(rule=obj_expression)\n",
    "\n",
    "# Constraints\n",
    "def ax_constraint_rule(model, i, t):\n",
    "    # return the expression for the constraint for i\n",
    "    return sum(model.a[i,j] * model.x[j, t] for j in model.J) >= model.b[i, t]\n",
    "\n",
    "model.AxbConstraint = Constraint(model.I, model.T, rule=ax_constraint_rule) # this creates one constraint for each member of the set model.I\n",
    "\n",
    "model.dual = Suffix(direction=Suffix.IMPORT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### <span style=\"color:blue\">Task</span>\n",
    "With pen and paper, determine the parameters a, b, and c to replicate the concrete model.\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By running the code, we create an abstract model. Now we need to create an instance of it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can create an instance without filling it with data\n",
    "instance = model.create_instance()\n",
    "instance.pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can load data from a file (written in AMPL format)\n",
    "data = DataPortal()\n",
    "data.load(filename='01_abstract.dat')\n",
    "# You can view the defined sets and parameters here\n",
    "list(data.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### <span style=\"color:blue\">Task</span>\n",
    "Compare your results for the parameters a, b, and c with the used values in _01_abstract.dat_.\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can create an instance that is filled with input data\n",
    "instance = model.create_instance(data)\n",
    "instance.pprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = SolverFactory('glpk')\n",
    "status = opt.solve(instance)\n",
    "status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "instance.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another way of reporting the results\n",
    "instance.solutions.store_to(status)\n",
    "status.write(filename='01_abstract_results.json', format='json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"Duals\")\n",
    "for c in instance.component_objects(Constraint, active=True):\n",
    "    print (\"   Constraint\",c)\n",
    "    for index in c:\n",
    "        try:\n",
    "            print (\"      \", index, instance.dual[c[index]])\n",
    "        except: # if there is an error, skip that constraint\n",
    "            continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### <span style=\"color:blue\">Task</span>\n",
    "1. Set the demand in the last time step to 140. What happens?\n",
    "2. How easy is it to add another power plant? Another time step?\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reporting into pandas DataFrame objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pandas is a package that allows you to organize your data in multidimensional \"tables\", so-called DataFrame objects. It is useful if you want to export your results into csv or Microsoft Excel format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the library first\n",
    "import pandas as pd\n",
    "\n",
    "# we will create a dictionary from the model instance variable x\n",
    "# the keys of the dictionary are the indices, the values are the x values\n",
    "supply_data = {(j, t): value(x) for (j, t), x in instance.x.items()}\n",
    "\n",
    "# create a DataFrame object from that dictionary\n",
    "df_supply = pd.DataFrame.from_dict(supply_data, orient=\"index\", columns=[\"SE [MWh]\"])\n",
    "df_supply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's make the index look better - through a multiindex\n",
    "df_supply.index = pd.MultiIndex.from_tuples(df_supply.index, names=('Technology', 't'))\n",
    "# Get rid of the letter t in the timesteps\n",
    "df_supply.index = df_supply.index.set_levels(df_supply.index.levels[1].str.replace('t', ''), level=1)\n",
    "# Show the DataFrame\n",
    "df_supply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This looks already good... but let's try to pivot the table, so that technologies appear as columns\n",
    "df_supply = df_supply.unstack(level=0)\n",
    "df_supply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The columns have two levels now, let's get rid of the first one\n",
    "df_supply = df_supply.droplevel(0, axis=1)\n",
    "df_supply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's add the unit to the names of the technologies\n",
    "df_supply.columns = [x + \" SE [MWh]\" for x in df_supply.columns]\n",
    "df_supply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We repeat this for the demand\n",
    "# Since I know that the first character of the timesteps is always t, I index from 1 in t[1:] to get rid of it\n",
    "demand_data = {t[1:]: value(x) for (i, t), x in instance.b.items() if i == \"Dem\"}\n",
    "\n",
    "# create a DataFrame object from that dictionary\n",
    "df_demand = pd.DataFrame.from_dict(demand_data, orient=\"index\", columns=[\"D [MWh]\"])\n",
    "df_demand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let' save the supply and demand into two different csv files\n",
    "# ! Adapt the decimal character and the column separator to your local settings !\n",
    "df_supply.to_csv(\"01_supply.csv\", sep=\";\", decimal=\",\")\n",
    "df_demand.to_csv(\"01_demand.csv\", sep=\";\", decimal=\",\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting with matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Most basic plot\n",
    "# pandas has in-built functions from matplotlib, so we can plot a DataFrame object\n",
    "%matplotlib inline\n",
    "df_supply.plot.area()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the library\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Create an empty figure\n",
    "fig = plt.figure()\n",
    "\n",
    "# Now let's make an interactive object that we can edit with a GUI\n",
    "%matplotlib\n",
    "plot_supply = df_supply.plot.area()\n",
    "# Add these options:\n",
    "# To change the colors: color=[\"darkgreen\", \"gray\"]\n",
    "# To add title: title = \"Example 1\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's add the demand to the same plot\n",
    "df_demand.plot.line(ax=plot_supply, color=\"k\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## <span style=\"color:red\">Homework</span>\n",
    "1. Expand the timesteps to cover one day in summer and one day in winter (total: 48h)\n",
    "2. Add a third technology, PV, with zero running costs and with varying upper bounds for every time step.\n",
    "3. Plot your results: <br>\n",
    "    a. separately for the summer day and the winter day (demand and supply curves, chronologically)<br>\n",
    "    b. together for the sorted load duration curve<br>\n",
    "    c. costs for each time step (same order as for the sorted load curve)\n",
    "***"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "urbs",
   "language": "python",
   "name": "urbs"
  },
  "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}