oemof-jordan.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# oemof jordan\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pkg_resources as pkg\n",
    "import pandas as pd\n",
    "\n",
    "import pyomo.environ as po \n",
    "from pyomo.opt import SolverFactory\n",
    "from oemof.solph import EnergySystem, Model, Bus\n",
    "from oemof.tools.economics import annuity as annuity\n",
    "from oemof.solph import constraints, NonConvex\n",
    "import oemof.tabular.tools.postprocessing as pp\n",
    "import oemof.tabular.facades as fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating and Setting the Datapaths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, the datapath for raw-data and results is set. Data handling looks more complex than it is. You can easily adapt this to a simple `pd.read_excel(filepath,...)` in the next block if your file is located somewhere else. Otherwise we will use data from the repository repository. \n",
    "\n",
    "In addition a results directory will be created in `home/user/oemof-results/results/<scenario-name>/output`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario = \"AUT\" # name of the scenario and the corresponding file, e.g. \"A.xls\"\n",
    "\n",
    "# datapath for input data from the oemof tabular pacakge\n",
    "datapath = os.path.join(\n",
    "    os.getcwd(),#pkg.resource_filename(\"oemof.jordan\", \"\"),\n",
    "    \"scenarios\", \n",
    "    scenario + \".xls\",\n",
    ")\n",
    "\n",
    "carrier_technology_path = os.path.join(\n",
    "    os.getcwd(), \n",
    "    \"scenarios\",\n",
    "    \"carrier-technology.xls\"\n",
    ")\n",
    "\n",
    "# results path\n",
    "results_path = os.path.join(\n",
    "    os.path.expanduser(\"~\"), \"oemof-results\", \"jordan\"\n",
    ")\n",
    "\n",
    "scenario_path = os.path.join(\n",
    "    results_path, scenario, \"output\"\n",
    ")\n",
    "\n",
    "if not os.path.exists(scenario_path):\n",
    "    os.makedirs(scenario_path)\n",
    "print(scenario_path)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the required input data will be read. The profiles index will be used for the `EnergySystem` object below. \n",
    "All generator data etc. will also be loaded. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = pd.ExcelFile(datapath)\n",
    "\n",
    "sheet_names = [typ for typ in file.sheet_names if typ in fc.TYPEMAP.keys()]\n",
    "\n",
    "data = {}\n",
    "\n",
    "for sheet in sheet_names: \n",
    "    data[sheet]= pd.read_excel(datapath, sheet_name=sheet, index_col=0)\n",
    "\n",
    "# profiles and tech data not be part of TYPEMAP and need to be read seperately\n",
    "profiles = pd.read_excel(\n",
    "    datapath,\n",
    "    sheet_name=\"profiles\",\n",
    "    index_col=[0],\n",
    "    parse_dates=True,\n",
    ")\n",
    "profiles.index.freq = \"1H\"\n",
    "\n",
    "technology = pd.read_excel(\n",
    "    carrier_technology_path, sheet_name=\"technology-data\", index_col=[0, 1, 2]\n",
    ")\n",
    "carrier = pd.read_excel(\n",
    "    carrier_technology_path, sheet_name=\"carrier\", index_col=[0, 1]\n",
    ")\n",
    "\n",
    "if \"co2-limit\" in file.sheet_names:\n",
    "    co2_limit = pd.read_excel(datapath, sheet_name=\"co2-limit\", index_col=0)\n",
    "else: \n",
    "    co2_limit = None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_components = pd.concat([v for k, v in data.items() if k not in [\"bus\", \"co2-limit\"]], sort=False)\n",
    "# Only be used for Latex export of tables \n",
    "#columns = ['profile', 'capacity_potential']\n",
    "#print(all_components.to_latex(columns=columns, na_rep=\"-\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _capacity_cost(g, scenario=scenario):\n",
    "    if bool(g[\"expandable\"]):\n",
    "        val = annuity(\n",
    "            technology.at[(g[\"carrier\"], g[\"tech\"], scenario), \"capex\"], \n",
    "            technology.at[(g[\"carrier\"], g[\"tech\"], scenario), \"lifetime\"],\n",
    "            technology.at[(g[\"carrier\"], g[\"tech\"], scenario), \"wacc\"]) * (\n",
    "            (1 + (technology.at[(g[\"carrier\"], g[\"tech\"], scenario), \"fom\"])\n",
    "            )) \n",
    "        return val * 1000 # $/kw -> $/MW\n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "def _marginal_cost(g, scenario):\n",
    "    if not pd.isna(g[\"marginal_cost\"]):\n",
    "        return g[\"marginal_cost\"]\n",
    "    else:\n",
    "        return (carrier.at[(g[\"carrier\"], scenario), \"cost\"] / g[\"efficiency\"])\n",
    "\n",
    "def _none(number):\n",
    "    if pd.isna(number):\n",
    "        return None\n",
    "    else:\n",
    "        return number\n",
    "    \n",
    "def _write_investment_cost(scenario_path):\n",
    "    invest_e = {}\n",
    "    invest_p = {}\n",
    "    for n in es.nodes: \n",
    "        if hasattr(n, \"storage_capacity_cost\"):\n",
    "            if n.storage_capacity_cost is not None:\n",
    "                invest_e[n.label] = (n.storage_capacity_cost, m.results[n, None][\"scalars\"][\"invest\"])\n",
    "                invest_p[n.label] = (n.capacity_cost, m.results[es.groups[\"JO-electricity\"], n][\"scalars\"][\"invest\"])\n",
    "        elif hasattr(n, \"capacity_cost\"):\n",
    "            if n.capacity_cost is not None:\n",
    "                invest_p[n.label] = (n.capacity_cost, m.results[n, es.groups[\"JO-electricity\"]][\"scalars\"][\"invest\"])\n",
    "\n",
    "        pd.DataFrame(invest_e, index=[\"$/MWha\", \"MWh\"]).to_csv(\n",
    "            os.path.join(scenario_path, \"investment_energy.csv\"))\n",
    "        pd.DataFrame(invest_p, index=[\"$/MWa\", \"MW\"]).to_csv(\n",
    "            os.path.join(scenario_path, \"investment_power.csv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the EnergySystem and its Nodes\n",
    "\n",
    "Firs, an `EnergySystem` object will be created. This holds all information (nodes, etc.) of the modelled energy system that will be added below. This is just the standard way of using the `oemof.solph` library for your modelling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "es = EnergySystem(timeindex=profiles.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Buses\n",
    "\n",
    "Before any components are added, all bus objects are added tothe energy system object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buses = {\n",
    "    name: Bus(label=name, balanced=bool(arg.balanced))\n",
    "    for name, arg in data[\"bus\"].iterrows()\n",
    "}\n",
    "es.add(*buses.values())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bus Constraints \n",
    "\n",
    "With the set of all Buses $B$ all inputs $x^{flow}_{i(b),b}$ to a bus $b$ must equal all its outputs $x^{flow}_{b,o(b)}$\n",
    "\n",
    "$$\\sum_i x^{flow}_{i(b), b}(t) - \\sum_o x^{flow}_{b, o(b)}(t) = 0 \\qquad \\forall t \\in T, \\forall b \\in B$$\n",
    "\n",
    "This equation will be build once the complete energy system is setup with its components. Every time a `Component` is created, the connected bus inputs/outputs will be updated. By this update every bus has all required information of its inputs and outputs available to construct the constraints. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load Constraint \n",
    "\n",
    "For the set of all Load denoted with $l \\in L$ the load $x_l$ at timestep t equals the exogenously defined  profile value $c^{profile}_l$ multiplied by the amount of this load $c^{amount}_l$\n",
    "\n",
    "$$ x^{flow}_{l}(t) = c^{profile}_{l}(t) \\cdot c^{amount}_{l} \\qquad \\forall t \\in T, \\forall l \\in L$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, l in data[\"load\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Load(\n",
    "            label=name,\n",
    "            bus=buses[\n",
    "                l.bus\n",
    "            ],  # reference the bus in the buses dictionary\n",
    "            amount=l.amount,  # amount column\n",
    "            profile=profiles[l.profile],\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dispatchable and Volatile component\n",
    "\n",
    "#### Dispatchable Generator Constraint\n",
    "\n",
    "A `Dispatachble` component can be used to model all types of dispatchble units in a energy system. This can include diesel generators oder coal fired power plants but also hot water boilers for heat. Every generator **must** be connected to an `Bus` object. \n",
    "\n",
    "This basic mathematical model for the component with the set of all dispatchable generators being $d \\in D$ looks as follows:\n",
    "\n",
    "$$x^{flow}_{d}(t) \\leq c^{capacity}_{d} \\qquad \\forall t \\in T,  \\forall d \\in D$$\n",
    "\n",
    "Meaning, the production of the generator $x^{flow}$ must be less than its maximum capacity $c^{capacity}_d$ in every timestep. *Note that this equation holds for the case where the `expandable` attribute is set to `False`*. For the investment case the following two equation hold:\n",
    "\n",
    "$$x^{flow}_{d}(t) \\leq x^{capacity}_{d} \\qquad \\forall t \\in T,  \\forall d \\in D$$\n",
    "$$c^{capacity}_d + c^{capacity\\_minium}_d \\leq x^{capacity}_{d} \\leq c^{capacity\\_potential}_{d} \\forall d \\in D$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, g in data[\"dispatchable\"].iterrows():\n",
    "    g = g.to_dict()\n",
    "    \n",
    "    if g.get(\"p_min\"):\n",
    "        milp_constraints = NonConvex()\n",
    "    else:\n",
    "        milp_constraints = None\n",
    "    es.add(\n",
    "        fc.Dispatchable(\n",
    "            label=name,\n",
    "            bus=buses[g[\"bus\"]],\n",
    "            carrier=g[\"carrier\"],\n",
    "            tech=g[\"tech\"],\n",
    "            marginal_cost=_marginal_cost(g, scenario),\n",
    "            expandable=g.get(\"expandable\", False),\n",
    "            capacity=g[\"capacity\"],\n",
    "            capacity_potential=None,\n",
    "            capacity_cost=_capacity_cost(g, scenario),\n",
    "            capacity_minimum=g.get(\"capacity_minimum\", 0),\n",
    "            output_parameters={\n",
    "                \"emission_factor\": (\n",
    "                    carrier.at[(g[\"carrier\"], scenario), \"emission_factor\"]\n",
    "                    / g[\"efficiency\"]\n",
    "                ),\n",
    "                \"summed_max\": g.get(\"max_flh\", 8760),\n",
    "                \"summed_min\": g.get(\"min_flh\", 0),\n",
    "                \"nonconvex\": milp_constraints,\n",
    "                \"min\": g.get(\"p_min\", 0)\n",
    "            },\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Volatile Generator Constraint\n",
    "\n",
    "For the `Volatile` $v \\in V$ generators that can be used to model PV and Wind units the flow will be fixed by a given profile multiplied with the (optimised) capacity $x^{capacity}_{v}$ in the investment case and the exougenously defined capacity $c^{capacity}_{v}$ in the non-investment case (i.e. `expandable` attribute `False`). \n",
    "\n",
    "$$ x^{flow}_{v}(t) = c^{profile}_{v}(t) \\cdot x^{capacity}_{v} \\qquad \\forall t \\in T, \\forall v \\in V$$\n",
    "\n",
    "The optimised capacity is limited by the following variable bounds:\n",
    "\n",
    "$$c^{capacity}_v + c^{capacity\\_minimum}_v \\leq x^{capacity}_{v}  \\leq c^{capacity\\_potential}_{v} \\qquad \\forall v \\in V$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, v in data[\"volatile\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Volatile(\n",
    "            label=name,\n",
    "            bus=buses[v.bus],\n",
    "            carrier=v.carrier,\n",
    "            tech=v.tech,\n",
    "            expandable=v.expandable,\n",
    "            capacity=v.capacity,\n",
    "            capacity_potential=v.capacity_potential,\n",
    "            capacity_cost=_capacity_cost(v, scenario),\n",
    "            capacity_minimum=v.capacity_minimum,\n",
    "            profile=profiles[v.profile],\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Storage components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, s in data[\"storage\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Storage(\n",
    "            label=name,\n",
    "            bus=buses[s.bus],\n",
    "            carrier=s.carrier,\n",
    "            tech=s.tech,\n",
    "            marginal_cost=s.marginal_cost,\n",
    "            capacity=s.capacity,\n",
    "            storage_capacity=s.storage_capacity,\n",
    "            storage_capacity_potential=s.storage_capacity_potential,\n",
    "            min_storage_level=s.min_storage_level,\n",
    "            invest_relation_output_capacity=_none(s.invest_relation_output_capacity),\n",
    "            invest_relation_input_capacity=_none(s.invest_relation_input_capacity),\n",
    "            expandable=s.expandable,\n",
    "            efficiency=s.efficiency,\n",
    "            loss_rate=s.loss_rate,\n",
    "            storage_capacity_cost=annuity(\n",
    "                technology.at[\n",
    "                    (s.carrier, s.tech, scenario), \"storage_capex\"\n",
    "                ],\n",
    "                technology.at[\n",
    "                    (s.carrier, s.tech, scenario), \"lifetime\"\n",
    "                ],\n",
    "                 technology.at[\n",
    "                    (s.carrier, s.tech, scenario), \"wacc\"\n",
    "                ],\n",
    "            )\n",
    "            * 1000,  # $/kW -> $/MW\n",
    "            capacity_cost=_capacity_cost(s, scenario)\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Storage Constraints \n",
    "\n",
    "The mathematical representation of the storage for all storages $s \\in S$ will include the flow into the storage, out of the storage and a storage level. The default efficiency for input/output is 1. Note that the provided efficiency is applied for charge and discharge. Hence, if you want to set the round trip efficiency you need to set $\\eta = \\sqrt{\\eta^{roundtrip}}$, where $\\eta$ is the efficiency provided to construct the component. \n",
    "\n",
    "Intertemporal energy balance of the storage:\n",
    "\n",
    "$$ x^{level}_{s}(t) = (1-\\eta^{loss\\_rate}) x^{level}_{s}(t) + \\eta_{in} x^{flow}_{s, in} - \\frac{x^{flow}_{s, out}(t)}{\\eta_{out}} \\qquad \\forall t \\in T,  \\forall s \\in S$$ \n",
    "\n",
    "Bounds of the storage level variable $x^{level}_s(t)$ with investment. The first and the last timestep are linked and set to a fixed value (e.g. 0.5 of the maximum storage capacity)\n",
    "\n",
    "$$ x^{level}_s(1) = x_s^{level}(t_{e}) = 0.5 \\cdot \\overline{x}_s^{level} \\qquad \\forall t \\in T,  \\forall s \\in S$$ \n",
    "\n",
    "Additionally, the level is bounded by the optimised storage capacity (energy) $\\overline{x}_s^{level}$:\n",
    "\n",
    "$$ x^{level}_s(t) \\leq \\overline{x}_s^{level} \\qquad \\forall t \\in T, \\forall s \\in S$$ \n",
    "\n",
    "The investment itself is limited by parameters of the existing storage capacity and its potential: \n",
    "\n",
    "$$ c_s^{storage\\_capacity} \\leq \\overline{x}_s^{level} \\leq c_s^{storage\\_capacity\\_potential} $$\n",
    "\n",
    "The inflow/outflow (i.e. power) of the storage also needs to be within the limit of the minimum and maximum power. \n",
    "\n",
    "$$ -x_s^{capacity} \\leq x^{flow}_s(t) \\leq x_s^{capacity} \\qquad \\forall t \\in T, \\forall s \\in S$$ \n",
    "\n",
    "The investment in capacity is again bounded. Similary to the bounds of the storage capacity it is limited by the existing capacity and the capacity potential.\n",
    "\n",
    "$$ c_s^{capacity} \\leq x_s^{capacity} \\leq c_s^{capacity\\_potential} $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conversion component\n",
    "\n",
    "A conversion unit will take from a bus and feed into another: \n",
    "\n",
    "$$x^{flow}_{c, to}(t) = c^{efficiencty}_{c} \\cdot x^{flow}_{c, from}(t), \\qquad \\forall c  \\in C, \\forall t \\in T$$ \n",
    "\n",
    "For the non-investment case the outflow, i.e. the flow **to** the bus where the conversion unit is connected with its output the following constraint applies:\n",
    "\n",
    "$$x^{flow}_{c, to} \\leq c^{capacity}_c$$\n",
    "\n",
    "In the case of investment compare constraints from the volatile generator units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, c in data[\"conversion\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Conversion(\n",
    "            label=name,\n",
    "            from_bus=buses[c.from_bus],\n",
    "            to_bus=buses[c.to_bus],\n",
    "            carrier=c.carrier,\n",
    "            tech=c.tech,\n",
    "            efficiency=c.efficiency,\n",
    "            marginal_cost=c.marginal_cost,\n",
    "            carrier_cost=c.carrier_cost,\n",
    "            expandable=c.expandable,\n",
    "            capacity=c.capacity,\n",
    "            capacity_minimum=c.capacity_minimum,\n",
    "            capacity_potential=c.capacity_potential,\n",
    "            capacity_cost=_capacity_cost(c, scenario),\n",
    "            output_parameters={\n",
    "                \"emission_factor\": (\n",
    "                    carrier.at[(c.carrier, scenario), \"emission_factor\"]\n",
    "                    / c.efficiency)\n",
    "            },\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Commodity components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the commodity components, the aggregated flow for the complete time horizon is limited by the user defined amount:\n",
    "$$ \\sum_t x^{flow}_k(t) \\leq c^{amount} \\qquad \\forall k \\in K$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, c in data[\"commodity\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Commodity(\n",
    "            label=name,\n",
    "            bus=buses[c.bus],\n",
    "            carrier=c.carrier,\n",
    "            tech=c.tech,\n",
    "            amount=c.amount,\n",
    "            output_parameters={\n",
    "                \"summed_min\": c.minimum_amount\n",
    "            }\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Link components\n",
    "\n",
    "The link components allows to model energy transfer from one bus to an other. For example transhipment modelling approach can be used for transmission in electricity systems. The equation for a line $n$ is considering a loss within the transfer process $from$ a bus $to$ the other bus.\n",
    "\n",
    "$$x^{flow}_{from, n}(t) = c^{loss}_{n} \\cdot x^{flow}_{n, to}(t), \\qquad \\forall n  \\in N, \\forall t \\in T$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, c in data[\"link\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Link(\n",
    "            label=name,\n",
    "            from_bus=buses[c.from_bus],\n",
    "            to_bus=buses[c.to_bus],\n",
    "            capacity=c.capacity,\n",
    "            expandable=c.expandable,\n",
    "             capacity_cost=annuity(\n",
    "                technology.at[(c.carrier, c.tech, scenario), \"capex\"],\n",
    "                technology.at[\n",
    "                    (c.carrier, c.tech, scenario), \"lifetime\"\n",
    "                ],\n",
    "                technology.at[\n",
    "                    (c.carrier, c.tech, scenario), \"wacc\"\n",
    "                ],\n",
    "            )\n",
    "            * 1000,  # $/kW -> $/MW\n",
    "            loss=c.loss\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective Function \n",
    "\n",
    "The objective function is created from all instantiated objects. It will use all operating costs (i.e. `marginal_cost` argument) and if set all investment costs (i.e. `capacity_cost` argument) for all generators (dispatchable and volatile, conversion units and storages). In addition the energy related storage costs are added.  \n",
    "\n",
    "$$ \\text{min:} \\sum_g \\sum_t \\overbrace{c^{marginal\\_cost}_g \\cdot x^{flow}_{g}(t)}^{\\text{operating cost}} \\\\ \n",
    "\\sum_g \\sum_t \\overbrace{c^{capacity\\_cost}_g \\cdot x^{capacity}_{g}(t)}^{\\text{investment cost}} + \\\\\n",
    "\\sum_s \\overbrace{c^{storage\\_capacity\\_cost}_s \\cdot \\overline{x}^{storage\\_capacity}_{s}(t)}^{\\text{storage energy cost}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add Shortage/Excess Slack Components\n",
    "\n",
    "If these are added, additional cost will occur inside the objective function. Shortage / Excess can be understood as a very expensive generator or load repectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, e in data[\"excess\"].iterrows():\n",
    "    es.add(fc.Excess(label=name, bus=buses[e.bus]))\n",
    "\n",
    "for name, s in data[\"shortage\"].iterrows():\n",
    "    es.add(\n",
    "        fc.Shortage(\n",
    "            label=name,\n",
    "            carrier=\"electricity\",\n",
    "            tech=\"shortage\",\n",
    "            bus=buses[s.bus],\n",
    "            marginal_cost=s.marginal_cost,\n",
    "        )\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the Mathematical Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create model based on energy system and its components\n",
    "m = Model(es)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add CO2 Constraint / RE share\n",
    "\n",
    "To add a CO2-constraint we will use the `oemof.solph.constraints` module which allows to add such a constraint in a easy way. \n",
    "\n",
    "$$ \\sum_t \\sum_f  x^{flow}_f(t) \\cdot c^{emission\\_factor}_f \\leq \\overline{L_{CO_2}} $$\n",
    "\n",
    "The constraint will sum all flows for the complete time horzion that have an attribute `emission_factor` and multiple the flow value with this factor. \n",
    "\n",
    "#### RE Share\n",
    "\n",
    "The renewable share $RE^{share}$ is implemented with an indiviual constraint by limiting the total production of conventional units $c \\in C$:\n",
    "\n",
    "$$ \\sum_{t \\in T} \\sum_{c \\in C}  x^{flow}_c(t) \\leq (1-RE^{share}) \\cdot  c^{amount}_{l} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if co2_limit is not None:\n",
    "    constraints.generic_integral_limit(m, keyword='emission_factor', limit=co2_limit.loc[\"JO-electricity\", \"value\"])\n",
    "#m.write(io_options={\"symbolic_solver_labels\":True})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#m.receive_duals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dispatchable = [(i,o) for i,o in m.flows if isinstance(i, fc.Dispatchable)]\n",
    "conversion = [(i,o) for i,o in m.flows if isinstance(o, fc.Conversion)]\n",
    "demand = [(i,o) for i,o in m.flows if isinstance(o, fc.Load)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the Model and Writing Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if cbc solver library is available\n",
    "cbc = SolverFactory('cbc').available()\n",
    "\n",
    "if cbc:\n",
    "    #  solve the model using cbc solver\n",
    "    m.solve(\"gurobi\", cmdline_options = {\"mipgap\":0.02})\n",
    "\n",
    "    # write results back to the model object\n",
    "    m.results = m.results()\n",
    "\n",
    "    # writing results with the standard oemof-tabular output formatt\n",
    "    pp.write_results(m, m.results, scenario_path)\n",
    "\n",
    "    _write_investment_cost(scenario_path)\n",
    "    \n",
    "    pd.Series(\n",
    "        [(m.integral_limit_emission_factor_constraint() / 1e6), \n",
    "         #m.dual[m.integral_limit_emission_factor_constraint], \n",
    "         m.objective()],\n",
    "        index=[\"CO2 (Mio. t)\", \n",
    "               #\"Shadow Price in $/t\", \n",
    "               \"Objective value\"]\n",
    "    ).to_csv(os.path.join(scenario_path, \"costs.csv\"))\n",
    "    \n",
    "    \n",
    "    \n",
    "    print(\n",
    "        \"Optimization done. Results are in {}.\".format(\n",
    "            results_path\n",
    "        )\n",
    "    )\n",
    "    \n",
    "    # write the lp-file\n",
    "    # m.write(io_options={'symbolic_solver_labels': True})\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"Emissions in Mio t: \" + str(m.integral_limit_emission_factor_constraint() / 1e6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(\"Shadow price of CO2: \" + str(m.dual[m.integral_limit_emission_factor_constraint]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if True:\n",
    "    # base scenario co2 emissions without co2 constraint as co2_base \n",
    "    co2_base = m.integral_limit_emission_factor_constraint() # take as base for emissions\n",
    "    sensitivities = [.9, .8, .7, .6, .5, .4, .3, .2, .1]\n",
    "    \n",
    "    for co2_reduction in sensitivities:\n",
    "        #print(co2_reduction * co2_base)\n",
    "        sensitivity_path = scenario_path.replace(scenario, scenario + \"-\" + str(int(co2_reduction * 100)))\n",
    "        \n",
    "        if not os.path.exists(sensitivity_path):\n",
    "            os.makedirs(sensitivity_path)\n",
    "        \n",
    "        m = Model(es)\n",
    "        \n",
    "        def _re_share(m):\n",
    "            lhs = sum(m.flow[i,o,t] for i,o in conversion for t in m.TIMESTEPS)\n",
    "            lhs += sum(m.flow[i, o, t] for i,o in dispatchable for t in m.TIMESTEPS)\n",
    "            rhs = sum(m.flow[i, o, t] for i,o in demand for t in m.TIMESTEPS)\n",
    "            return (lhs <= (1-co2_reduction) * rhs)\n",
    "        m.renewable_share = po.Constraint(rule=_re_share)  \n",
    "        \n",
    "        constraints.generic_integral_limit(\n",
    "            m, keyword='emission_factor', \n",
    "            limit=co2_limit.loc[\"JO-electricity\", \"value\"]) #co2_reduction * co2_base)\n",
    "        \n",
    "        #m.receive_duals()\n",
    "\n",
    "        m.solve(\"gurobi\", cmdline_options = {\"mipgap\":0.02})\n",
    "\n",
    "        m.results = m.results()\n",
    "\n",
    "        pp.write_results(m, m.results, sensitivity_path)\n",
    "        \n",
    "        _write_investment_cost(sensitivity_path)\n",
    "    \n",
    "        pd.Series(\n",
    "            [(m.integral_limit_emission_factor_constraint() / 1e6), \n",
    "             #m.dual[m.integral_limit_emission_factor_constraint], \n",
    "             m.objective()],\n",
    "            index=[\"CO2 (Mio. t)\", \n",
    "                   #\"Shadow Price in $/t\", \n",
    "                   \"Objective value\"]\n",
    "\n",
    "        ).to_csv(os.path.join(sensitivity_path, \"costs.csv\"))"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "jordan-env",
   "language": "python",
   "name": "jordan-env"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}