{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 7.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Integrated Energy Grids**\n",
    "\n",
    "\n",
    "**Problem 7.4**\n",
    "\n",
    "**In this problem, we will model a node that contains an electricity, a gas and a heat bus. There is a demand of 50 MWh of electricity and 40 MWh of heat. Electricity can be produced using an Open Cycle Gas Turbine (OCGT) with an efficiency of 0.35 or using a Combined Heat and Power (CHP) unit. Heat can be produced with the CHP unit or with a gas boiler. The gas boiler has an efficiency of 0.9. The OCGT, gas boiler and CHP unit have a marginal cost due to the fuel cost of 20 EUR/MWh$_{th}$**\n",
    "\n",
    "**(a) In this section, we assume that the CHP has a fixed efficiency of 0.3 when producing electricity and 0.3 when producing heat. Model the OCGT and gas boiler, using a link and the CHP unit using and multilink element in PyPSA. Add a gas store to the gas bus that represents an unlimited supply of gas. Optimize the system and calculate which technologies are supplying the electricity and heat demand.**\n",
    "\n",
    "**(b) In this section, we assume that the CHP unit can be operated either in condensing mode or in back-pressure mode. In practice, this means the feasible space for operating the CHP unit is determined by the iso-fuel lines (with constant $c_v$=-0.15) and the back-pressure line (with constant $c_m$=0.75). Optimize the system and calculate which technologies are supplying the electricity and heat demand.**\n",
    "\n",
    "_Note: This problem is based on the PyPSA examples [Multi-Link usage: CHP with fixed heat-power ratio](https://pypsa.readthedocs.io/en/stable/examples/chp-fixed-heat-power-ratio.html) and [Power to Gas with Heat Coupling](https://pypsa.readthedocs.io/en/stable/examples/power-to-gas-boiler-chp.html)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "If you have not yet set up Python on your computer, you can execute this tutorial in your browser via [Google Colab](https://colab.research.google.com/). Click on the rocket in the top right corner and launch \"Colab\". If that doesn't work download the `.ipynb` file and import it in [Google Colab](https://colab.research.google.com/).\n",
    "\n",
    "Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.\n",
    "\n",
    "```sh\n",
    "!pip install numpy pypsa\n",
    "```\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by creating the network and adding the buses and loads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pypsa\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['gas'], dtype='object')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network = pypsa.Network()\n",
    "\n",
    "network.add(\"Bus\", \"electricity\")\n",
    "network.add(\"Load\", \"electricity load\", bus=\"electricity\", p_set=50)\n",
    "\n",
    "network.add(\"Bus\", \"heat\")\n",
    "network.add(\"Load\", \"heat load\", bus=\"heat\", p_set=40)\n",
    "\n",
    "network.add(\"Bus\", \"gas\")\n",
    "\n",
    "# We add a gas store with energy capacity and an initial filling level much higher than the required gas consumption, \n",
    "# this way gas supply is unlimited\n",
    "network.add(\"Store\", \"gas\", e_initial=1e6, e_nom=1e6, bus=\"gas\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>v_nom</th>\n",
       "      <th>type</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>carrier</th>\n",
       "      <th>unit</th>\n",
       "      <th>v_mag_pu_set</th>\n",
       "      <th>v_mag_pu_min</th>\n",
       "      <th>v_mag_pu_max</th>\n",
       "      <th>control</th>\n",
       "      <th>generator</th>\n",
       "      <th>sub_network</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bus</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>electricity</th>\n",
       "      <td>1.0</td>\n",
       "      <td></td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>AC</td>\n",
       "      <td></td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>inf</td>\n",
       "      <td>PQ</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>heat</th>\n",
       "      <td>1.0</td>\n",
       "      <td></td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>AC</td>\n",
       "      <td></td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>inf</td>\n",
       "      <td>PQ</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gas</th>\n",
       "      <td>1.0</td>\n",
       "      <td></td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>AC</td>\n",
       "      <td></td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>inf</td>\n",
       "      <td>PQ</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             v_nom type    x    y carrier unit  v_mag_pu_set  v_mag_pu_min  \\\n",
       "Bus                                                                          \n",
       "electricity    1.0       0.0  0.0      AC                1.0           0.0   \n",
       "heat           1.0       0.0  0.0      AC                1.0           0.0   \n",
       "gas            1.0       0.0  0.0      AC                1.0           0.0   \n",
       "\n",
       "             v_mag_pu_max control generator sub_network  \n",
       "Bus                                                      \n",
       "electricity           inf      PQ                        \n",
       "heat                  inf      PQ                        \n",
       "gas                   inf      PQ                        "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.buses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We add the links representing the OCGT and the CHP. For the later, we can use a multilink, that is a link that outputs to two buses (electricity and heating)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>bus0</th>\n",
       "      <th>bus1</th>\n",
       "      <th>type</th>\n",
       "      <th>carrier</th>\n",
       "      <th>efficiency</th>\n",
       "      <th>active</th>\n",
       "      <th>build_year</th>\n",
       "      <th>lifetime</th>\n",
       "      <th>p_nom</th>\n",
       "      <th>p_nom_mod</th>\n",
       "      <th>...</th>\n",
       "      <th>min_down_time</th>\n",
       "      <th>up_time_before</th>\n",
       "      <th>down_time_before</th>\n",
       "      <th>ramp_limit_up</th>\n",
       "      <th>ramp_limit_down</th>\n",
       "      <th>ramp_limit_start_up</th>\n",
       "      <th>ramp_limit_shut_down</th>\n",
       "      <th>p_nom_opt</th>\n",
       "      <th>bus2</th>\n",
       "      <th>efficiency2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Link</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>OCGT</th>\n",
       "      <td>gas</td>\n",
       "      <td>electricity</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>0.35</td>\n",
       "      <td>True</td>\n",
       "      <td>0</td>\n",
       "      <td>inf</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td></td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gas boiler</th>\n",
       "      <td>gas</td>\n",
       "      <td>heat</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>0.90</td>\n",
       "      <td>True</td>\n",
       "      <td>0</td>\n",
       "      <td>inf</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td></td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CHP</th>\n",
       "      <td>gas</td>\n",
       "      <td>electricity</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>0.30</td>\n",
       "      <td>True</td>\n",
       "      <td>0</td>\n",
       "      <td>inf</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>heat</td>\n",
       "      <td>0.3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 36 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           bus0         bus1 type carrier  efficiency  active  build_year  \\\n",
       "Link                                                                        \n",
       "OCGT        gas  electricity                     0.35    True           0   \n",
       "gas boiler  gas         heat                     0.90    True           0   \n",
       "CHP         gas  electricity                     0.30    True           0   \n",
       "\n",
       "            lifetime   p_nom  p_nom_mod  ...  min_down_time  up_time_before  \\\n",
       "Link                                     ...                                  \n",
       "OCGT             inf  1000.0        0.0  ...              0               1   \n",
       "gas boiler       inf  1000.0        0.0  ...              0               1   \n",
       "CHP              inf  1000.0        0.0  ...              0               1   \n",
       "\n",
       "            down_time_before  ramp_limit_up  ramp_limit_down  \\\n",
       "Link                                                           \n",
       "OCGT                       0            NaN              NaN   \n",
       "gas boiler                 0            NaN              NaN   \n",
       "CHP                        0            NaN              NaN   \n",
       "\n",
       "            ramp_limit_start_up  ramp_limit_shut_down  p_nom_opt  bus2  \\\n",
       "Link                                                                     \n",
       "OCGT                        1.0                   1.0        0.0         \n",
       "gas boiler                  1.0                   1.0        0.0         \n",
       "CHP                         1.0                   1.0        0.0  heat   \n",
       "\n",
       "            efficiency2  \n",
       "Link                     \n",
       "OCGT                1.0  \n",
       "gas boiler          1.0  \n",
       "CHP                 0.3  \n",
       "\n",
       "[3 rows x 36 columns]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.add(\n",
    "    \"Link\",\n",
    "    \"OCGT\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"electricity\",\n",
    "    p_nom=1000,\n",
    "    marginal_cost=20,\n",
    "    efficiency=0.35,\n",
    ")\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"gas boiler\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"heat\",\n",
    "    marginal_cost=20,\n",
    "    p_nom=1000,\n",
    "    efficiency=0.9,\n",
    ")\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"CHP\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"electricity\",\n",
    "    bus2=\"heat\",\n",
    "    p_nom=1000,\n",
    "    marginal_cost=20,\n",
    "    efficiency=0.3,\n",
    "    efficiency2=0.3,\n",
    ")\n",
    "\n",
    "network.links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pypsa.consistency:The following stores have carriers which are not defined:\n",
      "Index(['gas'], dtype='object', name='Store')\n",
      "WARNING:pypsa.consistency:The following buses have carriers which are not defined:\n",
      "Index(['electricity', 'heat', 'gas'], dtype='object', name='Bus')\n",
      "WARNING:pypsa.consistency:The following links have carriers which are not defined:\n",
      "Index(['OCGT', 'gas boiler', 'CHP'], dtype='object', name='Link')\n",
      "WARNING:pypsa.consistency:The following stores have carriers which are not defined:\n",
      "Index(['gas'], dtype='object', name='Store')\n",
      "WARNING:pypsa.consistency:The following buses have carriers which are not defined:\n",
      "Index(['electricity', 'heat', 'gas'], dtype='object', name='Bus')\n",
      "WARNING:pypsa.consistency:The following links have carriers which are not defined:\n",
      "Index(['OCGT', 'gas boiler', 'CHP'], dtype='object', name='Link')\n",
      "INFO:linopy.model: Solve problem using Highs solver\n",
      "INFO:linopy.io: Writing time: 0.03s\n",
      "INFO:linopy.constants: Optimization successful: \n",
      "Status: ok\n",
      "Termination condition: optimal\n",
      "Solution: 5 primals, 12 duals\n",
      "Objective: 3.24e+03\n",
      "Solver model: available\n",
      "Solver message: optimal\n",
      "\n",
      "INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Link-fix-p-lower, Link-fix-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-energy_balance were not assigned to the network.\n"
     ]
    }
   ],
   "source": [
    "network.optimize();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Load</th>\n",
       "      <th>electricity load</th>\n",
       "      <th>heat load</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>snapshot</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>now</th>\n",
       "      <td>50.0</td>\n",
       "      <td>40.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Load      electricity load  heat load\n",
       "snapshot                             \n",
       "now                   50.0       40.0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.loads_t.p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Link</th>\n",
       "      <th>OCGT</th>\n",
       "      <th>gas boiler</th>\n",
       "      <th>CHP</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>snapshot</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>now</th>\n",
       "      <td>28.571429</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>133.333333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Link           OCGT  gas boiler         CHP\n",
       "snapshot                                   \n",
       "now       28.571429        -0.0  133.333333"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.links_t.p0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Link</th>\n",
       "      <th>OCGT</th>\n",
       "      <th>gas boiler</th>\n",
       "      <th>CHP</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>snapshot</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>now</th>\n",
       "      <td>-10.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-40.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Link      OCGT  gas boiler   CHP\n",
       "snapshot                        \n",
       "now      -10.0         0.0 -40.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.links_t.p1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Link</th>\n",
       "      <th>OCGT</th>\n",
       "      <th>gas boiler</th>\n",
       "      <th>CHP</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>snapshot</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>now</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-40.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Link      OCGT  gas boiler   CHP\n",
       "snapshot                        \n",
       "now        0.0         0.0 -40.0"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.links_t.p2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(b) In this section, we assume that the CHP unit can be operated either in condensing mode or in back-pressure mode. In practice, this means the feasible space for operating the CHP unit is determined by the iso-fuel lines (with constant $c_v$=-0.15) and the back-pressure line (with constant $c_m$=0.75). Optimize the system and calculate which technologies are supplying the electricity and heat demand.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combined-Heat-and-Power (CHP) parameterisation\n",
    "\n",
    "The parameters describing the operation of the CPH unit follow https://doi.org/10.1016/0301-4215(93)90282-K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# backpressure limit\n",
    "c_m = 0.75\n",
    "\n",
    "# marginal loss of electricity generation for each additional generation of heat\n",
    "c_v = 0.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot a figure that shows the operating space of the CHP units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAHpCAYAAAA/CfW/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAmVFJREFUeJzs3Xd8XPWd7//XmapRLzOSbXnkgjFgXLCtMS2EFhMwwYRmCbKBLLAbwi6E9W9zF7YlJGRJJc7eLC2XheQmF4sSWjABBwKhL6GF0IvBYNybbFmacs75/TEjedTPGY3a6P18PM4D64xm5siZePSec76ft2Hbto2IiIiIiIgUBM9oH4CIiIiIiIjkj0KeiIiIiIhIAVHIExERERERKSAKeSIiIiIiIgVEIU9ERERERKSAKOSJiIiIiIgUEIU8ERERERGRAuIb7QMYaZZl8emnn1JWVoZhGKN9OCIiIiIiIo7Yts2ePXuYMmUKHk//5+smXMj79NNPiUajo30YIiIiIiIiOfn444+ZOnVqv7dPuJBXVlYGwLp166iurh7lo5FCkEwmeeSRRzjppJPw+/2jfThSAPSaknzTa0ryTa8pyTe9ppxpbW0lGo12ZZr+TLiQ13mJZllZGeXl5aN8NFIIkskkxcXFlJeX6x8lyQu9piTf9JqSfNNrSvJNryl3Blt2psErIiIiIiIiBUQhT0REREREpIAo5ImIiIiIiBQQhTwREREREZECopAnIiIiIiJSQBTyRERERERECohCnoiIiIiISAFRyBMRERERESkgCnkiIiIiIiIFRCFPRERERESkgCjkiYiIiIiIFBCFPBERERERkQKikCciIiIiIlJAFPJEREREREQKiEKeiIiIiIhIAVHIExERERERKSAKeSIiIiIiIgVEIU9ERERERKSAKOSJiIiIiIgUEIU8ERERERGRAqKQJyIiIiIiUkAU8kRERERERAqIQp6IiIiIiEgBUcgTEREREREpIAp5IiIiIiIiBUQhT0REREREpIAo5ImIiIiIiBQQhTwREREREZECopAnIiIiIiJSQBTyRERERERECohCnoiIiIiISAFRyBMRERERESkgoxry/vjHP3LaaacxZcoUDMPg3nvvHfQ+TzzxBIsXL6aoqIiZM2dy4403Dv+BioiIiIiIjBOjGvLa2tpYsGABP/vZzxx9/7p161i2bBnHHHMML7/8Mv/8z//M5Zdfzt133z3MRyoiIiIiIjI++EbzyU855RROOeUUx99/44030tDQwKpVqwA45JBD+NOf/sSPfvQjzjrrLHdPvu0dCB4AoSrw+t3dV0REREREZIwa1ZDn1rPPPstJJ53Ubd/nP/95brnlFpLJJH5/77AWj8eJx+NdX7e2tgLgv+3zEDQAsIsqIFSNXVwDoWoorsEuroZQ53+77yNUCYaWM0paMpns9l+RodJrSvJNrynJN72mJN/0mnLG6d/PuAp5mzZtoq6urtu+uro6UqkU27ZtY/Lkyb3uc+2113L11Vf32h/3lGDTjoGN0bEbOnZj7Fzn6DhsDBK+UhLeUhK+MhK+MuK+zJ+9ZenbeuxPeUJgGLn94DIurF27drQPQQqMXlOSb3pNSb7pNSX5ptfUwPbt2+fo+8ZVyAMwegQl27b73N/pqquuYuXKlV1ft7a2Eo1GMY2LSBk+oANox6AdMptBO9jtGJnbsm83iGNgE0ztIZjaA/GNjo7b9vihOHNGMJT937721aT3+4vd/vXIKEgmk6xdu5alS5f2eTZZxC29piTf9JqSfNNrSvJNrylnOq9KHMy4CnmTJk1i06ZN3fZt2bIFn89HTU1Nn/cJBoMEg8Fe+/2Whd8GKMpsVQ6PwoSs8Nd76+u2FIaVhL2bYe9mHJ/P84X2B77imh5bP/t8vX9WGRl+v1//KEle6TUl+abXlOSbXlOSb3pNDczp3824CnlHHnkkDzzwQLd9jzzyCI2NjSP4YvACJZnNqSTp8LeP3iGwr1DYAZiQaofWT9KbU4EyF6GwJjN4Zly9DEREREREZACj+tv93r17ee+997q+XrduHa+88grV1dU0NDRw1VVXsWHDBn75y18CcMkll/Czn/2MlStX8jd/8zc8++yz3HLLLdx+++2j9SM45M9sZQ6/3yYdDJ2cJcy+zYbEnvS26yPnh1dUOXgoLAnv3x+sAI8Gz4iIiIiIjEWjGvL+9Kc/cfzxx3d93bl27oILLuC2225j48aNrF+/vuv2GTNmsGbNGv7hH/6B//qv/2LKlCn853/+p/v6hDHPAAKZrcLhfWwgTu/gN9DZw470XTt2pbcd7zs8PG+PENjPWcLs/YFSDZ4RERERERkBoxryjjvuuK7BKX257bbbeu079thjeemll4bxqMYrA/frCy26B8D+Lh3N/joBtgltW9ObU96Au1BYXAP+kPPHFxERERERYJytyZN88wDFmc0pE2dnCbsPnsFMwJ6N6c0pf7G7wTOhavAFXPwsIiIiIiKFRyFPXPICpZnNqc7BM06nkbYDFiT3we59sPtj508VLM8EwLCDM4eZYnuP18XPIiIiIiIytinkyQjI1+CZgcJhZvBMvDW97fzQ4XMZ6Qmj3cJfz697nEHU4BkRERERGcMU8mQMGsrgmb4uIe3v8tJ4+n7tO9Lb9ncdHl73wTPeoioWbGvD84eXoSyiwTMiIiIiMqoU8qRAZA+ecarn4JnBLiHtoK/BMx5gOsAzf+j/qRwPnsmqqtDgGRERERHJgUKeTGC5DJ5J0TMEmkacd+vKOXDTh3jZR+9waOY4eKbE3eCZ4mrw+l38LCIiIiJSiBTyRFzx0XPwjGV4eHvyAg7Y8ipey+rx/Z3rCwc6Y9jzMtIO0oNn2mB3G+xej2PBij4CoAbPiIiIiEwkCnkiwyp7fWG5w/vYQAJnl5Fm7wPiu9PbznXOj6/X4JlBegyLKrS+UERERGQMU8gTGXMMIJjZKh3exyI9SGagQTM9Q2GOg2c8vnQnYc/wVzJAbYW/WMFQREREZIQo5IkUBA8QymxOmewPff0Nm+kZDpNgpaBtS3pzylfk7Cxh9uYLuvhZRERERKSTQp7IhOUFSjKbUymcTyLNGjyT6oDWDenNqUCpi6EzNemzi179kyYiIiKi34hExAUf6VL7oRTbO+kxtCGxN73tcjF4pqjCQSjMqqkoqlSxvYiIiBQchTwRGUZDKbZ3cpYwex/QsTu97fjA4eF5+l5fONAlpcEyrS8UERGRMU0hT0TGmOxi+0qH9+kcPNPZU+ik5D4BtgX7tqU3pzx+55NIO7eAmy5GERERkaFRyBORApDr4Bmn9RSdWwqsJOzdlN6c8oXcFdv7nV4OKyIiItKbQp6ITFBeehbbD66/YvuBhtBYkGqH1k/SmwN+YJmnCN+6ut4BsKS/YvsqFduLiIgIoJAnIuKCP7MNZfCMk6mkNn6rA3Z9lN4cMSBU6bDYPnNbsEKDZ0RERAqQQp6IyLDJbfBM0kjwxMHTOe6tl/HZbQweCjtIF9vvTG/b33N4eN4eQdDBGsNAqQbPiIiIjHEKeSIiY4oBRoi2oknYxpT0cJhBWXQPggOtMeysqkiAbULb1vTmlDfgfvCM381aSRERERkqhTwRkXHPAxRnNqc6B8+4WWOYAjMBezamN6f8JZng53DwTKgafAEXP4uIiIhkU8gTEZmQhjJ4Jru4frBwaEGyDXa3we6PnT9VsNx5KNTgGRERkW4U8kRExKFcBs8kGPyMYXZozKwvjLemt50fOnwuDZ4RERHppJAnIiLDxACCmc0pm3SxfV+BsL+wGCf/g2fC/QyeKdHgGRERGfMU8kREZAwxgKLMVuXwPj0HzzgJh8kcB88EBzlT2Fe5fZHzxxcREckDhTwRERnnchk8k2LgSaR9rTs0wYzDnk/Tm1P+EnehsLgavH4XP4uIiEh3CnkiIjIB+UivLXSzvrAzGLopt88ePLPe+eEFK/oIgAOsMQxVavCMiIh0UcgTEREZlMH+wTPlDu+TPXhmsGCY/TUQ353edq5zeHgeKKpMB76ScK9QaAQrqdv9AcaGOiivTe8Plmt9oYhIgVLIExERGRbZg2cqHd7Hwt3gmY7099sWtO9Ib9vf7fWoPuAIgA+u27/T43N26Wi39YXFCoYiIuOAQp6IiMiY4QFCmc0pk/7XFqb/bNFOa8igon0HRufgGSsFezenN6d8RYNcOtrHfp+b6aoiIpIPCnkiIiLjmhcoyWx9Mz0enjh4ActefRW/ZdH3+sLBhtCYkOqA1g3pzalA6QBnCcO994eqwKtfT0REhmLC/iv6Z3sz9RRTjL9rC+LFQJehiIhIoctl8EwSZ8Nmsr+2IbE3ve1yMXimc32ho8Ez1envV7G9iEiXCRvyHuIdiuk+icyD0S30Odn8eBQMRUSkwBlAILNVOLxPX8X2g4XDjvRdO3altx3vOzw8D4T6K7bvJywGy7S+UEQK1oQNeTOowovNPpLsI0kCEwubvSTYS8Lx4/jwuA6GPvRpo4iIFLp8FNv39+c+Bs/s25benPL4nZ0lzP464KaLUURk9EzYkLfCmEsNga6vU1hdgc/plsIihUUrcVqJO37uAF5XoTCED6+CoYiIFLxciu0HHzzTe41hCqwk7N2U3pzyhfaHv5Ke6wn7CIqhavAFBn9cEZE8m7AhrycfHsoJUo6zKWA2NskewbCNBO2kBgyGFjYJTBKY7Oq8LMWBInyug6EuIxURkcI3+OCZ3vorth+oqsKEVDu0fpLenAqWO6unyB48o2J7ERkihbwcGRgE8BLASyVFju5jYxPH7DP89RcQ20liAx2k6CDFjs6i3EGPD0IuLyPV4BkREZkYch080zlp1OkaQxvirelt54cOn8uAUGWPNYYaPCMi7ijkjSADgyJ8FOGj2mEHkoVNxyBnB3sGxA5S2ND1tVMDDZ4p6XfwjD5tFBGRQpc9eMYpt4Nn9mW+34b2nekt34NnSrLXF5Zq8IxIAVPIG+Oyg5dTJlZX6OvvDGEbiW5fJ7FyGjzjdzl4JqTBMyIiMiHke/BMf/sSuQ2e8QbcD57xO/uAWkRGn0JeAfLioZQApQRwukYhidkrDA60xrCNBGZmXeJu4ux2MXgm6HrwTLqoQkREpLDlMngmRd9rCPsLh/sAE8wE7NmY3pzyF3eFP2+omkU7O/A88hSURjR4RmSMUcgTAPx48eN1NXgm0UcwHCwgWpl1iXFMdroYPBNyOXimSINnRERkQvABpZnNKafF9tlrEC1I7oPd+2D3x3iAKMALzw78VBo8IzIqFPIkJwYGQXwE8eU0eGaws4Sdt3UOnmknRTsptjscPOPBcB0MAxo8IyIiE4I/s5U7/H4bSJAdClNGB29NqeKQDe/j7TWQpvP7hjB4piv4hQe/nLSoQusLRXpQyJMRk4/BMz3XEva1xTPF9m0kaXMxeMY7wOCZ/jYNnhERkcJnAMHMlmYbHt6vXcBBG1/Fa1l93KevwTOD1VX0GDyz/T1nh+fxORs8k70/UKJgKAVNIU/GtFwHz7gttk9iYWKzhwR7XA6eCdoBPti2iR32ZsoHWW+owTMiIjIx5DJ4JrvYvr/LSHvelgQrBW1b0ptT3uAgobCPfX5nVy6JjAUKeVJwvHgoI0iZw/WFkB484zYYdg6e6SDOjuQOPmAHXuxBn0uDZ0RERPqSa7F99lAZJz2GJphx2PNpenPKX+IsFJaE968v9Dr/kFoknxTyREgPnqnAS4WL9YWJTDDcjclD1QcQ2/waSeJ9rDFMry3cl1lf6HbwTPqzUA2eERER6c3t4BmbdDAc7Axhz/0WJNtgdxvsXu/88IoqnA+dKa5Rsb3kjUKeSA6yB8+UGh6mBKewwNiK3+5rXUKa3U+xfXYI7Lm1Z4rtNXhGREQkHwxyGzwTx/klpJ1fAx2709uODxwenid9BtBRKMzsD5ZrfaH0opAnMkIMDEKZyy9rHN7Hwu43APYXDjV4RkREJJ+y1xdWOryPxf7BM04uI+1If79twb7t6c0pj89BKOyxL+DmklgZjxTyRMYwDwYlBCjBeZlsCov2Qc4QZt/WRpLUEAbPdAa+EgIO1hf68GrwjIiIFDwPEMps1Q7vkz14xkm5/T4glR48s3dzenPKV9Qj/IUHD4k+57MOZPQp5IkUGF8eBs+0ORg8Y2UGz+wmzm7ijp8rl/WFGjwjIiKFL9fBM30Fwr66CztvNyHVAa0b0ptTgTL3xfZeRY3Ror95Eclp8Ezc5UTSzmL7DlJ0kGKHw/WFBhByeRlpUOsLRURkQvABZZnNCRtIsj8EDjRsJvtrGxJ70tuuj5wfXlGl81AYKE9frip5oZAnIq7lo9jeydaRGTzT+bVTnpzWF3oUDEVEpMAZQCCzOdVz8IzTNYZAx670tuP9QZ/FDzRWLgG+4OLYpD8KeSIyInIttm93GQwTmcEze0mw18X6Ql/W+kKnm4rtRUSk8A118IzTqooECZ/TKgwZjEKeiIxZXjyUEqDUxSeOScxewbCNRL9hsY0EJjYpLFqJ0+pifWHAdbG9Bs+IiMhEkD14xpnkosN4PRVn6rAd08SikCciBcWPFz9eyh0OnrEzA2T6OzPYX0C0sElgksBkl8Nie3A2eMZvB2hNtbLPTlCu9YUiIjIReHyYXr3f5YtCnohMaAYGAbwE8FLpcvDMQGcIhzJ4xsTg3a2f8g6b8WF3DZ4pUbG9iIiIOKCQJyLiUvbgGafcDJ5pJUnAk75ENXvwzDaHz6ViexERkYlNIU9EZAS4GTyTNDysqVvA5ze+TIqE40tIOwfPDLXY3tn6Qg2eERERGasU8kRExiiv4aHIHvrgmYECYufgmVyK7YOuB8/4VWwvIiIyAhTyREQKyFAGzwy2xrAtc1axnRRWZl1iHJOdLgbPhBwMnsneivBpfaGIiIhLCnkiIhNYroNn3Bbbt5MCoJ0U7aTYPsDgme7Hh+v1hRo8IyIiE51CnoiIuGJgEMpcflnj8D4WNu0uQuE+ksQxsYE2krSRdHx8GjwjIiITnUKeiIgMOw8GJQQocbG+0Bygv7C/LYk1IoNnivGr2F5ERMYshTwRERmTvHgoI0iZw/WFkB484zQQtmX6CzV4RkRECo1CnoiIFAw/XirwUuFifWGij2DYGQD7C4g2uB48YwBFGjwjIiIjQCFPREQmLAODID6C+Kgi5Og+NnZmgMzAZwh7Dp6xcT94xoPhaiJpCYHM+UIFQxGRiUwhT0RExAUja7BLLoNnBjpL2HlbW6bY3sJ2PXjG18f6woAd4LU9HiL2J5T3ERpVbC8iUlgU8kRERIZZ9uCZiMP7pHIYPJPCIoVFK3Fas9YXmhi8uzdFnM14sXs9V8D1+kKfBs+IiIxhox7yrr/+en74wx+yceNGDj30UFatWsUxxxzT7/f/+te/5gc/+AHvvvsuFRUVnHzyyfzoRz+ipsbp56kiIiJjnw8P5QRzKrbvubWSwi6eyiFYJDKl9p2blVmXmMBkl4ti+1zWF2rwjIjIyBjVkNfS0sIVV1zB9ddfz9FHH81NN93EKaecwhtvvEFDQ0Ov73/qqac4//zz+clPfsJpp53Ghg0buOSSS7j44ou55557RuEnEBERGRsGKrZPGh7aKxawzAjgt62u/TY2cRcTSfdlLie1gQ5SdJBih4ti+5DLmoqgiu1FRHIyqiHvuuuu46KLLuLiiy8GYNWqVTz88MPccMMNXHvttb2+/7nnnmP69OlcfvnlAMyYMYOvfvWr/OAHP+j3OeLxOPH4/ktWWltbAXjIeJ86o7hrrHWvYlxDxbjiTNLj6fZfkaHSa0rybaDXlBcvZQQoc/hYlm3RQapX8OsrDLaRoD0TBgH2kGIPKXAxeKZnDcVAX3cW2xuGguFw079Tkm9JO30peTLpfA3yROT072fUQl4ikeDFF1/kyiuv7Lb/pJNO4plnnunzPkcddRT/8i//wpo1azjllFPYsmULd911F6eeemq/z3Pttddy9dVX99p/TziFP9ROf280PsNH0BMk6AkS8AS6/ty1Gd2/DngCeBUMJ7S18+aN9iFIgdFrSvJtJF5TfqAis3UybZOElSBuxYlbcRJ2+s8dVke3/V2bHSdlpfp5BjOz9X1pqdfw9nrPDngCvd63sze9f+dO/05J3ljpqwzWrl07ygcytu3bt8/R941ayNu2bRumaVJXV9dtf11dHZs2berzPkcddRS//vWvaWpqoqOjg1QqxfLly/nf//t/9/s8V111FStXruz6urW1lWg0ynnbQvgCZH3iuH+NgtnHonQgc3lKf28r3YtxQwT6LMIt6bk+wdAnYONd0uNh7bx5LH3tNfyWNfgdRAah15Tk29h+TXmB4sy2X9I2+6ii6PuMYfbgmf7EM1tfAt3ev/susu+5b6K/f4/t15SMR8mFC1lrWSxduhS/3z/ahzNmdV6VOJhRH7zS85IK27b7vczijTfe4PLLL+ff//3f+fznP8/GjRv5xje+wSWXXMItt9zS532CwSDBYO9F68fZDdTYgV77+yvGHawg1wZSpGglRaub9Qm21icUCr9l6Y1O8kqvKcm38fSa8mNQTIAaer9X92WgwTMDbRY2ZuYy0j0O378BQraK7WF8vaZkjMv8/u/3+xXyBuD072bUQl44HMbr9fY6a7dly5ZeZ/c6XXvttRx99NF84xvfAGD+/PmUlJRwzDHHcM011zB58uQhH9dQinHdvKl0ZIpxO792Knt9Qs+C3JI+zh6m1yeoGFdERArbQINn+tPX4JnOdYR9f9Cb6Hr/dltsb4DLmgp9sCsiuRu1kBcIBFi8eDFr167ljDPO6Nq/du1aTj/99D7vs2/fPny+7ofs9aavo7ftvi+xHAnZxbhOmVgDBsP2bm8q6a/jmWLcvSTYS8Lxc/VVjDvYpmJcEREpdAYGRfgowke1ww92Lexeg2cGC4idwdBtsX3PwTNONn2wKyIwypdrrly5ki9/+cs0NjZy5JFHcvPNN7N+/XouueQSIL2ebsOGDfzyl78E4LTTTuNv/uZvuOGGG7ou17ziiitYsmQJU6ZMGc0fxTUvHkoJUOrwMhTIbzHuYFSMKyIi0psnxw920x/g9v/hbs+AmBimD3ZLer1/64NdkUI0qiGvqamJ7du38+1vf5uNGzcyd+5c1qxZw7Rp0wDYuHEj69ev7/r+r3zlK+zZs4ef/exn/H//3/9HZWUlJ5xwAt///vdH60cYUfksxh1sfcJIFOOGCnR9goiISDYvHsoIUubw/RsgidkV+vo7Q9jWo9jexM7pg92g6w92tV5KZKwb9cErl156KZdeemmft91222299l122WVcdtllw3xUhSFf6xMG21SMKyIikl9+vPjxuvpgd6DBcX2/f6ewMu/7cUx2uvhgN2D7Wbf1YzbbWyhz8CFvoQ6eERmrRj3kydiS7/UJwzl4xunlKH7UfyQiIoUt18Fxbt+/2zPF9u2k2JPaw8e04u2neqr78bkbPFOMn4A+2BXJmUKeDFmu6xPcTiTNdX2CP2t9gtMzh1qfICIihc7AIJR5b6xxeB8Lm3aS7MbkwZpZHL75zySJD/j+Hc80ELsdPOPNafCMPtgVAYU8GSW5DJ5JZi5DcbpwvY0EZmZd4m7i7B6mwTN+O4BlqyNIREQKnweDEgIEDA+1gVoOMSL4B3kPTGH1OTW87/fx/YPjTGz2kGBPjh/sOt00OE4KkUKejBt+vFTgpcLh9/ccPDNQ91H27W4Hz5gYvLvpfd6wtztal6DBMyIiMpH4chw8018I7C8g5vrBbi6DZzx6/5YxTiFPClY+B88MFBD3ZtYndJAiSXJIg2d6jrbW+gQREZmI9n+w6+79u78A2F9AtMH14BmDwSeKlxDQ4DgZVQp5IllyGTwTB+6tPYTPbn6JZNY4656jrftan+B28IzWJ4iIiPSW/f6dy+AZJ5eQZg+OaydFOym2O/xgV8X2MtIU8kSGyGN4KPIWETZK8NvO3ljMHPoLkyO0PkHFuCIiMhHkMnhmJAfHDVZsr8FxMhCFPJFRkGsxrttgqPUJIiIi+ZPL4LhUDh/sprByKrZ3Mziucz6ABs8UJoU8kXEil/UJbotxh3N9Qs9NxbgiIjIR+PBQTtBVsX0yh2DodnBcJ7fv3xocNz4o5IkUqOEqxu25bqF9COsTQi7fWDR4RkRECl0+B8cNtLVnPtjtIEUHqSENjhts0+CZkaeQJyJdhlKMO9D0sp63xTPrE1SMKyIiMnS5DI6zHHyw23PrHDzjdnBc9uCZsziEOkq73Z60DR7ZBpucr2KRQSjkiciQdBbjlrhcn9BXEOwvILaR32LcnqOtsy9BKUbFuCIiUviyg5dT+Rg8k72Gf50VpMWMcJcZZttbBkfXebhwOH7YCUghT0RGXC7FuIkB+o86w2Fbj/3WMAye6SsgFuHT4BkRESl4+Rg8U2SHuMeqYXUqzPN2edf3hUsDVAScXTIqg1PIE5FxoXN9gtti3FzWJ+QyeEbrE0RERHrrHDzzsVVJixnhHrOGPZkI4jHguINqaYpFOeaAKtY+/LtRPtrCoZAnIgVpPK1P6DVgxg7wYXsF79nbqcgaTqNiXBERGU/22B7uM2toMSO8Zpd07a+vDNEUi3L24qlMqUy/RyeTzt9DZXAKeSIiGflYn+BkCE18kGJcE4N3d+1hG5vxYnftVzGuiIiMdbYNL9qlrDbDPGhW054ZgOb3Gpx06CSaY1GOPiCMx6MPLYeTQp6IyBAMRzFuKyn2BCcxiTbiJNiHinFFRGRs2277+I1Zw2ozwvv2/itoDqwtpSkW5cxFU6kucf5eKUOjkCciMsIGK8ZNGh5C1QtY9nE1ftsa88W4GjwjIjIxWTY8ZZWz2oyw1qokmfmQMOT3ctqCyTTFGljUUIlh6D1ipLkOeZs3b+Yf//EfefTRR9myZQu2bXe73TTNvB2ciIioGFdERMaWT+0Ad5ph7kiF2ZD1geWCaCVNjVFOWzCZsiLnSx8k/1yHvK985SusX7+ef/u3f2Py5MlK5iIiY9B4HjzTf7G9Bs+IiIyWpG3wqFXJajPME1YFdubf4/IiH2cumkpTLMohk8sHeRQZKa5D3lNPPcWTTz7JYYcdNgyHIyIio2Uog2faSAxYkJt9e89iXKc0eEZEZOS9bxVxhxnmbjPMtqz3hyNn1tC8JMrnD51Ekd87ikcofXEd8qLRaK9LNEX48EP4xS/gn/4JiorglVfgd7+DK68c7SPbb9cu+OlP4atfhUmTRvtoRApCLoNnkpjdAuFAAbEtM3jGxB7WwTMlXYNn/FpfKCITXrvtYY1VRUsqwv/YZV37I2VBzl48labGKNPDJQM8gow21yFv1apVXHnlldx0001Mnz59GA5JxqyPP4Zbb4WZM+Gv/mq0j2bsuvde6OiA5ub8P/bjj8Nbb8Ell+T/sUVGiB8vfrz9Dp7paaDBMwMFxFwHz4RyGDyjy0hFpBD8xSpmtRnhPrO6W2H5CQfX0hRr4PiDIvi8ukJiPHAd8pqamti3bx8HHHAAxcXF+P3dL+vZsWNH3g5OxpiXX4YlS+Cll2D3bqioGO0jEpEJINfBMx2kui4ldbO+sD1zv+0uBs+4vYw0oMEzIjJG7La93G9Ws9qM8HpWYXm0OkRzrIGzFk1lUoWzf3tl7MjpTJ5MQIkEvP46/M3fwN696csxjz3W3WO88AI8+2w6IFZVwTHHwIIF6dsefhh27IBzz01//dxz6X3nnguzZ6f3/exncMQR0NiY/vrll+GZZ2DnTqishMMPh1hs//Nt2AC//S1s3Qq1tennG0x7e/oy03fegVQKpk+Hk0+Gmpr07X2cSfM8+yxLf/Yz+Pu/T9/+6qvpG66+Ov3fCy5IH99PfwpnnQXPPw8bN0J1NSxbln4O6PsS17fegpYW+OY307c/8UT3xz79dND6WJFeDAxCmcsvcxk8k32GcKCQGMfEBtpI0uZi8Iw3p8EzWvMiIvlh2/CCXcrqVIQ1VjUdmfXLAa/ByXMn0xyLcsTMGhWWj2OuQ94FF1wwHMchY93rr0M4nN7mz4eHHoLPfhacTld98810gDn55PTlnu+8A/fdB+XlMGNGOui8/HL6Xx3DgI8+guLi9H9nz04Hy+3b9weiF19MB6ply9Lr6zZtggceAL8/HXoSCfh//y/92GeckV6P97vfDX6c992Xfp7mZggG4fe/Tz/OpZeC18EvWEcdBdu2QTyeDmAAoRDs2ZP+89q18PnPQySSDry33w5f/3r6Zx3MoYfCli3w3ntw/vnpfUFnl7uJyOCyB8+EcfD/SdKDZ9z2FyaxMLHZQ4I9LgbP+HMYPKNiexHJttX28RszTIsZ5oOswvKD6spYEYty5sJ6qlRYXhByKkM3TZN7772XN998E8MwmDNnDsuXL8fr5JdgGZ9efhnmzUv/edasdIhaty4d2Jx49tl0+Oo803bkkfDJJ+n9M2bAtGnpx9y4ESZPToe7o45Kh0NIP1dJSTpkAvzxj3DSSXDIIemvq6rSZ+xefDH9PK+9lg6Mp5+eDn61tdDaCg8+2P8xbt8Ob78NF14I0Wh635lnwk9+kj6jduihg/+cgQD4fOmzgKWlvW+PxWDOnPSfv/AFeP/99N/t0UcP/th+f/rxPZ6+H1tERpwXD2UEKXO4vhDSg2fcBkMzsy5xN3F2uxg8E8wMngnaAd7ZsYuEvZHyAYbRaPCMSOExbfijVUGLGeb3VhWpzP/HiwNeli+YQlMsymFRFZYXGtch77333mPZsmVs2LCBgw46CNu2eeedd4hGozz44IMccMABw3GcMpq2bUtf+rhiRfprjycdeF5+2XnI27oVFi3qvi8aTV+6COmJnJMmpad0er3ps3mLF6fP1sXj6f2dZ/Ha2tKB7f7702fvOllW+nE6n6+uLh2MOk2dOvjP6fFAff3+fcXF6WC5bZuzn3MwneER0s81ZUr+HltExgU/XirwUuFifWECk7ZMYX3/NRXdb7eBOCZxTEzifBr/lD+zGS/9T8g2gCINnhEpCB9b6cLyO80IG7MmIB8WreTcJVFOnT+F0mBO53tkHHD9v+zll1/OAQccwHPPPUd1dTUA27dv56/+6q+4/PLLeXCgMyUyPr38cjpAXXdd9/0eT3oNW8jZepc+L+3M3jdtWvoMntebDnShUPqyxo8/Tu8//PD093VWeJx2Wu/gNpRPofqrBsne39fjm2buz5ltOB9bRMYtA4MgPoL4wOH6QrtHsf0uUvy+YiaLNr9OgnifATHXwTMeDNcTSTV4RmR4JGyDtVYlq80IT1nlXYXllcV+zlhYT3OsgYMmlQ3yKFIIXIe8J554olvAA6ipqeF73/seRzu55EzGF8tKDxI56SToeZb2jjvSl0UuWTL440QisH79/kErkA5vnZdfwv51eR5P+hJOSAe/v/yl+3q80lIoK0sPXJk/v//n+/OfIZncfzbvk08GP0bLSp+17Dzjtm9f+rkjkfTXxcXp9YGdawcBY9Om7o/j9fYfGD/5JP0zQfq5Pv10/99fcXH6rGUikb4sE9JrDZ0+tohIRvbgmRpgkuHhg+KZHGXswW9bfd7Hwu7zTOFAZxDjmWL7XAfPlBBQsb1IHrxnFdFiRrjbrGFHVmH50bNqaIo1cNKcOhWWTzCuQ14wGGRP5xCJLHv37iUQGD8LNX9pv0wtoUHfVCb8+oR33kl3vi1cuP9SyE6HHLK/VmEwRx0Fd96ZXm83Y0b6cd98c/8AEdi/Lu/tt+H449P7pk9Ph8ni4v1BC+C449LDX4JBOPDA9Bq4Tz9NH+uRR6bXDz72WPqSzs9+Nj145dlnBz7Gmho46KD0JaBf+EI6aD36aHo4zEEH7T+eNWvg6afTa+veew/jvfe6r5GrrEyvtdu2LX3c2cNRXnghPVWzc/BK598tpM9K+v3p5zz88HTY7JzUmf3YO3emw195+f41gCIiQ+TBoIQAJS6K7VOZwTM9Q2DbAJeVpnIcPOO02H7/+7dPg2ekoO2zPTxoVtFiRvhTVmF5XXmQcxZHWdEYpaHG2RApKTyufzv8whe+wN/+7d9yyy23sCTzy/3zzz/PJZdcwvLly/N+gMNlI3vZ7fBSlAldjNu57q5nwIN0yHnqqfSwlMEcfHB6suYzz6TDWVVVeihK59k52L8ub/fu/YFu2rT0mavs74P0+j6/P/14v/99+s91dfsv6QwE0vULv/0t3HRT+vE+97l0YBzI6aenp3D+v/+XvlRy2jQ477z9kzUjETj1VHjyyfTwlzlzsI46Kn3WMPvYPvwQfv7zdGjtrFAAOPHEdEDctCn9d9DcvH+yZiiUHvSydm26i3DmzHRNxW9/u/+xDzkkHY5/8Yt0QFSFgoiMIh8eygnmpdh+oC3XYvtc1hdO6A92ZcyzbXjNTheW32/WsDdTreL1GBx/UC3NsSjHqbBcAMO23V37tWvXLi644AIeeOCBriL0VCrF8uXLue2226gY4wXZra2tVFRU8OyVXyUYtAd8U2knldNzqBh3Ykl6PKxZsIBlr76K3+r7Mih27Ur35H31q+kgKzIAR68pERfG82vKxiaemUjazuCXkGYPnnHLAEIu37+DE/T9ezy/psaj3baXe80aVpsR3rT3n52bVlPMisYoZy+eSl35+C4sTyaTrFmzhmXLlnVlDOmtM8vs3r2b8vLyfr/P9Zm8yspK7rvvPt59913eeustbNtmzpw5zJo1a0gHPNIONKqpGeSSlP7WJww02SyhYlwREZG8MTAowkfREAbPDLTGsPPS0s7BM537nfLk9P7tmZDBUNyxbXjeLqMlFWaNVU28s7Dc5+GUuZNoikU5YoYKy6VvOS/mOfDAAznwwAPzeSxjTq7rE/oLhn198tg2hPUJKsYVERHprefgGSdMLNr7CYb9bYnM4Jm9JNjr4v3bl8P7twbPTBxbbB93m2HuMCOss/efnTt4UhnNsShfXFhPZfH4mYMho8NRyFu5ciXf+c53KCkpYeXKlQN+73U9x+xPML4ci3GdXHqSz2Jc5wvXJ/jgmXyprIRvfnO0j0JERPrhxUMpAUpzGDzjZkthkcKilTitLt6/NXimsKUyheWrzQiPWpWYmd+9SgJelh9WT3MsyvypFSosF8cchbyXX36ZZDLZ9WfJLz9eKvFS6bIY18l46+zbsotxdzpcuK71CSIiIn0rtMEzoUIaHDdOfGwFuMOMcKcZZlPWBwyLGippXtLAqfMmU6LCcsmBo1fNH/7whz7/LKMjuxi3aojrEwYKiO1anyAiIpI3BgYBvARcfrAbz/pgt/f7eKLXZabtmQ92O0jRQYodDqeJu/lg128HSFgJXM7vEyBuGzxiVdFihnnK2j+wsKrYz5mLptIci3JgnQrLZWhcfzRw4YUX8tOf/pSysu4vvra2Ni677DL++7//O28HJ/mTy/qEwQbP9BUQ41qfICIikjfZg2eqHX6wa2V9sNtXCOwrILodPGNi8O7mD/kLWyhzecZwon6w+06msPw3Zg072T898pgDwzTFoiydU0fQp6F6kh+uQ94vfvELvve97/UKee3t7fzyl79UyCsgQxk8M9gaw87bOwfPaH2CiIhIfmRfURPGWRl2X4NnBgqIrZkgONwf7JZkPqAerx/sttkeHjSrWW1GeMku7dpfVx5kRWO6sDxarcJyyT/HIa+1tRXbtrFtmz179lCUVY5tmiZr1qyhtrZ2WA5Sxo9cB884vYS0c9S1inFFRETyx+3gmaTh4f5Jh3Lc5hdJkhj0MtK2zPeY2Dl9sDueBsfZNrxql9BihrnfrKEtq7D8xINraYpFOXa2CstleDkOeZWVlRiGgWEYzJ49u9fthmFw9dVX5/XgZGLw46UCLxV5Wp/Q9/rC4Vuf4LcDbEtsY4e9jwp8GjwjIiITgs/wUW4U4bedBcOeg+OcvX+nsDLv+24GxwGEcvhgdyjv37tsL/eYNbSYEd7KKiyfXlNMU6yBsxbXU1s2vgvLZfxwHPL+8Ic/YNs2J5xwAnfffTfV1dVdtwUCAaZNm8aUKVOG5SBFsg11fYLTzen6BBODd7dv4UM248XW4BkREZE+5HNw3GDBEKCdFO2k2O7ig92B3qtnU0NNj8tfLRues8poMSM8ZFWRyFxWGvR5WDZvMk2xKIfPqFb1gYw4xyHv2GOPBWDdunVEo1E8Hp1ilvEjO3g55bQYdw8pPvWWEMCLmfnEUYNnREREhm44Bsf1tcUxsYG2zNKQvlRQ1BXytth+7jTD3GGG+SirsPyQyeWcuyTK6QvqqSh2/juHSL65Hrwybdo0APbt28f69etJJLr/Ijt//vz8HJnIKHO6PiFpeJhSu4BlG1/FsPuqqeh/4fo+DZ4RERHJq6EMjhsoCJbbIX6fKSz/Q1ZheWnQx/LDpnBurIG59eU6aydjguuQt3XrVv76r/+ahx56qM/bTdMc8kGJjFe5FOMmMB2Nt87eVIwrIiKSPwMNjltvBWkxw3zfDLMlKzg2TquiKRbl1PmTKQ6osFzGFtevyCuuuIKdO3fy3HPPcfzxx3PPPfewefNmrrnmGn784x8PxzGKFKzs9Qlui3EHO0OYfftwF+OWZP05oMEzIiIyznVkFZY/nVVYXl0S4KxF9TTFosyqVWG5jF2uQ95jjz3GfffdRywWw+PxMG3aNJYuXUp5eTnXXnstp5566nAcp4hkZA+ecarn4BknAbFzfcJgg2d68uY0eEblryIiMvretkKsNsPcY4bZlXmfNQz4zKww5y5p4HOH1BHwaemDjH2uQ15bW1tXH151dTVbt25l9uzZzJs3j5deeinvBygiQ5fr4Bk3i9b3kSSJhYnNHhLscTF4xu9y8Mx4LsYVEZGxZa/t4beZwvJXsgrLJ1cUcU5jlBWNU5lapcJyGV9ch7yDDjqIt99+m+nTp3PYYYdx0003MX36dG688UYmT548HMcoIqPAm4dieyebiU0Si93E2V2gxbgiIjK22Da8bJfQYkZ4wKxmX+aKEp/H4HOH1NG0JMpnD4zg9eh9Q8annNbkbdy4EYBvfvObfP7zn+fXv/41gUCA2267Ld/HJyLjSC7F9m6Lcfdl1heOh2JcEREZW3baXn5jhmkxw7yTVVg+M1zCiliUsxZNJVLm/MNNkbHKdcj70pe+1PXnhQsX8uGHH/LWW2/R0NBAOBzO68GJSGEb68W4HgzXwVCDZ0RExhbLhmetMlabER7uUVh+6vzJNMcaiE2vUvWBFBRXIS+ZTHLQQQfx29/+ljlz5gBQXFzMokWLhuXgRER6yncxbltm+mhfxbgW9oDFuH3R4BkRkbFhk+3nrsxZu4+zCssPnVJO85IGli+YQkVIheVSmFyFPL/fTzwe1ycdIjKuDKUYt78Q2DMgtmWK7fMxeCZgB3m9NUWp/QnlPdYelhBQsb2ISD9SNvzBqqTFDPOYVYmVubKiLOjj9IVTaI41MLe+YpBHERn/XF+uedlll/H973+f//N//g8+n4ofRaQwDVSM2598DZ4xMXi3DWw248Xu87ncFtsX4dPgGREpWB9aQe4ww9zVo7B8yfRqmmJRls2bTCigqyZk4nCd0p5//nkeffRRHnnkEebNm0dJSUm323/zm9/k7eBERMaTfA2e2U0Kf2kD8zZ7SRDvc/DMcBXbZ29BrS8UkTGswzZ42KpitRnhWau8a39NSYCzF09lRSzKAZHSAR5BpHC5DnmVlZWcddZZw3EsIiITSn+DZ5KGh51lC1hmgN+2ut3Hxs4MkHG+xrCdVE7F9h6H6wuzh9No8IyIDLc3rRAtZoR7zBp2ZxWWHzs7QlNjlBNVWC7iPuTdeuutw3EcIiLigJEVvHIZPON0jWHn4Jm9JNjrYn3hQINnSgj0GRA1eEZEBrPH9vCAWUOLGebVrMLy+soQKxqjnN04lfpKZ1OaRSaCnBbVpVIpHn/8cd5//33OO+88ysrK+PTTTykvL6e0VKfFRUTGkuzBMxGH9+kcPNNXCGzrJyTmOngm4LrYXoNnRCYC24aX7FJWm2F+a1bTnvlAyO81OGnOJJpiUY6eFVZhuUgfXIe8jz76iJNPPpn169cTj8dZunQpZWVl/OAHP6Cjo4Mbb7xxOI5TRERGkNvBM3ZmgIzbwTNWZl1iApNdLortNXhGpHBtt33cY9aw2ozwnr3/7NwBkRLOXdLAGQvrqSlVYbnIQFyHvK9//es0Njby6quvUlOz/2KhM844g4svvjivByciIuODgUEALwG8VLoYPBN3OZG0XYNnRAqSZcNTVjktZoRHrEqSmbP1Ib83U1geZfE0FZaLOOU65D311FM8/fTTBALd+6amTZvGhg0b8nZgIiJS2AwMivBRhI9qnK2lsbDpIOUqGHYM8+CZ7sX2HgVDERc22n7uNMO0pCJsyLpyYP7UCppiUU5bMIXyIhWWi7jlOuRZloVpmr32f/LJJ5SVleXloERERPqSHbycMrFo7yMYtpHoc/8+kiRyHDzj61Fs72TzaX2hTDBJ2+BRq4IWM8ITVsX+wvIiH2csrKcpFuXQKSosFxkK1yFv6dKlrFq1iptvvhkAwzDYu3cv3/zmN1m2bFneD1BERGQovHgoJUApgcG/OSOJ2W8A7CsgtpHAxCaFRStxWjPF9k64HTxTnDlfKDLerLOCtJgR7jLDbMv6oGbJjGrOXRLllLmTKfJr2q5IPrgOeT/5yU84/vjjmTNnDh0dHZx33nm8++67hMNhbr/99uE4RhERkRHlx4sfL+U5DJ4Z6Axhz9uHY/BMSR/7vLbzgCuSTx22wUNWNatTYZ639xeWh0uD6cLyxqnMVGG5SN65DnlTpkzhlVde4fbbb+ell17CsiwuuugivvSlLxEKqZ9EREQmnnwNnskOgG0k+hg8kwLcD56xMFi/ZR0f2VspcziZVMX2MhSvZxWW78n8uukx4LiDammKRTnh4Fr8Xl2qLDJccurJC4VCXHjhhVx44YVDPoDrr7+eH/7wh2zcuJFDDz2UVatWccwxx/T7/fF4nG9/+9v86le/YtOmTUydOpV/+Zd/ycuxiIiIjJRcB8/0V2bf3xbHTE8kNTvYyj52YDt6roGK7fsfPKNL7SayVtvL/WY1LWaE1+ySrv1Tq0I0ZQrLJ1fohIDISMgp5L3zzjs8/vjjbNmyBcuyut327//+744fp6WlhSuuuILrr7+eo48+mptuuolTTjmFN954g4aGhj7vs2LFCjZv3swtt9zCrFmz2LJlC6lUKpcfQ0REZFzJLrZ3ysRiNyYPhGdz5OZXSBIfNBgmcyy29+cweEbF9uObbcMLVimrUxEetKroyAT9gNfgpEMn0Rxr4KgDavCosFxkRLkOeT//+c/52te+RjgcZtKkSd36SgzDcBXyrrvuOi666KKufr1Vq1bx8MMPc8MNN3Dttdf2+v7f/e53PPHEE3zwwQdUV1cDMH36dLc/AgC27eyTTBERkfHMi4cyw0elv5IZRhV+2xr0PkmX/YX7SGJm1iXuJs5uF4Nngi4Hz4Q0eGZM2G77eOxTg592zOGDrMLyA2tLac4UlleXaC2oyGhxHfKuueYavvvd7/JP//RPQ3riRCLBiy++yJVXXtlt/0knncQzzzzT533uv/9+Ghsb+cEPfsD//b//l5KSEpYvX853vvOdftcDxuNx4vH9bzatra0AfN94mmqjqNubRs83kd7FuD6VcEovSY+n239FhkqvKck396+pzjNyzth2eoBM7/DXfQhNe2Zf59fp/kKLfcTBYTA06D54JkRg0PfvIr1/54VpwzNWOXekwvzerCD1kQcIEfJ7OHXeZFY01nPY1Iquv+tk0nkvpUjn60Wvm4E5/ftxHfJ27tzJOeec4/qAetq2bRumaVJXV9dtf11dHZs2berzPh988AFPPfUURUVF3HPPPWzbto1LL72UHTt28N///d993ufaa6/l6quv7v1YkQgfh3r2LFkwwBuNx/AQ9AQJeAIEjSBBz+Cbz8jpilgZh9bOmzfahyAFRq8pybfRek15gNLM1sm2bRJ2grgVJ27FSVj7/xy34/v/3OP2vg3+/t3zvTvgCfR+3za6v38rGKbtiMPzWwye3+JhZ2L/38m0Upsjai0W1aQo8n3Extc+YuNro3igUhDWrl072ocwpu3bt8/R97lOIOeccw6PPPIIl1xyieuD6kvPf0Bt2+73H1XLsjAMg1//+tdUVKRLMq+77jrOPvts/uu//qvPs3lXXXUVK1eu7Pq6tbWVaDTKD7YdQFGAHp8uJmnL+nP24vYEvQvgOyUy254+bsten9Dz08bOTxxDPSad+QwtXB9Pkh4Pa+fNY+lrr+G3Br8MSmQwek1Jvo2/11T6fB09JpVadu9i+74G0bRnnS2MD/D+ncxse/u4LXvwTM/375J+rgDyF9D7d8I2eMys4M5UmCetcuzMJbIVIR/LF0zhzAV1rP/zsyxduhS/v+eH5iLuJZNJ1q5dq9fUIDqvShyM65A3a9Ys/u3f/o3nnnuOefPm9fof4fLLL3f0OOFwGK/X2+us3ZYtW3qd3es0efJk6uvruwIewCGHHIJt23zyyScceOCBve4TDAYJBnv3HFXaPmpc9AZ1FuMO1H/Uc/S1iY2FyV5M9rroPwraWp8wHvkta5z88iTjhV5Tkm+F8JoK4qMSHzicSGpm+gudbm0kSWEBNvuIZy4ldSbg+v3bN+YGz7xvFXGHGebuHoXlRx1QQ1MsyucPnUSR30symWT9n8Hv9+sXcskrvaYG5vTvxnXIu/nmmyktLeWJJ57giSee6HabYRiOQ14gEGDx4sWsXbuWM844o2v/2rVrOf300/u8z9FHH82dd97J3r17KS1NX/Txzjvv4PF4mDp1qtsfxZWhFOO6Kce1Mr1JcUx2ugiGPc8GDrYV4VP/kYiIFDwvHsoIUubw/Rv6HzyTfbVPX4Nn8lls39/7+nB8sNtue1hjVdGSivA/dlnX/khZkHMWT2VFY5Tp4ZIBHkFExhrXIW/dunV5e/KVK1fy5S9/mcbGRo488khuvvlm1q9f33Up6FVXXcWGDRv45S9/CcB5553Hd77zHf76r/+aq6++mm3btvGNb3yDCy+8cMwVsedajNtBqlcYHKgct7MYtz1zv+0Oi3ENcD3mWsW4IiIyEfjxUoGXChfv3wnMAUNgX7fZuC+27zl4pq+tpMflpcF+3r//YhWz2oxwn1ndrbD8hINraYo1cPxBEXwqLBcZl3KeCpJIJFi3bh0HHHAAPl9uD9PU1MT27dv59re/zcaNG5k7dy5r1qxh2rRpAGzcuJH169d3fX9paSlr167lsssuo7GxkZqaGlasWME111yT648xphgYmTV6/hEpxm3LvOk4pWJcERGR3gwMgvgIuriMtPODXSdnCTtvayeFjfsPdo9jOscxHYDdmcLy1WaE17MKy6PVmcLyxVEmVTgLtyIydrlOZ/v27eOyyy7jF7/4BZC+XHLmzJlcfvnlTJkypVclwmAuvfRSLr300j5vu+2223rtO/jggzV1J0uuxbhu+49UjCsiIpI/2R/s1ji8T/YHu20O378TmIRsP8/bpbSkIjxoVRPPvM92Fpafu6SBI2eqsFykkLgOeVdddRWvvvoqjz/+OCeffHLX/s997nN885vfdB3yZOTlc32CinFFRERGRvYHuxEH37/V9nGnWcU/mTV8mLXWbnZdKc2xdGF5lQrLRQqS65B377330tLSwhFHHNGt6mDOnDm8//77eT04GTtyXZ/gNhja4HrwjJP1CRo8IyIiE4Fpwx+tClrMML+3qkhl3uuKA16WL5hCUyzKYdFKdQCKFDjXIW/r1q3U1tb22t/W1qZ/MKRL9vqEqhzXJzjZcl2f4MFwPZFUg2dERGSs+sQOcEcqzJ1mhI1ZSzgWRCs5NxblCwumUBrMeRSDiIwzrv/fHovFePDBB7nsssuA/WXmP//5zznyyCPze3QyoQx1fYLThetxTCzsvA2eCdhB3mkLEbU3U9HjUlMNnhERkeGSsA1+b1VyuxnhqazC8spiP2csrKcpFuXgSeWjfJQiMhpch7xrr72Wk08+mTfeeINUKsVPf/pTXn/9dZ599tlevXkiwy2XwTMprD6DYH8BsbMYt7/BMyYG77a208pmvNjdbnMzeKaEwJgsxhURkbHlPauIFjPC3WYNO9hfjHz0rBqaYg2cNKeOIr8+ZBSZyFyHvKOOOoqnn36aH/3oRxxwwAE88sgjLFq0iGeffZZ58+YNxzGK5JUvh8EzCcx+zxC2YtJRNJnpJEgQH/LgmVzWF2rwjIhIYdtne3jQrKLFjPCnrCEqdeVBzl48labGBhpqikfxCEVkLMnp4ux58+Z1VSiITASdxfZ9DZ5JGh6oWsCy9WX4bQtIry+Muxw8007uxbghlzUV/RXjiojI2GHb8JqdLiy/36xhb2YJgNdjcPxBtTTHohynwnIR6YPrkNfa2trnfsMwCAaDBAIaxStiYFCEjyJ8rort3Q6e6cgMnun82ilPTsX2HgVDEZERsNv2cq9Zw2ozwpv2/rNz02qKWdEY5ezFU6krV2G5iPTPdcirrBx47O7UqVP5yle+wje/+U08Hn2yJOJUdvByysSi3WUwTGQGz+wlwV4Xxfa+HIrtfVpfKCLiiG3D83YZLakwa7ILy30eTpk7iaZYlCNmqLBcRJxxHfJuu+02/uVf/oWvfOUrLFmyBNu2eeGFF/jFL37Bv/7rv7J161Z+9KMfEQwG+ed//ufhOGYRyfDioZQApS4Hz7jtL0xhkcKilTitLtYXBlwW2xer2F5EJpgtto+7zTAtZoQP7f1n5w6eVEZzLMoXF9ZTWayrpETEHdch7xe/+AU//vGPWbFiRde+5cuXM2/ePG666SYeffRRGhoa+O53v6uQJzIG+fBQTpByh4Nn7MwAmf3TRhMDnj3svN3CJoFJApNdDovtwf3gmZCK7UVknEllCstXmxEetSoxM/+GlQS8LD8sXX2wYGqF+odFJGeuQ96zzz7LjTfe2Gv/woULefbZZwH4zGc+w/r164d+dCIy6gyMrsEzlX0MnulLf4NnBgqIGjwjIoXuYyvAHWaEO80wm7KuwFg8rYqmWJRT502mRIXlIpIHrv8lmTp1Krfccgvf+973uu2/5ZZbiEajAGzfvp2qqqr8HKGIjDu5Dp5pJ9lvCOwrIOY6eKa/YvuBB8+oc0pE3IvbBo9YVbSYYZ6yKrr2VxX7OWvRVJpiUQ6sKxvgEURE3HMd8n70ox9xzjnn8NBDDxGLxTAMgxdeeIG33nqLu+66C4AXXniBpqamvB+siBSuXIrtOwfPdAbANhKDri9MDlBsPxA3xfbpy0g1eEZkInvHKmK1GeEes4admYFahgGfmRWmOdbA5+bUEvTpwyMRGR6uQ97y5ct55513uPHGG3n77bexbZtTTjmFe++9l+nTpwPwta99Ld/HKSLSSy6DZ5Iu+wuHUmwfdDl4JqTBMyLjWpvt4bdmNavNCC/bpV37J5UXsaJxKuc0RolWq7BcRIZfThd+T5s2jWuvvTbfxyIiMuz8eKnop9i+L3ZmgIy7Yvv04Jk4JnFMdroYPBPCR9AO8N72LbTamykfJCgWafCMyKiybXjVLmG1GeEBs5q2zKXdPo/BiYfU0hxr4LOzI3hVfSAiI0ire0VEBmBgEMRHEB9VDtcXdg6eGWjQTM9LSzvXF7aTYi8m2xLbeIdteLEHfC4PBiGXE0kDGjwjMmS7bC/3mDW0mBHeyiosnxEuoSkW5cxF9dSWqbBcREaHQp6ISJ5lD55xqnPwzD6S7Mbkd1UziW3+Cwni/Z4xjGeK7dtI0qbBMyLDzrLhOauM1WaE31lVJDLrboM+D6fOm0xTLMqSGdWqPhCRUaeQJyIyBmQPnqk0PESLoiwyduC3rX7vY+ZQbD9Sg2eK8ePV4BkpEJttP3eZYVrMMOuzCsvnTC7n3CVRlh9WT0XIP4pHKCLSnUKeiMg45cVDGUHKHBbbgwbPiDiVsuHxTGH5H7IKy8uCPpYfNoVzlzQwt75ikEcRERkdrkPet771Lf76r/+aadOmDcfxiIjIMBquwTNtmUL7zq9tcD14xgCKXK4v1OAZybePrCB3mGHuNMNsyZrcG5teRVOsgWXzJlEc0GfkIjK2uf5X6oEHHuCaa67h2GOP5aKLLuLMM8+kqEgLi0VEClGug2c6egyc6RkCe96WPXimnRTbaXf0XBo8I/nQYRs8bFXRYkZ4xirv2l9dEuCsRfU0xRqYVVs6wCOIiIwtrkPeiy++yJ///GduvfVW/uEf/oG/+7u/o7m5mQsvvJBYLDYcxygiIuOIgUEoc/lljcP7ZA+ecRoOh3PwTAmBrMtIfRo8U6DeskK0mGHuMcPsyvxKZBhwzIERmmNRPndIHQGf1paKyPiT0/UG8+fP5yc/+Qk//OEPeeCBB7j11ls5+uijOeigg7j44ov5yle+QkWFrlMXERFnsgfPOJXKDJ4Z6Axhz32pHAfPBFyvL/Rp8MwYtTersPyVrMLyKRVFnNMY5ZzGqUytUmG5iIxvQ7qo3LIsEokE8Xgc27aprq7mhhtu4N/+7d/4+c9/TlNTU76OU0REpBsfHsoJUu5w8IydGSDjdvCMlVmXmMBkl4ti+1zWF2rwzPCwbXjZLqElU1i+L6uwfOmcOppiUY45UIXlIlI4cgp5L774Irfeeiu33347wWCQ888/n//6r/9i1qxZAPz4xz/m8ssvV8gTEZExw8AggJcAXipdDJ6Ju5xI2p4ZPNNBig5S7HC4vtAAQi5rKoJaXzigHbYvU1ge5p2swvKZkRKaY1HOXDSVcKnz6bQiIuOF65A3f/583nzzTU466SRuueUWTjvtNLze7msVzj//fL7xjW/k7SBFRERGQ3axfbXDwTNWH4NnBts6B890fu2Ux2Gxvd8OsDe1l4SdyswjLdxgaNnwjFXOajPMI1mF5UV+D6fOm0LzkiiN06pUWC4iBc11yDvnnHO48MILqa+v7/d7IpEIltV/ga+IiEihyg5eTplYtA8SDNtIdPueRGbwzF4S7B1kfaGJwbtbP+EtNhN0GAyzN984WF+4sauwPMIn9v6zc/PqK2iKRVl+2BTKi1RYLiITg+uQZ9s2VVVVvfa3t7fzwx/+kH//93/Py4GJiIhMFF48lBKg1MXgmSRmt9DXMwRmb60k8RjpoJbCopU4rS6K7d0OnikeoWL7pG3wB6uCFjPCH6wKrM7C8iIfXzysnqZYVIXlIjIhuQ55V199NZdccgnFxd0nT+3bt4+rr75aIU9ERGQE+PHix+to8EzS8PBg3XyWbnqJJIlBzxKO1OCZkhyL7T+0grSYEe4ya9iaFYyXzKimORbllLmTCQVUeyEiE1dOZ/L6uo791Vdfpbq6Oi8HNRI22nvwUpJeq4CnoNcniIiIGIZBwPBRYntcDZ7pyBTUO7mMtI1E1/rCXAbP9AyBB1BNI1OAdGH576wqVpsRnssqLA+XBjhr8VRWNEY5IKLCchERcBHyqqrSi5QNw2D27Nndgp5pmuzdu5dLLrlkWA5yOPySVyjuHKGMpyDXJ4iIiAxFdrG9m8Ez/XUX9rfFMbGhV7F9CQHetEK0mBHuMWvYnfm1xWPAsbMjNMUaOPGQWvxevSeLiGRzHPJWrVqFbdtceOGFXH311d3KzgOBANOnT+fII48cloMcDmUEMbBIZbbhXp+gYlwREZkIcim2N7uK7dNnBLfaKR62illlRnnPmtn1ffWVIVZkCsunVDoLnSIiE5HjkHfBBRcAMGPGDI466ij8/vE9oepSYwnV+Md0MW7IxfoEERGR8cqLh1I7yDt2DavNML81q2nPXG3j9xqcNGcSTbEon5kVxqPCchGRQTkKea2trZSXp69/X7hwIe3t7bS3932Nfef3jQcqxhURERld2zOF5avNCO/Z+8/OzaotpTkW5YyF9dSosFxExBVHIa+qqoqNGzdSW1tLZWVln4NXOgeymKaZ94McSwqlGLdbSa4Gz4iIyAiybHjKKqfFjPCIVUkys5wh5PfyhfmTaV4SZVGDCstFRHLlKOQ99thjXZMzH3vsMf2j69JwFeP2nGzmphg3m5PBM9ljrkMaPCMiIjn41A5wpxnmjlSYDVnVD/OnpgvLT1ugwnIRkXxwFPKOPfbYrj8fd9xxw3UskiUfxbiDdR/tI5nz4Jmg68EzI1OMKyIiY0vSNng0U1j+RFZheXmRjy8uTBeWHzpFheUiIvnkuifv1ltvpbS0lHPOOafb/jvvvJN9+/Z1DWiRkeemGBfS6wsTfQTDwQKilVmXGMdkp4vBMyGXg2fcFOOKiMjY8kGmsPxuM8y2rCtZjphZTXOsgZPnTqLIr8JyEZHh4Drkfe973+PGG2/stb+2tpa//du/VcgbRwwMgvgI4stp8MxAZwmzb+scPNOeKdTdPoRi3MG2gAbPiIiMmnbbw0NWFatTYf7Hzi4sD3JOY7qwfEa4ZBSPUERkYnAd8j766CNmzJjRa/+0adNYv359Xg5Kxq58DZ4ZLCD2V4w7GG9Og2f0SbKIyFC8boVYbUa416xhT1Zh+fEH1bIiFuWEg1VYLiIyklyHvNraWv785z8zffr0bvtfffVVampq8nVcUkDyNXimjcSAA2iSWJjY7CHBHheDZ/wOBs/03FRsLyITXavt5X6zmhYzwmv2/rNz0eoQTY1Rzl4cZVKFs6tEREQkv1yHvObmZi6//HLKysr47Gc/C8ATTzzB17/+dZqbm/N+gDIx5Tp4xm2xvYlNEovdxNmd4+CZgB3krV3t+OwNlPczkEaDZ0SkENg2/MkuZXUqwoNWFR2ZKyECXoOTDp1Ec6yBow6oUWG5iMgocx3yrrnmGj766CNOPPFEfL703S3L4vzzz+c//uM/8n6AIk758VKBlwoX6wsTOQRDG7oNnjHZy7p2Hz4248Xu9/k0eEZExqttto/fmGFWm2E+yCosn11XyorGKGcumkp1ifMP5UREZHi5DnmBQICWlha+853v8OqrrxIKhZg3bx7Tpk0bjuMTGTbZg2eqHK4vtPtYX7ibFI+VTWfB5jdJEO8VCttJAe4Hz3gwXAdDDZ4RkXwxbXgyU1i+1qoilfm3pTjg5bT5U2haEmVhtFLduSIiY5DrkNdp9uzZzJ49O5/HIjLmGRiEMpdfdq5ATRoeNpYewueMBH7b6nUfC5v2Acvsk71uj2eK7TV4RkRG2id2gDtTYe40w3yaVclzWLSS5liULyyYQmkw518fRERkBDj6V3rlypV85zvfoaSkhJUrVw74vdddd11eDmwi+z0fcD0vcC0nchDh0T4cGSIPBiUEKHGxvjCFRXs/AbCvcNiWKbbX4BkRyUXCNvi9VclqM8KTVjl25qxdRcjPGZnC8kMmlw/yKCIiMlY4Cnkvv/wyyWT6bMJLL73U76UZumQjPxYyiZUcyWTKgPSkSQNDgzsmEB8eyghS5rDYHnoPnnESEPMxeKZzKyEw4PpCvX5Fxp73rCLuMMPcbYbZnjUB+agDamiKRfn8oSosFxEZjxyFvD/84Q9df3788ceH61gko4ZiPkND19cDnTWxsbtGfeiX6Iktn4Nn2voJh53F9tmDZ5wwgJDLs4VBrS8UGRb7bA9rrCpaUhFesMu69teW7S8sn1ajwnIRkfHM1UX1qVSKoqIiXnnlFebOnTtcxzThnc89zCHClXyGdezkeTZQSREHEyZKeVfos7Excvg1eDcdXMMfeYqPmUElqziZWvSGPtHkMnimr2L7wbYOUtjQ9bVTnpzWF3oUDEX68fFe+PdEA79NVXUrLD/h4FqaYg0cf1AEnwrLRUQKgquQ5/P5mDZtGqZpDtfxTHg2Nm+yjUmUso6d/BX3kMRka6bV7WgaaOFsIP1L+qN8wDtsZyrlLGYKUygb8LENDP6Fx3iK9TQzFwubUO7zd2SCyVex/WBbIjN4Zi8J9rpYX+jLYX2hT+sLpYDtzhSW327W8sZrPiACQEN1MU2xKGctmqrCchGRAuT6t/t//dd/5aqrruJXv/oV1dXVw3FME1oCEwOopYT7eJu32MZ2/hcA77CdJ/kIgDgpbuZFVvE8kyhlI3uIUMJNfIHDmNTnYxsY7KSd59lAM3O5ks8MeCwpLFe/AP+ZzUynknIX68ik8OVSbJ/Cct1fmMIihUUrcVpdrC8M9FNg398WwqfBMzKm2Tb8j11KSyrCg1Y18czr1WvYnDJ3MucdPo0jZqqwXESkkLkOef/5n//Je++9x5QpU5g2bRolJd0v83vppZfydnAT0V4SJDCJUo6FjQeDh3mPgwkzmxpmZwb3v8CnfJs/8vfE+CbHsZN2/pr7+Ar38gqXYGVW6nX+14eHj9nNNFbhx8tbbON+3uZf+SwHUcPbbGcZB3Y7loECXudKwOxL4z7Df/NPHM2/8Nk+72NisYE9NFCR+1+QTAg+PJQTdPyBgZ0ZINN/TUWiz7OJVmZdYgKTXQ7XFwIUuewvDKnYXkbAVtvH3WaYFjPCOnv/2bmD6so4e/EUSra8zjmnz8fvd34mXkRExifXIe/000/XFM1htI19+PFSTYgjifI4H/Jv/IHZ1HA2c/giBwPwIO9QTxnf5DgAqghxGUu4kPt5ho85iijQfRhLlAoe4wIu5n7OZg7zqKWWEm7nL3ybJ0jwb12XdN7Dm/ySP3M3K7oeI4XFXhJUUtT1C2tnEN1BO9WEiA4Q4N5iG/O4gb9wKXOI8Awf8488wqOcT8jF5X8iPRkYBPASwEuli8Ezccx+A2BfAbFz8EwHKTpIscNhsb3bwTN+O4Bt24M+rohpwx+tClabER61KrsVli9fMIWmWJTDopWkUinWrHl9lI9WRERGiuuQ961vfWsYDkM6fcoeivBRQzGlBLiJ03ieT2jhdZq5i+/xOa7gCF5jS9dlmQlMAnjpIMUkSmkjwTp2cjVPEMJHAxUs40AWMIkdtFNCgFOYxbFMB+A2XukKhZ1eYiN/4tOugPcBO7mFl1jN62xiL4uZzL/yWU7iAAA+ZjcmNtUDDPDYxF5KCVCT+Z56yjiDgxXwZFQYGBTho8jFP4MjNXjGxOD9ze/ztr2dsqyzhiUqtpeMj60Ad5ph7jQjbMy6FHphQ7qw/NT5KiwXEZnIXL8DzJw5kxdeeIGamppu+3ft2sWiRYv44IMP8nZwE9E6dhGmmC208QIbmE0NhzOVw5nK62zlWT7hCtJnEoKZX+i8mSD2Ma348VBCgDDFzKeO19jC7fyFKkIsYBLvs6Pb4Ix2knzAzl5n4D5hD3OpBdLB82qeYB07uZcmSgnwU57nmzxOhGIWMpl17KIIH1WZsyidZwSzvc9OKijq6n6bRiXf4Og+v9fGxsrsVzWEjBUjNXimHQvLTp85b3exvtBtsX1Ig2fGlXhWYflTWYXllcV+zlw4laZYlIMm9T98S0REJg7XIe/DDz/sc7pmPB7nk08+yctBTWRvspUUFjY2P+E5klg0MpkUFi+wge/zOQCOZCp38gYvs5GFTAbgDl5nMmXMoJIygqzkyK7HNbEA2MheKgh2jcyPY/IJrRzB1K7vtbFZx04Oyqz/u5+3eYVN/C+O4lBq8WCwipP5PL/iXt5iIZN5jx1UEOzqaDOxMTJBzYsHDwbvsYNJlFKMHxub7/BHUlh8m+OJk+oKuGGKMTC6wmt/OtcbKgTKWJbL4Jl92NxXezCf2fwSSRKDrjFsI5HXYvvBgqH+Pzey3rWKaDEj/MasYUfWBwyfmRWmKRblpEPrCPp0FldERPZzHPLuv//+rj8//PDDVFTsP/NjmiaPPvooM2bMyO/RTUDf5URaiVOEDxObtbzPw7yPFw83cxqncxAAl3M4z7OBf2Qtx9DAq2zmNbbwO77E5D5qFDrPlH3KHqoJUZE5m2Zjs5V9XQNdDAzaSPAJrV2DWF5lE+vYyRU8zNd4EA8GdZTyLtuZlznb9xG7qMp63M6zA9m/dnzATqZlzhgmMHmY97sqHzawhy/xG8IUcwqz+BaPM5daruEEjmN6t7N9JlZXcByIlamK1yREGW/8hpdibzGTjDL8tjXo92cPnnGzxtDKrEt0U2wPEOpj8Ex/aw5L8FOkwTOu7bM9/NaspsUM82JWYXldeZAVjVHOWRyloaZ4FI9QRETGMsch74tf/CIAhmFwwQUXdLvN7/czffp0fvzjH+f14CaiAF7CpN+4v8BsvsDsPr8vQgk/4fP8ij/zDB8zmTKe4CvMyXQg9dT569Vm2phHbddZhSpCdJBiG/u6vvcVNrGBPRxIuiJjF3E+x0xu4FT8ePmY3XxCKx+ws+sM4CfsoYYQkyjl//Iq1/IUJQQoJ8j/5hTmEGEduzghsw4w/UtlO8szP18CEz8eXmUT5zCHp7iQf+L3fJ+nOYCqrstJ/8hHXM8L/J4PqKGYr7KYS4l1rasysUhgZp1t0C+WUviGMnjGyRTS/ZeRpgfPtJOinRTbXQyecdtfGMA74YKhbcNrdjGrzQj3mzXs7bwk32Nw4sG1NC+J8tkDVVguIiKDcxzyLCv9afKMGTN44YUXCIfDw3ZQ4syh1HJt5vLNwXT+stRKnJLMp+6dzuBg/puXmUElXjz8K48RJ8WsTMhbyky+xePsJcEBVFOdWd+XbSN7OIooJjZfYDZ1lPI+O3iVzZRknmsDrcykCkhXRWymres59pLgDbZyPafyV8wH4BIWcylreJvtRKngDl7nuzzJScxkDV/iTbbyI57Fj4fLOJwUFjfyJ27iRdazmwjFNHEol3E4kygdwt+0SOHJHjwz0MCkbBY27S77C+OY2EAbSdocDp6B9Fpjt8FwvA6e2W17udesYbUZ4U17/9m5aTXpwvKzF02ltlyF5SIi4pzrNXnr1q0bjuOQEfIif0t7j1+0ruEE/oGHuZgHWMRkLmYRl/IgNZkzil/kYG7gT1zFo/wDRzCDKjayh+f4hEtoxMBgFx1MopQgXorwZaZuHtD1HK2ZdUIzMiFvNx20Eu8KfTtpx8LmcOq77tNABTto7zoz8V2e5Him8y2Oo4QAS6hnLwke5F0u4DAe5j1+yDNcxWc4jum8xEbeZGuvn1dEcuPBoIQAJS7WF5o5FNsnsTCx2UOCPSQcP5fbwTPF+Eftcm7bhuesMlrMCGusKhKZ4wj4PCybO4mmWANHzKxWZZGIiOTEdci7/PLLmTVrFpdffnm3/T/72c947733WLVqVb6OTYZJz8qCKZTRwtldX6ewOJookzNnv6oJcRfncAUP83l+hQeDA6nheKZjZDryNrCnW39eT1tpI4VFfWYN3jb2YWJ1FaNvpx0vHqZS3nWfHbTTTpIZVJLEZD27uZkX+W9eJoiPSZTiweA1NtNOkgQmVRRxEgcwkyoOJsxm9nYNmRGRkefFQxnBrqm6TiR7XEbqZBsvg2e22H7uMmu4w4zwYVZh+cGTyjh3SQNfPKyeimLVyoiIyNC4Dnl33313tyEsnY466ii+973vKeQVAB8e5lHXbd80KrmHJgB20cF29nWt64uTYjkHcUDmrFxf1rELE6vr7OCn7KGcIDUUY2KxlTZKCXQLoJvYiw8PNRSzg3ZaiXM/zcyllg3sYR07Wc9uljGLCCWcxAH8mGc5kV9yNodwAYd11UCIyPjhx0sF3q5pvYOxsUnkEAxtcD14xgCK+hg803OroKjrMvGUDU9kCssfsyoxMyGxNOhj+WFTaI5FmVdfobN2IiKSN65D3vbt27tN1uxUXl7Otm3b8nJQMrZVUtRtuMNkyvg1Zw56v2OY1nW/N9hKBykg3fn3NtsJZV6O6ZmY6SBYSwmQvuRrNjW8ymZOZTbTqOxV4B6hhJf4KnfyOnfzJifzK/6Vz3IJjfn4sUVkjDIwCOIjiM/xmXsbOzNApr8ppMlet7Vniu2dDJ5poILPWUdwhxnmTjPM5qxLXBc1VNK8pIEvzJ9McUCF5SIikn+u311mzZrF7373O/7+7/++2/6HHnqImTNn5u3ApLB8jpl8jv2vj8s4nOWZOggPBgdQxcnM6ro9hcUrbO66/LOcIH/LIr7P0xxKhBOZSTtJ3mQbFQSZRx0b2UMdpZzDoZzCgfwHT/ITnuMSGvssXBeRicvIGuxS4/A+2YNn+gqB+0iyy07xou3nD9ZUvp2a33XfqmI/Zy1KF5YfWKfCchERGV6uQ97KlSv5+7//e7Zu3coJJ5wAwKOPPsqPf/xjXaopjnWWnkN6jeA/ZBW3d46Dv5KjuYiFAATx8fcsIYnFJTzIXhLUU8Z86riUGFto4/s8zUImcSIzaSPBS2zsqlYQERmq7MEzPctq3rZCrDbD3GOG2ZX5d8cw0oXlzbEGPjenVoXlIiIyYlz/BnzhhRcSj8f57ne/y3e+8x0Apk+fzg033MD555+f9wOUievQzHq6zrNwfrz8L47mHzmKj9nNBvYQxMuh1NJBikmU8n2e5u95iHKCmeqH40b3hxCRgtWWKSxfbUZ42d5f0zKpvIgVjVM5pzFKtFqF5SIiMvJyOs3xta99ja997Wts3bqVUChEaWnuHWTXX389P/zhD9m4cSOHHnooq1at4phjjhn0fk8//TTHHnssc+fO5ZVXXsn5+WXs63mZpQeDaVQyjcqufUX4+GeO4Z9Jv3a20EYwa3CDLtUUkXywbXjFLqHFjPCAWU1bppvP5zE48ZBammMNfHZ2BK9H/+aIiMjoySnkpVIpHn/8cd5//33OO+88AD799FPKy8tdBb6WlhauuOIKrr/+eo4++mhuuukmTjnlFN544w0aGhr6vd/u3bs5//zzOfHEE9m8eXMuP4IUuM6BLSIi+bDL9vIbM0yLGebtrMLyGeESmmJRzlxUT22ZCstFRGRscB3yPvroI04++WTWr19PPB5n6dKllJWV8YMf/ICOjg5uvPFGx4913XXXcdFFF3HxxRcDsGrVKh5++GFuuOEGrr322n7v99WvfpXzzjsPr9fLvffe6/ZHAOBleyOTCY2ZYlwRERlbrExh+Wozwu+yCsuDPg+nzJ1E85IGDp+hwnIRERl7XIe8r3/96zQ2NvLqq69SU7N/JtkZZ5zRFdacSCQSvPjii1x55ZXd9p900kk888wz/d7v1ltv5f333+dXv/oV11xzzaDPE4/Hicf3F+O2trYC8DvjfUJG70Xw2cW4IQJ9FuF2/9qHx1AwnMiSHk+3/4oMlV5To2uz5ec3Zg13psJ8bO8vcT94UhlNjfWcNn8yFaF0p2cqlRqtw3QlmUx2+6/IUOk1Jfmm15QzTv9+XIe8p556iqeffppAINBt/7Rp09iwYYPjx9m2bRumaVJX1710u66ujk2bNvV5n3fffZcrr7ySJ598Ep/P2aFfe+21XH311b32J6KHkQhC3IoTt+PErTi2bffxCDaQyGy9GYZBwAgQ8AQIeoLdNyPYe58niN/w65PfArR23rzRPgQpMHpNjRzThjd2Gjy3xeD1nQZ2Zh1vkddmcdjmyFqLqSU7Mbbv5Ok//GWUjzZ3a9euHe1DkAKj15Tkm15TA9u3b5+j73Md8izLwjTNXvs/+eQTysrcd//0DDu2bfcZgEzT5LzzzuPqq69m9uzZjh//qquuYuXKlV1ft7a2Eo1G+ecN5dRkBVXbtukgRRsJ2kllld/27EFKdO1rp/9PcFOZra2P2zxZ/Ux9nx309zqD6MerYDhGJT0e1s6bx9LXXsNvWaN9OFIA9JoaOR9ZAe5OhbnbrGGLvf89oXFaJWcvqueUuXUFUVieTCZZu3YtS5cuxe/3j/bhSAHQa0ryTa8pZzqvShyM63eupUuXsmrVKm6++WYgHdL27t3LN7/5TZYtW+b4ccLhMF6vt9dZuy1btvQ6uwewZ88e/vSnP/Hyyy93FbFbloVt2/h8Ph555JGu3r5swWCQYDDYa7/fsnr98hTASzkhxz9Dz2LcnqW4fW0JTMCmnTjtxNnu8Ll8eCi2e4e/gTaf1heOqL5eUyJDodfU8OiwDR62qmgxIzxjlXftrykJcNbiqaxojDKrNvep0WOZ3+/XL0+SV3pNSb7pNTUwp383rkPeT37yE44//njmzJlDR0cH5513Hu+++y7hcJjbb7/d8eMEAgEWL17M2rVrOeOMM7r2r127ltNPP73X95eXl/Paa69123f99dfz2GOPcddddzFjxgy3P8qQDVSM258UlqMwmL2lsEhh0UqcVuKDP0lGIGt9oZMthE+DZ0SkYL1lhVhtRrjHrGF3VmH5Zw+M0ByLcuIhdQR8+jdQRETGP9chb8qUKbzyyivcfvvtvPTSS1iWxUUXXcSXvvQlQiHnZ8EAVq5cyZe//GUaGxs58sgjufnmm1m/fj2XXHIJkL7UcsOGDfzyl7/E4/Ewd+7cbvevra2lqKio1/6xzIeHcoKU0/vsYl9sbJI5BEMLmwQmCUx20eH4+IrwuQ6G6qATkbFqr+3hgUxh+atZheX1lSHOXjyVFbEo9ZXu3rtERETGupwWGoRCIS688EIuvPDCIT15U1MT27dv59vf/jYbN25k7ty5rFmzhmnTpgGwceNG1q9fP6TnGO8MDAJ4CeClEmcdTDY2cUxXobCdJDbQQYoOUuyg3eHx0ecawoG2IF4FQxEZNrYNL2UKy39rVrMvU1ju9xosnVNHU6yBz8wKq7BcREQKlqOQd//99zt+wOXLl7s6gEsvvZRLL720z9tuu+22Ae/7rW99i29961uunm8iMDAowkcRPqodrjG0SA+ecRMMO0hhQ9fXTvUcPONk8+NRMBSRAe2wffzGrKHFjPCuvf/fvpmREs6NNXDGonrCpc6uohARERnPHIW8L37xi44ezDCMPidvytiXHbycMrG6TSJ1OnjGwmYvCfb2U0vRFx8e18FQg2dECp9lw9NWOavNCGutyq7C8iK/h1PnTaEpFiU2vUrTiUVEZEJxFPIsTXeTPnjxUEqAUgKDf3OGBs+ISD5stP3cZYZpMSN8klVYPq++gqZYlOWHTaG8SNPZRERkYhr/5T8yrgx18EzPHsOeW+ftIzl4RkRGRtI2eMyqoMWM8LhVgZW5hLusyMcZC+tZ0Rhlbn3FKB+liIjI6HP8G+qyZcu4/fbbqahIv4F+97vf5e/+7u+orKwEYPv27RxzzDG88cYbw3KgMjHlc/DMQAFxKINnAnaAdVs/YYu9hTIHIVGDZ0TcWWcFaTEj3GWG2ZZ1SfmSGdU0NUZZNm8yoYB3FI9QRERkbHEc8h5++GHi8f2Xyn3/+9/n3HPP7Qp5qVSKt99+O+8HKOJWroNn2kn2GwL7CojZg2daU62sZzde7EGfa6DBMyX9Dp7RL7AysXTYBr+zqrg9FeF5e39hebg0XVje1BhlZqQwC8tFRESGynHIs217wK9FxrPsYnunOgfP7CbFg9WzOHzzn0mS6BYMe4bFJFZOg2f8LgfPhDR4RsapN6wQLZnC8tbMW5THgGNnR2iKNXDiIbX4vXpti4iIDEQLikRy1Dl4JmgUUResY45Ri98eeEhR0mV/4T6SmJl1ibuJs9vF4Jmg68Ez6aIKkZG2x/Zwf6b64M92Sdf++soQTbEoZy+eyhQVlouIiDjmOOQZhtFrBLVGUou448dLBV4qXKwvTLgutk8PnoljEsdkp4vBMyGXg2eK8Gl9oeQkXVheymozzG/NatqzCstPmjOJpliUz8wK41FhuYiIiGuuLtf8yle+QjCYnorY0dHBJZdcQklJ+lPX7PV6IpIfBgZBfATxUeVwfaGdQ7F9OykA2knRTortLgbPuO0vDGjwzIS23fZxj1nDajPCe1mF5bNqS2mORTljYT01KiwXEREZEsch74ILLuj29V/91V/1+p7zzz9/6EckIkNiYBDKXH5Z4/A+nYNnnATCtsw00jgmNtCW2eeUd4DBM/1tGjwzvlk2PGWV02JGeMSqJJlZLxrye/nC/Mk0xaIsnqbCchERkXxxHPJuvfXW4TwOERlFuQyeSWH1GQw7Q2Bf+1NYmNjsIcGeYRw8U4xfxfZjwKd2gDvNMHekwmzI6sacPzVTWL5gCmUqLBcREck7DV4RkZz48FBGkDKHxfYACcwBzxj2FRCHY/BMCYE+1xdq8MzQJW2DRzOF5U9kFZaXdxaWx6IcOkWF5SIiIsNJIU9ERkxnsX0ug2f6O0OYfVvnf21wPXjGAEIuzxaq2H6/DzKF5Xf3KCw/YmY1zbEGTp47iSK/LrsVEREZCQp5IjJm5TJ4xsph8Ex2sf0+F+sLByq27399oadggmG77eEhq4rVqTD/k1VYHikLcvbiqaxojDIjXDLAI4iIiMhwUMgTkYKSHbyc6iy2dz6NND14Jpdie5+D9YV+O8DO5E5a7Q4q8I25Yvu/WMW0mGHuNWvYk1VYftxBtTTHohx/sArLRURERpNCnohMeJ3F9qUuB890Bj4naww7B8+ksGglTusA6wtNDN7dton32YwXm4DrYntf3gfP7La93G9W02JG+EtWYXm0OkRTY5SzF0eZVOHsMlwREREZXgp5IiI58OGhnCDlLgbPJLOK7dsGCId7SPGJN5QZBJNel5jAZJeLYvuiHIrtew6esW14wS5ldSrCGquKjkyVRcBrcNKhkzh3SQNHzqxRYbmIiMgYo5AnIjJC/HipcDB4Jml4WFO7gFM+fQXLxdrCzjOLNtBBig5S7HBRbB/CTxVFfNE+nN+YYVrMMB9kFZYfWFtK85IGzlhYT3WJ87OeIiIiMrIU8kRExijDMCiyfRTho3qYB8+YwOu2h/esSr6fPIxU5qxeccDLafOn0LQkysJopQrLRURExgGFPBGRAuJ28MwndoCWVDV3mGV8CthY+DA4LFpJcyzKFxZMoTSotwoREZHxRO/cIiITTMI2+L1VyWozwpNWOXbmrF11yM8ZC+tpikU5ZHL5II8iIiIiY5VCnojIBPGeVZQpLK9hR9aZvqMOqKEpFuXzh6qwXEREpBAo5ImIFLB9tocHzSpazAh/ssu69teWBTmnMV1YPq1GheUiIiKFRCFPRKTA2Db8xS5mtRnhfrO6q7Dc6zE4/qAIzbEGjjsogk+F5SIiIgVJIU9EpEDstr3cZ9aw2gzzRlZheUN1MU2xKGcvnkpduQrLRURECp1CnojIOGbb8LxdRksqzBqrmjjps3MBr8HJcyfTHItyhArLRUREJhSFPBGRcWir7eNuM0yLGWGdvf/s3EF1ZTTFopyxsJ4qFZaLiIhMSAp5IiLjhGnDH60KVpsRHrUquwrLSwJelh82haZYAwumVqiwXEREZIJTyBMRGeM+tgLcaYa504ywkf1n5xY1VNIca+DU+ZMpUWG5iIiIZOi3AhGRMShlwUOpSu5MhXkqq7C8qtjPGQun0rwkyuy6skEeRURERCaiCRvy1tjvMokiivH3uRXhw0CXPInIyHrXKuL2VC0tL3ppSx3Qtf+YA8M0xaIsnVNH0KfCchEREenfhA15r7GJ9+n/FyUD+g2A/W0BvAqGIuJam+3hQbOaFjPMi1mF5XXlQZoao5zTGCVaXTyKRygiIiLjyYQNeZ9lOn5s9pHstcUxsYE2krSRdPyYXgzXwdA/QNAUkcJl2/Bnu4TVZpgHzBr2Zv4t8HoMTjgownR7IyubjyFUFBzlIxUREZHxZsKGvCONKDX0PV7cxOoz/A20JbEwsdlDgj0kHB+HH4/rYOjN9GCJyPizy/Zyr1nDajPCW/b+s3PTa4ppijVw1uJ6qoq8rFnzKT6v/r8uIiIi7k3YkDcQLx7KCFKG80/Qk5i0kaS9nxDY120mNkksdhNnN3HHzxXE6yoUhvDj0WWkIqPGsuE5q4wWM8JDVhWJzAc1QZ+HZfMm0xSLcviM6q7qg2TS+RUEIiIiIj0p5OWJHy+VeKmkaPBvBmxsEpiOzhJmB0QbiGMSx2QnHY6eywCK8LkKhho8IzJ0W2w/d5ph7jDDfJRVWH7wpDLOXdLAFw+rp6LYP4pHKCIiIoVIIW+UGBgE8RHERxUhR/exsekg5egsYefWTgobaCdFOym20+7ouTwYhLKCYQkBB+sLPQqGMuGlbHgiU1j+mFWJmfn/RGnQx/LDptAcizKvXoXlIiIiMnwU8sYRA4NQ5vLLGof3sbD7DYD9BcQ4Jha268EzvhzWF/q0vlAKxHoryB1mmDvNMJuz1vs2TquiKRbl1PmTKQ7on1wREREZfvqNo8B5MCghQEk/Q2b6ksKifZAzhNkBsY0kKSxSWLQSp9XF+sKA6/WFPg2ekTGjwzZ4xKqixQzztFXRtb+6JMCZC+tpikU5UIXlIiIiMsIU8qQXX46DZ5xeQtp5u5VZl5jAZJfD9YWQ2/pCDZ6RfHrbCtFihvmNGWZX5p9Rw4DPzArTHGvgc3NqVVguIiIio0YhT/LCj5cKvFTkYfDMQAHRBjpI0UGKHQ7XFxpAyOVlpEEV20sPbbaH35rVrDYjvGyXdu2fXFHEOY1Rzlk8VYXlIiIiMiYo5MmoyHXwTHsfg2cG2joyg2c6v3bK46LY3m8HSFpJbNvO8W9DxirbhlfsElrMCA+Y1bRlCst9HoPPHVJH05Ionz0wgtejDwRERERk7FDIk3HDyApeTplYroNhIjN4Zi8J9jootjcxeHfzOl5nK2UuLyXV4Jmxaaft5R4zTIsZ5u2swvIZ4RKaYlHOXFRPbZmzs9YiIiIiI00hTwqaFw+lBCh1OXjGTSjcQ6rrfho8M35ZNjxrlbHajPBwj8LyUzOF5UuyCstFRERExiqFPJEefHgoJ0i5w8EzCQzurzuUYze/SJJEH2sME32eTRypwTMhFdsPaLPt567MWbv1WYXlcyaXc+6SKMsPq6cipMJyERERGT8U8kSGyDAM/B4/lUYIv+0sGNrYxDODZ/oLgT1va9fgmbxJ2fAHq5IWM8xjViVW5mctC/o4feEUmmMNzK2vGORRRERERMYmhTyRUWBgUISPInxUOxw8Y2HTkRUGBwqHnbfHMYd98ExJ5wAaxn5lwEdWkBYzzF1mmC1Zl/DGplfRHGtg2bzJhAJj/+cQERERGYhCnsg44cnj4JmBAqLbwTOd/HhcXkY6MoNnOmyDh60qVpsRnrXKu/bXlAQ4e/FUVsSiHBApHeARRERERMYXhTyRApbL4JmexfZONhObJBa7ibPbxeCZoOvBM37HxfZvWiFazAj3mDXsziosP3Z2hOZYlBMOriPg0xAbERERKTwKeSLSTT6L7fvb2klhZdYlxjHZ6WLwTGiAwTOH2FF+Z9bSYoZ5NauwfEpFEStiUc5pjFJf6ezyWBEREZHxSiFPRIYk12L7Dpf9he2Zqop2UrSTYntm8IwNbLQD/MUuZkfiNDoyZy39XoOlc+poijXwmVlhFZaLiIjIhKGQJyIjzsAglLn8ssbhfSxs2rNC36e2xRqzhIesYj7GxiZFMQEOiJTQHGvgjEX1hEudTTsVERERKSQKeSIyLngwCNkBXrLCrDYjPGJVkswMbqnyezh13hTOXRJl8bQqFZaLiIjIhKaQJyJj3kbbz51mmDvMCJ9kdRHOq6+geUmU0xZMobxIheUiIiIioJAnImNU0jZ4zKpgtRnhCatif2F5kY8zFtazojGqwnIRERGRPijkiciYss4K0mJGuMsMsy2rE3DJjGrOXRLllLmTKfKrsFxERESkPwp5IjLqOmyDh6xqVqfCPG/vLywPlwbTheWNU5mpwnIRERERRxTyRGTUvJFVWN6a+efIkyksb4o1cOIhtfi9KiwXERERcUMhT0RG1B7bw/1mDS1mhD/bJV376ytDrGiMck7jVKaosFxEREQkZwp5IjLsbBtetEtZbYZ50KymnfSaOr/X4KRDJ9Eci3L0AWE8KiwXERERGTKFPBEZNtttH78xa1htRnjf3n92blZtKc2xKGcsrKdGheUiIiIieaWQJyJ5ZdnwlFVOS4/C8pDfyxfmT6YppsJyERERkeGkkCciefGpHUgXlqfCbGD/2bkFUytoijVw2oLJlKmwXERERGTYKeSJSM4StsFjViWrzTBPWBXYmcLy8iIfZy6ayorGKHOmlA/yKCIiIiKSTwp5IuLa+1YRd5hh7u5RWH7kzBqaYlFOnjtJheUiIiIio0QhT0Qcabc9rLGqaElF+B+7rGt/pCzIOYvTZ+2mh0sGeAQRERERGQkKeSIyoL9Yxaw2I9xnVrMnq7D8hINraYo1cPxBEXwqLBcREREZM0b9N7Prr7+eGTNmUFRUxOLFi3nyySf7/d7f/OY3LF26lEgkQnl5OUceeSQPP/zwCB6tyMSw2/byf1MRTo3P4QuJQ/mVWcsefESrQ/zjSbN55soT+T8XxFg6p04BT0RERGSMGdUzeS0tLVxxxRVcf/31HH300dx0002ccsopvPHGGzQ0NPT6/j/+8Y8sXbqU//iP/6CyspJbb72V0047jeeff56FCxeOwk8gUjhsG16wS1mdirDGqqYj8xlQwGvw+bmTaY5FOXJmjQrLRURERMa4UQ151113HRdddBEXX3wxAKtWreLhhx/mhhtu4Nprr+31/atWrer29X/8x39w33338cADD/Qb8uLxOPF4vOvr1tZWAH7OS0SMECH8FGdtPb8uxo8frzq9pF9Jj6fbf8ebbbaPe1I13JkKs84u6tp/YG0JKxqncvqCyVQVBwAwzRSmOVpHOnEkk8lu/xUZKr2mJN/0mpJ802vKGad/P6MW8hKJBC+++CJXXnllt/0nnXQSzzzzjKPHsCyLPXv2UF1d3e/3XHvttVx99dW99r8QKcEf6pwKaGa2jj4fw2t4CXqC3baAJ0DQCPbaX+QpIuAJ4DU0WXCiWTtv3mgfgmOWDW/tMnh2i8FfdhpYdvpDjIDHZlHY5shai2mluzF27ubZx18f5aOduNauXTvahyAFRq8pyTe9piTf9Joa2L59+xx936iFvG3btmGaJnV1dd3219XVsWnTJkeP8eMf/5i2tjZWrFjR7/dcddVVrFy5suvr1tZWotEo/2vbJAIBaCfJvqytjSQdXX9OYGL3+9jxzNaXIN6ss4OBXmcHe59B9OExxueZoIku6fGwdt48lr72Gn7LGu3DGdAGK8DdqRruMsNstANd+xdMrWDF4nqWzZtEaVDzmEZbMplk7dq1LF26FL9fBfIydHpNSb7pNSX5pteUM51XJQ5m1H+b63kZpG3bji6NvP322/nWt77FfffdR21tbb/fFwwGCQaDvfbPsiuoyfolty82NkmsHiEwQTupXvs6/9xOCgubFClaSdFK+6A/S6eQ7esVBgfaivBhoMtIxwq/ZY3JkJewDX5vVbLajPCkVd5VWF5Z7OeMhfU0xaIcPEmF5WOR3+/XG53klV5Tkm96TUm+6TU1MKd/N6MW8sLhMF6vt9dZuy1btvQ6u9dTS0sLF110EXfeeSef+9znhu0YDQwCeAngpZKiwe9AOhh29AiBg23tpABoJ0U7KbY7DIYGuAqFxfgJ4FUwnCDes4poMSPcbdawI6uw/KgD0oXlnz9UheUiIiIihWjUQl4gEGDx4sWsXbuWM844o2v/2rVrOf300/u93+23386FF17I7bffzqmnnjoSh+qKgUEoczlmjcP7WNi9LhsdbItnLiRty1xi6pQXw3Uw9KMgMF7ssz08aFbRYkb4U1ZheW1ZkHMa04Xl02pUWC4iIiJSyEb1cs2VK1fy5S9/mcbGRo488khuvvlm1q9fzyWXXAKk19Nt2LCBX/7yl0A64J1//vn89Kc/5Ygjjug6CxgKhaioqBi1n2OoPBiUEKCEgS8fzWZmLiNtJ9XtctGBtiQWJjZ7SLCHhOPn8uNxHQy9o1/BOGHYNvzFLuZ2M8L9Zg17M6Hc6zE4/qBammNRjlNhuYiIiMiEMaohr6mpie3bt/Ptb3+bjRs3MnfuXNasWcO0adMA2LhxI+vXr+/6/ptuuolUKsXf/d3f8Xd/93dd+y+44AJuu+22kT78UeXFQxlByggCzs7MJDFdnS3cRxIzsy5xN3F29ztmprfswTNOthB+PLqM1JXdtpf7zBpWm2HesPe/Bhqqi2mKRTl78VTqyp1dZiwiIiIihWPUB69ceumlXHrppX3e1jO4Pf7448N/QAXMj5cKvFS4WF+YcBkMOwfPxDGJY7Kzn1qKngygCA2eGYxtw/N2GS2pMGusauKdheU+D6fMnURTLMoRM1RYLiIiIjKRjXrIk7HLwCCIjyA+qgg5uk+ug2ds3A+e8WAQchkMx+vgmS22j7vNMHeYkW6F5QdPKqM5FuWLC+upLHZ+ua+IiIiIFC6FPMmrkRw8Y2EX9OAZ04Y/WhWsNiM8alWSyoTTkoCX5YfV0xyLMn9qhaPKERERERGZOBTyZNTlMngmhdVnMGwj2e/+VI6DZwKDrC/02wE2xzezxd5LBT5C+IY0eOZjK8CdmbN2m7L+ThZPq6IpFuXUeZMpUWG5iIiIiPRDvynKuOTrNnjGmQTmgGcM+wqIZmZdYgKTXf2sLzQxeHfHNj5mM15sIL2+cLBLSUsIdP3ZY/tZa1XTYoZ5yto/Kbaq2M+Zi6bSHItyYF1Zn88vIiIiIpJNIU8mjM5i+1wGz/QVANsy/91Dis2+corZQYIENtBBig5Sgw6e2W77+ItdzNbEF9lLXdf+Yw4M0xSLsnROHUGfegpFRERExDmFPJF+OB08kzQ81EYWsOzTV/Ha5qCDZ3baKZ6zAjxnB1gP2CQop4Qp5UWsaJzKOY1RotXFI/eDioiIiEhBUcgTySNP1mCXbLYNr9oltJhh7jdraMsMc6nxGBx3UA3nLZnOsQfV4lX1gYiIiIgMkUKeyDDaZXu5x6yhxYzwlr3/7Nz0mmKaYg2ctbie2jIVlouIiIhI/ijkieSZZcNzVhktZoSHrCoSmUmbQZ+HZfMm0xSLcviMalUfiIiIiMiwUMgTyZMtto97UzXcYYb5KKuw/JDJ5enC8sPqqSj2D/AIIiIiIiJDp5AnMgQpGx5NVfDztzysbJ+PmSksLw36WH7YFJpjUebVq7BcREREREaOQp5IDtZbQVrMMHeaYbYQoLNbvbGzsHz+ZIoD+r+XiIiIiIw8/RYq4lCHbfCIVUWLGebpHoXlCyri/NPZn+GQ+qpRPEIREREREYU8kUG9bYVYbYa5xwyzK/N/GcOAYw6M0ByLcuysan7/yO+YVVs6ykcqIiIiIqKQJ9KnNtvDA2Y1q80Ir9j7w9vkiiLOaYxyzuKpXYXlyWRytA5TRERERKQXhTyRDNuGV+wSWswID5jVXYXlPo/B5w6poykW5bOzIyosFxEREZExTSFPJrydtpd7zDAtZpi3swrLZ4ZLaIpFOXPRVCJlwVE8QhERERER5xTyZEKybHjWKmO1GeHhrMLyIn+6sLw51kBsepWqD0RERERk3FHIkwlls+3nrsxZu/VZheVz68tpaoyy/LB6KkIqLBcRERGR8UshTwpeyoY/WJW0mGEesyqxMoXlZUU+vnhYPU2xKHPrKwZ5FBERERGR8UEhTwrWh1aQO8wwd3UWlmcsmVFNcyzKKXMnEwp4R/EIRURERETyTyFPCkqHbfCwVcVqM8KzVnnX/pqSAGcvnsqKWJQDIuqzExEREZHCpZAnBeFNK0SLGeEes4bdWYXlx85OF5afcHAdAZ9nlI9SRERERGT4KeTJuLXH9vCAWUOLGebVrMLy+soQKxqjnN04lfrK0CgeoYiIiIjIyFPIk3HFtuElu5TVZpjfmtW0ZwrL/V6DpXPqWNEY5ZgDVVguIiIiIhOXQp6MC9ttH/eYNaw2I7xn7z87d0CkhOZYA2csqidcqsJyERERERGFPBmzLBuetspZbUZ4xKokmVVY/oX5U2iORVk8TYXlIiIiIiLZFPJkzNlo+7nTDNOSirCB/Wfn5tVX0BSLsvywKZQXqbBcRERERKQvCnkyJiRtg0etClrMCE9YFd0Ky89YWM+KRhWWi4iIiIg4MWFDXtK2RvsQBFhnBVltRrjbDLON/WfnDp9RTfOSdGF5kV+F5SIiIiIiTk3YkHcdT1NJgGL8jrcQPryoa22oOmyDh6xqVqfCPG/vLywPlwbTheWNU5mpwnIRERERkZxM2JAHkMAkgckuOhzfpwifq2BYhA8PGgwC8HpWYfmezEvPY8BxB9XSFItywsG1+L0K0SIiIiIiQzFhQ97XOZIQBvtIOtraSWIDHaToIMUO2h09jwGEXITCYvwE8WIUSDBstb3cb1bTYkZ4zS7p2l9fGaIpFuXsxVOZosJyEREREZG8mbAhr8jwUU2AapwFDAubDlIDBsE2ErRnfU8HKWzo+topD4arUFiMHz+eMRMMbRtetEu5PRXhQauKjqzC8pMOnURzLMrRB4TxqLBcRERERCTvJmzIcys7eDllYnULfT1DYF8BMYGJhc1eEuwl4fi5fHhcB0NfntcXbssqLH8/q7D8wNpSmmJRzlw0leqSQF6fU0REREREulPIG0ZePJQSoBTnwSaJ2W8Q7C8cmtiksGglTitxx88VwOs6GPZcX2ja8JRVTosZYW1WYXnI7+UL8yfTvKSBRQ2VKiwXERERERkhCnljjB8vfryUZ5WAD8TGJuEyGO4jiZW5X66DZ+J2EX+2angzFWObPbXr9gXRSppjUb4wfzJlKiwXERERERlxCnnjnIFBEB9BfFRS5Og+NjZxTMeXkHYOnkkBr9k+/mKF+MguAvZRgkE45OeMhfU0xaIcMrl8sKcXEREREZFhpJA3ARkYFOGjCJ+jwTPvWUWsNqu516xgGya2EaeEOIdOLeKvlxzBGYcdqMJyEREREZExQiFP+tRue1hjVbE6FeEFu6xr/6SyIOcsnsqKxijTwyUDPIKIiIiIiIwGhTzp5i9WMavNCPeZ1d0Ky084uJamWAPHHxTBp8JyEREREZExSyFP2J0pLF9tRng9q7A8Wh2iqTHK2YujTKpwtt5PRERERERGl0LeBGXb8D92KS2pCA9a1cQz1QcBr8Hn506mORblyJk1KiwXERERERlnFPImmK22j9+YYVrMMB9kFZbPriulOdbAGQvrqVJhuYiIiIjIuKWQNwGYNvzRqqDFDPN7q4pUptC8OODltPlTaF4S5bCoCstFRERERAqBQl4B+8QOcEcqzJ1mhI3sPzu3sCFdWH7q/CmUBvUSEBEREREpJPoNv8AkbIO1ViWrzQhPWeXYmbN2lcV+zlw4laZYlIMmlQ3yKCIiIiIiMl4p5BWI96wiWswId5s17MDftf/oWTU0xRo4aU6dCstFRERERCYAhbxxbJ/t4UGzihYzwp+yCsvryoOcszjKisYoDTXFo3iEIiIiIiIy0hTyxhnbhtfsdGH5/WYNe0mfnfN6DE44uJbmWJRjZ6uwXERERERkolLIGyd2217uNWtYbUZ4095/dm5aTTFNsShnL5pKbbkKy0VEREREJjqFvDHMtuF5u4yWVJg12YXlPg/L5k6iKdbA4TOqVVguIiIiIiJdFPLGoC22n7vMGu4wI3xo7z87d/CkMs5d0sAXD6unotg/wCOIiIiIiMhEpZA3RqRseMKqYLUZ4TGrEjNTfVAS8LL8sHrOXRJlXn2FCstFRERERGRACnmj7GMrwB1mhDvNMJuyCssXT6uiKRbl1HmTKVFhuYiIiIiIOKT0MAritsEjVhUtZpinrIqu/VXFfs5cNJXmWJQD61RYLiIiIiIi7inkjaB3rCJWmxF+Y4bZlfVXf8yBYZpiUZbOqSPoU2G5iIiIiIjkTiFvmLXZHn5rVrPajPCyXdq1f1J5ESsap3JOY5RotQrLRUREREQkPxTyhoFtw6t2CavNCA+Y1bRlCst9mcLyc5c08NnZEbyqPhARERERkTxTyMujXbaXe8waWswIb2UVls8Il9AUi3Lmonpqy1RYLiIiIiIiw0chb4gsG56zylhtRvidVUUiU1ge9HlYNm8yTbEoh8+oVvWBiIiIiIiMCIW8HG22/dxlhmkxw6zPKiw/ZHI55y6JcvoCFZaLiIiIiMjIU8hzIWXD45nC8j9kFZaXBn2cftgUmmMNzK0v11k7EREREREZNQp5DnxkBbnDDHOnGWZLVmF5bHoVTbEGls2bRHFAf5UiIiIiIjL6lEz60WEbPGxV0WJGeMYq79pfUxLgzEX1NMUamFVbOsAjiIiIiIiIjDyFvB7eskKsNiPcY9awO/PXYxhwzIERzo1FOfGQOgI+zygfpYiIiIiISN9GPa1cf/31zJgxg6KiIhYvXsyTTz454Pc/8cQTLF68mKKiImbOnMmNN9445GPYa3tYnQpzevwQTk7M5Tazjt34mFJRxNdPPJAn/9fx/PLCJZwyb7ICnoiIiIiIjGmjeiavpaWFK664guuvv56jjz6am266iVNOOYU33niDhoaGXt+/bt06li1bxt/8zd/wq1/9iqeffppLL72USCTCWWed5eq5bRteskpoyRSW78sqLF86p44VsSifPVCF5SIiIiIiMr6Masi77rrruOiii7j44osBWLVqFQ8//DA33HAD1157ba/vv/HGG2loaGDVqlUAHHLIIfzpT3/iRz/6keuQ9/+3d+9BUZX/H8DfCwsiq+CNEIVATEPGEAFvOMjX4jLhWOr0g/Jumu54QWXUKG+Q5S0l07xNopii4SUdYkBhTFHQCSFIAyJCFP2FGpC4gIrA8/vDr/tr3Q1ZXHaX7f2a2RnOc55z9vPAZ5bz2edc/ueRG0pFd+Wyq50M7w5xwgQvR/To1KH1gyIiIiIiIjIggxV59fX1yMnJQWRkpEp7UFAQLl68qHGbS5cuISgoSKUtODgYsbGxePz4MSws1J9L9+jRIzx69Ei5XF1dDQAofijQUfoQgW52GD+4NwY72T559MGjGlQ+qnnR4dG/yOPHj1FXV4fKykqNOUikLeYU6RpzinSNOUW6xpxqGYVCAQAQQjTbz2BFXkVFBRobG2Fvb6/Sbm9vj9u3b2vc5vbt2xr7NzQ0oKKiAg4ODmrbrFu3DtHR0Wrt/7tzOgDgdwA7WzcEIiIiIiIivVMoFLC1tf3H9Qa/u+azDw4XQjT7MHFN/TW1P/XRRx8hIiJCuXzv3j04OzujrKys2V8MUUvdv38fTk5OuHnzJmxsbJ6/AdFzMKdI15hTpGvMKdI15lTLCCGgUCjQq1evZvsZrMjr0aMHzM3N1Wbt7t69qzZb91TPnj019pdKpejevbvGbTp06IAOHdSvsbO1tWUCkU7Z2Ngwp0inmFOka8wp0jXmFOkac+r5WjJRZbDnAVhaWsLb2xtpaWkq7WlpafD19dW4zYgRI9T6p6amwsfHh+fuEhERERERwcDPyYuIiMCePXuwd+9eFBYWYvHixSgrK4NcLgfw5FTLqVOnKvvL5XLcuHEDERERKCwsxN69exEbG4slS5YYaghERERERERGxaDX5IWFhaGyshKffPIJysvLMXDgQCQnJ8PZ2RkAUF5ejrKyMmX/Pn36IDk5GYsXL8b27dvRq1cvbN26VavHJ3To0AGrV6/WeAonUWswp0jXmFOka8wp0jXmFOkac0q3JOJ5998kIiIiIiKidsOgp2sSERERERGRbrHIIyIiIiIiMiEs8oiIiIiIiEwIizwiIiIiIiITYpJF3o4dO9CnTx9YWVnB29sbFy5caLZ/eno6vL29YWVlBVdXV+zatUtPkVJ7oU1OfffddwgMDISdnR1sbGwwYsQInD59Wo/RUnug7efUU5mZmZBKpfD09GzbAKnd0TanHj16hOXLl8PZ2RkdOnRA3759sXfvXj1FS+2BtjkVHx+PQYMGwdraGg4ODpgxYwYqKyv1FC0Zs/Pnz2Ps2LHo1asXJBIJTp48+dxteHz+YkyuyEtISMCiRYuwfPly5Obmws/PD2+++abKoxj+rrS0FCEhIfDz80Nubi4+/vhjhIeH4/jx43qOnIyVtjl1/vx5BAYGIjk5GTk5ORg9ejTGjh2L3NxcPUdOxkrbnHqquroaU6dOxRtvvKGnSKm9aE1OhYaG4syZM4iNjUVRUREOHz4MNzc3PUZNxkzbnMrIyMDUqVMxc+ZM5Ofn4+jRo7h8+TJmzZql58jJGNXW1mLQoEH46quvWtSfx+c6IEzM0KFDhVwuV2lzc3MTkZGRGvsvW7ZMuLm5qbTNmTNHDB8+vM1ipPZF25zSxN3dXURHR+s6NGqnWptTYWFhYsWKFWL16tVi0KBBbRghtTfa5lRKSoqwtbUVlZWV+giP2iFtc+rzzz8Xrq6uKm1bt24Vjo6ObRYjtU8AxIkTJ5rtw+PzF2dSM3n19fXIyclBUFCQSntQUBAuXryocZtLly6p9Q8ODkZ2djYeP37cZrFS+9CanHpWU1MTFAoFunXr1hYhUjvT2pzat28fSkpKsHr16rYOkdqZ1uRUYmIifHx8sHHjRvTu3Rv9+/fHkiVL8ODBA32ETEauNTnl6+uLW7duITk5GUII3LlzB8eOHcOYMWP0ETKZGB6fvzipoQPQpYqKCjQ2NsLe3l6l3d7eHrdv39a4ze3btzX2b2hoQEVFBRwcHNosXjJ+rcmpZ23evBm1tbUIDQ1tixCpnWlNThUXFyMyMhIXLlyAVGpSH9ukA63JqWvXriEjIwNWVlY4ceIEKioqMHfuXFRVVfG6PGpVTvn6+iI+Ph5hYWF4+PAhGhoa8NZbb2Hbtm36CJlMDI/PX5xJzeQ9JZFIVJaFEGptz+uvqZ3+vbTNqacOHz6MqKgoJCQk4KWXXmqr8KgdamlONTY2YuLEiYiOjkb//v31FR61Q9p8TjU1NUEikSA+Ph5Dhw5FSEgIYmJiEBcXx9k8UtImpwoKChAeHo5Vq1YhJycHp06dQmlpKeRyuT5CJRPE4/MXY1JfCffo0QPm5uZq3zLdvXtX7duAp3r27Kmxv1QqRffu3dssVmofWpNTTyUkJGDmzJk4evQoAgIC2jJMake0zSmFQoHs7Gzk5uZi/vz5AJ4coAshIJVKkZqaitdff10vsZNxas3nlIODA3r37g1bW1tl24ABAyCEwK1bt9CvX782jZmMW2tyat26dRg5ciSWLl0KAPDw8IBMJoOfnx8+/fRTzryQVnh8/uJMaibP0tIS3t7eSEtLU2lPS0uDr6+vxm1GjBih1j81NRU+Pj6wsLBos1ipfWhNTgFPZvCmT5+OQ4cO8XoEUqFtTtnY2ODq1avIy8tTvuRyOV599VXk5eVh2LBh+gqdjFRrPqdGjhyJP/74AzU1Ncq23377DWZmZnB0dGzTeMn4tSan6urqYGamelhpbm4O4P9nYIhaisfnOmCgG760mW+//VZYWFiI2NhYUVBQIBYtWiRkMpm4fv26EEKIyMhIMWXKFGX/a9euCWtra7F48WJRUFAgYmNjhYWFhTh27JihhkBGRtucOnTokJBKpWL79u2ivLxc+bp3756hhkBGRtucehbvrknP0janFAqFcHR0FO+8847Iz88X6enpol+/fmLWrFmGGgIZGW1zat++fUIqlYodO3aIkpISkZGRIXx8fMTQoUMNNQQyIgqFQuTm5orc3FwBQMTExIjc3Fxx48YNIQSPz9uCyRV5Qgixfft24ezsLCwtLYWXl5dIT09Xrps2bZrw9/dX6X/u3DkxePBgYWlpKVxcXMTOnTv1HDEZO21yyt/fXwBQe02bNk3/gZPR0vZz6u9Y5JEm2uZUYWGhCAgIEB07dhSOjo4iIiJC1NXV6TlqMmba5tTWrVuFu7u76Nixo3BwcBCTJk0St27d0nPUZIzOnj3b7LERj891TyIE59CJiIiIiIhMhUldk0dERERERPRvxyKPiIiIiIjIhLDIIyIiIiIiMiEs8oiIiIiIiEwIizwiIiIiIiITwiKPiIiIiIjIhLDIIyIiIiIiMiEs8oiIiIiIiEwIizwiIiIT5+Ligi1bthg6DCIi0hMWeUREpHfTp0/HuHHj1NrPnTsHiUSCe/fu6ey92mKfxiouLg5dunRRa798+TJmz56t/4CIiMggpIYOgIiIiJpXX18PS0vLVm9vZ2enw2iIiMjYcSaPiIiM2sWLFzFq1Ch07NgRTk5OCA8PR21trXL9wYMH4ePjg86dO6Nnz56YOHEi7t69CwC4fv06Ro8eDQDo2rUrJBIJpk+f/o/v9fXXX8PJyQnW1tYYP348YmJi1GbGvv/+e3h7e8PKygqurq6Ijo5GQ0ODcr1EIsGePXswfvx4WFtbo1+/fkhMTFTZR0FBAUJCQtCpUyfY29tjypQpqKioUK7/z3/+g/nz5yMiIgI9evRAYGAgACAmJgavvfYaZDIZnJycMHfuXNTU1AB4MmM5Y8YMVFdXQyKRQCKRICoqCoD66ZplZWV4++230alTJ9jY2CA0NBR37txRro+KioKnpycOHDgAFxcX2Nra4t1334VCoXjOX4uIiIwBizwiIjJaV69eRXBwMCZMmIArV64gISEBGRkZmD9/vrJPfX091qxZg59//hknT55EaWmpspBzcnLC8ePHAQBFRUUoLy/Hl19+qfG9MjMzIZfLsXDhQuTl5SEwMBCfffaZSp/Tp09j8uTJCA8PR0FBAXbv3o24uDi1ftHR0QgNDcWVK1cQEhKCSZMmoaqqCgBQXl4Of39/eHp6Ijs7G6dOncKdO3cQGhqqso/9+/dDKpUiMzMTu3fvBgCYmZlh69at+OWXX7B//3788MMPWLZsGQDA19cXW7ZsgY2NDcrLy1FeXo4lS5aojVMIgXHjxqGqqgrp6elIS0tDSUkJwsLCVPqVlJTg5MmTSEpKQlJSEtLT07F+/fpm/15ERGQkBBERkZ5NmzZNmJubC5lMpvKysrISAMRff/0lhBBiypQpYvbs2SrbXrhwQZiZmYkHDx5o3HdWVpYAIBQKhRBCiLNnz6rs85+EhYWJMWPGqLRNmjRJ2NraKpf9/PzE2rVrVfocOHBAODg4KJcBiBUrViiXa2pqhEQiESkpKUIIIVauXCmCgoJU9nHz5k0BQBQVFQkhhPD39xeenp7NxiuEEEeOHBHdu3dXLu/bt08l3qecnZ3FF198IYQQIjU1VZibm4uysjLl+vz8fAFAZGVlCSGEWL16tbC2thb3799X9lm6dKkYNmzYc2MiIiLD4zV5RERkEKNHj8bOnTtV2n788UdMnjxZuZyTk4Pff/8d8fHxyjYhBJqamlBaWooBAwYgNzcXUVFRyMvLQ1VVFZqamgA8OSXR3d29xfEUFRVh/PjxKm1Dhw5FUlKSSjyXL19WmblrbGzEw4cPUVdXB2trawCAh4eHcr1MJkPnzp2Vp5Dm5OTg7Nmz6NSpk1oMJSUl6N+/PwDAx8dHbf3Zs2exdu1aFBQU4P79+2hoaMDDhw9RW1sLmUzWonEWFhbCyckJTk5OyjZ3d3d06dIFhYWFGDJkCIAnp3h27txZ2cfBwUE5BiIiMm4s8oiIyCBkMhleeeUVlbZbt26pLDc1NWHOnDkIDw9X2/7ll19GbW0tgoKCEBQUhIMHD8LOzg5lZWUIDg5GfX29VvEIISCRSNTano0nOjoaEyZMUNveyspK+bOFhYXKOolEoiw+m5qaMHbsWGzYsEFtHw4ODsqfny3abty4gZCQEMjlcqxZswbdunVDRkYGZs6cicePH7dwlJrHqam9uTEQEZFxY5FHRERGy8vLC/n5+WrF4FNXr15FRUUF1q9fr5yZys7OVunz9K6UjY2Nzb6Xm5sbsrKyVNqe3ZeXlxeKior+MZ6W8PLywvHjx+Hi4gKptOX/hrOzs9HQ0IDNmzfDzOzJJfVHjhxR6WNpafnccbq7u6OsrAw3b95U/s4KCgpQXV2NAQMGaDkaIiIyRrzxChERGa0PP/wQly5dwrx585CXl4fi4mIkJiZiwYIFAJ7M5llaWmLbtm24du0aEhMTsWbNGpV9ODs7QyKRICkpCX/++afybpTPWrBgAZKTkxETE4Pi4mLs3r0bKSkpKrNbq1atwjfffIOoqCjk5+ejsLAQCQkJWLFiRYvHNG/ePFRVVeG9995DVlYWrl27htTUVLz//vvNFmh9+/ZFQ0ODcqwHDhzArl27VPq4uLigpqYGZ86cQUVFBerq6tT2ExAQAA8PD0yaNAk//fQTsrKyMHXqVPj7+2s8RZSIiNofFnlERGS0PDw8kJ6ejuLiYvj5+WHw4MFYuXKl8rRGOzs7xMXF4ejRo3B3d8f69euxadMmlX307t0b0dHRiIyMhL29vcqdOf9u5MiR2LVrF2JiYjBo0CCcOnUKixcvVjkNMzg4GElJSUhLS8OQIUMwfPhwxMTEwNnZucVj6tWrFzIzM9HY2Ijg4GAMHDgQCxcuhK2trXKGThNPT0/ExMRgw4YNGDhwIOLj47Fu3TqVPr6+vpDL5QgLC4OdnR02btyoth+JRIKTJ0+ia9euGDVqFAICAuDq6oqEhIQWj4GIiIybRDx7wQEREREBAD744AP8+uuvuHDhgqFDISIiajFek0dERPRfmzZtQmBgIGQyGVJSUrB//37s2LHD0GERERFphTN5RERE/xUaGopz585BoVDA1dUVCxYsgFwuN3RYREREWmGRR0REREREZEJ44xUiIiIiIiITwiKPiIiIiIjIhLDIIyIiIiIiMiEs8oiIiIiIiEwIizwiIiIiIiITwiKPiIiIiIjIhLDIIyIiIiIiMiEs8oiIiIiIiEzI/wEHl3ut7k+MZAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 900x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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(\"Heat generation\")\n",
    "ax.set_ylabel(\"Electricity generation\")\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": [
    "We can build the network and add the buses, loads and links representing the OCGT and CHP unit. We use two idependent link to represent the generation of electricity and heat in the CHP unit and later we will impose constraints limiting the feasible space (as shown in the figure above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['gas boiler'], dtype='object')"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network = pypsa.Network()\n",
    "\n",
    "network.add(\"Bus\", \"electricity\")\n",
    "network.add(\"Load\", \"electricity load\", bus=\"electricity\", p_set=50)\n",
    "\n",
    "network.add(\"Bus\", \"heat\")\n",
    "network.add(\"Load\", \"heat load\", bus=\"heat\", p_set=40)\n",
    "\n",
    "network.add(\"Bus\", \"gas\")\n",
    "\n",
    "# We add a gas store with energy capacity and an initial filling level much higher than the required gas consumption, \n",
    "# this way gas supply is unlimited\n",
    "network.add(\"Store\", \"gas\", e_initial=1e6, e_nom=1e6, bus=\"gas\")\n",
    "\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"OCGT\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"electricity\",\n",
    "    p_nom=1000,\n",
    "    marginal_cost=20,\n",
    "    efficiency=0.35,\n",
    ")\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"CHP generator\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"electricity\",\n",
    "    efficiency=0.3,\n",
    "    p_nom=1000,\n",
    "    marginal_cost=20,\n",
    ")\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"CHP boiler\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"heat\",\n",
    "    p_nom=1000,\n",
    ")\n",
    "\n",
    "\n",
    "network.add(\n",
    "    \"Link\",\n",
    "    \"gas boiler\",\n",
    "    bus0=\"gas\",\n",
    "    bus1=\"heat\",\n",
    "    marginal_cost=20,\n",
    "    p_nom=1000,\n",
    "    efficiency=0.9,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>bus</th>\n",
       "      <th>carrier</th>\n",
       "      <th>type</th>\n",
       "      <th>p_set</th>\n",
       "      <th>q_set</th>\n",
       "      <th>sign</th>\n",
       "      <th>active</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Load</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>electricity load</th>\n",
       "      <td>electricity</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>50.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>heat load</th>\n",
       "      <td>heat</td>\n",
       "      <td></td>\n",
       "      <td></td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-1.0</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                          bus carrier type  p_set  q_set  sign  active\n",
       "Load                                                                  \n",
       "electricity load  electricity                50.0    0.0  -1.0    True\n",
       "heat load                heat                40.0    0.0  -1.0    True"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.loads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We add the CHP constraints to limite the feasible operation space\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pypsa.consistency:The following stores have carriers which are not defined:\n",
      "Index(['gas'], dtype='object', name='Store')\n",
      "WARNING:pypsa.consistency:The following buses have carriers which are not defined:\n",
      "Index(['electricity', 'heat', 'gas'], dtype='object', name='Bus')\n",
      "WARNING:pypsa.consistency:The following links have carriers which are not defined:\n",
      "Index(['OCGT', 'CHP generator', 'CHP boiler', 'gas boiler'], dtype='object', name='Link')\n",
      "INFO:linopy.model: Solve problem using Highs solver\n",
      "INFO:linopy.io: Writing time: 0.04s\n",
      "INFO:linopy.constants: Optimization successful: \n",
      "Status: ok\n",
      "Termination condition: optimal\n",
      "Solution: 6 primals, 16 duals\n",
      "Objective: 3.14e+03\n",
      "Solver model: available\n",
      "Solver message: optimal\n",
      "\n",
      "INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Link-fix-p-lower, Link-fix-p-upper, Store-fix-e-lower, Store-fix-e-upper, Store-energy_balance, backpressure, top_iso_fuel_line were not assigned to the network.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "('ok', 'optimal')"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Guarantees isofuel lines, i.e. increase in heat generation is proportional to decrease in electricity generation\n",
    "network.links.at[\"CHP boiler\", \"efficiency\"] = (network.links.at[\"CHP generator\", \"efficiency\"] / c_v)\n",
    "\n",
    "model = network.optimize.create_model()\n",
    "\n",
    "link_p = model.variables[\"Link-p\"]\n",
    "\n",
    "# Guarantees back-pressure line\n",
    "model.add_constraints(\n",
    "    c_m * network.links.at[\"CHP boiler\", \"efficiency\"] * link_p.sel(Link=\"CHP boiler\")\n",
    "    - network.links.at[\"CHP generator\", \"efficiency\"] * link_p.sel(Link=\"CHP generator\")\n",
    "    <= 0,\n",
    "    name=\"backpressure\",)\n",
    "\n",
    "\n",
    "# Guarantees top iso fuel line\n",
    "model.add_constraints(\n",
    "    link_p.sel(Link=\"CHP boiler\") + link_p.sel(Link=\"CHP generator\")\n",
    "    - network.links.at[\"CHP generator\", \"p_nom\"]\n",
    "    <= 0,\n",
    "    name=\"top_iso_fuel_line\",\n",
    ")\n",
    "\n",
    "network.optimize.solve_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Link</th>\n",
       "      <th>OCGT</th>\n",
       "      <th>CHP generator</th>\n",
       "      <th>CHP boiler</th>\n",
       "      <th>gas boiler</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>snapshot</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>now</th>\n",
       "      <td>-20.0</td>\n",
       "      <td>-30.0</td>\n",
       "      <td>-40.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Link      OCGT  CHP generator  CHP boiler  gas boiler\n",
       "snapshot                                             \n",
       "now      -20.0          -30.0       -40.0         0.0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.links_t.p1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}