examples/notebooks/power-to-gas-boiler-chp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Power to Gas with Heat Coupling\n",
    "\n",
    "This is an example for power to gas with optional coupling to heat sector (via boiler OR Combined-Heat-and-Power (CHP))\n",
    "\n",
    "A location has an electric, gas and heat bus. The primary source is wind power, which can be converted to gas. The gas can be stored to convert into electricity or heat (with either a boiler or a CHP)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import pypsa\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "logging.basicConfig(level=\"INFO\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combined-Heat-and-Power (CHP) parameterisation\n",
    "\n",
    "This setup follows http://www.ea-energianalyse.dk/reports/student-reports/integration_of_50_percent_wind%20power.pdf pages 35-6 which follows http://www.sciencedirect.com/science/article/pii/030142159390282K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ratio between max heat output and max electric output\n",
    "nom_r = 1.0\n",
    "\n",
    "# backpressure limit\n",
    "c_m = 0.75\n",
    "\n",
    "# marginal loss for each additional generation of heat\n",
    "c_v = 0.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graph for the case that max heat output equals max electric output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(9, 5))\n",
    "\n",
    "t = 0.01\n",
    "ph = np.arange(0, 1.0001, t)\n",
    "\n",
    "ax.plot(ph, c_m * ph)\n",
    "ax.set_xlabel(\"P_heat_out\")\n",
    "ax.set_ylabel(\"P_elec_out\")\n",
    "ax.grid(True)\n",
    "\n",
    "ax.set_xlim([0, 1.1])\n",
    "ax.set_ylim([0, 1.1])\n",
    "ax.text(0.1, 0.7, \"Allowed output\", color=\"r\")\n",
    "ax.plot(ph, 1 - c_v * ph)\n",
    "\n",
    "for i in range(1, 10):\n",
    "    k = 0.1 * i\n",
    "    x = np.arange(0, k / (c_m + c_v), t)\n",
    "    ax.plot(x, k - c_v * x, color=\"g\", alpha=0.5)\n",
    "\n",
    "ax.text(0.05, 0.41, \"iso-fuel-lines\", color=\"g\", rotation=-7)\n",
    "ax.fill_between(ph, c_m * ph, 1 - c_v * ph, facecolor=\"r\", alpha=0.5)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network = pypsa.Network()\n",
    "network.set_snapshots(pd.date_range(\"2016-01-01 00:00\", \"2016-01-01 03:00\", freq=\"H\"))\n",
    "\n",
    "network.add(\"Bus\", \"0\", carrier=\"AC\")\n",
    "network.add(\"Bus\", \"0 gas\", carrier=\"gas\")\n",
    "\n",
    "network.add(\"Carrier\", \"wind\")\n",
    "network.add(\"Carrier\", \"gas\", co2_emissions=0.2)\n",
    "\n",
    "network.add(\"GlobalConstraint\", \"co2_limit\", sense=\"<=\", constant=0.0)\n",
    "\n",
    "network.add(\n",
    "    \"Generator\",\n",
    "    \"wind turbine\",\n",
    "    bus=\"0\",\n",
    "    carrier=\"wind\",\n",
    "    p_nom_extendable=True,\n",
    "    p_max_pu=[0.0, 0.2, 0.7, 0.4],\n",
    "    capital_cost=1000,\n",
    ")\n",
    "\n",
    "network.add(\"Load\", \"load\", bus=\"0\", p_set=5.0)\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"P2G\",\n",
    "    bus0=\"0\",\n",
    "    bus1=\"0 gas\",\n",
    "    efficiency=0.6,\n",
    "    capital_cost=1000,\n",
    "    p_nom_extendable=True,\n",
    ")\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"generator\",\n",
    "    bus0=\"0 gas\",\n",
    "    bus1=\"0\",\n",
    "    efficiency=0.468,\n",
    "    capital_cost=400,\n",
    "    p_nom_extendable=True,\n",
    ")\n",
    "\n",
    "network.add(\"Store\", \"gas depot\", bus=\"0 gas\", e_cyclic=True, e_nom_extendable=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add heat sector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.add(\"Bus\", \"0 heat\", carrier=\"heat\")\n",
    "\n",
    "network.add(\"Carrier\", \"heat\")\n",
    "\n",
    "network.add(\"Load\", \"heat load\", bus=\"0 heat\", p_set=10.0)\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"boiler\",\n",
    "    bus0=\"0 gas\",\n",
    "    bus1=\"0 heat\",\n",
    "    efficiency=0.9,\n",
    "    capital_cost=300,\n",
    "    p_nom_extendable=True,\n",
    ")\n",
    "\n",
    "network.add(\"Store\", \"water tank\", bus=\"0 heat\", e_cyclic=True, e_nom_extendable=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add CHP constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Guarantees ISO fuel lines, i.e. fuel consumption p_b0 + p_g0 = constant along p_g1 + c_v p_b1 = constant\n",
    "network.links.at[\"boiler\", \"efficiency\"] = (\n",
    "    network.links.at[\"generator\", \"efficiency\"] / c_v\n",
    ")\n",
    "\n",
    "model = network.optimize.create_model()\n",
    "\n",
    "link_p = model.variables[\"Link-p\"]\n",
    "link_p_nom = model.variables[\"Link-p_nom\"]\n",
    "\n",
    "# Guarantees heat output and electric output nominal powers are proportional\n",
    "model.add_constraints(\n",
    "    network.links.at[\"generator\", \"efficiency\"] * nom_r * link_p_nom[\"generator\"]\n",
    "    - network.links.at[\"boiler\", \"efficiency\"] * link_p_nom[\"boiler\"]\n",
    "    == 0,\n",
    "    name=\"heat-power output proportionality\",\n",
    ")\n",
    "\n",
    "# Guarantees c_m p_b1  \\leq p_g1\n",
    "model.add_constraints(\n",
    "    c_m * network.links.at[\"boiler\", \"efficiency\"] * link_p.sel(Link=\"boiler\")\n",
    "    - network.links.at[\"generator\", \"efficiency\"] * link_p.sel(Link=\"generator\")\n",
    "    <= 0,\n",
    "    name=\"backpressure\",\n",
    ")\n",
    "\n",
    "# Guarantees p_g1 +c_v p_b1 \\leq p_g1_nom\n",
    "model.add_constraints(\n",
    "    link_p.sel(Link=\"boiler\")\n",
    "    + link_p.sel(Link=\"generator\")\n",
    "    - link_p_nom.sel({\"Link-ext\": \"generator\"})\n",
    "    <= 0,\n",
    "    name=\"top_iso_fuel_line\",\n",
    ")\n",
    "\n",
    "network.optimize.solve_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.loads_t.p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.links.p_nom_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CHP is dimensioned by the heat demand met in three hours when no wind\n",
    "4 * 10.0 / 3 / network.links.at[\"boiler\", \"efficiency\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# elec is set by the heat demand\n",
    "28.490028 * 0.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.links_t.p0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.links_t.p1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame({attr: network.stores_t[attr][\"gas depot\"] for attr in [\"p\", \"e\"]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame({attr: network.stores_t[attr][\"water tank\"] for attr in [\"p\", \"e\"]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame({attr: network.links_t[attr][\"boiler\"] for attr in [\"p0\", \"p1\"]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.stores.loc[\"gas depot\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.generators.loc[\"wind turbine\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.links.p_nom_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the overall efficiency of the CHP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta_elec = network.links.at[\"generator\", \"efficiency\"]\n",
    "\n",
    "r = 1 / c_m\n",
    "\n",
    "# P_h = r*P_e\n",
    "(1 + r) / ((1 / eta_elec) * (1 + c_v * r))"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "",
   "language": "python",
   "name": ""
  },
  "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "9da88e3b8ddecaa5eb87bf19fcbb093fc705eb8087ff652e7fa90136336dde1c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}