electrical_power_generation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cO-7xILPmxfp"
      },
      "source": [
        "# Electrical Power Generation\n",
        "**2021-03-14**\n",
        "\n",
        "Notebook created by **Arthur Verrez** for the course of *Optimisation Combinatoire et Stochastique* at MINES ParisTech, PSL University.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0jITx0l6Itl5"
      },
      "source": [
        "# 1: Basic problem description & notation\n",
        "\n",
        "#### Basic definition of the problem\n",
        "**Periods**\n",
        "\n",
        "Let's define $p\\in \\{0,...,P\\}$ the different time period indices and $\\Delta p \\in \\mathbb{N}$ the duration of each of those periods in hours.\n",
        "\n",
        "For each of those periods, we define a demand in power noted $D_p \\in \\mathbb{R}_{+}$ (in MW).\n",
        "\n",
        "\n",
        "**Power generation units**\n",
        "\n",
        "All the power generation units are grouped by types noted $t\\in T$. Let's note $N_t \\in \\mathbb{N}$ the total number of units by type.\n",
        "\n",
        "If a unit of type $t$ is on during the period $p$, then its power output is contant during a period and has a value between $\\underline{L_t}\\in \\mathbb{R}_{+}$ and $\\overline{L_t}\\in ]\\underline{L_t}, +\\infty [$ (both in MW).\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "**Costs**\n",
        "\n",
        "Each power generation unit type has two types of costs:\n",
        "$C_{t}^{b} \\in \\mathbb{R}_{+}$ (in euro/hour) to operate the unit at its minimum power output, and $C_{t}^{r} \\in \\mathbb{R}_{+}$ to generate additional MWh (in euro/MWh).\n",
        "\n",
        "#### Definition of the model 2.1.1\n",
        "Let's note $x_{t,i}^{p}\\in \\{0,1\\}$ the binary variable that defines if the $i^{th}$ unit of type $t$ is on or off at the period $p$ and $x_{t}^{p}:=\\sum\\limits_{i=0}^{N_{t}} x_{t,i}^{p}$, the number of units of type $t$ on at the period $p$.\n",
        "\n",
        "Let's note $l_{t,i}^{p}$ the additional MW generated by said unit at the period $p$ and $l_{t}^{p}:=\\sum\\limits_{i=0}^{N_{t}} l_{t,i}^{p}$, the total additional MW generated by all the units of type $t$ at the period $p$.\n",
        "\n",
        "We can therefore define the output power by type $t$ and period $p$:\n",
        "$o_t^p:= \\Delta_p(\\underline{L_t}\\times x_{t}^{p}+ l_{t}^{p})$ and $o^p:=\\sum\\limits_{t\\in T} \\Delta_p(\\underline{L_t}\\times x_{t}^{p}+ l_{t}^{p})$\n",
        "\n",
        "\n",
        "Let's note $c_{t,i}^{p}$ the cost of running the $i^{th}$ unit of type $t$ during the period $p$, we have\n",
        "$c_{t,i}^{p}=\\Delta_p( x_{t,i}^{p}\\times C_{t}^{b} + l_{t,i}^{p}\\times C_{t}^{r})$\n",
        "and\n",
        "$c_t^{p}:=\\sum\\limits_{i=0}^{N_t}\\Delta_p( x_{p}^{t,i}\\times C_{b}^{t} + l_{p}^{t,i}\\times C_{r}^{t})\n",
        "=\\Delta_p(x_{p}^{t}\\times C_{b}^{t} + l_{p}^{t}\\times C_{r}^{t})$\n",
        "the total cost of the units of type $t$ at the period $p$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hUskOCH9myzl"
      },
      "source": [
        "#2: Discrete Optimization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pekX3vTzmy2Q"
      },
      "source": [
        "## 2.1: Modelling (un)differentiated units\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L8UbKZOvDldW"
      },
      "source": [
        "### Question 2.1.1\n",
        "\n",
        "\n",
        "**Constraints**\n",
        "\n",
        "We have the following constraints, based on binary and continuous variables:\n",
        "\n",
        "\n",
        "(C1) $\\forall t \\in T, \\forall (i,p)\\in N_{t}\\times [\\![ 0, P ]\\!]:  0 \\leq l_{t,i}^{p} \\leq (\\overline{L_t}-\\underline{L_t})x_{t,i}^{p}$\n",
        "\n",
        "\n",
        "(C2) $\\forall p \\in [\\![ 0, P ]\\!]:\\sum\\limits_{i=0}^{N_{t}}\\underline{L_t}\\times  x_{t,i}^{p}+ l_{t}^{p}$\n",
        "\n",
        "(C3) $\\forall (t,p) \\in T\\times [\\![ 0, P ]\\!]:  0 \\leq \\sum\\limits_{i=0}^{N_{t}} x_{t,i}^{p} \\leq N_t$\n",
        "\n",
        "\n",
        "\n",
        "**The objective**\n",
        "\n",
        "The objective is to minimize the following expression:\n",
        "\n",
        "(O1) $\\sum\\limits_{t\\in T}\\sum\\limits_{p=0}^P c_t^{p}=\\sum\\limits_{t\\in T}\\sum\\limits_{p=0}^P\\Delta_p(\\sum\\limits_{i=0}^{N_t}(x_{p}^{t,i})\\times C_{b}^{t} + l_{p}^{t}\\times C_{r}^{t})$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dZRKC4-3ZROd"
      },
      "source": [
        "### Question 2.1.2\n",
        "\n",
        "This model is extremely similar to  the first one, with an added assumption, using the same conventions as in the question 2.1.1, we can express it as:\n",
        "$ \\forall (p,t)\\in [\\![ 0, P ]\\!]\\times T, \\exists K_t^{p}\\in \\mathbb{R}_{+} / \\forall i \\in [\\![ 0, N_t ]\\!]$, if $ x_{t,i}^{p}=1 $ then $l_{t,i}^{p}=K $\n",
        "\n",
        "\n",
        "We can therefore rewrite the constraints as following:\n",
        "\n",
        "**Constraints**\n",
        "\n",
        "(C1.2)\n",
        "$\\forall (t,p) \\in T \\times [\\![ 0, P ]\\!]:  0 \\leq l_{t}^{p} \\leq (\\overline{L_t}-\\underline{L_t})x_{t}^{p}$\n",
        "\n",
        "(C2) $ \\forall p \\in [\\![ 0, P ]\\!]: o_p \\geq \\Delta_p D_p$\n",
        "\n",
        "(C3) $\\forall (t,p) \\in T\\times [\\![ 0, P ]\\!]:  0 \\leq x_{t}^{p} \\leq N_t$\n",
        "\n",
        "\n",
        "\n",
        "**The objective**\n",
        "\n",
        "The objective is to minimize the following expression:\n",
        "\n",
        "(O1) $\\sum\\limits_{t\\in T}\\sum\\limits_{p=0}^P c_t^{p}=\\sum\\limits_{t\\in T}\\sum\\limits_{p=0}^P\\Delta_p(x_{p}^{t}\\times C_{b}^{t} + l_{p}^{t}\\times C_{r}^{t})$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UXfItSBcZXyE"
      },
      "source": [
        "### Question 2.1.3\n",
        "\n",
        "The model 2.1.1 is a relaxation of the model 2.1.2 because the model 2.1.2 adds the constraint that all the committed units of a same type operate with the same power output, which is not a constraint present in the model 2.1.1.\n",
        "\n",
        "But in reality, both models are identical in the sense that both model will produce the same objective cost.\n",
        "\n",
        "#### Proof\n",
        "\n",
        "The model 2.1.1 being a relaxation of the model 2.1.2, the space of solutions of the model 2.1.2 is included in the space of solutions of the model 2.1.1.\n",
        "\n",
        "That means that the optimal solution for the model 2.1.2 will be a solution for the model 2.1.1, which directly implies that the optimal solution for the model 2.1.1 is better or as good as the optimal solution for the model 2.1.2.\n",
        "\n",
        "Let's take, for a given problem, the optimal solution of the model 2.1.1 and create a solution that satisfies the conditions of the model 2.1.2, the solution will then be optimal for the model 2.1.2 and we will have shown the fact that both model are identical.\n",
        "\n",
        "The optimal solution for the model 2.1.1 is caracterised by $l_{t,i}^{p}$ and $x_{t,i}^{p}$ for all relevant $p$, $t$ and $i$.\n",
        "\n",
        "If $ \\forall (p,t)\\in [\\![ 0, P ]\\!]\\times T, \\exists K_t^{p}\\in \\mathbb{R}_{+} / \\forall i \\in [\\![ 0, N_t ]\\!]$, if $ x_{t,i}^{p}=1 $ then $l_{t,i}^{p}=K $, then this solution is already valid for the model 2.1.2.\n",
        "\n",
        "Else, we create a solution with the same $ x_{t,i}^{p} $ and for all $p$ and $t$, we define:\n",
        "\n",
        "$\\widetilde{l_t^p}=  \\begin{cases}\n",
        "               0 \\text{ if } \\sum\\limits_{i=0}^{N_t} x_{t,i}^{p} = 0\\\\\n",
        "               \\frac{\\sum\\limits_{i=0}^{N_t} l_{t,i}^{p}}{\\sum\\limits_{i=0}^{N_t} x_{t,i}^{p}}  \\text{ otherwise }\n",
        "               \n",
        "            \\end{cases}\n",
        "$\n",
        "\n",
        "It is trivial to verify that a solution caracterised by the given $ x_{t,i}^{p} $ and $\\widetilde{l_t^p}$ outputs the same power for each $t$ and each $p$ as the optimal solution of the model 2.1.1, and it is therefore also optimal. But this solution also verifies the conditions of the model 2.1.2.\n",
        "\n",
        "Which concludes the proof.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EFzI-_HSn2pd"
      },
      "source": [
        "##2.2: Additional Models"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BhDZ_H-ImQ1-"
      },
      "source": [
        "### Question 2.2.1\n",
        "\n",
        "With a reserve factor $F\\in [0,1]$ denoting the **additional** possible demand, it simply adds a constraint that is the following:\n",
        "\n",
        "(C4) $\\forall p \\in [\\![ 0, P ]\\!]:  \\sum\\limits_{t\\in T}\\overline{L_t}\\times x_{t}^{p} \\geq (1+F) \\times D_p$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6-RbFBi-oKll"
      },
      "source": [
        "### Question 2.2.2\n",
        "\n",
        "Considering for a given type $t$, $A_t \\in \\mathbb{N}$ the number of units of this type that is already on at time 0 and a startup cost $C_t^s \\in \\mathbb{R}^{+}$ (in euros) of transitioning a unit from off to on, the only thing that changes is the objective.\n",
        "\n",
        "We have to add the cost of starting up machines.\n",
        "This one is equal to:\n",
        "$C_t^s\\sum\\limits_{t\\in T}(max(x_t^0-A_t,0) + \\sum\\limits_{p=1}^P  max(x_t^{p}-x_t^{p-1},0))$\n",
        "\n",
        "\n",
        "The final objective to minimize therefore becomes:\n",
        "\n",
        "(O2) $\\sum\\limits_{t\\in T}\\sum\\limits_{p=0}^P c_t^{p}+C_t^s\\sum\\limits_{t\\in T}(max(x_t^0-A_t,0) + \\sum\\limits_{p=1}^P max(x_t^{p}-x_t^{p-1},0))$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JFD8zmk4oKtq"
      },
      "source": [
        "\n",
        "### Question 2.2.3\n",
        "\n",
        "The hydropower generation changes a lot of parts of the models.\n",
        "We'll consider that for our use case, the available reservoirs of the hydroelectric will never be entirely depleted.\n",
        "\n",
        "Let's note $y_h^p$ the binary variable being equal to 1 if the $h$ hydro power unit is on at the period $p$ and otherwise 0.\n",
        "\n",
        "**Regarding the constraints**\n",
        "\n",
        "The power output can be re-written as $o^p:=\\sum\\limits_{t\\in T} \\Delta_p(\\underline{L_t}\\times x_{t}^{p}+ l_{t}^{p})+\\sum\\limits_{h\\in H} \\Delta_p(L_h\\times y_{h}^{p})$\n",
        "\n",
        "We now have to take into account the power consumption of pumping back the water into the reservoirs. We will consider that pumping back the water into the reservoir can be done during the same period as the usage for simplification reasons.\n",
        "\n",
        "The depletion of water at period $p$ will be equal to:\n",
        "$\\sum\\limits_{h\\in H} R_h\\times y_{h}^{p}$ meters/hour.\n",
        "\n",
        "Therefore, the power needed to pump back the water per hour will be:\n",
        "$w_p:=E\\sum\\limits_{h\\in H} R_h\\times y_{h}^{p}$.\n",
        "\n",
        "The constraint 2, with the updated $o_p$, becomes:\n",
        "\n",
        "(C2.3) $ \\forall p \\in [\\![ 0, P ]\\!]: o_p \\geq \\Delta_p(D_p + w_p)$\n",
        "\n",
        "Furthermore, all hydroelectric units participate to the reserve, therefore, constraint 4 becomes:\n",
        "\n",
        "(C4.3) $\\forall (t,p) \\in T\\times [\\![ 0, P ]\\!]:  \\sum\\limits_{t\\in T}\\overline{L_t}\\times x_{t}^{p} + \\sum\\limits_{h\\in H} L_h\\geq (1+F) \\times D_p$\n",
        "\n",
        "**Regarding the objectives**\n",
        "\n",
        "We have to take into account the costs of the hydroelectric units.\n",
        "\n",
        "The operation cost at each period $p$:\n",
        "\n",
        "$\\Delta_p\\sum\\limits_{h\\in H}  y_{h}^{p}\\times C_h^b$\n",
        "\n",
        "\n",
        "The startup costs on all periods will be (considering that no hyroelectric plant is turned on in the beginning):\n",
        "\n",
        "$C_h^s\\sum\\limits_{h\\in H}\\sum\\limits_{p=1}^P max(y_h^{p}-y_h^{p-1},0)$\n",
        "\n",
        "The objective will therefore be to minimize the following expression:\n",
        "\n",
        "\n",
        "(O3) $\\sum\\limits_{t\\in T}\\sum\\limits_{p=0}^P c_t^{p}+C_t^s\\sum\\limits_{t\\in T}(max(x_t^1-A_t,0) + \\sum\\limits_{p=1}^P \\times max(x_t^{p}-x_t^{p-1},0))+C_h^s\\sum\\limits_{h\\in H}\\sum\\limits_{p=1}^P (max(y_h^{p}-y_h^{p-1},0)+\\Delta_p  y_{h}^{p}\\times C_h^b)$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f8UiL3OnoUTE"
      },
      "source": [
        "##2.3: Implementation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KsnJo8Lgs0Uz"
      },
      "source": [
        "### 2.3.0 General set up"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rvFOSjd6rIDw",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "d1c45520-42e4-41c5-f151-a9e02d695d19"
      },
      "source": [
        "# Installation of Gurobi\n",
        "!pip install -i https://pypi.gurobi.com gurobipy"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Looking in indexes: https://pypi.gurobi.com\n",
            "Requirement already satisfied: gurobipy in /usr/local/lib/python3.7/dist-packages (9.1.1)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s_djxor7rJMQ"
      },
      "source": [
        "# Importing the libraries we'll use\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "import gurobipy as gp\n",
        "from gurobipy import GRB, max_"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tx7gw1XwrdIw"
      },
      "source": [
        "# Defining the inputs of the problem\n",
        "# Names are self explanatory\n",
        "\n",
        "\n",
        "n_units_by_type= [12, 10, 5]\n",
        "\n",
        "n_power_types= len(n_units_by_type)\n",
        "\n",
        "min_power_by_type= [850, 1250, 1500] # in MW\n",
        "\n",
        "max_power_by_type= [2000, 1750, 4000] # in MW\n",
        "\n",
        "base_cost_by_type= [1000, 2600, 3000] # in euro/hour\n",
        "\n",
        "per_mw_cost_by_type= [2, 1.3, 3] # in euro/MW\n",
        "\n",
        "max_units_by_type=max(n_units_by_type)\n",
        "\n",
        "\n",
        "power_demand_by_period= [15000, 30000, 25000, 40000, 27000] # in MW\n",
        "\n",
        "period_durations= [6, 3, 6, 3, 6] # in hours\n",
        "\n",
        "n_periods= len(period_durations)\n",
        "\n",
        "#For section 2.3\n",
        "min_down_time = [9,9,9]\n",
        "min_up_time = [9,9,9]\n",
        "max_ramp_down = [800,1200,1700]\n",
        "max_ramp_up = [800,1200,1700]\n",
        "max_switch_on = [1000,1500,2000]\n",
        "max_switch_off = [1000,1500,2000]\n",
        "\n",
        "# For question 2.2.1\n",
        "reserve_factor = 0.15\n",
        "\n",
        "#For question 2.2.1\n",
        "\n",
        "startup_cost_by_type = [2000, 1000, 500]\n",
        "initial_units_on_by_type=[0,0,0]\n",
        "\n",
        "#For question 2.2.3\n",
        "n_hydro_units = 2\n",
        "hydro_load = [900, 1400]\n",
        "hydro_cost = [90, 150]\n",
        "hydro_height_reduction = [0.31, 0.47]\n",
        "hydro_startup_cost = [1500, 1200]\n",
        "pumping_consumption=3000 #in MWh/meter of height\n",
        "\n",
        "\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WRdsfoWbrM3m"
      },
      "source": [
        "### Question 2.3.1 (Model 1)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ibRrS6airM3t"
      },
      "source": [
        "# Definition of the model 1\n",
        "\n",
        "model_1 = gp.Model('PowerGenerationModel1')\n",
        "\n",
        "\n",
        "# Definition of the variables\n",
        "\n",
        "# A 2 dimensional array that gives by type and period the integer correponding to the number of units on\n",
        "max_units_every_type=max(n_units_by_type)\n",
        "\n",
        "units_on=model_1.addVars(n_power_types, n_periods, max_units_every_type, vtype=GRB.BINARY, name=\"units_on\")\n",
        "\n",
        "# A 2 dimensional array that gives the power generated at given period\n",
        "gen_power = model_1.addVars(n_power_types, n_periods, vtype=GRB.CONTINUOUS, name=\"gen_power\")\n",
        "\n",
        "\n",
        "# Definition of the constraints\n",
        "\n",
        "# The maximum number of units on at a given time (C3)\n",
        "max_units_on_contr = model_1.addConstrs(gp.quicksum(units_on[type_id, period_id, unit_id] for unit_id in range(max_units_every_type)) <= n_units_by_type[type_id]\n",
        "                        for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# Additional constraint due to the definition of gen_power to insure it produces at least the base power per unit\n",
        "min_gen_power_constr = model_1.addConstrs((gen_power[type_id, period_id] >= min_power_by_type[type_id] * gp.quicksum(units_on[type_id, period_id, unit_id] for unit_id in range(max_units_every_type)))\n",
        "                              for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# The maximum power output given the units on at a given period, equivalent to (C1.2)\n",
        "max_gen_power_constr = model_1.addConstrs((gen_power[type_id, period_id] <= max_power_by_type[type_id] * gp.quicksum(units_on[type_id, period_id, unit_id] for unit_id in range(max_units_every_type)))\n",
        "                              for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# Constraint to meet the demand in power (C2)\n",
        "meet_demand_constr = model_1.addConstrs(gp.quicksum(gen_power[type_id, period_id]\n",
        "                                                  for type_id in range(n_power_types))\n",
        "                                                  >= power_demand_by_period[period_id]\n",
        "                                                  for period_id in range(n_periods))\n",
        "\n",
        "# Definition of the costs\n",
        "base_cost = gp.quicksum(base_cost_by_type[type_id]*period_durations[period_id]*gp.quicksum(units_on[type_id, period_id, unit_id] for unit_id in range(max_units_every_type))\n",
        "                    for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "additional_cost = gp.quicksum(per_mw_cost_by_type[type_id]*(gen_power[type_id,period_id] - min_power_by_type[type_id]*gp.quicksum(units_on[type_id, period_id, unit_id] for unit_id in range(max_units_every_type)))*period_durations[period_id]\n",
        "                      for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# Setting the objective for the model\n",
        "model_1.setObjective(base_cost+additional_cost)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9MjQf4ffrM3t",
        "outputId": "4110096a-eb99-44ce-8215-007e9cfb87d1"
      },
      "source": [
        "model_1.optimize()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "Optimize a model with 50 rows, 195 columns and 585 nonzeros\n",
            "Model fingerprint: 0xfb5006c2\n",
            "Variable types: 15 continuous, 180 integer (180 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 4e+03]\n",
            "  Objective range  [4e+00, 9e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [5e+00, 4e+04]\n",
            "Found heuristic solution: objective 1294050.0000\n",
            "Presolve removed 43 rows and 189 columns\n",
            "Presolve time: 0.01s\n",
            "Presolved: 7 rows, 6 columns, 15 nonzeros\n",
            "Found heuristic solution: objective 1011000.0000\n",
            "Variable types: 0 continuous, 6 integer (0 binary)\n",
            "\n",
            "Root relaxation: objective 9.780000e+05, 3 iterations, 0.00 seconds\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 978000.000    0    1 1011000.00 978000.000  3.26%     -    0s\n",
            "H    0     0                    980340.00000 978000.000  0.24%     -    0s\n",
            "*    0     0               0    978900.00000 978900.000  0.00%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  MIR: 1\n",
            "\n",
            "Explored 1 nodes (4 simplex iterations) in 0.04 seconds\n",
            "Thread count was 2 (of 2 available processors)\n",
            "\n",
            "Solution count 4: 978900 980340 1.011e+06 1.29405e+06 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 9.789000000000e+05, best bound 9.789000000000e+05, gap 0.0000%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 195
        },
        "id": "5xxZXXA7rM3u",
        "outputId": "64fe5c50-083e-44a2-a214-b5c9fa8e8c54"
      },
      "source": [
        "periods_info=[]\n",
        "for p in range(n_periods):\n",
        "  periods_info.append({})\n",
        "  cost_additional=0\n",
        "  cost_base=0\n",
        "  total_gen_power=0\n",
        "  for t in range(n_power_types):\n",
        "    periods_info[p][\"Power type \"+str(t)]=gen_power[t,p].x\n",
        "    periods_info[p][\"N units on \"+str(t)]=sum(units_on[t,p,i].x for i in range(max_units_every_type))\n",
        "    total_gen_power+=gen_power[t,p].x\n",
        "    cost_additional+=per_mw_cost_by_type[t]*(gen_power[t,p].x - min_power_by_type[t]*sum(units_on[t,p,i].x for i in range(max_units_every_type)))*period_durations[p]\n",
        "    cost_base+=base_cost_by_type[t]*period_durations[p]*sum(units_on[t,p,i].x for i in range(max_units_every_type))\n",
        "\n",
        "  periods_info[p][\"Total gen power\"]=total_gen_power\n",
        "  periods_info[p][\"Demand\"]=power_demand_by_period[p]\n",
        "  periods_info[p][\"Cost Base\"]=cost_base\n",
        "  periods_info[p][\"Cost Additional\"]=cost_additional\n",
        "\n",
        "\n",
        "pd.DataFrame(periods_info)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "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>Power type 0</th>\n",
              "      <th>N units on 0</th>\n",
              "      <th>Power type 1</th>\n",
              "      <th>N units on 1</th>\n",
              "      <th>Power type 2</th>\n",
              "      <th>N units on 2</th>\n",
              "      <th>Total gen power</th>\n",
              "      <th>Demand</th>\n",
              "      <th>Cost Base</th>\n",
              "      <th>Cost Additional</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>11500.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>3500.0</td>\n",
              "      <td>2.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>15000.0</td>\n",
              "      <td>15000</td>\n",
              "      <td>103200.0</td>\n",
              "      <td>23400.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>12500.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>30000.0</td>\n",
              "      <td>30000</td>\n",
              "      <td>114000.0</td>\n",
              "      <td>33300.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>11000.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>14000.0</td>\n",
              "      <td>8.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>25000.0</td>\n",
              "      <td>25000</td>\n",
              "      <td>196800.0</td>\n",
              "      <td>40800.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>22500.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>40000.0</td>\n",
              "      <td>40000</td>\n",
              "      <td>114000.0</td>\n",
              "      <td>93300.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>11250.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>15750.0</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>27000.0</td>\n",
              "      <td>27000</td>\n",
              "      <td>212400.0</td>\n",
              "      <td>47700.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Power type 0  N units on 0  Power type 1  ...  Demand  Cost Base  Cost Additional\n",
              "0       11500.0          12.0        3500.0  ...   15000   103200.0          23400.0\n",
              "1       12500.0          12.0       17500.0  ...   30000   114000.0          33300.0\n",
              "2       11000.0          12.0       14000.0  ...   25000   196800.0          40800.0\n",
              "3       22500.0          12.0       17500.0  ...   40000   114000.0          93300.0\n",
              "4       11250.0          12.0       15750.0  ...   27000   212400.0          47700.0\n",
              "\n",
              "[5 rows x 10 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 304
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E3xQTKG9CNCA"
      },
      "source": [
        "# Definition of the model\n",
        "\n",
        "model_2 = gp.Model('PowerGenerationModel2')\n",
        "\n",
        "\n",
        "# Definition of the variables\n",
        "\n",
        "# A 2 dimensional array that gives by type and period the integer correponding to the number of units on\n",
        "n_units_on=model_2.addVars(n_power_types, n_periods, vtype=GRB.INTEGER, name=\"n_units_on\")\n",
        "\n",
        "# A 2 dimensional array that gives the power generated at given period\n",
        "gen_power = model_2.addVars(n_power_types, n_periods, vtype=GRB.CONTINUOUS, name=\"gen_power\")\n",
        "\n",
        "\n",
        "# Definition of the constraints\n",
        "\n",
        "# The maximum number of units on at a given time (C3)\n",
        "max_units_on_contr = model_2.addConstrs(n_units_on[type_id, period_id] <= n_units_by_type[type_id]\n",
        "                        for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# Additional constraint due to the definition of gen_power to insure it produces at least the base power per unit\n",
        "min_gen_power_constr = model_2.addConstrs((gen_power[type_id, period_id] >= min_power_by_type[type_id] * n_units_on[type_id, period_id])\n",
        "                              for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# The maximum power output given the units on at a given period, equivalent to (C1.2)\n",
        "max_gen_power_constr = model_2.addConstrs((gen_power[type_id, period_id] <= max_power_by_type[type_id] * n_units_on[type_id, period_id])\n",
        "                              for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# Constraint to meet the demand in power (C2)\n",
        "meet_demand_constr = model_2.addConstrs(gp.quicksum(gen_power[type_id, period_id]\n",
        "                                                  for type_id in range(n_power_types))\n",
        "                                                  >= power_demand_by_period[period_id]\n",
        "                                                  for period_id in range(n_periods))\n",
        "\n",
        "# Definition of the costs\n",
        "base_cost = gp.quicksum(base_cost_by_type[type_id]*period_durations[period_id]*n_units_on[type_id, period_id]\n",
        "                    for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "additional_cost = gp.quicksum(per_mw_cost_by_type[type_id]*(gen_power[type_id,period_id] - min_power_by_type[type_id]*n_units_on[type_id, period_id])*period_durations[period_id]\n",
        "                      for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "# Setting the objective for the model\n",
        "model_2.setObjective(base_cost+additional_cost)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_bRR_1KUDKLI",
        "outputId": "b8197853-42b0-45dc-b994-e3a12e598c7f"
      },
      "source": [
        "model_2.optimize()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "Optimize a model with 50 rows, 30 columns and 90 nonzeros\n",
            "Model fingerprint: 0x1d9a580e\n",
            "Variable types: 15 continuous, 15 integer (0 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 4e+03]\n",
            "  Objective range  [4e+00, 9e+03]\n",
            "  Bounds range     [0e+00, 0e+00]\n",
            "  RHS range        [5e+00, 4e+04]\n",
            "Presolve removed 43 rows and 24 columns\n",
            "Presolve time: 0.00s\n",
            "Presolved: 7 rows, 6 columns, 15 nonzeros\n",
            "Variable types: 0 continuous, 6 integer (0 binary)\n",
            "Found heuristic solution: objective 992400.00000\n",
            "\n",
            "Root relaxation: objective 9.780000e+05, 3 iterations, 0.00 seconds\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 978000.000    0    1 992400.000 978000.000  1.45%     -    0s\n",
            "H    0     0                    980340.00000 978000.000  0.24%     -    0s\n",
            "*    0     0               0    978900.00000 978900.000  0.00%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  MIR: 1\n",
            "\n",
            "Explored 1 nodes (4 simplex iterations) in 0.03 seconds\n",
            "Thread count was 2 (of 2 available processors)\n",
            "\n",
            "Solution count 3: 978900 980340 992400 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 9.789000000000e+05, best bound 9.789000000000e+05, gap 0.0000%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 195
        },
        "id": "XC-77-Gph-5a",
        "outputId": "01f57c63-8773-4eb5-979d-34ef5ebd0601"
      },
      "source": [
        "periods_info=[]\n",
        "for p in range(n_periods):\n",
        "  periods_info.append({})\n",
        "  cost_additional=0\n",
        "  cost_base=0\n",
        "  total_gen_power=0\n",
        "  for t in range(n_power_types):\n",
        "    periods_info[p][\"Power type \"+str(t)]=gen_power[t,p].x\n",
        "    periods_info[p][\"N units on \"+str(t)]=n_units_on[t,p].x\n",
        "    total_gen_power+=gen_power[t,p].x\n",
        "    cost_additional+=per_mw_cost_by_type[t]*(gen_power[t,p].x - min_power_by_type[t]*n_units_on[t, p].x)*period_durations[p]\n",
        "    cost_base+=base_cost_by_type[t]*period_durations[p]*n_units_on[t, p].x\n",
        "\n",
        "  periods_info[p][\"Total gen power\"]=total_gen_power\n",
        "  periods_info[p][\"Demand\"]=power_demand_by_period[p]\n",
        "  periods_info[p][\"Cost Base\"]=cost_base\n",
        "  periods_info[p][\"Cost Additional\"]=cost_additional\n",
        "\n",
        "\n",
        "pd.DataFrame(periods_info)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "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>Power type 0</th>\n",
              "      <th>N units on 0</th>\n",
              "      <th>Power type 1</th>\n",
              "      <th>N units on 1</th>\n",
              "      <th>Power type 2</th>\n",
              "      <th>N units on 2</th>\n",
              "      <th>Total gen power</th>\n",
              "      <th>Demand</th>\n",
              "      <th>Cost Base</th>\n",
              "      <th>Cost Additional</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>11500.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>3500.0</td>\n",
              "      <td>2.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>15000.0</td>\n",
              "      <td>15000</td>\n",
              "      <td>103200.0</td>\n",
              "      <td>23400.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>12500.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>30000.0</td>\n",
              "      <td>30000</td>\n",
              "      <td>114000.0</td>\n",
              "      <td>33300.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>11000.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>14000.0</td>\n",
              "      <td>8.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>25000.0</td>\n",
              "      <td>25000</td>\n",
              "      <td>196800.0</td>\n",
              "      <td>40800.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>22500.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>40000.0</td>\n",
              "      <td>40000</td>\n",
              "      <td>114000.0</td>\n",
              "      <td>93300.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>11250.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>15750.0</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>27000.0</td>\n",
              "      <td>27000</td>\n",
              "      <td>212400.0</td>\n",
              "      <td>47700.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Power type 0  N units on 0  Power type 1  ...  Demand  Cost Base  Cost Additional\n",
              "0       11500.0          12.0        3500.0  ...   15000   103200.0          23400.0\n",
              "1       12500.0          12.0       17500.0  ...   30000   114000.0          33300.0\n",
              "2       11000.0          12.0       14000.0  ...   25000   196800.0          40800.0\n",
              "3       22500.0          12.0       17500.0  ...   40000   114000.0          93300.0\n",
              "4       11250.0          12.0       15750.0  ...   27000   212400.0          47700.0\n",
              "\n",
              "[5 rows x 10 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 307
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mMCoX8gcDpsr"
      },
      "source": [
        "### Question 2.3.2 (Model 2)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zv5VESyZvHd8"
      },
      "source": [
        "Since we proved at question 2.1.3 that both models were the same, we will use the model 2 which is easier and faster to implement."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CAxBJHlDug9K"
      },
      "source": [
        "model_2_2 = gp.Model('PowerGeneration2')\n",
        "\n",
        "####Defining the variables\n",
        "\n",
        "#The number of thermal generators on at a given period\n",
        "n_units_on = model_2_2.addVars(n_power_types, n_periods, vtype=GRB.INTEGER, name=\"n_units_on\")\n",
        "\n",
        "#The number of thermal generators that need to be start at a period p\n",
        "n_start = model_2_2.addVars(n_power_types, n_periods, vtype=GRB.INTEGER, name=\"n_start\")\n",
        "\n",
        "#The total generated power for the thermal generators (in MW)\n",
        "gen_power = model_2_2.addVars(n_power_types, n_periods, vtype=GRB.CONTINUOUS, name=\"gen_power\")\n",
        "\n",
        "#The hydro units that are on or off at a given time p\n",
        "hydro_on = model_2_2.addVars(n_hydro_units, n_periods, vtype=GRB.BINARY, name=\"hydro_on\")\n",
        "\n",
        "#To count the number of starts of hydro generators over the periods\n",
        "hydro_start = model_2_2.addVars(n_hydro_units, n_periods, vtype=GRB.BINARY, name=\"hydro_start\")\n",
        "\n",
        "#The pumping needed to keep the same level of water\n",
        "pumping = model_2_2.addVars(n_periods, vtype=GRB.CONTINUOUS, name=\"pumping\")\n",
        "\n",
        "#The height of the water total by period\n",
        "height = model_2_2.addVars(n_periods, vtype=GRB.CONTINUOUS, name=\"height\")\n",
        "\n",
        "###### Constraints\n",
        "\n",
        "#The maximum number of thermal generators that can be on at a given period\n",
        "max_units_on_contr = model_2_2.addConstrs(n_units_on[type_id, period_id] <= n_units_by_type[type_id]\n",
        "                         for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "\n",
        "\n",
        "# The minimum and maximum output of the thermal generators\n",
        "min_gen_power_constr = model_2_2.addConstrs((gen_power[type_id, period_id] >= min_power_by_type[type_id] * n_units_on[type_id, period_id])\n",
        "                              for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "max_gen_power_constr = model_2_2.addConstrs((gen_power[type_id, period_id] <= max_power_by_type[type_id] * n_units_on[type_id, period_id])\n",
        "                              for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "\n",
        "# The constraint to satisfy the demand\n",
        "meet_demand_constr = model_2_2.addConstrs(gp.quicksum(gen_power[type_id, period_id] for type_id in range(n_power_types)) +\n",
        "                               gp.quicksum(hydro_load[unit]*hydro_on[unit,period_id] for unit in range(n_hydro_units))\n",
        "                               >= power_demand_by_period[period_id] + pumping[period_id]\n",
        "                               for period_id in range(n_periods))\n",
        "\n",
        "# Fixing the height of the water\n",
        "reservoir_constr = model_2_2.addConstrs(height[period_id] == height[period_id-1] + period_durations[period_id]*pumping[period_id]/pumping_consumption -\n",
        "                             gp.quicksum(hydro_height_reduction[unit]*period_durations[period_id]*hydro_on[unit,period_id] for unit in range(n_hydro_units))\n",
        "                             for period_id in range(1,n_periods))\n",
        "\n",
        "#Adding the constraint that the height of the water must be the same at the beginning and end considering the pumping\n",
        "reservoir_same_constr = model_2_2.addConstr(height[0] == height[n_periods-1] + period_durations[0]*pumping[0]/pumping_consumption -\n",
        "                             gp.quicksum(hydro_height_reduction[unit]*period_durations[0]*hydro_on[unit,0]\n",
        "                             for unit in range(n_hydro_units)))\n",
        "\n",
        "# Keeping a reserve of power for the reserve_factor\n",
        "reserve_constr = model_2_2.addConstrs(gp.quicksum(max_power_by_type[type_id]*n_units_on[type_id, period_id] for type_id in range(n_power_types)) >=\n",
        "                           (1+reserve_factor)*power_demand_by_period[period_id] - sum(hydro_load)\n",
        "                           for period_id in range(n_periods))\n",
        "\n",
        "# The initialization of the startup of the thermal generators\n",
        "startup_at_0_constr = model_2_2.addConstrs((n_units_on[type_id,0] <= initial_units_on_by_type[type_id] + n_start[type_id,0])\n",
        "                            for type_id in range(n_power_types))\n",
        "\n",
        "# Fixing the maximum number of startup for the thermal generators\n",
        "startup_constr = model_2_2.addConstrs((n_units_on[type_id,period_id] <= n_units_on[type_id,period_id-1] + n_start[type_id,period_id])\n",
        "                           for type_id in range(n_power_types) for period_id in range(1,n_periods))\n",
        "\n",
        "# Hydro generator startup constraints (basically the same as above for the thermal generators)\n",
        "hydro_startup_at_0_constr = model_2_2.addConstrs((hydro_on[unit,0] <= hydro_start[unit,0])\n",
        "                                    for unit in range(n_hydro_units))\n",
        "\n",
        "hydro_startup_constr = model_2_2.addConstrs((hydro_on[unit,period_id] <= hydro_on[unit,period_id-1] + hydro_start[unit,period_id])\n",
        "                           for unit in range(n_hydro_units) for period_id in range(1,n_periods))\n",
        "\n",
        "####### The costs (self explanatory names)\n",
        "\n",
        "base_cost = gp.quicksum(base_cost_by_type[type_id]*period_durations[period_id]*n_units_on[type_id,period_id]\n",
        "                    for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "additional_cost = gp.quicksum(per_mw_cost_by_type[type_id]*period_durations[period_id]*(gen_power[type_id,period_id] - min_power_by_type[type_id]*n_units_on[type_id,period_id])\n",
        "                       for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "startup_cost = gp.quicksum(startup_cost_by_type[type_id]*n_start[type_id,period_id]\n",
        "                         for type_id in range(n_power_types) for period_id in range(n_periods))\n",
        "\n",
        "hydro_base_cost = gp.quicksum(hydro_cost[unit]*period_durations[period_id]*hydro_on[unit,period_id]\n",
        "                           for unit in range(n_hydro_units) for period_id in range(n_periods))\n",
        "\n",
        "hydro_startup = gp.quicksum(hydro_startup_cost[unit]*hydro_start[unit,period_id]\n",
        "                            for unit in range(n_hydro_units) for period_id in range(n_periods))\n",
        "\n",
        "\n",
        "# Setting the objective which is to minize the costs\n",
        "model_2_2.setObjective(base_cost + additional_cost + startup_cost + hydro_base_cost + hydro_startup)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cFQVee3ewkva",
        "outputId": "3c389fea-e249-49d8-8d63-153f087b6a27"
      },
      "source": [
        "model_2_2.optimize()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "Optimize a model with 85 rows, 75 columns and 215 nonzeros\n",
            "Model fingerprint: 0x47ce49dc\n",
            "Variable types: 25 continuous, 50 integer (20 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-03, 4e+03]\n",
            "  Objective range  [4e+00, 9e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [5e+00, 4e+04]\n",
            "Presolve removed 20 rows and 5 columns\n",
            "Presolve time: 0.00s\n",
            "Presolved: 65 rows, 70 columns, 190 nonzeros\n",
            "Variable types: 25 continuous, 45 integer (18 binary)\n",
            "\n",
            "Root relaxation: objective 1.011757e+06, 32 iterations, 0.00 seconds\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 1011757.14    0    7          - 1011757.14      -     -    0s\n",
            "H    0     0                    1121190.0000 1011757.14  9.76%     -    0s\n",
            "H    0     0                    1021935.0000 1011757.14  1.00%     -    0s\n",
            "H    0     0                    1019220.0000 1011757.14  0.73%     -    0s\n",
            "H    0     0                    1018470.0000 1011757.14  0.66%     -    0s\n",
            "     0     0 1012876.91    0   10 1018470.00 1012876.91  0.55%     -    0s\n",
            "     0     0 1012876.91    0   10 1018470.00 1012876.91  0.55%     -    0s\n",
            "H    0     0                    1014620.0000 1012876.91  0.17%     -    0s\n",
            "     0     2 1013040.67    0   10 1014620.00 1013040.67  0.16%     -    0s\n",
            "H    3     2                    1013630.0000 1013253.65  0.04%   2.3    0s\n",
            "\n",
            "Cutting planes:\n",
            "  MIR: 3\n",
            "\n",
            "Explored 10 nodes (75 simplex iterations) in 0.12 seconds\n",
            "Thread count was 2 (of 2 available processors)\n",
            "\n",
            "Solution count 5: 1.01363e+06 1.01462e+06 1.01847e+06 ... 1.12119e+06\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 1.013630000000e+06, best bound 1.013567142857e+06, gap 0.0062%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 230
        },
        "id": "VqpNWxW4wkxl",
        "outputId": "22721dc8-3906-4355-80cd-956229f056b7"
      },
      "source": [
        "def describe_model(model):\n",
        "  periods_info=[]\n",
        "  for p in range(n_periods):\n",
        "    periods_info.append({})\n",
        "    cost_additional=0\n",
        "    cost_base=0\n",
        "    total_gen_power=0\n",
        "    for t in range(n_power_types):\n",
        "      periods_info[p][\"Power type \"+str(t)]=gen_power[t,p].x\n",
        "      periods_info[p][\"N units on \"+str(t)]=n_units_on[t,p].x\n",
        "      total_gen_power+=gen_power[t,p].x\n",
        "      cost_additional+=per_mw_cost_by_type[t]*(gen_power[t,p].x - min_power_by_type[t]*n_units_on[t, p].x)*period_durations[p]\n",
        "      cost_base+=base_cost_by_type[t]*period_durations[p]*n_units_on[t, p].x\n",
        "\n",
        "    for h in range(n_hydro_units):\n",
        "      periods_info[p][\"Hydro unit \"+str(h)]=hydro_on[h,p].x\n",
        "      total_gen_power+=hydro_on[h,p].x*hydro_load[h]\n",
        "\n",
        "    periods_info[p][\"Total gen power\"]=total_gen_power\n",
        "    periods_info[p][\"Demand\"]=power_demand_by_period[p]\n",
        "    periods_info[p][\"Cost Base\"]=cost_base\n",
        "    periods_info[p][\"Cost Additional\"]=cost_additional\n",
        "\n",
        "    periods_info[p][\"Pumping\"]=pumping[p].x\n",
        "\n",
        "  print(\"Costs:\", model.getObjective().getValue())\n",
        "  return pd.DataFrame(periods_info)\n",
        "\n",
        "describe_model(model_2_2)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Costs: 1013630.0\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "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>Power type 0</th>\n",
              "      <th>N units on 0</th>\n",
              "      <th>Power type 1</th>\n",
              "      <th>N units on 1</th>\n",
              "      <th>Power type 2</th>\n",
              "      <th>N units on 2</th>\n",
              "      <th>Hydro unit 0</th>\n",
              "      <th>Hydro unit 1</th>\n",
              "      <th>Total gen power</th>\n",
              "      <th>Demand</th>\n",
              "      <th>Cost Base</th>\n",
              "      <th>Cost Additional</th>\n",
              "      <th>Pumping</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>10565.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>5250.0</td>\n",
              "      <td>3.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>15815.0</td>\n",
              "      <td>15000</td>\n",
              "      <td>118800.0</td>\n",
              "      <td>16080.0</td>\n",
              "      <td>815.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>14250.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>15750.0</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>30000.0</td>\n",
              "      <td>30000</td>\n",
              "      <td>106200.0</td>\n",
              "      <td>41850.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>10200.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>15750.0</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>25950.0</td>\n",
              "      <td>25000</td>\n",
              "      <td>212400.0</td>\n",
              "      <td>35100.0</td>\n",
              "      <td>950.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>21350.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>15750.0</td>\n",
              "      <td>9.0</td>\n",
              "      <td>1500.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>40000.0</td>\n",
              "      <td>40000</td>\n",
              "      <td>115200.0</td>\n",
              "      <td>84450.0</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>10200.0</td>\n",
              "      <td>12.0</td>\n",
              "      <td>15750.0</td>\n",
              "      <td>9.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>27350.0</td>\n",
              "      <td>27000</td>\n",
              "      <td>212400.0</td>\n",
              "      <td>35100.0</td>\n",
              "      <td>350.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Power type 0  N units on 0  ...  Cost Additional  Pumping\n",
              "0       10565.0          12.0  ...          16080.0    815.0\n",
              "1       14250.0          12.0  ...          41850.0      0.0\n",
              "2       10200.0          12.0  ...          35100.0    950.0\n",
              "3       21350.0          12.0  ...          84450.0      0.0\n",
              "4       10200.0          12.0  ...          35100.0    350.0\n",
              "\n",
              "[5 rows x 13 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 353
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0edbukiSo4uG"
      },
      "source": [
        "### Question 2.3.3\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TwkFESlwGcZt",
        "outputId": "55528765-8624-475b-8082-061d07d4768c"
      },
      "source": [
        "for period_id in range(n_periods):\n",
        "    demand_reserve = (1+reserve_factor)*power_demand_by_period[period_id] - sum(hydro_load)\n",
        "    max_potential = sum(hydro_load)\n",
        "    for type_id in range(n_power_types):\n",
        "        max_potential += max_power_by_type[type_id]*n_units_on[type_id, period_id].x\n",
        "\n",
        "\n",
        "    print(f\"Period {period_id} : max_potential= {max_potential}, demand={power_demand_by_period[period_id]}, effective reserve factor={(max_potential-power_demand_by_period[period_id])/power_demand_by_period[period_id]}\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Period 0 : max_potential= 31550.0, demand=15000, effective reserve factor=1.1033333333333333\n",
            "Period 1 : max_potential= 42050.0, demand=30000, effective reserve factor=0.40166666666666667\n",
            "Period 2 : max_potential= 42050.0, demand=25000, effective reserve factor=0.682\n",
            "Period 3 : max_potential= 46050.0, demand=40000, effective reserve factor=0.15125\n",
            "Period 4 : max_potential= 42050.0, demand=27000, effective reserve factor=0.5574074074074075\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CRrbZfqbIy64"
      },
      "source": [
        "We never have a period where \"Diff reserve - demand\" is at 0.15 (the initial value of reserve_factor), so there is no reserve constraint active.\n",
        "\n",
        "The maximum reserve factor satisfied by the current optimal solution is $0.15125$.\n",
        "\n",
        "After rerunning the model with reserve_factor$=0.20$, we find that the cost went from $1.01363\\times 10^6$ to $1.01413\\times 10^6$ which is a difference of $500$€."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XKlgBsKdo41R"
      },
      "source": [
        "\n",
        "\n",
        "### Question 2.3.4\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FmPmGvmxLvAn",
        "outputId": "5c6502fd-5705-4e6a-a4eb-3c6c220c33d1"
      },
      "source": [
        "# Careful, running this will overwrite the run of the model_2_2 and the following questions relying on the model_2_2 will not carry the same result\n",
        "model_2_2_relax=model_2_2.relax()\n",
        "model_2_2_relax.optimize()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "Optimize a model with 85 rows, 75 columns and 215 nonzeros\n",
            "Model fingerprint: 0x716f752c\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-03, 4e+03]\n",
            "  Objective range  [4e+00, 9e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [5e+00, 4e+04]\n",
            "Presolve removed 20 rows and 5 columns\n",
            "Presolve time: 0.01s\n",
            "Presolved: 65 rows, 70 columns, 190 nonzeros\n",
            "\n",
            "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
            "       0    9.2877000e+03   6.422496e+03   0.000000e+00      0s\n",
            "      30    1.0116196e+06   0.000000e+00   0.000000e+00      0s\n",
            "\n",
            "Solved in 30 iterations and 0.01 seconds\n",
            "Optimal objective  1.011619643e+06\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1otum96YL792",
        "outputId": "dc713df2-d785-4201-ea4e-3059e331cd12"
      },
      "source": [
        "meet_demand_constr"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{0: <gurobi.Constr R45>,\n",
              " 1: <gurobi.Constr R46>,\n",
              " 2: <gurobi.Constr R47>,\n",
              " 3: <gurobi.Constr R48>,\n",
              " 4: <gurobi.Constr R49>}"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 356
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VpsxSWHfL6L2",
        "outputId": "2f69a4d7-c295-42ff-d750-bc65532d90d5"
      },
      "source": [
        "for i in range(45,50):\n",
        "    print(f\"Period {i-45}: Dual values demand satisfaction: {model_2_2_relax.getConstrs()[i].getAttr('Pi') }\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Period 0: Dual values demand satisfaction: 11.142857142857144\n",
            "Period 1: Dual values demand satisfaction: 6.0\n",
            "Period 2: Dual values demand satisfaction: 10.714285714285715\n",
            "Period 3: Dual values demand satisfaction: 6.0\n",
            "Period 4: Dual values demand satisfaction: 11.160714285714286\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NBJPDO7Go46m"
      },
      "source": [
        "### Question 2.3.5\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N_JHgwJZeUb-"
      },
      "source": [
        "Not done."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2fsx_13AobYS"
      },
      "source": [
        "#3: Stochastic Optimization Model: handling the random power demand"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bKA6C4R8ok8n"
      },
      "source": [
        "##Question 3.1\n",
        "\n",
        "Let's reformulate the probability constraint statistically with $N_{scen}$ scenarios $demand_p^i\\in \\mathbb{R}^5$ of respective probabilities $(Pr_i)_{i \\in [\\![ 1, N_{scen} ]\\!]}$.\n",
        "\n",
        "We have $\\sum\\limits_{i=1}^{N_{scen}}Pr_i = 1$\n",
        "\n",
        "Given a solution caracterised by:\n",
        "\n",
        "$\\text{output}_{t,p}, \\text{hydro}_{load h}, \\text{hydro}_{h,p}, \\text{pumping}_p \\text{ for }(t,p) \\in \\text{Types} \\times \\text{Periods} \\text{ and } h \\in \\text{HydroUnits}$\n",
        "\n",
        "\n",
        "For each $i \\in [\\![ 1, N_{scen} ]\\!]$, let's define:\n",
        "\n",
        "\n",
        "$x_i=  \\begin{cases}\n",
        "               Pr_i \\text{  if  } \\forall p \\in \\text{Periods},\\sum\\limits_{t \\in \\text{Types}}\\text{output}_{t,p}+\n",
        "\\sum\\limits_{h \\in \\text{HydroUnits}}\\text{hydro}_{load h} \\times \\text{hydro}_{h,p} \\geq \\text{demand}_p^i + \\text{pumping}_p \\\\\n",
        "               0 \\text{ otherwise }\n",
        "               \n",
        "            \\end{cases}\n",
        "$\n",
        "\n",
        "The given solution must respect the following condition:\n",
        "\n",
        "$\\sum\\limits_{i =0 }^{N_{scen}}x_i \\geq p_{level}$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YnipwvMiooeD"
      },
      "source": [
        "## Question 3.2"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "P0lu-XglEvRc",
        "outputId": "a46e8247-b462-4428-e592-ce253d2c1c99"
      },
      "source": [
        "n_scen=10\n",
        "cov=[[265000, 475000, 650000, 700000],\n",
        "[475000, 1300000, 750000, 1650000],\n",
        "[650000, 750000, 2500000, 2000000],\n",
        "[700000, 1650000, 2000000, 4700000]]\n",
        "\n",
        "np.random.seed(0)\n",
        "\n",
        "def generate_scenarios(demand, cov, n_scen):\n",
        "  random_demands=list(np.random.multivariate_normal(demand[1:], cov, n_scen))\n",
        "\n",
        "  # Here we fix the first demand value for all of the outputs\n",
        "  for i in range(len(random_demands)):\n",
        "    random_demands[i]=[demand[0]]+list(random_demands[i])\n",
        "  return np.array(random_demands)\n",
        "\n",
        "\n",
        "generate_scenarios(power_demand_by_period, cov, n_scen)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[15000.        , 29171.49941369, 24404.01661527, 37488.65957728,\n",
              "        23207.4126291 ],\n",
              "       [15000.        , 29625.14387839, 23857.55996075, 38702.19021692,\n",
              "        22299.54764945],\n",
              "       [15000.        , 29837.46683819, 25349.71768921, 39741.23249569,\n",
              "        27372.48601415],\n",
              "       [15000.        , 29731.84136744, 24721.84903329, 38952.02249572,\n",
              "        25333.98381587],\n",
              "       [15000.        , 29592.52254608, 23881.38500525, 38368.75144285,\n",
              "        23689.6250807 ],\n",
              "       [15000.        , 31244.82135745, 27916.84212642, 42392.20287483,\n",
              "        32369.99112224],\n",
              "       [15000.        , 29324.21841524, 22742.46956666, 38719.20723368,\n",
              "        21510.84566532],\n",
              "       [15000.        , 29129.2430066 , 24204.48739933, 36656.89621975,\n",
              "        24524.04852357],\n",
              "       [15000.        , 30574.61575014, 25002.35649817, 43098.20794966,\n",
              "        27928.90218964],\n",
              "       [15000.        , 29249.4016992 , 23974.74105569, 37290.24121146,\n",
              "        25195.28798648]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 321
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_J-u9Moforkk"
      },
      "source": [
        "## Question 3.3"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TKfVKDeLN3wH"
      },
      "source": [
        "scenarios=generate_scenarios(power_demand_by_period, cov, n_scen=300)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 623
        },
        "id": "VlM6OOgpOgdl",
        "outputId": "08f376b4-78af-48e1-835a-301fe4c18149"
      },
      "source": [
        "# Plot\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.interpolate import make_interp_spline, BSpline\n",
        "plt.rcParams['figure.figsize'] = [15, 5]\n",
        "\n",
        "def plot_scenario(ids_scen_to_plot):\n",
        "  x=list(range(len(demand)))\n",
        "  xnew = np.linspace(x[0], x[-1], 300)\n",
        "  spl = make_interp_spline(x, scenarios[id_scen_to_plot], k=3)  # type: BSpline\n",
        "  power_smooth = spl(xnew)\n",
        "  plt.plot(xnew, power_smooth)\n",
        "  plt.show()\n",
        "\n",
        "plot_scenario(1)\n",
        "plot_scenario(100)\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAEvCAYAAAD4uAgWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3iUVfrG8e9J7wmpJCQhdAgdQldEWRQriqKggh277uqubdfV1W2u7lrXgh0bIGBXbIiC1ARClV5SICEhpJE+c35/ZNZl/SFSEt6U+3Ndc5k5887MPQTGeeac9znGWouIiIiIiIi0DF5OBxAREREREZGGoyJPRERERESkBVGRJyIiIiIi0oKoyBMREREREWlBVOSJiIiIiIi0ICryREREREREWhCfIz3QGOMNpAO51tpzjDGvAacAJZ5DrrTWZhpjDPAkcBZQ4Rlf6XmMK4A/eI7/s7X2dc/4QOA1IBD4FLjd/sLeDtHR0TYlJeVI44uIiIiIiLQoGRkZhdbamJ+OH3GRB9wO/ACEHTT2O2vt7J8cdybQxXMZAjwHDDHGRAIPAGmABTKMMR9aa/d7jrkOWEZ9kTcW+OxwYVJSUkhPTz+K+CIiIiIiIi2HMWbXocaPaLmmMSYROBt46QgOHwdMt/WWAhHGmHjgDOBLa22Rp7D7EhjruS3MWrvUM3s3HTj/SHKJiIiIiIjI/zrSc/KeAO4C3D8Z/4sxZo0x5nFjjL9nrB2QfdAxOZ6xw43nHGJcREREREREjtIvFnnGmHOAvdbajJ/cdC/QHRgERAJ3N3y8/5dlqjEm3RiTXlBQ0NhPJyIiIiIi0uwcyUzeCOA8Y8xOYAZwmjHmTWvtHs+SzGrgVWCw5/hcIOmg+yd6xg43nniI8f/HWjvNWptmrU2Lifl/5xeKiIiIiIi0er9Y5Flr77XWJlprU4CJwHxr7eWec+nwdNM8H1jnucuHwBRTbyhQYq3dA3wOnG6MaWOMaQOcDnzuua3UGDPU81hTgA8a+HWKiIiIiIi0CkfTXfOn3jLGxAAGyARu8Ix/Sv32CVup30LhKgBrbZEx5mFghee4h6y1RZ6fb+K/Wyh8xi901hQREREREZFDM7+wHV2TlZaWZrWFgoiIiIiItFbGmAxrbdpPx4+0u6aIiIiIiIg0AyryREREREREWpDjOSdPRERERKRZqap1UXSg5sdLeXUdLrfFbS0ud/3FGEOIvw9hAT6EBvgSGuBDRJAv4YG+1PcJFGnaVOSJiIiISIuSX1rFprwydu07wI7CCnbuO8DOfQfYU1xFZa3rmB832M+bpMggkiKDSI4Mon1UEKnxYfRMCCfQz7sBX4HI8VGRJyIiIiLNVnFFDauyi1mbU8KanGLW5JSwt6z6x9sDfb1pHxVEt7hQTusWS2SIH5FBfrQJ9iMy2I/QAB98vAzGGLyNwdvL4LaWsqo6z6WW8uo6ig7UkLO/kpz9Fezad4BFWwp/LBi9vQxdYkPokxhO36QIhneKpkN0sFN/JCIq8kRERESk+aiqdbFiZxGLthayeOs+1u0uwVowBjpGB3NS52h6J4bTIz6MDtHBxIb6N8oSS2st+aXVrM39b3H55YZ8ZqXnANA+KohTusZwStcYhnWKIshPH7vlxNEWCiIiIiLSpO0treLzDfl8sT6PZTuKqKlz4+tt6J/chhGdohncIZJe7cIIDfB1NKe1ll37Kli4pYAFmwpYvG0flbUu/Ly9GNUthvED2nFq91j8fbS0UxrGz22hoCJPRERERJqc7KIK5q3LY976PFZm7cda6BgTzOjusYzoXF/YNfXZseo6Fyt27Ofrjfl8tHoPheXVhAf6cnafeMb3b8fA9m3UyEWOi4o8EREREWnSDlTX8enaPczOyGHZjiIAeiaEMbZnW8b2akuXuFCHEx67OpebRVsLeW9VLp+vz6Oq1k3XuBCuPakj5/VLIMBXs3ty9FTkiYiIiEiTY61l6fYi3s3I5rO1eVTWuugQHcxFAxM5r28CSZFBTkdscOXVdXy6Zg+vfL+DjXllRIf4MWVYCpcNSSYqxN/peNKMqMgTERERkSajoqaO91bl8vrinWzOLyfU34dz+sZz0cBEBiS3jmWM1loWb9vHSwu3882mAvx9vJg4KIlbTutCTKiKPfllKvJERERExHHZRRVMX7KTmSuyKa2qo2dCGFcMT+HcPgmteq+5rXvLePG7HcxemYO/jxfXnNSB60Z2JMzhZjLStKnIExERERHHbNhdyrMLtvLp2j0YYzizV1uuHJ6i5iM/sb2gnH99uZmP1+whIsiXm0Z1YsqwFJ2zJ4ekIk9ERERETrj0nUU8u2Ab8zfuJcTfh8uGJnPV8A60DQ9wOlqTti63hH98vonvNheQFBnIQ+N6cWq3WKdjSROjIk9ERERETpil2/fxry83s3xHEZHBflw9IoXJw1IID9Tyw6Px/dZC/vjBOrYVHOCs3m154NyexIWpQJZ6KvJEREREpNGtzi7msS82sXBLIXFh/lw/shMTByc1+T3tmrLqOhcvfredp+dvxdfbi9+e3pXJw1Lw9tIy19ZORZ6IiIiINJrN+WX884tNfL4+nzZBvtx8amcuH9pe55I1oF37DvCH99excEshfZMieOKSfnSIDnY6ljhIRZ6IiIiINLi9pVU89sUm3s3IIdjPh+tO7sjVJ6UQqq6QjcJay0dr9nD/++uodbl54NxULk5LUvOaVurnijzNm4uIiIjIUauqdfHSwu08u2AbtS43V4/owC2ndqZNsJ/T0Vo0Ywzn9U1gUEob7py1mrvnrOWbjQX8bXxv/dnLjzSTJyIiIiJHzFrLh6t38495m8gtruSMnnHcc2YPLRt0gNtteWnRdh79fBORwX786+J+jOgc7XQsOYF+bibPy4kwIiIiItL8bMor45IXlnL7jEwignx557qhvDA5TQWeQ7y8DFNHduL9m0cQGuDL5JeX8dyCbTTXSRxpOFquKSIiIiKHVV5dx5NfbeaV73cSFuDD38f3ZkJakro7NhE9E8L58JYR3DV7DY/M28i63SU8elEfdTRtxfSbFxEREZFDstby6do8Hvp4Pfml1UwanMRdZ3TXuV9NUJCfD09P6k/vduE8Mm8j2/aWM21yGslRQU5HEwdouaaIiIiI/D+5xZVc9doKbn57JVHB/sy9aTh/G99HBV4TZozh+lM68dpVg9lTUsW5zyxi4ZYCp2OJA1TkiYiIiMiP3G7L64t3cvq/vmXZ9iLuPyeVD28ZwYDkNk5HkyM0smsMH94ygvjwAK58dQWzVmQ7HUlOMC3XFBEREREAtu4t4+45a8nYtZ+Tu0Tz1wt6kxSp5X7NUfuoYObcOJwb3szgrjlryC+t4pbTOms/vVZCRZ6IiIhIK+dyW174bhtPfLmFQD9vHpvQlwsHtFNB0MwF+/vw8hWDuHvOGv755Wbyy6r403m91DCnFVCRJyIiItKK7dp3gDtmrSZj137O7NWWh8b1IibU3+lY0kD8fLz454S+xIb588K32yksq+GJif0I8PV2Opo0IhV5IiIiIq2QtZZ3lmfz50824O1lePySvpzfT7N3LZGXl+HeM3sQGxrAwx9vYMrLy3n5yjRCA3ydjiaNREWeiIiISCuzt7SKu+es4ZtNBYzoHMWjF/UlISLQ6VjSyK45qQMxof7cMTOTK15ZzutXD1ah10KpyBMRERFpRT5bu4f73ltLRY2LB85N5YphKXjpHK1W47y+Cfh5G255exVTPIVemAq9FkdbKIiIiIi0AmVVtdwxM5Mb31pJUmQQn9x2MleN6KACrxUa2yueZy4dwNqcEqa8vJzSqlqnI0kDU5EnIiIi0sKtySnm7KcW8cHq3dw+ugtzbhxO59gQp2OJg8b2asu/LxvAutwSJr+8nJJKFXotiYo8ERERkRbKWsvLi3Zw4XOLqXO5mTl1KL8Z0xVfb30EFDijZ1uevWwAG3aXMOXlZZRpRq/FOOJ/4cYYb2PMKmPMx57rHYwxy4wxW40xM40xfp5xf8/1rZ7bUw56jHs945uMMWccND7WM7bVGHNPw708ERERkdZp/4EarpuezsMfb2BUt1g+vf1k0lIinY4lTczpPdvy7GUDWb+7lKnTM6iqdTkdSRrA0XyNczvww0HXHwEet9Z2BvYD13jGrwH2e8Yf9xyHMSYVmAj0BMYCz3oKR2/g38CZQCowyXOsiIiIiByDFTuLOOuphXy3uZAHz01l2uSBRAT5OR1LmqgxqXH846I+LNm+j9/MzMTltk5HkuN0REWeMSYROBt4yXPdAKcBsz2HvA6c7/l5nOc6nttHe44fB8yw1lZba3cAW4HBnstWa+12a20NMMNzrIiIiIgcBZfb8sz8LUycthR/Hy/m3jScK0d00N538ovGD0jkD2f34LN1edz/wTqsVaHXnB3pFgpPAHcBoZ7rUUCxtbbOcz0HaOf5uR2QDWCtrTPGlHiObwcsPegxD75P9k/GhxwqhDFmKjAVIDk5+Qiji4iIiLR8+8qr+fXMTBZuKWRcvwT+ckFvQvy1W5YcuWtP7khheQ3Pf7uN6BB/7hjT1elIcox+8V++MeYcYK+1NsMYM6rxI/08a+00YBpAWlqavl4QERERATKzi7npzQwKD9TwyIW9uTgtSbN3ckzuHtuNfeXVPPX1FqJD/JgyLMXpSHIMjuTrnRHAecaYs4AAIAx4Eogwxvh4ZvMSgVzP8blAEpBjjPEBwoF9B43/x8H3+blxEREREfkZ1lreWZ7Ngx+uJzbMn7k3DqdXu3CnY0kzZozhb+N7s7+ilgc+XE9cWABn9GzrdCw5Sr94Tp619l5rbaK1NoX6xinzrbWXAd8AF3kOuwL4wPPzh57reG6fb+sX9X4ITPR03+wAdAGWAyuALp5unX6e5/iwQV6diIiISAtVVevirtlruO+9tQzrFMVHt5ykAk8ahI+3F89c2p8+7cL5zcxM1u8ucTqSHKXj2STlbuAOY8xW6s+5e9kz/jIQ5Rm/A7gHwFq7HpgFbADmATdba12emcBbgM+p7945y3OsiIiIiBxCdlEFFz63mHczcrhtdBdeuXIQbYLVPVMaToCvNy9OSSMswJfrXk9nb1mV05HkKJjm2jknLS3NpqenOx1DRERE5IT6ZtNefj0jE2stj1/Sj9E94pyOJC3YutwSLnp+MT3iw3jnuqEE+Ho7HUkOYozJsNam/XT8eGbyREREROQEsbZ+e4SrX1tBQkQgH916kgo8aXS92oXz+MX9WJVVzD1z1mhrhWZCRZ6IiIhIE1dZ4+KWd1bx2BebOa9vAnNvHE77qGCnY0krcWbveO4c05X3M3fz7IJtTseRI6DNU0RERESasNziSqZOT2fDnlLuObM714/sqO0R5IS75bTObNlbzqOfb6JbXCi/StUsclOmmTwRERGRJmrFziLGPbOIrH0VvHLFIG44pZMKPHGEMYZ/XNSHXu3CuGNWJtlFFU5HksNQkSciIiLSBM1YnsWlLy4lNMCX924ewandY52OJK1cgK83z146EIAb38qgqtblcCL5OSryRERERJqQOpebBz9czz1z1zKsUzTv3zSCzrEhTscSASA5Koh/XtyPdbmlPPTxBqfjyM9QkSciIiLSRJRU1HLFq8t5bfFOrju5A69ckUZ4kK/TsUT+x5jUOG44pRNvL8ti7socp+PIIajxioiIiEgTsGvfAa56bQU5RZU8NqEvFw1MdDqSyM/67eldWZm1n9+/t46eCeF0axvqdCQ5iGbyRERERBy2YmcR5//7e/YfqOHNa4eowJMmz8fbi2cm9SfY34cb38qgvLrO6UhyEBV5IiIiIg76IDOXy15cRpsgP967aQSDO0Q6HUnkiMSGBfD0pP7sLDzAHz9Y53QcOYiKPBEREREHWGt54qvN3D4jkwHtI5h703BSorXBuTQvwzpFccupnZm7MpeP1+x2Oo54qMgTEREROcGq61z8ZmYmT3y1hQsHJDL96iFEBPk5HUvkmNw6ugt9kyK4b+5adhdXOh1HUJEnIiIickIVHajh8peW8X7mbn53Rjcem9AHPx99JJPmy9fbiycv6Ued23LnrNW43dbpSK2e3lFERERETpBtBeVc8Oz3rM4p4elJ/bn51M4YY5yOJXLcUqKDefDcnizZvo8XF253Ok6rpyJPRERE5ARYvK2QC/79PeVVdbxz3VDO7ZvgdCSRBjUhLZGxPdvy2BebWJdb4nScVk1FnoiIiEgjm7syhykvLyc2LID3bx7BwPZtnI4k0uCMMfxtfG8ig/24fcYqKmtcTkdqtVTkiYiIiDQSay3//mYrd8xazaCUSObcOJykyCCnY4k0mjbBfjw2oS/bCg7wyLyNTsdptVTkiYiIiDSCOpebP7y/jkc/38S4fgm8fvVgwgN9nY4l0uhO7hLDFcPa8/qSnaTvLHI6TqukIk9ERESkgVXWuLjhzQzeWpbFDad04vGL+6mDprQqd43tTkJ4IHfNXkNVrZZtnmh6txERERFpQPvKq5n04lK+3riXh8b15J4zu+PlpQ6a0roE+/vwyIV92F54gCe+2uJ0nFZHRZ6IiIhIA9m17wAXPreYH/aU8vzlA5kyLMXpSCKOOalLNBMHJTHtu22szi52Ok6roiJPREREpAFkZhcz/tnFFFfW8vZ1QzijZ1unI4k47r6zexAbGsDdc9ZQU+d2Ok6roSJPRERE5Dh9/UM+k6YtJcjfmzk3Dmdg+0inI4k0CWEBvvzlgl5szCvj2QVbnY7TaqjIExERETkOby/L4rrp6XSODWHOjcPpFBPidCSRJmV0jzgu6N+OZ+Zv5Yc9pU7HaRVU5ImIiIgcA2st//pyM/e9t5aRXWOYMXUosaEBTscSaZL+eE4qEUG+3Dt3LW63dTpOi6ciT0REROQoudyWP7y/jqe+3sJFAxN5cUoawf4+TscSabLaBPvx+7N7kJldzIwV2U7HafFU5ImIiIgcheo6F7e+s/LHPfAevagPvt76SCXyS87v144hHSJ5ZN5G9pVXOx2nRdM7koiIiMgRKquq5cpXVvDp2jz+cHYP7jmzO8ZoDzyRI2GM4c/n9+JAdR1//2yj03FaNBV5IiIiIkegoKyaidOWsmJnEf+6uC/XntzR6UgizU6XuFCuG9mRdzNyWLGzyOk4LZaKPBEREZFfkF1UwYTnF7OtoJwXr0hj/IBEpyOJNFu3ntaZdhGB/OG9ddS6tHdeY1CRJyIiInIYG3aXMv65xeyvqOWta4dyardYpyOJNGtBfj48cG4qm/LLeO37nU7HaZFU5ImIiIj8jGXb93HJC0vwNoZ3bxjGwPZtnI4k0iKMSY1jdPdYHv9qM7uLK52O0+KoyBMRERE5hC/W5zH5leXEhvkz56bhdI0LdTqSSIthjOHB83riclv+8skPTsdpcX6xyDPGBBhjlhtjVhtj1htj/uQZf80Ys8MYk+m59POMG2PMU8aYrcaYNcaYAQc91hXGmC2eyxUHjQ80xqz13OcpozZVIiIi4qCZK7K44c0MUuPDePeG4bSLCHQ6kkiLkxQZxI2jOvHJ2j0s36EmLA3pSGbyqoHTrLV9gX7AWGPMUM9tv7PW9vNcMj1jZwJdPJepwHMAxphI4AFgCDAYeMAY8581D88B1x10v7HH/cpEREREjpK1lmcXbOXuOWs5qUsMb107hMhgP6djibRY14/sRHx4AA99vB632zodp8X4xSLP1iv3XPX1XA73GxgHTPfcbykQYYyJB84AvrTWFllr9wNfUl8wxgNh1tql1loLTAfOP47XJCIiInLU3G7Lwx//wD/mbeK8vgm8NCWNYH8fp2OJtGiBft7cPbY763JLmbMyx+k4LcYRnZNnjPE2xmQCe6kv1JZ5bvqLZ0nm48YYf89YOyD7oLvneMYON55ziHERERGRE6Kmzs0dszJ55fsdXDUihScu6Yefj1oXiJwI5/VNoF9SBP/4fBMHquucjtMiHNG7l7XWZa3tByQCg40xvYB7ge7AICASuLvRUnoYY6YaY9KNMekFBQWN/XQiIiLSClTU1HHd9HTez9zN787oxh/PScXLS+0BRE4ULy/DH89NpaCsmucWbHM6TotwVF9RWWuLgW+AsdbaPZ4lmdXAq9SfZweQCyQddLdEz9jhxhMPMX6o559mrU2z1qbFxMQcTXQRERGR/2f/gRoufXEZC7cU8Pfxvbn51M6o/5vIiTcguQ3j+iUwbeF2cvZXOB2n2TuS7poxxpgIz8+BwBhgo+dcOjydMM8H1nnu8iEwxdNlcyhQYq3dA3wOnG6MaeNpuHI68LnntlJjzFDPY00BPmjYlykiIiLyv3YXVzLhhSVs2FPKc5cPZOLgZKcjibRqd4/tjpeBv3+20ekozd6RnE0cD7xujPGmviicZa392Bgz3xgTAxggE7jBc/ynwFnAVqACuArAWltkjHkYWOE57iFr7X96pd4EvAYEAp95LiIiIiKNYkt+GVNeWU55VR3Trx7M0I5RTkcSafUSIgKZOrITT329hSuHF5GWEul0pGbL1De0bH7S0tJsenq60zFERESkmVmZtZ+rX1uBr7cXr181mNSEMKcjiYhHRU0dpz32LXHhAbx/03Atn/4FxpgMa23aT8fVNkpERERajW827eWyF5cRHujLnBuGq8ATaWKC/Hy44/SurM4uZt66PKfjNFsq8kRERKRVeH9VLte9nk7HmGBm3zCc5KggpyOJyCGM79+OzrEhPPrFJupcbqfjNEsq8kRERKTFe3nRDn49M5NBKZHMmDqUmFD/X76TiDjCx9uL353Rje0FB5idoQ3Sj4WKPBEREWmxrLU8Mm8jD3+8gTN7teXVqwYRGuDrdCwR+QWnp8bRPzmCJ77aQlWty+k4zY6KPBEREWmR6lxu7p6zhucWbOPSIck8c+kAAny9nY4lIkfAGMPdY7uTV1rF64t3Oh2n2VGRJyIiIi1OVa2LG99ayaz0HG4b3YW/nN8Lby916RNpToZ2jGJUtxieXbCNkspap+M0KyryREREpEUpqahl8svL+OqHfP50Xk/uGNNVbdhFmqnfndGNkspanv92m9NRmhUVeSIiItJi5JdWcfELS8jMLubpSf25YniK05FE5Dj0TAhnXL8EXv1+B/mlVU7HaTZU5ImIiEiLsL2gnPHPLiZnfwWvXTWYc/okOB1JRBrAnWO6UeeyPPn1FqejNBsq8kRERKTZW51dzEXPL6Gq1sWMqcMY0Tna6Ugi0kCSo4KYNDiZWSuyyS6qcDpOs6AiT0RERJq17zYXMOnFpQT7ezP7xuH0Tgx3OpKINLCbT+2Ml5fh399sdTpKs6AiT0RERJqtDzJzufq1FbSPCmbODcPpEB3sdCQRaQRtwwO4dHAyszNyNJt3BFTkiYiISLP0yqId3D4jk4Ht2zDz+qHEhgU4HUlEGtENp3TCy8vwzHzN5v0SFXkiIiLSrFhr+ce8jTz08QbG9mzL61cPJizA1+lYItLI/jObN2dlDln7NJt3OCryREREpNmoc7m5e84anl2wjUmDk/n3ZQMI8PV2OpaInCA3jvLM5n2jTpuHoyJPREREmoWqWhc3vLmSWek53Da6C3+9oBfeXtrkXKQ1iQv7z2xermbzDkNFnoiIiDR5JRW1TH55GV9vzOfhcT25Y0xXjFGBJ9Ia3TSqEz6azTssFXkiIiLSpOWVVHHxC0tYnV3CM5MGMHlYitORRMRBsWEBXDqkfjZv174DTsdpklTkiYiISJO1dW85Fz63mNziSl67ahBn94l3OpKINAE3nuKZzVOnzUNSkSciIiJN0qqs/Ux4fjHVdS5mTB3K8M7RTkcSkSbiP7N5c1flat+8Q1CRJyIiIk3OlxvymfTiUsICfZl9w3B6tQt3OpKINDFTR3bEy8CLC7c7HaXJUZEnIiIiTcobS3dx/RvpdIsLZc6Nw0mJDnY6kog0QfHhgYzvn8jMFdkUlFU7HadJUZEnIiIiTcJ/Njm///11nNotlnemDiU6xN/pWCLShE09pSM1Ljevfr/D6ShNioo8ERERcVxNnZs7Z63+cZPzFyYPJMjPx+lYItLEdYoJ4cxebXljyS5Kq2qdjtNkqMgTERERR5VV1XL1ayuYuyqXO8d05a8X9MLHWx9RROTI3DSqM2XVdby1NMvpKE2G3kFFRETEMfmlVUx4fglLt+/j0Yv6cOvoLtrkXESOSq924ZzcJZqXF+2gqtbldJwmQUWeiIiIOGJzfhkX/Pt7sosqeOXKQUxIS3I6kog0UzeN6kxheTXvZuQ4HaVJ0GJ3kSbAWkuNy02ty1Jb56bW5abObfHxMnh7GXy8vfD1Nvh6e+GrJUwi0gIs3b6PqdPT8ff1Zub1w7RFgogcl6EdI+mfHMG077YxaVBSq1/yrSJPpJG43ZaC8mpyiyvZ/eOlir1lVRRX1FJSWfvjf8ur6474cf19vAgN8CUswIfQAB/CAn2JDvEnNtSfmFB/YsMCiAv1JzEyiLZhAXh7admTiDQtH67ezW9nrSYpMpDXrhpMUmSQ05FEpJkzxnDjKZ2Y+kYGH6/Zw/n92zkdyVEq8kSOk7WW3SVV/LC7lC17y9mSX8aWveVs3VtO5U/WhYf6+xAb5k+bID/ahgXQrW0o4YG+hAb44u9TP1vn5+2Fr48X3sbgspY6l/1xZq+2zk15dR2lVXWUVdVSVlVHSWUt2wsOUFBWTY3L/T/P5+ttSGwTRFJkEO0jg+gYE0zXuFC6xIUQE+Kv815E5ISy1vLM/K3888vNDE6JZNqUgUQE+TkdS0RaiF/1iKNLbAjPLdjGuH4Jrfpzjoo8kaNUUVPHmpwSVmUVk5m9n1VZxew9aAPOtmEBdIkLYdLgZDrEBJMYEUhCRCDxEQGEBfg2Wi5rLaWVdewtqyKvtIqc/ZVkFVWQta+CrKIKVmXtp6zqvzOGbYJ86RIXSmp8GL3bhdM7MZxOMSGa+RORRlFd5+LeOWuZuyqXC/q34+8X9sbfx9vpWCLSgnh5GW44pRN3vruabzcXMKpbrNORHGOstU5nOCZpaWk2PT3d6RjSCrjclvW7S1i4pZCFWwrI2LWfWlf9v5uUqCD6J7ehX1IEvdqF0yUupFELueNhraWgrJrN+WX09UIAACAASURBVOVszi9jy94yNuWV8cOesh9nHAN9vUlNCKN/UgSDOkSS1r4NUdqIWESO0/4DNVz/ZgbLdxRxx5iu3Hpa51b9DbuINJ6aOjcnPTKfbm1DeeOaIU7HaXTGmAxrbdpPxzWTJ3IIFTV1LNhUwGfr8li4pYDiivrNNVPjw7h6RAeGdIykX1IbIoObzzIjYwyxYQHEhgVwUpfoH8ddbsv2gnLW5pbUX3JKmL50Fy8t2gFAp5hgBneIJK19JIM7RJLYJlAfzkTkiO0oPMDVr60gd38lT07sx7h+rfs8GRFpXH4+XlwxPIVHP9/ExrxSurcNczqSIzSTJ+JRVlXL/I17+WxtHgs276Wq1k1UsB+jusUysms0IzpHE91KZrWq61yszSlhxc79rNhZRPrOIko9Sz3bhgUwonM0I7tGc3KXmGZV6IrIibVs+z6ufzMDL2OYNnkgaSmRTkcSkVaguKKGYX+bzzl94nl0Ql+n4zSqn5vJ+8UizxgTAHwH+FM/8zfbWvuAMaYDMAOIAjKAydbaGmOMPzAdGAjsAy6x1u70PNa9wDWAC7jNWvu5Z3ws8CTgDbxkrf37L70gFXnSENxuy9Id+5i1IpvP1uVRXecmNtSfsb3acmaveAZ3iNQ5atT/OW3KLyN9ZxFLtxexaGshJZW1GAO924UzsksMI7vG0D85Qls8iAgAczJyuGfuGpIjg3jlykG0jwp2OpKItCL3v7+OmSuyWXTPqcSGBjgdp9EcT5FngGBrbbkxxhdYBNwO3AHMtdbOMMY8D6y21j5njLkJ6GOtvcEYMxG4wFp7iTEmFXgHGAwkAF8BXT1PsxkYA+QAK4BJ1toNh8ulIk+Ox+7iSuZk5PBuRg5ZRRWEBvgwrl8C5/drx4DkNnipsDssl9uyJqeYhVsK+W5zAauyi3G5LaH+PozoHM2Y1DhO6x5LG83yibQ61loe/3IzT83fyvBOUTx32UDCg5rmucoi0nJtLyhn9L++5dZTO3PH6d2cjtNojvmcPFtfBZZ7rvp6LhY4DbjUM/468CDwHDDO8zPAbOAZT6E4Dphhra0GdhhjtlJf8AFstdZu9wSd4Tn2sEWeyLFYmbWflxfu4LN1e3BbGNYxijvGdGVsr7YE+KrL25Hy9jL0T25D/+Q23Da6CyWVtSzZVsi3mwv5ZuNe5q3Pw9vLMDglktN7xjEmNY7ENtoHS6Slq6p1cdfsNXy4ejeXpCXx5wt6aXZfRBzRMSaE0d3jeGPpLm46tXOr+5x3RI1XjDHe1C/J7Az8G9gGFFtr/9OPPQf4z5nU7YBsAGttnTGmhPolne2ApQc97MH3yf7JeMtvhSMnTJ3LzRcb8nlp4XZWZhUTGuDDdSd35LIh7UmOUuHREMIDfRnbK56xveKx1rI2t4Qv1ufzxYY8/vTRBv700QZS48M4vWccZ/WOp2tcqNORRaSB7S2t4vo3M1iVVcw9Z3bn+pEd1aRJRBx17ckd+OqHfOauzOXSIclOxzmhjqjIs9a6gH7GmAjgPaB7o6b6GcaYqcBUgOTk1vWLkqNX63IzJyOHfy/YSnZRJe2jgnjw3FQmpCUR7K/Gso3FGEOfxAj6JEbw2zO6saPwAF9uyOPLDfk8+fUWnvhqC13jQji3TwLn9E2gQ7TO0xFp7lZnF3P9GxmUVtXy/OUDGNsr3ulIIiIM6RBJr3ZhvLxoOxMHJbWq03GO6pOutbbYGPMNMAyIMMb4eGbzEoFcz2G5QBKQY4zxAcKpb8Dyn/H/OPg+Pzf+0+efBkyD+nPyjia7tB51LjdzV+Xy9PwtZBdV0i8pgj+cncqvesSpiYoDOkQHM3VkJ6aO7ERBWTWfrdvDx6v38M8vN/PPLzfTq10Y5/RJ4Oze8SRFamZVpLn5IDOXu2avISbUnzk3DqdHfOtsVy4iTY8xhmtP6sivZ2by7eYCTu3eejZHP5LGKzFArafACwS+AB4BrgDmHNR4ZY219lljzM1A74Mar4y31l5sjOkJvM1/G698DXQBDPWNV0ZTX9ytAC611q4/XC41XpGfcrktH2Tm8tTXW9i5r4I+ieH8ZkxXRnWN0ZKhJmhPSSWfrNnDR2v2sDq7GID+yRGM79+Oc/smEBGkpi0iTZnLbfnH5xt54dvtDOkQyXOXD9SWKiLS5NTUuRn5j2/oFBvMW9cOdTpOgzue7pp9qG+s4g14AbOstQ8ZYzpSv4VCJLAKuNxaW+3ZcuENoD9QBEw8qKnK74GrgTrg19bazzzjZwFPeJ7jFWvtX37pBanIk4Mt2baPhz/ewIY9paTGh3HHmK6M7hGr4q6ZyNpXwcdrd/PBqt1syi/Dz9uL0T1iuXBAIqd0i1HjBpEmprSqltvfWcU3mwq4fGgyD5zbU/9ORaTJem7BNh6Zt5HPbj+5xa02OOYir6lSkSdQXxz89dMfmLc+j3YRgdxzZnfO6ROv4q6Zstayfncpc1bm8EHmbooO1BAd4sd5fdtx4cB29EwIdzqiSKu3o/AA176+gl37KnjwvJ5cPrS905FERA6ruKKGIX/9mvEDEvnb+N5Ox2lQKvKkRTlQXcfT87fyyqId+HgbbhrViWtP7tjq2uO2ZLUuNws2FTAnI4evN+ZT67J0bxvKRQMTGdevHTGh/k5HFGl1vtm4l9tnrMLby/Dc5QMZ2jHK6UgiIkfkrtmr+Wj1HpbeN5rwwJazd6eKPGkx5m/M5w/vrWN3SRUXDkjkrrHdiAsLcDqWNKL9B2r4aM1u5mTksDqnBG8vw2ndY7l0cDIju8aooY5II3O7LU9+vYWn5m+he9swpk0eqEZJItKsrMst4ZynF3H/Oalcc1IHp+M0GBV50uwVlFXz0Mcb+Gj1brrGhfC38X0Y2L6N07HkBNuSX8bsjBzmrMyhsLyGhPAAJqQlcfGgJNpFBDodT6TFKa6o4TczM/lmUwEXDkjkz+f3ItBPqyZEpPm58LnF7CuvZv6do1rMdgoq8qTZstbybkYOf/nkByprXNxyWmduOKUTfj46yb81q6lz89UP+byzPItFWwsBGNU1homDkzmte6yaQIg0gHW5Jdz4VgZ5JVU8cG5PLhuSrHOeRaTZ+iAzl9tnZPLaVYMY1a1lbKegIk+apb1lVdw1ew0LNhUwKKUNfxvfh86xIU7HkiYmu6iCmSuymZWezd6yamJD/ZmQlsjEQclaUiZyjGZn5PD799YSGezHs5cNoH+yVk6ISPNWU+dm+N/n0ycxnFeuHOR0nAahIk+ana825HPXnDUcqK7jvrN6MHlo+xYztS6No87l5ptNBcxYnsU3m/bitjCyawyTh7bntO6xOndP5AhU17l46KMNvLUsi2Edo3j60v5Eh6jRkYi0DP/6cjNPz9/Ct789leSo5v9FsIo8aTYqaur48yc/8PayLFLjw3hyYj+6xIU6HUuamT0llcxckc07y7PIL62mXUQglw5J5pJBSfrAKvIzsosquOWdVazOLub6Uzryu9O74aOlzyLSguSVVDHikflcc1IH7jurh9NxjpuKPGkW1uWWcNs7q9ix7wBTT+7IHad3xd9HJ/jLsat1ufn6h3zeWLqL77fuw9fbcFbveCYPbc/A9m10fpGIx6dr93D3nDUAPHpRH8b2inc4kYhI47j5rZUs2lrI0ntHN/tGUj9X5Pk4EUbkUGatyOYPH6wjMsiPt64dwvBO0U5HkhbA19uLsb3iGdsrnq17y3lz6S7mZNRvtt69bSiTh7Xn/H7tCPbX26G0TlW1Lh7+uH55Zr+kCJ6e1F/nsopIizZlWHs+WbuHD1fncsmgZKfjNArN5InjqmpdPPjhemasyOakztE8ObEfUVpOJ42ooqaODzJ3M33JLn7YU0qovw8XDkzk8qHJdI7V0mBpPbbuLeeWt1eyMa+M60d25LdndFNnWhFp8ay1nPnkQryM4ZPbTmrWq3q0XFOapOyiCm58K4N1uaXccmpnfjOmq5pjyAljrWVl1n7eWLKLT9fmUeNyM6xjFJOHtWdMapw+7EqLNjsjh/vfX0egnzf/vLgvp7aQduIiIkfi7WVZ3PfeWubcOIyB7SOdjnPMVORJk/Pd5gJum7EKl9vy+MX9+FVqnNORpBUrLK9mVno2by3NIre4krgwfyYOSubSIcnEhQU4HU+kwZRX1/HH99cxd1UuQztG8uTE/vo7LiKtzoHqOob89WvO6NmWf17c1+k4x0xFnjQpbyzZyQMfrqdrXCjPXz6QlOhgpyOJAOByW77ZuJfpS3fx3eYCfLwMZ/Rsy+Rh7RnSIbJZL+kQWbGziDtnrSZnfwW3je7Crad10eoJEWm1fv/eWmZn5LD8vl8RHuTrdJxjosYr0iS43JY/f7KBV7/fyejusTw1qb8aXkiT4u1l+FVqHL9KjWNn4QHeXLqLWenZfLJ2D13jQpg8tD0XDEgkRH9vpRmprnPx+JdbeOG7bSS2CWTG1GEM7tB8lyeJiDSESYOTeWtZFu9n5nLF8BSn4zQozeTJCVNeXcft76zi6417uXpEB35/dg99gyzNQmWNi49W72b60p2syy0lxN+H8QPaMXloe+3hKE3exrxSfj0jk415ZUwclMQfzknVlxQiIh7nPbOI6lo38359crNcraOZPHHU7uJKrnk9nc35ZTw8rieTh6U4HUnkiAX6eXPxoCQmpCWyKruYN5fsYsbybKYv2cWwjlFMGdaeX6lRizQxLrflpYXb+ecXmwkL9OGlKWk691lE5CcmDU7m3rlrWZlVzMD2bZyO02A0kyeNbkt+GZNfXk55dR3PXNqfUergJi3AvvJqZv6kUculg9szaXASsWpiIQ7LLqrgzlmrWb6ziDN6xvHXC3praxoRkUM4UF3H4L98xZm943lsQvNrwKLGK+KI1dnFXPnqcry9vJh+9WBSE8KcjiTSoA7VqGVsr7ZMGZbCoJQ2zXLphzRfbrflzWW7eOSzjXgZw4Pn9WT8gHb6eygichj3vbeWuStzWHbfrwgPbF4NWLRcU064xVsLuW56OpEhfrx5zRDaR6mDprQ8Bzdq2eFp1PJuejYfr9lD97ahXD60PRf0b6cGQ9Lotu4t5965a1ixcz8nd4nmb+N7k9gmyOlYIiJN3qWDk3l7WRbvr2o5DVg0kyeNYt66PG57ZxUp0UG8cc0Q7cEkrUpljYsPMnOZvmQXG/aUEurvw4UDE7l8aHs6x4Y4HU9amFqXm2nfbefJr7YQ6OfN/eekcqFm70REjsq5Ty+i1uXms9ubVwMWzeTJCfNuejZ3z1lD36QIXr1yEBFBfk5HEjmhAv28mTg4mUsGJbEyq5g3luzk7WVZvLZ4JyM6RzF5aHt+1SMOHzVqkeO0Mms/v39vHT/sKeXs3vE8eF5PYkJ17p2IyNGaNDiZ+95by6rsYgYkN/8GLJrJkwY1c0UWd89Zy8ldonlh8kCC/PQ9gghAYXk1M1dk8/ay/zZquSQtiYsHJWlJnRy1/QdqeGTeRmasyKZtWAB/GteTM3q2dTqWiEizVV5dx5C/fMVZveN5tBk1YFHjFWl0s1Zkc/fcNYzsEsMLkwcS4OvtdCSRJsfltnz9Qz7vLM9iweYCAEZ2iWHS4GRG94jVNgxyWG63ZVZ6No/M20hpVR1Xj0jh9l911b53IiIN4N65a3lvVfNqwKLlmtKoZqXXF3gnq8ATOSxvL8PpPdtyes+25BZXMmtFNrPSs7nhzQxiQv2ZMDCRiYOSSY7S7J78r9XZxTz40XpWZRUzKKUND5/fi+5t1bFYRKShTBqcxDvLs/hw9W4mD23vdJzjopk8OW7vpmdz15w1nNQ5mhenpKnAEzlKdS43324u4J3l2czfmI/bwkmdo5k0OJkxqXH4+Wh2rzXbU1LJo/M2MXdVLtEhftxzZg81VhERaQTWWs58ciH+Pl58cMtJTsc5IprJk0YxOyNHBZ7IcfLx9mJ0jzhG94hjT0kl76bnMHNFNje/vZLIYD/O65vARQMT6ZkQpg/2rUhFTR0vfLudF77bhtvCjaM6cdOoToQGNI8lRCIizY0xhosGJvLnT35gc34ZXeNCnY50zDSTJ8fs07V7uOXtlQzvFM1LV6jAE2lILrdl4ZYC3k3P4csN+dS43HRvG8qFAxIZ1z+B2FBtS9JS1bnczFmZw+NfbiGvtIqz+8Rzz9juJEVqCa+ISGPbV17NkL9+zVUjUvj92alOx/lFarwiDWrRlkKufm0FvRPDefOaIQT6qcATaSwlFbV8tGY3szNyyMwuxtvLcErXGC4ckMjoHrH6gqWFcLstH6/dw+NfbmZH4QH6JUXwh7N7kJYS6XQ0EZFWZer0dFZm7WfJvaObfEM0LdeUBpOZXczUN9LpGBPMK1cMUoEn0sjCg3y5fGh7Lh/anq17y5m7Moe5K3OZv3ElYQE+nNcvgfEDEumfFKHlnM2QtZavf9jLY19sYmNeGd3bhvLSlDRG94jV71NExAET0pL4YkM+CzYVMCY1zuk4x0QzeXJUtu4t46LnlxAa4MOcG4YTG6YlYyJOcLkti7cVMjsjh8/X51FV6yYpMpBz+yRwbt8EurcNVYHQxLndls/X5/HvBVtZl1tKSlQQvxnTlXP7JODlpd+diIhTal1uhv3tawYkt2HalP83SdakaCZPjlvO/gouf2k5Pl5evHnNEBV4Ig7y9jKc3CWGk7vEUFZVy7x1eXy0Zg8vfLedZxdso3NsCOf2SeCcvvF0iglxOq4cpNbl5oPM3Ty3YCvbCg6QEhXEIxf2ZvyAxCa/LEhEpDXw9fbigv7tePX7nRSWVxMd4u90pKOmmTw5IvsP1HDh84spKKtm5tRhpCZobyaRpmhfeTWfrcvjo9W7Wb6zCGshNT6Ms/vEc0bPODrHNt9OYc1dSWUt76Zn8+r3O8ktrqRHfBg3jerEWb3j8dbMnYhIk7I5v4zTH/+O+89J5ZqTOjgd52ep8Yocs+o6F5NfXk5mVjFvXjuEwR3UBECkOcgrqeKTtXv4aPVuMrOLAegYE8zpqW05vWcc/RIjtCzwBNiSX8brS3Yyd2UuFTUuBqW04aZRnRnVLUZLakVEmrBxzyyius7NZ7ef3GTfr1XkyTGx1nLHrNW8tyqXJyf2Y1y/dk5HEpFjsKekkq825PP5+nyWbt9HndsSG+rPmNQ4Tu/ZlqEdI/H3UROlhlJT52b+xnzeXJrFoq2F+Pl4cV7fBK4cnkKvduFOxxMRkSPwxtJd3P/+Oj6+9aQm+959zEWeMSYJmA7EARaYZq190hjzIHAdUOA59D5r7aee+9wLXAO4gNustZ97xscCTwLewEvW2r97xjsAM4AoIAOYbK2tOVwuFXknxhNfbeaJr7Zw55iu3Dq6i9NxRKQBlFTU8s2mvXy+Po8FmwqorHUR6OvN0I6RnNI1hpFdY+gQHdxkv7Vsqqy1rN9dyuyMHD7IzGV/RS1twwKYPKw9EwclEdUMz+kQEWnNSipqGfTXr5g0KIk/jevldJxDOp4iLx6It9auNMaEUl+EnQ9cDJRbax/7yfGpwDvAYCAB+Aro6rl5MzAGyAFWAJOstRuMMbOAudbaGcaY54HV1trnDpdLRV7je29VDr+ZuZoLByTy2IQ++sAn0gJV1bpYvK2QbzcV8N2WQnYUHgAgKTKQkV3qC77hnaIIDfB1OGnTtbPwAPPW5/H+qlw25pXh5+3FmJ5xXDQwkZM7R+OjZioiIs3Wre+sYuGWApbdN7pJrng55u6a1to9wB7Pz2XGmB+Aw63ZGwfMsNZWAzuMMVupL/gAtlprt3sCzQDGeR7vNOBSzzGvAw8Chy3ypHEt31HE3bPXMrRjJH8b31sFnkgLFeDrzWnd4zite/0+QLv2HeC7zQV8u7mQ91bl8tayLLy9DKnxYQzuEMmglEgGpbRp1bNS1lo25Zcxb10e89blsTGvDIC+SRE8fH4vzuuTQHiQimIRkZZgwsBEPlq9m6827OXsPvFOxzliR7WFgjEmBegPLANGALcYY6YA6cCd1tr91BeASw+6Ww7/LQqzfzI+hPolmsXW2rpDHP/T558KTAVITk4+muhyFLL2VTD1jXQSIwN54fI0/Hz0LbRIa9E+KpjJw4KZPCyFmjo3Gbv2s3hbIct3FPHG0l28vGgHAJ1jQxiUEsmQDpEMSG5DUmRgi/4yqOhADYu3FfL91kIWbikkZ38lxsCg9pHcf04qp6fGkRQZ5HRMERFpYCM6RzNhYCJtw5vX1mFHXOQZY0KAOcCvrbWlxpjngIepP0/vYeCfwNWNktLDWjsNmAb1yzUb87laq4qaOqa+kY618OqVg/RttEgr5ufjxbBOUQzrFAXUd9pdm1PC8p1FLN9RxMerd/PO8iwAwgN96ZMYTu929ZeeCeEktglslt07rbXkFleSmV1MZlYxS7bvY/3uUgBC/X0Y2imKG0d1YkxqHLGhzet/+iIicnS8vQyPTujrdIyjdkRFnjHGl/oC7y1r7VwAa23+Qbe/CHzsuZoLJB1090TPGD8zvg+IMMb4eGbzDj5eTiBrLb+bvYbN+WW8dtVg2kcFOx1JRJoQfx9v0lIiSUuJ5KZR4HJbNuaVsjq7hLW5xazJKWHad9upc9d/Bxfo603XuBC6xIXSLS6UzrEhJEcFkdgmsMmc11Bd52JnYQVb95azZW8Za3NKWJ1TTGF5fe8vPx8v+idF8NvTuzK8czR92oXrHDsREWnyfrHIM/Xrb14GfrDW/uug8XjP+XoAFwDrPD9/CLxtjPkX9Y1XugDLAQN08XTSzAUmApdaa60x5hvgIuo7bF4BfNAQL06OzgvfbeeTNXu458zujOwa43QcEWnivL0MPRPqZ+2gfgl9Va2LjXllbNxTyqb8Mjbnl7FgUwGzM3J+vJ8xkBAeSHJkEEmRgcSFBRAbFkBcqD9xYQFEh/oTFuBDiL/PcS0BdbstZdV1lFTUkldaxZ6SSnYXV5FXUklucRXbC8rZVVSBy1OUGgOdYkI4pWss/ZIj6JcYQbe2oVqyLiIizc6RzOSNACYDa40xmZ6x+4BJxph+1C/X3AlcD2CtXe/plrkBqANutta6AIwxtwCfU7+FwivW2vWex7sbmGGM+TOwivqiUk6g7zYX8I95GzmnTzzXj+zodBwRaaYCfL3plxRBv6SI/xkvOlDDjsJydhZWsKuogqx9B9hVVMGCTQUUllfjPsQCfC8DYYG+hAX4EuLvg6+PF37eBl9vL3y8vfDxMtS63NS5LHVuN7UuS02dm7LqWkoqaimrruNQDaTDAnxIiAika1woZ/eJp3NsCJ1i6i+Bfk1jhlFEROR4aDN0Yde+A5z3zPfEhwcw96bhBPkdVT8eEZHj4nJb9pVXk19aTX5pFYXl1ZRV1VFaVUtpZS2lVXWUVdVR63IfdKkv7Hy8vPA9qPDz8zaEBfjWF4eBvoQF+BAe6Evb8ADiwwOJDw8g2F/vcSIi0jIc8xYK0rJV1NRx/RsZAEybnKYCT0ROOG8vQ6xnyWZvwp2OIyIi0uzpE30rd//769mUX8brVw0mOUrtv0VEREREmjudTd6KvZuezZyVOdx6Whc1WhERERERaSFU5LVSm/PLuP+DdQztGMnto7s4HUdERERERBqIirxWqKKmjpvfWkmIvw9PTeyPdzPcrFhERERERA5N5+S1Qn/8YD1bC8p54+ohxIYFOB1HREREREQakGbyWpnZGTnMzsjh1lM7c1KXaKfj/F97dx9tVV0mcPz7cLmCKAgKKsOLkECKL5mSYeZkSIYvRbNyJqeVL2VZkzZazSqrVTQ1q5kax5rKdDnFqE0rLV+SHFyE5Uu10kR8AxG8WAqEgKAio6Jwn/njbuouvMhBOGffs8/3s9ZZ95zf/p1zn/OcH/fsh9/evy1JkiRpF7PIayEdqzfwxZ8V5+FNnVB2OJIkSZLqwCKvRby0qZMLr72P/u19+E/Pw5MkSZIqy3PyWsS3bl3CghXrufwDR7Gf5+FJkiRJleVMXgv4/R/WcdkdS3nfpFFMO3T/ssORJEmSVEcWeRW3/sWX+eS19zN67wF86V0Tyw5HkiRJUp15uGbFfXnWQp5c/yI//dgx7NHPj1uSJEmqOmfyKuzmB//EDfNXcP7bx3Hk6CFlhyNJkiSpASzyKurJZ1/kCzcu4IhRgzl/yriyw5EkSZLUIBZ5FZSZXHTDg7y0qZNvvu8I2tv8mCVJkqRW4d5/BV0/fwW3L17DZ6a9nrFD9yg7HEmSJEkNZJFXMavWv8hXfr6QN40ZwlnHjCk7HEmSJEkNZpFXIZnJF25cwMZNnXzjtDfQp0+UHZIkSZKkBrPIq5BZD/yJWxet4p9O9DBNSZIkqVVZ5FXEmuc2MmPWQt44ejAfeuvYssORJEmSVBKLvIr40k0LeP6lzfz7aYfT5mGakiRJUsuyyKuA2Q+t5JYFT/LJqRMYt+/AssORJEmSVCKLvCa3/sWXmTFrIYeOGMRHjvMwTUmSJKnV9S07AO2ci+csZu2Gjcw860309aLnkiRJUsuzKmhi9y97hh/e9ThnHjOGw0buVXY4kiRJknoBi7wmtWlzJ5+74SH2HdiPT584oexwJEmSJPUSHq7ZpP77t39k0cr1XP6BIxnYv73scCRJkiT1Es7kNaEVz7zAJXOXcMJB+/LOQ/YvOxxJkiRJvYhFXpPJTGbctACAf55+CBFeE0+SJEnSX1jkNZk5C1dx66LVXDh1PCOHDCg7HEmSJEm9jEVeE3nhpc189eaHef1+A/nQW70mniRJkqRXcuGVJnL5HUtZ8cwLXHPuZNq9Jp4kSZKkHmy3UoiIURFxW0Q8HBELI+KCon3viJgbEY8WP4cU7RER346Ijoh4MCKO7PZaZxX9H42Is7q1HxURDxXP+XZ4otkrLFv3PJffsZRTDx/O5NftU3Y4kiRJknqpWqaDNgGfzsyJwGTgvIiYCFwE/DIzxwO/LB4DnASML27nApdBV1EIzADeDBwNzNhSGBZ9M3Qn5gAADKdJREFUPtLtedN2/q1Vy9dmL6JPBJ8/+eCyQ5EkSZLUi223yMvMlZk5v7j/HLAIGAFMB64qul0FvKe4Px24OrvcBQyOiOHAO4G5mbkuM58G5gLTim2DMvOuzEzg6m6vJeC3HU9xy4InOe/tB/JXg3cvOxxJkiRJvdgOndgVEWOANwJ3A/tl5spi05PAfsX9EcCybk9bXrS9WvvyHtoFvLy5ky/PWsjovQfw4eNeV3Y4kiRJknq5mou8iNgTuB64MDPXd99WzMDlLo6tpxjOjYh5ETFvzZo19f51vcLVv3ucR1dv4IunTqR/e1vZ4UiSJEnq5Woq8iKina4C70eZeUPRvKo41JLi5+qifQUwqtvTRxZtr9Y+sof2V8jMKzJzUmZOGjZsWC2hN7WnNmzkW3OX8NcThjH14H3LDkeSJElSE6hldc0AfgAsysxLum2aBWxZIfMs4KZu7WcWq2xOBp4tDuucA5wYEUOKBVdOBOYU29ZHxOTid53Z7bVa2sVzFvPCy5uZ8a6JuOCoJEmSpFrUcp28Y4EzgIci4v6i7fPAvwE/iYhzgMeBvyu2zQZOBjqA54EPAmTmuoj4KnBP0e8rmbmuuP9x4Epgd+CW4tbSFq1cz7XzlvGhY8dy4LA9yw5HkiRJUpPYbpGXmb8BtjWNdEIP/RM4bxuvNROY2UP7PODQ7cXSSv71lkcY1L+dT0wZV3YokiRJkprIDq2uqca4Y8ka7lyyhk9MGcfgAbuVHY4kSZKkJmKR18ts7ky+9r+LGL33AM445oCyw5EkSZLUZCzyepnr7l3G4lXPcdFJB9Gvr5dMkCRJkrRjLPJ6kf/buImLf7GEI0cP5qRD9y87HEmSJElNyCKvF7nizsdY89xGvnCKl0yQJEmS9NpY5PUSq9a/yBV3PsYphw3nqAOGlB2OJEmSpCZlkddLXPKLJWzq7OSz0w4qOxRJkiRJTcwirxfoWL2Bn967jDMmj2H0PgPKDkeSJElSE7PI6wUumbuY3dvbOO/tB5YdiiRJkqQmZ5FXsoeWP8vsh57knONexz579is7HEmSJElNziKvZN+Y8whDBrTzkePGlh2KJEmSpAqwyCvR75au5dePPsXHjx/HwP7tZYcjSZIkqQIs8kqSmXxjziPsP6g/ZxxzQNnhSJIkSaoIi7yS3LpoNfc98QwXTB1P//a2ssORJEmSVBEWeSXY3JlcPGcxY4fuwd8eNbLscCRJkiRViEVeCWY9sILFq57jU++YQN82PwJJkiRJu44VRoO9vLmTb859lInDB3HKYcPLDkeSJElSxVjkNdiN963giXXP86l3TKBPnyg7HEmSJEkVY5HXQC9v7uS7v+rg0BGDOOHgfcsOR5IkSVIFWeQ10JZZvAtPmECEs3iSJEmSdj2LvAbZMot32Ii9nMWTJEmSVDcWeQ3y51m8qeOdxZMkSZJUNxZ5DdB9Fm/KQc7iSZIkSaofi7wGuHG+s3iSJEmSGsMir85e3tzJd2571Fk8SZIkSQ1hkVdnN85fwbJ1LziLJ0mSJKkhLPLqaNPmTr57m+fiSZIkSWoci7w6uvnBlTyx7nk+MWWcs3iSJEmSGsIir046O5Pv3d7BhP32ZOrB+5UdjiRJkqQWYZFXJ7cuWsWSVRv4+PHj6NPHWTxJkiRJjWGRVweZyaW3L2XU3rtz6uHDyw5HkiRJUguxyKuD3y1dywPLnuFjbzuQvm2mWJIkSVLjWIHUwaW3d7DvwH6898iRZYciSZIkqcVY5O1i9z3xNL/tWMuHjxtL//a2ssORJEmS1GK2W+RFxMyIWB0RC7q1fTkiVkTE/cXt5G7bPhcRHRGxOCLe2a19WtHWEREXdWsfGxF3F+3XRsRuu/INNtr3bl/KXru38/43H1B2KJIkSZJaUC0zeVcC03po/2ZmHlHcZgNExETgdOCQ4jnfi4i2iGgDLgVOAiYCf1/0Bfh68VrjgKeBc3bmDZVpyarnmPvwKs5+yxj27Ne37HAkSZIktaDtFnmZeSewrsbXmw5ck5kbM/MPQAdwdHHryMzHMvMl4BpgenRdIXwKcF3x/KuA9+zge+g1Lrt9KQN2a+Pst4wpOxRJkiRJLWpnzsk7PyIeLA7nHFK0jQCWdeuzvGjbVvs+wDOZuWmr9h5FxLkRMS8i5q1Zs2YnQt/1ntqwkZ8/8Cfef/RohuzR1EecSpIkSWpir7XIuww4EDgCWAn8xy6L6FVk5hWZOSkzJw0bNqwRv7JmQ/fsx+wLjuOjbzuw7FAkSZIktbDXdOJYZq7acj8i/gu4uXi4AhjVrevIoo1ttK8FBkdE32I2r3v/pjNhv4FlhyBJkiSpxb2mmbyIGN7t4d8AW1benAWcHhH9ImIsMB74PXAPML5YSXM3uhZnmZWZCdwGnFY8/yzgptcSkyRJkiSphpm8iPgxcDwwNCKWAzOA4yPiCCCBPwIfBcjMhRHxE+BhYBNwXmZuLl7nfGAO0AbMzMyFxa/4LHBNRPwLcB/wg1327iRJkiSpxUTXZFrzmTRpUs6bN6/sMCRJkiSpFBFxb2ZO2rp9Z1bXlCRJkiT1MhZ5kiRJklQhFnmSJEmSVCEWeZIkSZJUIRZ5kiRJklQhFnmSJEmSVCEWeZIkSZJUIU17nbyIWAM8XnYcPRgKPFV2EC3K3JfH3JfH3JfH3JfH3JfH3JfL/Jent+b+gMwctnVj0xZ5vVVEzOvpgoSqP3NfHnNfHnNfHnNfHnNfHnNfLvNfnmbLvYdrSpIkSVKFWORJkiRJUoVY5O16V5QdQAsz9+Ux9+Ux9+Ux9+Ux9+Ux9+Uy/+Vpqtx7Tp4kSZIkVYgzeZIkSZJUIRZ5r0FETIuIxRHREREX9bC9X0RcW2y/OyLGND7K6qoh/2dHxJqIuL+4fbiMOKsmImZGxOqIWLCN7RER3y4+lwcj4shGx1hVNeT++Ih4ttuY/1KjY6yqiBgVEbdFxMMRsTAiLuihj2O/DmrMvWO/DiKif0T8PiIeKHL/zz30cV+nDmrMvfs5dRQRbRFxX0Tc3MO2phn3fcsOoNlERBtwKfAOYDlwT0TMysyHu3U7B3g6M8dFxOnA14H3NT7a6qkx/wDXZub5DQ+w2q4EvgtcvY3tJwHji9ubgcuKn9p5V/LquQf4dWae2phwWsom4NOZOT8iBgL3RsTcrf7mOPbro5bcg2O/HjYCUzJzQ0S0A7+JiFsy865ufdzXqY9acg/u59TTBcAiYFAP25pm3DuTt+OOBjoy87HMfAm4Bpi+VZ/pwFXF/euAEyIiGhhjldWSf9VBZt4JrHuVLtOBq7PLXcDgiBjemOiqrYbcq04yc2Vmzi/uP0fXF/+Irbo59uugxtyrDoqxvKF42F7ctl7EwX2dOqgx96qTiBgJnAJ8fxtdmmbcW+TtuBHAsm6Pl/PKL50/98nMTcCzwD4Nia76ask/wHuLw6aui4hRjQmt5dX62ag+jikO77klIg4pO5gqKg7LeSNw91abHPt19iq5B8d+XRSHrN0PrAbmZuY2x737OrtWDbkH93Pq5VvAZ4DObWxvmnFvkacq+jkwJjMPB+byl/9xkapqPnBAZr4B+A7ws5LjqZyI2BO4HrgwM9eXHU8r2U7uHft1kpmbM/MIYCRwdEQcWnZMraKG3LufUwcRcSqwOjPvLTuWXcEib8etALr/j8nIoq3HPhHRF9gLWNuQ6Kpvu/nPzLWZubF4+H3gqAbF1upq+behOsjM9VsO78nM2UB7RAwtOazKKM6LuR74UWbe0EMXx36dbC/3jv36y8xngNuAaVttcl+nzraVe/dz6uZY4N0R8Ue6TgeaEhH/s1Wfphn3Fnk77h5gfESMjYjdgNOBWVv1mQWcVdw/DfhVekHCXWW7+d/qXJh303Ueh+pvFnBmsdLgZODZzFxZdlCtICL233JOQEQcTdff9l75pdNsirz+AFiUmZdso5tjvw5qyb1jvz4iYlhEDC7u707XYmePbNXNfZ06qCX37ufUR2Z+LjNHZuYYuvYvf5WZH9iqW9OMe1fX3EGZuSkizgfmAG3AzMxcGBFfAeZl5iy6vpR+GBEddC2WcHp5EVdLjfn/x4h4N10rs60Dzi4t4AqJiB8DxwNDI2I5MIOuE8LJzMuB2cDJQAfwPPDBciKtnhpyfxrwDxGxCXgBOL23fuk0oWOBM4CHinNkAD4PjAbHfp3VknvHfn0MB64qVrTuA/wkM292X6chasm9+zkN1KzjPvxbKEmSJEnV4eGakiRJklQhFnmSJEmSVCEWeZIkSZJUIRZ5kiRJklQhFnmSJEmSVCEWeZIkSZJUIRZ5kiRJklQhFnmSJEmSVCH/D5Br2+c203+tAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 1080x360 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAEvCAYAAAD4uAgWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3iUVfrG8e9J7wmpJCQhdAgdQldEWRQriqKggh277uqubdfV1W2u7lrXgh0bIGBXbIiC1ARClV5SICEhpJE+c35/ZNZl/SFSEt6U+3Ndc5k5887MPQTGeeac9znGWouIiIiIiIi0DF5OBxAREREREZGGoyJPRERERESkBVGRJyIiIiIi0oKoyBMREREREWlBVOSJiIiIiIi0ICryREREREREWhCfIz3QGOMNpAO51tpzjDGvAacAJZ5DrrTWZhpjDPAkcBZQ4Rlf6XmMK4A/eI7/s7X2dc/4QOA1IBD4FLjd/sLeDtHR0TYlJeVI44uIiIiIiLQoGRkZhdbamJ+OH3GRB9wO/ACEHTT2O2vt7J8cdybQxXMZAjwHDDHGRAIPAGmABTKMMR9aa/d7jrkOWEZ9kTcW+OxwYVJSUkhPTz+K+CIiIiIiIi2HMWbXocaPaLmmMSYROBt46QgOHwdMt/WWAhHGmHjgDOBLa22Rp7D7EhjruS3MWrvUM3s3HTj/SHKJiIiIiIjI/zrSc/KeAO4C3D8Z/4sxZo0x5nFjjL9nrB2QfdAxOZ6xw43nHGJcREREREREjtIvFnnGmHOAvdbajJ/cdC/QHRgERAJ3N3y8/5dlqjEm3RiTXlBQ0NhPJyIiIiIi0uwcyUzeCOA8Y8xOYAZwmjHmTWvtHs+SzGrgVWCw5/hcIOmg+yd6xg43nniI8f/HWjvNWptmrU2Lifl/5xeKiIiIiIi0er9Y5Flr77XWJlprU4CJwHxr7eWec+nwdNM8H1jnucuHwBRTbyhQYq3dA3wOnG6MaWOMaQOcDnzuua3UGDPU81hTgA8a+HWKiIiIiIi0CkfTXfOn3jLGxAAGyARu8Ix/Sv32CVup30LhKgBrbZEx5mFghee4h6y1RZ6fb+K/Wyh8xi901hQREREREZFDM7+wHV2TlZaWZrWFgoiIiIiItFbGmAxrbdpPx4+0u6aIiIiIiIg0AyryREREREREWpDjOSdPRERERKRZqap1UXSg5sdLeXUdLrfFbS0ud/3FGEOIvw9hAT6EBvgSGuBDRJAv4YG+1PcJFGnaVOSJiIiISIuSX1rFprwydu07wI7CCnbuO8DOfQfYU1xFZa3rmB832M+bpMggkiKDSI4Mon1UEKnxYfRMCCfQz7sBX4HI8VGRJyIiIiLNVnFFDauyi1mbU8KanGLW5JSwt6z6x9sDfb1pHxVEt7hQTusWS2SIH5FBfrQJ9iMy2I/QAB98vAzGGLyNwdvL4LaWsqo6z6WW8uo6ig7UkLO/kpz9Fezad4BFWwp/LBi9vQxdYkPokxhO36QIhneKpkN0sFN/JCIq8kRERESk+aiqdbFiZxGLthayeOs+1u0uwVowBjpGB3NS52h6J4bTIz6MDtHBxIb6N8oSS2st+aXVrM39b3H55YZ8ZqXnANA+KohTusZwStcYhnWKIshPH7vlxNEWCiIiIiLSpO0treLzDfl8sT6PZTuKqKlz4+tt6J/chhGdohncIZJe7cIIDfB1NKe1ll37Kli4pYAFmwpYvG0flbUu/Ly9GNUthvED2nFq91j8fbS0UxrGz22hoCJPRERERJqc7KIK5q3LY976PFZm7cda6BgTzOjusYzoXF/YNfXZseo6Fyt27Ofrjfl8tHoPheXVhAf6cnafeMb3b8fA9m3UyEWOi4o8EREREWnSDlTX8enaPczOyGHZjiIAeiaEMbZnW8b2akuXuFCHEx67OpebRVsLeW9VLp+vz6Oq1k3XuBCuPakj5/VLIMBXs3ty9FTkiYiIiEiTY61l6fYi3s3I5rO1eVTWuugQHcxFAxM5r28CSZFBTkdscOXVdXy6Zg+vfL+DjXllRIf4MWVYCpcNSSYqxN/peNKMqMgTERERkSajoqaO91bl8vrinWzOLyfU34dz+sZz0cBEBiS3jmWM1loWb9vHSwu3882mAvx9vJg4KIlbTutCTKiKPfllKvJERERExHHZRRVMX7KTmSuyKa2qo2dCGFcMT+HcPgmteq+5rXvLePG7HcxemYO/jxfXnNSB60Z2JMzhZjLStKnIExERERHHbNhdyrMLtvLp2j0YYzizV1uuHJ6i5iM/sb2gnH99uZmP1+whIsiXm0Z1YsqwFJ2zJ4ekIk9ERERETrj0nUU8u2Ab8zfuJcTfh8uGJnPV8A60DQ9wOlqTti63hH98vonvNheQFBnIQ+N6cWq3WKdjSROjIk9ERERETpil2/fxry83s3xHEZHBflw9IoXJw1IID9Tyw6Px/dZC/vjBOrYVHOCs3m154NyexIWpQJZ6KvJEREREpNGtzi7msS82sXBLIXFh/lw/shMTByc1+T3tmrLqOhcvfredp+dvxdfbi9+e3pXJw1Lw9tIy19ZORZ6IiIiINJrN+WX884tNfL4+nzZBvtx8amcuH9pe55I1oF37DvCH99excEshfZMieOKSfnSIDnY6ljhIRZ6IiIiINLi9pVU89sUm3s3IIdjPh+tO7sjVJ6UQqq6QjcJay0dr9nD/++uodbl54NxULk5LUvOaVurnijzNm4uIiIjIUauqdfHSwu08u2AbtS43V4/owC2ndqZNsJ/T0Vo0Ywzn9U1gUEob7py1mrvnrOWbjQX8bXxv/dnLjzSTJyIiIiJHzFrLh6t38495m8gtruSMnnHcc2YPLRt0gNtteWnRdh79fBORwX786+J+jOgc7XQsOYF+bibPy4kwIiIiItL8bMor45IXlnL7jEwignx557qhvDA5TQWeQ7y8DFNHduL9m0cQGuDL5JeX8dyCbTTXSRxpOFquKSIiIiKHVV5dx5NfbeaV73cSFuDD38f3ZkJakro7NhE9E8L58JYR3DV7DY/M28i63SU8elEfdTRtxfSbFxEREZFDstby6do8Hvp4Pfml1UwanMRdZ3TXuV9NUJCfD09P6k/vduE8Mm8j2/aWM21yGslRQU5HEwdouaaIiIiI/D+5xZVc9doKbn57JVHB/sy9aTh/G99HBV4TZozh+lM68dpVg9lTUsW5zyxi4ZYCp2OJA1TkiYiIiMiP3G7L64t3cvq/vmXZ9iLuPyeVD28ZwYDkNk5HkyM0smsMH94ygvjwAK58dQWzVmQ7HUlOMC3XFBEREREAtu4t4+45a8nYtZ+Tu0Tz1wt6kxSp5X7NUfuoYObcOJwb3szgrjlryC+t4pbTOms/vVZCRZ6IiIhIK+dyW174bhtPfLmFQD9vHpvQlwsHtFNB0MwF+/vw8hWDuHvOGv755Wbyy6r403m91DCnFVCRJyIiItKK7dp3gDtmrSZj137O7NWWh8b1IibU3+lY0kD8fLz454S+xIb588K32yksq+GJif0I8PV2Opo0IhV5IiIiIq2QtZZ3lmfz50824O1lePySvpzfT7N3LZGXl+HeM3sQGxrAwx9vYMrLy3n5yjRCA3ydjiaNREWeiIiISCuzt7SKu+es4ZtNBYzoHMWjF/UlISLQ6VjSyK45qQMxof7cMTOTK15ZzutXD1ah10KpyBMRERFpRT5bu4f73ltLRY2LB85N5YphKXjpHK1W47y+Cfh5G255exVTPIVemAq9FkdbKIiIiIi0AmVVtdwxM5Mb31pJUmQQn9x2MleN6KACrxUa2yueZy4dwNqcEqa8vJzSqlqnI0kDU5EnIiIi0sKtySnm7KcW8cHq3dw+ugtzbhxO59gQp2OJg8b2asu/LxvAutwSJr+8nJJKFXotiYo8ERERkRbKWsvLi3Zw4XOLqXO5mTl1KL8Z0xVfb30EFDijZ1uevWwAG3aXMOXlZZRpRq/FOOJ/4cYYb2PMKmPMx57rHYwxy4wxW40xM40xfp5xf8/1rZ7bUw56jHs945uMMWccND7WM7bVGHNPw708ERERkdZp/4EarpuezsMfb2BUt1g+vf1k0lIinY4lTczpPdvy7GUDWb+7lKnTM6iqdTkdSRrA0XyNczvww0HXHwEet9Z2BvYD13jGrwH2e8Yf9xyHMSYVmAj0BMYCz3oKR2/g38CZQCowyXOsiIiIiByDFTuLOOuphXy3uZAHz01l2uSBRAT5OR1LmqgxqXH846I+LNm+j9/MzMTltk5HkuN0REWeMSYROBt4yXPdAKcBsz2HvA6c7/l5nOc6nttHe44fB8yw1lZba3cAW4HBnstWa+12a20NMMNzrIiIiIgcBZfb8sz8LUycthR/Hy/m3jScK0d00N538ovGD0jkD2f34LN1edz/wTqsVaHXnB3pFgpPAHcBoZ7rUUCxtbbOcz0HaOf5uR2QDWCtrTPGlHiObwcsPegxD75P9k/GhxwqhDFmKjAVIDk5+Qiji4iIiLR8+8qr+fXMTBZuKWRcvwT+ckFvQvy1W5YcuWtP7khheQ3Pf7uN6BB/7hjT1elIcox+8V++MeYcYK+1NsMYM6rxI/08a+00YBpAWlqavl4QERERATKzi7npzQwKD9TwyIW9uTgtSbN3ckzuHtuNfeXVPPX1FqJD/JgyLMXpSHIMjuTrnRHAecaYs4AAIAx4Eogwxvh4ZvMSgVzP8blAEpBjjPEBwoF9B43/x8H3+blxEREREfkZ1lreWZ7Ngx+uJzbMn7k3DqdXu3CnY0kzZozhb+N7s7+ilgc+XE9cWABn9GzrdCw5Sr94Tp619l5rbaK1NoX6xinzrbWXAd8AF3kOuwL4wPPzh57reG6fb+sX9X4ITPR03+wAdAGWAyuALp5unX6e5/iwQV6diIiISAtVVevirtlruO+9tQzrFMVHt5ykAk8ahI+3F89c2p8+7cL5zcxM1u8ucTqSHKXj2STlbuAOY8xW6s+5e9kz/jIQ5Rm/A7gHwFq7HpgFbADmATdba12emcBbgM+p7945y3OsiIiIiBxCdlEFFz63mHczcrhtdBdeuXIQbYLVPVMaToCvNy9OSSMswJfrXk9nb1mV05HkKJjm2jknLS3NpqenOx1DRERE5IT6ZtNefj0jE2stj1/Sj9E94pyOJC3YutwSLnp+MT3iw3jnuqEE+Ho7HUkOYozJsNam/XT8eGbyREREROQEsbZ+e4SrX1tBQkQgH916kgo8aXS92oXz+MX9WJVVzD1z1mhrhWZCRZ6IiIhIE1dZ4+KWd1bx2BebOa9vAnNvHE77qGCnY0krcWbveO4c05X3M3fz7IJtTseRI6DNU0RERESasNziSqZOT2fDnlLuObM714/sqO0R5IS75bTObNlbzqOfb6JbXCi/StUsclOmmTwRERGRJmrFziLGPbOIrH0VvHLFIG44pZMKPHGEMYZ/XNSHXu3CuGNWJtlFFU5HksNQkSciIiLSBM1YnsWlLy4lNMCX924ewandY52OJK1cgK83z146EIAb38qgqtblcCL5OSryRERERJqQOpebBz9czz1z1zKsUzTv3zSCzrEhTscSASA5Koh/XtyPdbmlPPTxBqfjyM9QkSciIiLSRJRU1HLFq8t5bfFOrju5A69ckUZ4kK/TsUT+x5jUOG44pRNvL8ti7socp+PIIajxioiIiEgTsGvfAa56bQU5RZU8NqEvFw1MdDqSyM/67eldWZm1n9+/t46eCeF0axvqdCQ5iGbyRERERBy2YmcR5//7e/YfqOHNa4eowJMmz8fbi2cm9SfY34cb38qgvLrO6UhyEBV5IiIiIg76IDOXy15cRpsgP967aQSDO0Q6HUnkiMSGBfD0pP7sLDzAHz9Y53QcOYiKPBEREREHWGt54qvN3D4jkwHtI5h703BSorXBuTQvwzpFccupnZm7MpeP1+x2Oo54qMgTEREROcGq61z8ZmYmT3y1hQsHJDL96iFEBPk5HUvkmNw6ugt9kyK4b+5adhdXOh1HUJEnIiIickIVHajh8peW8X7mbn53Rjcem9AHPx99JJPmy9fbiycv6Ued23LnrNW43dbpSK2e3lFERERETpBtBeVc8Oz3rM4p4elJ/bn51M4YY5yOJXLcUqKDefDcnizZvo8XF253Ok6rpyJPRERE5ARYvK2QC/79PeVVdbxz3VDO7ZvgdCSRBjUhLZGxPdvy2BebWJdb4nScVk1FnoiIiEgjm7syhykvLyc2LID3bx7BwPZtnI4k0uCMMfxtfG8ig/24fcYqKmtcTkdqtVTkiYiIiDQSay3//mYrd8xazaCUSObcOJykyCCnY4k0mjbBfjw2oS/bCg7wyLyNTsdptVTkiYiIiDSCOpebP7y/jkc/38S4fgm8fvVgwgN9nY4l0uhO7hLDFcPa8/qSnaTvLHI6TqukIk9ERESkgVXWuLjhzQzeWpbFDad04vGL+6mDprQqd43tTkJ4IHfNXkNVrZZtnmh6txERERFpQPvKq5n04lK+3riXh8b15J4zu+PlpQ6a0roE+/vwyIV92F54gCe+2uJ0nFZHRZ6IiIhIA9m17wAXPreYH/aU8vzlA5kyLMXpSCKOOalLNBMHJTHtu22szi52Ok6roiJPREREpAFkZhcz/tnFFFfW8vZ1QzijZ1unI4k47r6zexAbGsDdc9ZQU+d2Ok6roSJPRERE5Dh9/UM+k6YtJcjfmzk3Dmdg+0inI4k0CWEBvvzlgl5szCvj2QVbnY7TaqjIExERETkOby/L4rrp6XSODWHOjcPpFBPidCSRJmV0jzgu6N+OZ+Zv5Yc9pU7HaRVU5ImIiIgcA2st//pyM/e9t5aRXWOYMXUosaEBTscSaZL+eE4qEUG+3Dt3LW63dTpOi6ciT0REROQoudyWP7y/jqe+3sJFAxN5cUoawf4+TscSabLaBPvx+7N7kJldzIwV2U7HafFU5ImIiIgcheo6F7e+s/LHPfAevagPvt76SCXyS87v144hHSJ5ZN5G9pVXOx2nRdM7koiIiMgRKquq5cpXVvDp2jz+cHYP7jmzO8ZoDzyRI2GM4c/n9+JAdR1//2yj03FaNBV5IiIiIkegoKyaidOWsmJnEf+6uC/XntzR6UgizU6XuFCuG9mRdzNyWLGzyOk4LZaKPBEREZFfkF1UwYTnF7OtoJwXr0hj/IBEpyOJNFu3ntaZdhGB/OG9ddS6tHdeY1CRJyIiInIYG3aXMv65xeyvqOWta4dyardYpyOJNGtBfj48cG4qm/LLeO37nU7HaZFU5ImIiIj8jGXb93HJC0vwNoZ3bxjGwPZtnI4k0iKMSY1jdPdYHv9qM7uLK52O0+KoyBMRERE5hC/W5zH5leXEhvkz56bhdI0LdTqSSIthjOHB83riclv+8skPTsdpcX6xyDPGBBhjlhtjVhtj1htj/uQZf80Ys8MYk+m59POMG2PMU8aYrcaYNcaYAQc91hXGmC2eyxUHjQ80xqz13OcpozZVIiIi4qCZK7K44c0MUuPDePeG4bSLCHQ6kkiLkxQZxI2jOvHJ2j0s36EmLA3pSGbyqoHTrLV9gX7AWGPMUM9tv7PW9vNcMj1jZwJdPJepwHMAxphI4AFgCDAYeMAY8581D88B1x10v7HH/cpEREREjpK1lmcXbOXuOWs5qUsMb107hMhgP6djibRY14/sRHx4AA99vB632zodp8X4xSLP1iv3XPX1XA73GxgHTPfcbykQYYyJB84AvrTWFllr9wNfUl8wxgNh1tql1loLTAfOP47XJCIiInLU3G7Lwx//wD/mbeK8vgm8NCWNYH8fp2OJtGiBft7cPbY763JLmbMyx+k4LcYRnZNnjPE2xmQCe6kv1JZ5bvqLZ0nm48YYf89YOyD7oLvneMYON55ziHERERGRE6Kmzs0dszJ55fsdXDUihScu6Yefj1oXiJwI5/VNoF9SBP/4fBMHquucjtMiHNG7l7XWZa3tByQCg40xvYB7ge7AICASuLvRUnoYY6YaY9KNMekFBQWN/XQiIiLSClTU1HHd9HTez9zN787oxh/PScXLS+0BRE4ULy/DH89NpaCsmucWbHM6TotwVF9RWWuLgW+AsdbaPZ4lmdXAq9SfZweQCyQddLdEz9jhxhMPMX6o559mrU2z1qbFxMQcTXQRERGR/2f/gRoufXEZC7cU8Pfxvbn51M6o/5vIiTcguQ3j+iUwbeF2cvZXOB2n2TuS7poxxpgIz8+BwBhgo+dcOjydMM8H1nnu8iEwxdNlcyhQYq3dA3wOnG6MaeNpuHI68LnntlJjzFDPY00BPmjYlykiIiLyv3YXVzLhhSVs2FPKc5cPZOLgZKcjibRqd4/tjpeBv3+20ekozd6RnE0cD7xujPGmviicZa392Bgz3xgTAxggE7jBc/ynwFnAVqACuArAWltkjHkYWOE57iFr7X96pd4EvAYEAp95LiIiIiKNYkt+GVNeWU55VR3Trx7M0I5RTkcSafUSIgKZOrITT329hSuHF5GWEul0pGbL1De0bH7S0tJsenq60zFERESkmVmZtZ+rX1uBr7cXr181mNSEMKcjiYhHRU0dpz32LXHhAbx/03Atn/4FxpgMa23aT8fVNkpERERajW827eWyF5cRHujLnBuGq8ATaWKC/Hy44/SurM4uZt66PKfjNFsq8kRERKRVeH9VLte9nk7HmGBm3zCc5KggpyOJyCGM79+OzrEhPPrFJupcbqfjNEsq8kRERKTFe3nRDn49M5NBKZHMmDqUmFD/X76TiDjCx9uL353Rje0FB5idoQ3Sj4WKPBEREWmxrLU8Mm8jD3+8gTN7teXVqwYRGuDrdCwR+QWnp8bRPzmCJ77aQlWty+k4zY6KPBEREWmR6lxu7p6zhucWbOPSIck8c+kAAny9nY4lIkfAGMPdY7uTV1rF64t3Oh2n2VGRJyIiIi1OVa2LG99ayaz0HG4b3YW/nN8Lby916RNpToZ2jGJUtxieXbCNkspap+M0KyryREREpEUpqahl8svL+OqHfP50Xk/uGNNVbdhFmqnfndGNkspanv92m9NRmhUVeSIiItJi5JdWcfELS8jMLubpSf25YniK05FE5Dj0TAhnXL8EXv1+B/mlVU7HaTZU5ImIiEiLsL2gnPHPLiZnfwWvXTWYc/okOB1JRBrAnWO6UeeyPPn1FqejNBsq8kRERKTZW51dzEXPL6Gq1sWMqcMY0Tna6Ugi0kCSo4KYNDiZWSuyyS6qcDpOs6AiT0RERJq17zYXMOnFpQT7ezP7xuH0Tgx3OpKINLCbT+2Ml5fh399sdTpKs6AiT0RERJqtDzJzufq1FbSPCmbODcPpEB3sdCQRaQRtwwO4dHAyszNyNJt3BFTkiYiISLP0yqId3D4jk4Ht2zDz+qHEhgU4HUlEGtENp3TCy8vwzHzN5v0SFXkiIiLSrFhr+ce8jTz08QbG9mzL61cPJizA1+lYItLI/jObN2dlDln7NJt3OCryREREpNmoc7m5e84anl2wjUmDk/n3ZQMI8PV2OpaInCA3jvLM5n2jTpuHoyJPREREmoWqWhc3vLmSWek53Da6C3+9oBfeXtrkXKQ1iQv7z2xermbzDkNFnoiIiDR5JRW1TH55GV9vzOfhcT25Y0xXjFGBJ9Ia3TSqEz6azTssFXkiIiLSpOWVVHHxC0tYnV3CM5MGMHlYitORRMRBsWEBXDqkfjZv174DTsdpklTkiYiISJO1dW85Fz63mNziSl67ahBn94l3OpKINAE3nuKZzVOnzUNSkSciIiJN0qqs/Ux4fjHVdS5mTB3K8M7RTkcSkSbiP7N5c1flat+8Q1CRJyIiIk3OlxvymfTiUsICfZl9w3B6tQt3OpKINDFTR3bEy8CLC7c7HaXJUZEnIiIiTcobS3dx/RvpdIsLZc6Nw0mJDnY6kog0QfHhgYzvn8jMFdkUlFU7HadJUZEnIiIiTcJ/Njm///11nNotlnemDiU6xN/pWCLShE09pSM1Ljevfr/D6ShNioo8ERERcVxNnZs7Z63+cZPzFyYPJMjPx+lYItLEdYoJ4cxebXljyS5Kq2qdjtNkqMgTERERR5VV1XL1ayuYuyqXO8d05a8X9MLHWx9RROTI3DSqM2XVdby1NMvpKE2G3kFFRETEMfmlVUx4fglLt+/j0Yv6cOvoLtrkXESOSq924ZzcJZqXF+2gqtbldJwmQUWeiIiIOGJzfhkX/Pt7sosqeOXKQUxIS3I6kog0UzeN6kxheTXvZuQ4HaVJ0GJ3kSbAWkuNy02ty1Jb56bW5abObfHxMnh7GXy8vfD1Nvh6e+GrJUwi0gIs3b6PqdPT8ff1Zub1w7RFgogcl6EdI+mfHMG077YxaVBSq1/yrSJPpJG43ZaC8mpyiyvZ/eOlir1lVRRX1FJSWfvjf8ur6474cf19vAgN8CUswIfQAB/CAn2JDvEnNtSfmFB/YsMCiAv1JzEyiLZhAXh7admTiDQtH67ezW9nrSYpMpDXrhpMUmSQ05FEpJkzxnDjKZ2Y+kYGH6/Zw/n92zkdyVEq8kSOk7WW3SVV/LC7lC17y9mSX8aWveVs3VtO5U/WhYf6+xAb5k+bID/ahgXQrW0o4YG+hAb44u9TP1vn5+2Fr48X3sbgspY6l/1xZq+2zk15dR2lVXWUVdVSVlVHSWUt2wsOUFBWTY3L/T/P5+ttSGwTRFJkEO0jg+gYE0zXuFC6xIUQE+Kv815E5ISy1vLM/K3888vNDE6JZNqUgUQE+TkdS0RaiF/1iKNLbAjPLdjGuH4Jrfpzjoo8kaNUUVPHmpwSVmUVk5m9n1VZxew9aAPOtmEBdIkLYdLgZDrEBJMYEUhCRCDxEQGEBfg2Wi5rLaWVdewtqyKvtIqc/ZVkFVWQta+CrKIKVmXtp6zqvzOGbYJ86RIXSmp8GL3bhdM7MZxOMSGa+RORRlFd5+LeOWuZuyqXC/q34+8X9sbfx9vpWCLSgnh5GW44pRN3vruabzcXMKpbrNORHGOstU5nOCZpaWk2PT3d6RjSCrjclvW7S1i4pZCFWwrI2LWfWlf9v5uUqCD6J7ehX1IEvdqF0yUupFELueNhraWgrJrN+WX09UIAACAASURBVOVszi9jy94yNuWV8cOesh9nHAN9vUlNCKN/UgSDOkSS1r4NUdqIWESO0/4DNVz/ZgbLdxRxx5iu3Hpa51b9DbuINJ6aOjcnPTKfbm1DeeOaIU7HaXTGmAxrbdpPxzWTJ3IIFTV1LNhUwGfr8li4pYDiivrNNVPjw7h6RAeGdIykX1IbIoObzzIjYwyxYQHEhgVwUpfoH8ddbsv2gnLW5pbUX3JKmL50Fy8t2gFAp5hgBneIJK19JIM7RJLYJlAfzkTkiO0oPMDVr60gd38lT07sx7h+rfs8GRFpXH4+XlwxPIVHP9/ExrxSurcNczqSIzSTJ+JRVlXL/I17+WxtHgs276Wq1k1UsB+jusUysms0IzpHE91KZrWq61yszSlhxc79rNhZRPrOIko9Sz3bhgUwonM0I7tGc3KXmGZV6IrIibVs+z6ufzMDL2OYNnkgaSmRTkcSkVaguKKGYX+bzzl94nl0Ql+n4zSqn5vJ+8UizxgTAHwH+FM/8zfbWvuAMaYDMAOIAjKAydbaGmOMPzAdGAjsAy6x1u70PNa9wDWAC7jNWvu5Z3ws8CTgDbxkrf37L70gFXnSENxuy9Id+5i1IpvP1uVRXecmNtSfsb3acmaveAZ3iNQ5atT/OW3KLyN9ZxFLtxexaGshJZW1GAO924UzsksMI7vG0D85Qls8iAgAczJyuGfuGpIjg3jlykG0jwp2OpKItCL3v7+OmSuyWXTPqcSGBjgdp9EcT5FngGBrbbkxxhdYBNwO3AHMtdbOMMY8D6y21j5njLkJ6GOtvcEYMxG4wFp7iTEmFXgHGAwkAF8BXT1PsxkYA+QAK4BJ1toNh8ulIk+Ox+7iSuZk5PBuRg5ZRRWEBvgwrl8C5/drx4DkNnipsDssl9uyJqeYhVsK+W5zAauyi3G5LaH+PozoHM2Y1DhO6x5LG83yibQ61loe/3IzT83fyvBOUTx32UDCg5rmucoi0nJtLyhn9L++5dZTO3PH6d2cjtNojvmcPFtfBZZ7rvp6LhY4DbjUM/468CDwHDDO8zPAbOAZT6E4Dphhra0GdhhjtlJf8AFstdZu9wSd4Tn2sEWeyLFYmbWflxfu4LN1e3BbGNYxijvGdGVsr7YE+KrL25Hy9jL0T25D/+Q23Da6CyWVtSzZVsi3mwv5ZuNe5q3Pw9vLMDglktN7xjEmNY7ENtoHS6Slq6p1cdfsNXy4ejeXpCXx5wt6aXZfRBzRMSaE0d3jeGPpLm46tXOr+5x3RI1XjDHe1C/J7Az8G9gGFFtr/9OPPQf4z5nU7YBsAGttnTGmhPolne2ApQc97MH3yf7JeMtvhSMnTJ3LzRcb8nlp4XZWZhUTGuDDdSd35LIh7UmOUuHREMIDfRnbK56xveKx1rI2t4Qv1ufzxYY8/vTRBv700QZS48M4vWccZ/WOp2tcqNORRaSB7S2t4vo3M1iVVcw9Z3bn+pEd1aRJRBx17ckd+OqHfOauzOXSIclOxzmhjqjIs9a6gH7GmAjgPaB7o6b6GcaYqcBUgOTk1vWLkqNX63IzJyOHfy/YSnZRJe2jgnjw3FQmpCUR7K/Gso3FGEOfxAj6JEbw2zO6saPwAF9uyOPLDfk8+fUWnvhqC13jQji3TwLn9E2gQ7TO0xFp7lZnF3P9GxmUVtXy/OUDGNsr3ulIIiIM6RBJr3ZhvLxoOxMHJbWq03GO6pOutbbYGPMNMAyIMMb4eGbzEoFcz2G5QBKQY4zxAcKpb8Dyn/H/OPg+Pzf+0+efBkyD+nPyjia7tB51LjdzV+Xy9PwtZBdV0i8pgj+cncqvesSpiYoDOkQHM3VkJ6aO7ERBWTWfrdvDx6v38M8vN/PPLzfTq10Y5/RJ4Oze8SRFamZVpLn5IDOXu2avISbUnzk3DqdHfOtsVy4iTY8xhmtP6sivZ2by7eYCTu3eejZHP5LGKzFArafACwS+AB4BrgDmHNR4ZY219lljzM1A74Mar4y31l5sjOkJvM1/G698DXQBDPWNV0ZTX9ytAC611q4/XC41XpGfcrktH2Tm8tTXW9i5r4I+ieH8ZkxXRnWN0ZKhJmhPSSWfrNnDR2v2sDq7GID+yRGM79+Oc/smEBGkpi0iTZnLbfnH5xt54dvtDOkQyXOXD9SWKiLS5NTUuRn5j2/oFBvMW9cOdTpOgzue7pp9qG+s4g14AbOstQ8ZYzpSv4VCJLAKuNxaW+3ZcuENoD9QBEw8qKnK74GrgTrg19bazzzjZwFPeJ7jFWvtX37pBanIk4Mt2baPhz/ewIY9paTGh3HHmK6M7hGr4q6ZyNpXwcdrd/PBqt1syi/Dz9uL0T1iuXBAIqd0i1HjBpEmprSqltvfWcU3mwq4fGgyD5zbU/9ORaTJem7BNh6Zt5HPbj+5xa02OOYir6lSkSdQXxz89dMfmLc+j3YRgdxzZnfO6ROv4q6Zstayfncpc1bm8EHmbooO1BAd4sd5fdtx4cB29EwIdzqiSKu3o/AA176+gl37KnjwvJ5cPrS905FERA6ruKKGIX/9mvEDEvnb+N5Ox2lQKvKkRTlQXcfT87fyyqId+HgbbhrViWtP7tjq2uO2ZLUuNws2FTAnI4evN+ZT67J0bxvKRQMTGdevHTGh/k5HFGl1vtm4l9tnrMLby/Dc5QMZ2jHK6UgiIkfkrtmr+Wj1HpbeN5rwwJazd6eKPGkx5m/M5w/vrWN3SRUXDkjkrrHdiAsLcDqWNKL9B2r4aM1u5mTksDqnBG8vw2ndY7l0cDIju8aooY5II3O7LU9+vYWn5m+he9swpk0eqEZJItKsrMst4ZynF3H/Oalcc1IHp+M0GBV50uwVlFXz0Mcb+Gj1brrGhfC38X0Y2L6N07HkBNuSX8bsjBzmrMyhsLyGhPAAJqQlcfGgJNpFBDodT6TFKa6o4TczM/lmUwEXDkjkz+f3ItBPqyZEpPm58LnF7CuvZv6do1rMdgoq8qTZstbybkYOf/nkByprXNxyWmduOKUTfj46yb81q6lz89UP+byzPItFWwsBGNU1homDkzmte6yaQIg0gHW5Jdz4VgZ5JVU8cG5PLhuSrHOeRaTZ+iAzl9tnZPLaVYMY1a1lbKegIk+apb1lVdw1ew0LNhUwKKUNfxvfh86xIU7HkiYmu6iCmSuymZWezd6yamJD/ZmQlsjEQclaUiZyjGZn5PD799YSGezHs5cNoH+yVk6ISPNWU+dm+N/n0ycxnFeuHOR0nAahIk+ana825HPXnDUcqK7jvrN6MHlo+xYztS6No87l5ptNBcxYnsU3m/bitjCyawyTh7bntO6xOndP5AhU17l46KMNvLUsi2Edo3j60v5Eh6jRkYi0DP/6cjNPz9/Ct789leSo5v9FsIo8aTYqaur48yc/8PayLFLjw3hyYj+6xIU6HUuamT0llcxckc07y7PIL62mXUQglw5J5pJBSfrAKvIzsosquOWdVazOLub6Uzryu9O74aOlzyLSguSVVDHikflcc1IH7jurh9NxjpuKPGkW1uWWcNs7q9ix7wBTT+7IHad3xd9HJ/jLsat1ufn6h3zeWLqL77fuw9fbcFbveCYPbc/A9m10fpGIx6dr93D3nDUAPHpRH8b2inc4kYhI47j5rZUs2lrI0ntHN/tGUj9X5Pk4EUbkUGatyOYPH6wjMsiPt64dwvBO0U5HkhbA19uLsb3iGdsrnq17y3lz6S7mZNRvtt69bSiTh7Xn/H7tCPbX26G0TlW1Lh7+uH55Zr+kCJ6e1F/nsopIizZlWHs+WbuHD1fncsmgZKfjNArN5InjqmpdPPjhemasyOakztE8ObEfUVpOJ42ooqaODzJ3M33JLn7YU0qovw8XDkzk8qHJdI7V0mBpPbbuLeeWt1eyMa+M60d25LdndFNnWhFp8ay1nPnkQryM4ZPbTmrWq3q0XFOapOyiCm58K4N1uaXccmpnfjOmq5pjyAljrWVl1n7eWLKLT9fmUeNyM6xjFJOHtWdMapw+7EqLNjsjh/vfX0egnzf/vLgvp7aQduIiIkfi7WVZ3PfeWubcOIyB7SOdjnPMVORJk/Pd5gJum7EKl9vy+MX9+FVqnNORpBUrLK9mVno2by3NIre4krgwfyYOSubSIcnEhQU4HU+kwZRX1/HH99cxd1UuQztG8uTE/vo7LiKtzoHqOob89WvO6NmWf17c1+k4x0xFnjQpbyzZyQMfrqdrXCjPXz6QlOhgpyOJAOByW77ZuJfpS3fx3eYCfLwMZ/Rsy+Rh7RnSIbJZL+kQWbGziDtnrSZnfwW3je7Crad10eoJEWm1fv/eWmZn5LD8vl8RHuTrdJxjosYr0iS43JY/f7KBV7/fyejusTw1qb8aXkiT4u1l+FVqHL9KjWNn4QHeXLqLWenZfLJ2D13jQpg8tD0XDEgkRH9vpRmprnPx+JdbeOG7bSS2CWTG1GEM7tB8lyeJiDSESYOTeWtZFu9n5nLF8BSn4zQozeTJCVNeXcft76zi6417uXpEB35/dg99gyzNQmWNi49W72b60p2syy0lxN+H8QPaMXloe+3hKE3exrxSfj0jk415ZUwclMQfzknVlxQiIh7nPbOI6lo38359crNcraOZPHHU7uJKrnk9nc35ZTw8rieTh6U4HUnkiAX6eXPxoCQmpCWyKruYN5fsYsbybKYv2cWwjlFMGdaeX6lRizQxLrflpYXb+ecXmwkL9OGlKWk691lE5CcmDU7m3rlrWZlVzMD2bZyO02A0kyeNbkt+GZNfXk55dR3PXNqfUergJi3AvvJqZv6kUculg9szaXASsWpiIQ7LLqrgzlmrWb6ziDN6xvHXC3praxoRkUM4UF3H4L98xZm943lsQvNrwKLGK+KI1dnFXPnqcry9vJh+9WBSE8KcjiTSoA7VqGVsr7ZMGZbCoJQ2zXLphzRfbrflzWW7eOSzjXgZw4Pn9WT8gHb6eygichj3vbeWuStzWHbfrwgPbF4NWLRcU064xVsLuW56OpEhfrx5zRDaR6mDprQ8Bzdq2eFp1PJuejYfr9lD97ahXD60PRf0b6cGQ9Lotu4t5965a1ixcz8nd4nmb+N7k9gmyOlYIiJN3qWDk3l7WRbvr2o5DVg0kyeNYt66PG57ZxUp0UG8cc0Q7cEkrUpljYsPMnOZvmQXG/aUEurvw4UDE7l8aHs6x4Y4HU9amFqXm2nfbefJr7YQ6OfN/eekcqFm70REjsq5Ty+i1uXms9ubVwMWzeTJCfNuejZ3z1lD36QIXr1yEBFBfk5HEjmhAv28mTg4mUsGJbEyq5g3luzk7WVZvLZ4JyM6RzF5aHt+1SMOHzVqkeO0Mms/v39vHT/sKeXs3vE8eF5PYkJ17p2IyNGaNDiZ+95by6rsYgYkN/8GLJrJkwY1c0UWd89Zy8ldonlh8kCC/PQ9gghAYXk1M1dk8/ay/zZquSQtiYsHJWlJnRy1/QdqeGTeRmasyKZtWAB/GteTM3q2dTqWiEizVV5dx5C/fMVZveN5tBk1YFHjFWl0s1Zkc/fcNYzsEsMLkwcS4OvtdCSRJsfltnz9Qz7vLM9iweYCAEZ2iWHS4GRG94jVNgxyWG63ZVZ6No/M20hpVR1Xj0jh9l911b53IiIN4N65a3lvVfNqwKLlmtKoZqXXF3gnq8ATOSxvL8PpPdtyes+25BZXMmtFNrPSs7nhzQxiQv2ZMDCRiYOSSY7S7J78r9XZxTz40XpWZRUzKKUND5/fi+5t1bFYRKShTBqcxDvLs/hw9W4mD23vdJzjopk8OW7vpmdz15w1nNQ5mhenpKnAEzlKdS43324u4J3l2czfmI/bwkmdo5k0OJkxqXH4+Wh2rzXbU1LJo/M2MXdVLtEhftxzZg81VhERaQTWWs58ciH+Pl58cMtJTsc5IprJk0YxOyNHBZ7IcfLx9mJ0jzhG94hjT0kl76bnMHNFNje/vZLIYD/O65vARQMT6ZkQpg/2rUhFTR0vfLudF77bhtvCjaM6cdOoToQGNI8lRCIizY0xhosGJvLnT35gc34ZXeNCnY50zDSTJ8fs07V7uOXtlQzvFM1LV6jAE2lILrdl4ZYC3k3P4csN+dS43HRvG8qFAxIZ1z+B2FBtS9JS1bnczFmZw+NfbiGvtIqz+8Rzz9juJEVqCa+ISGPbV17NkL9+zVUjUvj92alOx/lFarwiDWrRlkKufm0FvRPDefOaIQT6qcATaSwlFbV8tGY3szNyyMwuxtvLcErXGC4ckMjoHrH6gqWFcLstH6/dw+NfbmZH4QH6JUXwh7N7kJYS6XQ0EZFWZer0dFZm7WfJvaObfEM0LdeUBpOZXczUN9LpGBPMK1cMUoEn0sjCg3y5fGh7Lh/anq17y5m7Moe5K3OZv3ElYQE+nNcvgfEDEumfFKHlnM2QtZavf9jLY19sYmNeGd3bhvLSlDRG94jV71NExAET0pL4YkM+CzYVMCY1zuk4x0QzeXJUtu4t46LnlxAa4MOcG4YTG6YlYyJOcLkti7cVMjsjh8/X51FV6yYpMpBz+yRwbt8EurcNVYHQxLndls/X5/HvBVtZl1tKSlQQvxnTlXP7JODlpd+diIhTal1uhv3tawYkt2HalP83SdakaCZPjlvO/gouf2k5Pl5evHnNEBV4Ig7y9jKc3CWGk7vEUFZVy7x1eXy0Zg8vfLedZxdso3NsCOf2SeCcvvF0iglxOq4cpNbl5oPM3Ty3YCvbCg6QEhXEIxf2ZvyAxCa/LEhEpDXw9fbigv7tePX7nRSWVxMd4u90pKOmmTw5IvsP1HDh84spKKtm5tRhpCZobyaRpmhfeTWfrcvjo9W7Wb6zCGshNT6Ms/vEc0bPODrHNt9OYc1dSWUt76Zn8+r3O8ktrqRHfBg3jerEWb3j8dbMnYhIk7I5v4zTH/+O+89J5ZqTOjgd52ep8Yocs+o6F5NfXk5mVjFvXjuEwR3UBECkOcgrqeKTtXv4aPVuMrOLAegYE8zpqW05vWcc/RIjtCzwBNiSX8brS3Yyd2UuFTUuBqW04aZRnRnVLUZLakVEmrBxzyyius7NZ7ef3GTfr1XkyTGx1nLHrNW8tyqXJyf2Y1y/dk5HEpFjsKekkq825PP5+nyWbt9HndsSG+rPmNQ4Tu/ZlqEdI/H3UROlhlJT52b+xnzeXJrFoq2F+Pl4cV7fBK4cnkKvduFOxxMRkSPwxtJd3P/+Oj6+9aQm+959zEWeMSYJmA7EARaYZq190hjzIHAdUOA59D5r7aee+9wLXAO4gNustZ97xscCTwLewEvW2r97xjsAM4AoIAOYbK2tOVwuFXknxhNfbeaJr7Zw55iu3Dq6i9NxRKQBlFTU8s2mvXy+Po8FmwqorHUR6OvN0I6RnNI1hpFdY+gQHdxkv7Vsqqy1rN9dyuyMHD7IzGV/RS1twwKYPKw9EwclEdUMz+kQEWnNSipqGfTXr5g0KIk/jevldJxDOp4iLx6It9auNMaEUl+EnQ9cDJRbax/7yfGpwDvAYCAB+Aro6rl5MzAGyAFWAJOstRuMMbOAudbaGcaY54HV1trnDpdLRV7je29VDr+ZuZoLByTy2IQ++sAn0gJV1bpYvK2QbzcV8N2WQnYUHgAgKTKQkV3qC77hnaIIDfB1OGnTtbPwAPPW5/H+qlw25pXh5+3FmJ5xXDQwkZM7R+OjZioiIs3Wre+sYuGWApbdN7pJrng55u6a1to9wB7Pz2XGmB+Aw63ZGwfMsNZWAzuMMVupL/gAtlprt3sCzQDGeR7vNOBSzzGvAw8Chy3ypHEt31HE3bPXMrRjJH8b31sFnkgLFeDrzWnd4zite/0+QLv2HeC7zQV8u7mQ91bl8tayLLy9DKnxYQzuEMmglEgGpbRp1bNS1lo25Zcxb10e89blsTGvDIC+SRE8fH4vzuuTQHiQimIRkZZgwsBEPlq9m6827OXsPvFOxzliR7WFgjEmBegPLANGALcYY6YA6cCd1tr91BeASw+6Ww7/LQqzfzI+hPolmsXW2rpDHP/T558KTAVITk4+muhyFLL2VTD1jXQSIwN54fI0/Hz0LbRIa9E+KpjJw4KZPCyFmjo3Gbv2s3hbIct3FPHG0l28vGgHAJ1jQxiUEsmQDpEMSG5DUmRgi/4yqOhADYu3FfL91kIWbikkZ38lxsCg9pHcf04qp6fGkRQZ5HRMERFpYCM6RzNhYCJtw5vX1mFHXOQZY0KAOcCvrbWlxpjngIepP0/vYeCfwNWNktLDWjsNmAb1yzUb87laq4qaOqa+kY618OqVg/RttEgr5ufjxbBOUQzrFAXUd9pdm1PC8p1FLN9RxMerd/PO8iwAwgN96ZMYTu929ZeeCeEktglslt07rbXkFleSmV1MZlYxS7bvY/3uUgBC/X0Y2imKG0d1YkxqHLGhzet/+iIicnS8vQyPTujrdIyjdkRFnjHGl/oC7y1r7VwAa23+Qbe/CHzsuZoLJB1090TPGD8zvg+IMMb4eGbzDj5eTiBrLb+bvYbN+WW8dtVg2kcFOx1JRJoQfx9v0lIiSUuJ5KZR4HJbNuaVsjq7hLW5xazJKWHad9upc9d/Bxfo603XuBC6xIXSLS6UzrEhJEcFkdgmsMmc11Bd52JnYQVb95azZW8Za3NKWJ1TTGF5fe8vPx8v+idF8NvTuzK8czR92oXrHDsREWnyfrHIM/Xrb14GfrDW/uug8XjP+XoAFwDrPD9/CLxtjPkX9Y1XugDLAQN08XTSzAUmApdaa60x5hvgIuo7bF4BfNAQL06OzgvfbeeTNXu458zujOwa43QcEWnivL0MPRPqZ+2gfgl9Va2LjXllbNxTyqb8Mjbnl7FgUwGzM3J+vJ8xkBAeSHJkEEmRgcSFBRAbFkBcqD9xYQFEh/oTFuBDiL/PcS0BdbstZdV1lFTUkldaxZ6SSnYXV5FXUklucRXbC8rZVVSBy1OUGgOdYkI4pWss/ZIj6JcYQbe2oVqyLiIizc6RzOSNACYDa40xmZ6x+4BJxph+1C/X3AlcD2CtXe/plrkBqANutta6AIwxtwCfU7+FwivW2vWex7sbmGGM+TOwivqiUk6g7zYX8I95GzmnTzzXj+zodBwRaaYCfL3plxRBv6SI/xkvOlDDjsJydhZWsKuogqx9B9hVVMGCTQUUllfjPsQCfC8DYYG+hAX4EuLvg6+PF37eBl9vL3y8vfDxMtS63NS5LHVuN7UuS02dm7LqWkoqaimrruNQDaTDAnxIiAika1woZ/eJp3NsCJ1i6i+Bfk1jhlFEROR4aDN0Yde+A5z3zPfEhwcw96bhBPkdVT8eEZHj4nJb9pVXk19aTX5pFYXl1ZRV1VFaVUtpZS2lVXWUVdVR63IfdKkv7Hy8vPA9qPDz8zaEBfjWF4eBvoQF+BAe6Evb8ADiwwOJDw8g2F/vcSIi0jIc8xYK0rJV1NRx/RsZAEybnKYCT0ROOG8vQ6xnyWZvwp2OIyIi0uzpE30rd//769mUX8brVw0mOUrtv0VEREREmjudTd6KvZuezZyVOdx6Whc1WhERERERaSFU5LVSm/PLuP+DdQztGMnto7s4HUdERERERBqIirxWqKKmjpvfWkmIvw9PTeyPdzPcrFhERERERA5N5+S1Qn/8YD1bC8p54+ohxIYFOB1HREREREQakGbyWpnZGTnMzsjh1lM7c1KXaKfj/F97dx9tVV0mcPz7cLmCKAgKKsOLkECKL5mSYeZkSIYvRbNyJqeVL2VZkzZazSqrVTQ1q5kax5rKdDnFqE0rLV+SHFyE5Uu10kR8AxG8WAqEgKAio6Jwn/njbuouvMhBOGffs8/3s9ZZ95zf/p1zn/OcH/fsh9/evy1JkiRpF7PIayEdqzfwxZ8V5+FNnVB2OJIkSZLqwCKvRby0qZMLr72P/u19+E/Pw5MkSZIqy3PyWsS3bl3CghXrufwDR7Gf5+FJkiRJleVMXgv4/R/WcdkdS3nfpFFMO3T/ssORJEmSVEcWeRW3/sWX+eS19zN67wF86V0Tyw5HkiRJUp15uGbFfXnWQp5c/yI//dgx7NHPj1uSJEmqOmfyKuzmB//EDfNXcP7bx3Hk6CFlhyNJkiSpASzyKurJZ1/kCzcu4IhRgzl/yriyw5EkSZLUIBZ5FZSZXHTDg7y0qZNvvu8I2tv8mCVJkqRW4d5/BV0/fwW3L17DZ6a9nrFD9yg7HEmSJEkNZJFXMavWv8hXfr6QN40ZwlnHjCk7HEmSJEkNZpFXIZnJF25cwMZNnXzjtDfQp0+UHZIkSZKkBrPIq5BZD/yJWxet4p9O9DBNSZIkqVVZ5FXEmuc2MmPWQt44ejAfeuvYssORJEmSVBKLvIr40k0LeP6lzfz7aYfT5mGakiRJUsuyyKuA2Q+t5JYFT/LJqRMYt+/AssORJEmSVCKLvCa3/sWXmTFrIYeOGMRHjvMwTUmSJKnV9S07AO2ci+csZu2Gjcw860309aLnkiRJUsuzKmhi9y97hh/e9ThnHjOGw0buVXY4kiRJknoBi7wmtWlzJ5+74SH2HdiPT584oexwJEmSJPUSHq7ZpP77t39k0cr1XP6BIxnYv73scCRJkiT1Es7kNaEVz7zAJXOXcMJB+/LOQ/YvOxxJkiRJvYhFXpPJTGbctACAf55+CBFeE0+SJEnSX1jkNZk5C1dx66LVXDh1PCOHDCg7HEmSJEm9jEVeE3nhpc189eaHef1+A/nQW70mniRJkqRXcuGVJnL5HUtZ8cwLXHPuZNq9Jp4kSZKkHmy3UoiIURFxW0Q8HBELI+KCon3viJgbEY8WP4cU7RER346Ijoh4MCKO7PZaZxX9H42Is7q1HxURDxXP+XZ4otkrLFv3PJffsZRTDx/O5NftU3Y4kiRJknqpWqaDNgGfzsyJwGTgvIiYCFwE/DIzxwO/LB4DnASML27nApdBV1EIzADeDBwNzNhSGBZ9M3Qn5gAADKdJREFUPtLtedN2/q1Vy9dmL6JPBJ8/+eCyQ5EkSZLUi223yMvMlZk5v7j/HLAIGAFMB64qul0FvKe4Px24OrvcBQyOiOHAO4G5mbkuM58G5gLTim2DMvOuzEzg6m6vJeC3HU9xy4InOe/tB/JXg3cvOxxJkiRJvdgOndgVEWOANwJ3A/tl5spi05PAfsX9EcCybk9bXrS9WvvyHtoFvLy5ky/PWsjovQfw4eNeV3Y4kiRJknq5mou8iNgTuB64MDPXd99WzMDlLo6tpxjOjYh5ETFvzZo19f51vcLVv3ucR1dv4IunTqR/e1vZ4UiSJEnq5Woq8iKina4C70eZeUPRvKo41JLi5+qifQUwqtvTRxZtr9Y+sof2V8jMKzJzUmZOGjZsWC2hN7WnNmzkW3OX8NcThjH14H3LDkeSJElSE6hldc0AfgAsysxLum2aBWxZIfMs4KZu7WcWq2xOBp4tDuucA5wYEUOKBVdOBOYU29ZHxOTid53Z7bVa2sVzFvPCy5uZ8a6JuOCoJEmSpFrUcp28Y4EzgIci4v6i7fPAvwE/iYhzgMeBvyu2zQZOBjqA54EPAmTmuoj4KnBP0e8rmbmuuP9x4Epgd+CW4tbSFq1cz7XzlvGhY8dy4LA9yw5HkiRJUpPYbpGXmb8BtjWNdEIP/RM4bxuvNROY2UP7PODQ7cXSSv71lkcY1L+dT0wZV3YokiRJkprIDq2uqca4Y8ka7lyyhk9MGcfgAbuVHY4kSZKkJmKR18ts7ky+9r+LGL33AM445oCyw5EkSZLUZCzyepnr7l3G4lXPcdFJB9Gvr5dMkCRJkrRjLPJ6kf/buImLf7GEI0cP5qRD9y87HEmSJElNyCKvF7nizsdY89xGvnCKl0yQJEmS9NpY5PUSq9a/yBV3PsYphw3nqAOGlB2OJEmSpCZlkddLXPKLJWzq7OSz0w4qOxRJkiRJTcwirxfoWL2Bn967jDMmj2H0PgPKDkeSJElSE7PI6wUumbuY3dvbOO/tB5YdiiRJkqQmZ5FXsoeWP8vsh57knONexz579is7HEmSJElNziKvZN+Y8whDBrTzkePGlh2KJEmSpAqwyCvR75au5dePPsXHjx/HwP7tZYcjSZIkqQIs8kqSmXxjziPsP6g/ZxxzQNnhSJIkSaoIi7yS3LpoNfc98QwXTB1P//a2ssORJEmSVBEWeSXY3JlcPGcxY4fuwd8eNbLscCRJkiRViEVeCWY9sILFq57jU++YQN82PwJJkiRJu44VRoO9vLmTb859lInDB3HKYcPLDkeSJElSxVjkNdiN963giXXP86l3TKBPnyg7HEmSJEkVY5HXQC9v7uS7v+rg0BGDOOHgfcsOR5IkSVIFWeQ10JZZvAtPmECEs3iSJEmSdj2LvAbZMot32Ii9nMWTJEmSVDcWeQ3y51m8qeOdxZMkSZJUNxZ5DdB9Fm/KQc7iSZIkSaofi7wGuHG+s3iSJEmSGsMir85e3tzJd2571Fk8SZIkSQ1hkVdnN85fwbJ1LziLJ0mSJKkhLPLqaNPmTr57m+fiSZIkSWoci7w6uvnBlTyx7nk+MWWcs3iSJEmSGsIir046O5Pv3d7BhP32ZOrB+5UdjiRJkqQWYZFXJ7cuWsWSVRv4+PHj6NPHWTxJkiRJjWGRVweZyaW3L2XU3rtz6uHDyw5HkiRJUguxyKuD3y1dywPLnuFjbzuQvm2mWJIkSVLjWIHUwaW3d7DvwH6898iRZYciSZIkqcVY5O1i9z3xNL/tWMuHjxtL//a2ssORJEmS1GK2W+RFxMyIWB0RC7q1fTkiVkTE/cXt5G7bPhcRHRGxOCLe2a19WtHWEREXdWsfGxF3F+3XRsRuu/INNtr3bl/KXru38/43H1B2KJIkSZJaUC0zeVcC03po/2ZmHlHcZgNExETgdOCQ4jnfi4i2iGgDLgVOAiYCf1/0Bfh68VrjgKeBc3bmDZVpyarnmPvwKs5+yxj27Ne37HAkSZIktaDtFnmZeSewrsbXmw5ck5kbM/MPQAdwdHHryMzHMvMl4BpgenRdIXwKcF3x/KuA9+zge+g1Lrt9KQN2a+Pst4wpOxRJkiRJLWpnzsk7PyIeLA7nHFK0jQCWdeuzvGjbVvs+wDOZuWmr9h5FxLkRMS8i5q1Zs2YnQt/1ntqwkZ8/8Cfef/RohuzR1EecSpIkSWpir7XIuww4EDgCWAn8xy6L6FVk5hWZOSkzJw0bNqwRv7JmQ/fsx+wLjuOjbzuw7FAkSZIktbDXdOJYZq7acj8i/gu4uXi4AhjVrevIoo1ttK8FBkdE32I2r3v/pjNhv4FlhyBJkiSpxb2mmbyIGN7t4d8AW1benAWcHhH9ImIsMB74PXAPML5YSXM3uhZnmZWZCdwGnFY8/yzgptcSkyRJkiSphpm8iPgxcDwwNCKWAzOA4yPiCCCBPwIfBcjMhRHxE+BhYBNwXmZuLl7nfGAO0AbMzMyFxa/4LHBNRPwLcB/wg1327iRJkiSpxUTXZFrzmTRpUs6bN6/sMCRJkiSpFBFxb2ZO2rp9Z1bXlCRJkiT1MhZ5kiRJklQhFnmSJEmSVCEWeZIkSZJUIRZ5kiRJklQhFnmSJEmSVCEWeZIkSZJUIU17nbyIWAM8XnYcPRgKPFV2EC3K3JfH3JfH3JfH3JfH3JfH3JfL/Jent+b+gMwctnVj0xZ5vVVEzOvpgoSqP3NfHnNfHnNfHnNfHnNfHnNfLvNfnmbLvYdrSpIkSVKFWORJkiRJUoVY5O16V5QdQAsz9+Ux9+Ux9+Ux9+Ux9+Ux9+Uy/+Vpqtx7Tp4kSZIkVYgzeZIkSZJUIRZ5r0FETIuIxRHREREX9bC9X0RcW2y/OyLGND7K6qoh/2dHxJqIuL+4fbiMOKsmImZGxOqIWLCN7RER3y4+lwcj4shGx1hVNeT++Ih4ttuY/1KjY6yqiBgVEbdFxMMRsTAiLuihj2O/DmrMvWO/DiKif0T8PiIeKHL/zz30cV+nDmrMvfs5dRQRbRFxX0Tc3MO2phn3fcsOoNlERBtwKfAOYDlwT0TMysyHu3U7B3g6M8dFxOnA14H3NT7a6qkx/wDXZub5DQ+w2q4EvgtcvY3tJwHji9ubgcuKn9p5V/LquQf4dWae2phwWsom4NOZOT8iBgL3RsTcrf7mOPbro5bcg2O/HjYCUzJzQ0S0A7+JiFsy865ufdzXqY9acg/u59TTBcAiYFAP25pm3DuTt+OOBjoy87HMfAm4Bpi+VZ/pwFXF/euAEyIiGhhjldWSf9VBZt4JrHuVLtOBq7PLXcDgiBjemOiqrYbcq04yc2Vmzi/uP0fXF/+Irbo59uugxtyrDoqxvKF42F7ctl7EwX2dOqgx96qTiBgJnAJ8fxtdmmbcW+TtuBHAsm6Pl/PKL50/98nMTcCzwD4Nia76ask/wHuLw6aui4hRjQmt5dX62ag+jikO77klIg4pO5gqKg7LeSNw91abHPt19iq5B8d+XRSHrN0PrAbmZuY2x737OrtWDbkH93Pq5VvAZ4DObWxvmnFvkacq+jkwJjMPB+byl/9xkapqPnBAZr4B+A7ws5LjqZyI2BO4HrgwM9eXHU8r2U7uHft1kpmbM/MIYCRwdEQcWnZMraKG3LufUwcRcSqwOjPvLTuWXcEib8etALr/j8nIoq3HPhHRF9gLWNuQ6Kpvu/nPzLWZubF4+H3gqAbF1upq+behOsjM9VsO78nM2UB7RAwtOazKKM6LuR74UWbe0EMXx36dbC/3jv36y8xngNuAaVttcl+nzraVe/dz6uZY4N0R8Ue6TgeaEhH/s1Wfphn3Fnk77h5gfESMjYjdgNOBWVv1mQWcVdw/DfhVekHCXWW7+d/qXJh303Ueh+pvFnBmsdLgZODZzFxZdlCtICL233JOQEQcTdff9l75pdNsirz+AFiUmZdso5tjvw5qyb1jvz4iYlhEDC7u707XYmePbNXNfZ06qCX37ufUR2Z+LjNHZuYYuvYvf5WZH9iqW9OMe1fX3EGZuSkizgfmAG3AzMxcGBFfAeZl5iy6vpR+GBEddC2WcHp5EVdLjfn/x4h4N10rs60Dzi4t4AqJiB8DxwNDI2I5MIOuE8LJzMuB2cDJQAfwPPDBciKtnhpyfxrwDxGxCXgBOL23fuk0oWOBM4CHinNkAD4PjAbHfp3VknvHfn0MB64qVrTuA/wkM292X6chasm9+zkN1KzjPvxbKEmSJEnV4eGakiRJklQhFnmSJEmSVCEWeZIkSZJUIRZ5kiRJklQhFnmSJEmSVCEWeZIkSZJUIRZ5kiRJklQhFnmSJEmSVCH/D5Br2+c203+tAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 1080x360 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "25DypWxsOAOl"
      },
      "source": [
        "output=np.zeros(shape=(n_power_types, n_periods))\n",
        "for t in range(n_power_types):\n",
        "  for p in range(n_periods):\n",
        "    output[t,p]=gen_power[t,p].x\n",
        "\n",
        "hydro=np.zeros(shape=(n_hydro_units, n_periods))\n",
        "\n",
        "for h in range(n_hydro_units):\n",
        "  for p in range(n_periods):\n",
        "    hydro[h,p]=hydro_on[h,p].x\n",
        "\n",
        "pumping_hydro=np.zeros(shape=(n_periods))\n",
        "for p in range(n_periods):\n",
        "  pumping_hydro[p]=pumping[p].x"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kfn5Ti7KQcbY"
      },
      "source": [
        "def satisfies_demand(demandp, output, hydro, pumping_hydro):\n",
        "  for p in range(n_periods):\n",
        "    sum_output=0\n",
        "    for t in range(n_power_types):\n",
        "      sum_output+=output[t, p]\n",
        "\n",
        "    sum_hydro=0\n",
        "    for h in range(n_hydro_units):\n",
        "      sum_hydro+=hydro[h,p]*hydro_load[h]\n",
        "\n",
        "    # print(sum_hydro)\n",
        "    # print(sum_output)\n",
        "    # print(demandp[p])\n",
        "    # print(pumping_hydro[p])\n",
        "    if sum_output+sum_hydro<demandp[p]+pumping_hydro[p]:\n",
        "      return False\n",
        "  return True"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2Yu6POUiQM-U",
        "outputId": "d7fc777f-53b0-46f8-93e1-cd9591eca830"
      },
      "source": [
        "satisfied=0\n",
        "for i in range(len(scenarios)):\n",
        "  if(satisfies_demand(scenarios[i], output, hydro, pumping_hydro)):\n",
        "    satisfied+=1\n",
        "\n",
        "print(f\"The percentage of scenarios satisfied by the determinist solution is {(satisfied/len(scenarios))*100}%\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "The percentage of scenarios satisfied by the determinist solution is 26.0%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7JU_PjGEorry"
      },
      "source": [
        "## Question 3.4\n",
        "\n",
        "For this, we just have to change the meet demand constraint and apply it to 90% of the cases."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BWhmvuglXu8L"
      },
      "source": [
        "sample_train_size=int(0.9*len(scenarios))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cxVnlT-4XcZr",
        "outputId": "33577a27-0f55-4ea9-ba39-f1c3f20417be"
      },
      "source": [
        "meet_demand_constr = model_2_2.addConstrs(gp.quicksum(gen_power[type_id, period_id] for type_id in range(n_power_types)) +\n",
        "                               gp.quicksum(hydro_load[unit]*hydro_on[unit,period_id] for unit in range(n_hydro_units))\n",
        "                               >= scenarios[s][period_id] + pumping[period_id]\n",
        "                               for period_id in range(n_periods) for s in range(sample_train_size))\n",
        "model_2_2.optimize()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "Optimize a model with 1435 rows, 75 columns and 8315 nonzeros\n",
            "Model fingerprint: 0x1c3ed18e\n",
            "Variable types: 25 continuous, 50 integer (20 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-03, 4e+03]\n",
            "  Objective range  [4e+00, 9e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [5e+00, 4e+04]\n",
            "\n",
            "MIP start from previous solve produced solution with objective 1.15113e+06 (0.02s)\n",
            "Loaded MIP start from previous solve with objective 1.15113e+06\n",
            "\n",
            "Presolve removed 1370 rows and 5 columns\n",
            "Presolve time: 0.01s\n",
            "Presolved: 65 rows, 70 columns, 190 nonzeros\n",
            "Variable types: 25 continuous, 45 integer (18 binary)\n",
            "\n",
            "Root relaxation: objective 1.141064e+06, 34 iterations, 0.00 seconds\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 1141064.18    0    7 1151131.02 1141064.18  0.87%     -    0s\n",
            "H    0     0                    1142671.8969 1141064.18  0.14%     -    0s\n",
            "     0     0 1141256.37    0    5 1142671.90 1141256.37  0.12%     -    0s\n",
            "     0     0 1141389.71    0    9 1142671.90 1141389.71  0.11%     -    0s\n",
            "     0     0 1141676.48    0    6 1142671.90 1141676.48  0.09%     -    0s\n",
            "H    0     0                    1141868.8969 1141676.48  0.02%     -    0s\n",
            "     0     0 1141693.54    0    6 1141868.90 1141693.54  0.02%     -    0s\n",
            "     0     0 1141693.54    0    5 1141868.90 1141693.54  0.02%     -    0s\n",
            "     0     0 1141693.54    0    3 1141868.90 1141693.54  0.02%     -    0s\n",
            "     0     0 1141693.54    0    3 1141868.90 1141693.54  0.02%     -    0s\n",
            "     0     0 1141766.90    0    1 1141868.90 1141766.90  0.01%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 1\n",
            "  MIR: 2\n",
            "\n",
            "Explored 1 nodes (73 simplex iterations) in 0.08 seconds\n",
            "Thread count was 2 (of 2 available processors)\n",
            "\n",
            "Solution count 3: 1.14187e+06 1.14267e+06 1.15113e+06 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 1.141868896929e+06, best bound 1.141766896929e+06, gap 0.0089%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 230
        },
        "id": "8JJhbKwzYEBn",
        "outputId": "cc1a756e-5db4-49b6-90e1-7d05a82b4046"
      },
      "source": [
        "describe_model(model_2_2)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Costs: 1141868.8969292203\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "execute_result",
          "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>Power type 0</th>\n",
              "      <th>N units on 0</th>\n",
              "      <th>Power type 1</th>\n",
              "      <th>N units on 1</th>\n",
              "      <th>Power type 2</th>\n",
              "      <th>N units on 2</th>\n",
              "      <th>Hydro unit 0</th>\n",
              "      <th>Hydro unit 1</th>\n",
              "      <th>Total gen power</th>\n",
              "      <th>Demand</th>\n",
              "      <th>Cost Base</th>\n",
              "      <th>Cost Additional</th>\n",
              "      <th>Pumping</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>10200.000000</td>\n",
              "      <td>12.0</td>\n",
              "      <td>6915.0</td>\n",
              "      <td>4.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>17115.000000</td>\n",
              "      <td>15000</td>\n",
              "      <td>134400.0</td>\n",
              "      <td>14937.000000</td>\n",
              "      <td>2115.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>13776.931773</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>31276.931773</td>\n",
              "      <td>30000</td>\n",
              "      <td>114000.0</td>\n",
              "      <td>40961.590636</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>10897.145393</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>28397.145393</td>\n",
              "      <td>25000</td>\n",
              "      <td>228000.0</td>\n",
              "      <td>47365.744717</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>23789.709221</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>1500.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>44189.709221</td>\n",
              "      <td>40000</td>\n",
              "      <td>123000.0</td>\n",
              "      <td>101038.255327</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>13043.025521</td>\n",
              "      <td>12.0</td>\n",
              "      <td>17500.0</td>\n",
              "      <td>10.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>-0.0</td>\n",
              "      <td>1.0</td>\n",
              "      <td>31943.025521</td>\n",
              "      <td>27000</td>\n",
              "      <td>228000.0</td>\n",
              "      <td>73116.306250</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Power type 0  N units on 0  ...  Cost Additional  Pumping\n",
              "0  10200.000000          12.0  ...     14937.000000   2115.0\n",
              "1  13776.931773          12.0  ...     40961.590636      0.0\n",
              "2  10897.145393          12.0  ...     47365.744717      0.0\n",
              "3  23789.709221          12.0  ...    101038.255327      0.0\n",
              "4  13043.025521          12.0  ...     73116.306250      0.0\n",
              "\n",
              "[5 rows x 13 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 365
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7CMO349dgTaW"
      },
      "source": [
        "**a) Costs**\n",
        "\n",
        "About 1,146B€ vs 1.014B€, so around 132k€ of difference.\n",
        "\n",
        "\n",
        "**b) Thermal production**\n",
        "\n",
        "30GW average vs 27GW average over the periods.\n",
        "\n",
        "**c) Hydro production**\n",
        "\n",
        "Same generation.\n",
        "\n",
        "**d) Pumping**\n",
        "\n",
        "Exactly the same total of pumping (2115) but in the last model, it is all done in the beginning."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "LdjT_11wezW1",
        "outputId": "eae1734f-5726-4421-8f49-1b801d5455c3"
      },
      "source": [
        "#Let's reload new scenarios to not test the model on situations it trained for\n",
        "scenarios=generate_scenarios(power_demand_by_period, cov, n_scen=300)\n",
        "\n",
        "output=np.zeros(shape=(n_power_types, n_periods))\n",
        "for t in range(n_power_types):\n",
        "  for p in range(n_periods):\n",
        "    output[t,p]=gen_power[t,p].x\n",
        "\n",
        "hydro=np.zeros(shape=(n_hydro_units, n_periods))\n",
        "\n",
        "for h in range(n_hydro_units):\n",
        "  for p in range(n_periods):\n",
        "    hydro[h,p]=hydro_on[h,p].x\n",
        "\n",
        "pumping_hydro=np.zeros(shape=(n_periods))\n",
        "for p in range(n_periods):\n",
        "  pumping_hydro[p]=pumping[p].x\n",
        "\n",
        "satisfied=0\n",
        "for i in range(len(scenarios)):\n",
        "  if(satisfies_demand(scenarios[i], output, hydro, pumping_hydro)):\n",
        "    satisfied+=1\n",
        "\n",
        "print(f\"The percentage of scenarios satisfied by the determinist solution is {(satisfied/len(scenarios))*100}%\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "The percentage of scenarios satisfied by the determinist solution is 98.33333333333333%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Yz2oe28KfcaX",
        "outputId": "898f2af5-dd58-44f0-e8f2-26e3ed2771e1"
      },
      "source": [
        "sample_train_size=int(0.95*len(scenarios))\n",
        "meet_demand_constr = model_2_2.addConstrs(gp.quicksum(gen_power[type_id, period_id] for type_id in range(n_power_types)) +\n",
        "                               gp.quicksum(hydro_load[unit]*hydro_on[unit,period_id] for unit in range(n_hydro_units))\n",
        "                               >= scenarios[s][period_id] + pumping[period_id]\n",
        "                               for period_id in range(n_periods) for s in range(sample_train_size))\n",
        "model_2_2.optimize()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "Optimize a model with 1510 rows, 75 columns and 8765 nonzeros\n",
            "Model fingerprint: 0xa4a8f7e4\n",
            "Variable types: 25 continuous, 50 integer (20 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-03, 4e+03]\n",
            "  Objective range  [4e+00, 9e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [5e+00, 4e+04]\n",
            "Presolve removed 1445 rows and 5 columns\n",
            "Presolve time: 0.00s\n",
            "Presolved: 65 rows, 70 columns, 190 nonzeros\n",
            "Variable types: 25 continuous, 45 integer (18 binary)\n",
            "\n",
            "Root relaxation: objective 1.150983e+06, 34 iterations, 0.00 seconds\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 1150983.45    0    7          - 1150983.45      -     -    0s\n",
            "H    0     0                    1295989.5438 1150983.45  11.2%     -    0s\n",
            "H    0     0                    1154009.2659 1150983.45  0.26%     -    0s\n",
            "     0     0 1151382.01    0    5 1154009.27 1151382.01  0.23%     -    0s\n",
            "H    0     0                    1152797.5429 1151382.01  0.12%     -    0s\n",
            "     0     0 1151516.82    0    9 1152797.54 1151516.82  0.11%     -    0s\n",
            "     0     0 1151832.94    0    6 1152797.54 1151832.94  0.08%     -    0s\n",
            "H    0     0                    1151994.5429 1151832.94  0.01%     -    0s\n",
            "     0     0 1151849.01    0    5 1151994.54 1151849.01  0.01%     -    0s\n",
            "     0     0 infeasible    0      1151994.54 1151994.54  0.00%     -    0s\n",
            "\n",
            "Explored 1 nodes (49 simplex iterations) in 0.06 seconds\n",
            "Thread count was 2 (of 2 available processors)\n",
            "\n",
            "Solution count 5: 1.15199e+06 1.15199e+06 1.1528e+06 ... 1.29599e+06\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 1.151994542852e+06, best bound 1.151994542852e+06, gap 0.0000%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "lPVA59sxfq7p",
        "outputId": "ca74933e-77f7-4106-8022-ded420c3d38a"
      },
      "source": [
        "#Let's reload new scenarios to not test the model on situations it trained for\n",
        "scenarios=generate_scenarios(power_demand_by_period, cov, n_scen=300)\n",
        "\n",
        "output=np.zeros(shape=(n_power_types, n_periods))\n",
        "for t in range(n_power_types):\n",
        "  for p in range(n_periods):\n",
        "    output[t,p]=gen_power[t,p].x\n",
        "\n",
        "hydro=np.zeros(shape=(n_hydro_units, n_periods))\n",
        "\n",
        "for h in range(n_hydro_units):\n",
        "  for p in range(n_periods):\n",
        "    hydro[h,p]=hydro_on[h,p].x\n",
        "\n",
        "pumping_hydro=np.zeros(shape=(n_periods))\n",
        "for p in range(n_periods):\n",
        "  pumping_hydro[p]=pumping[p].x\n",
        "\n",
        "satisfied=0\n",
        "for i in range(len(scenarios)):\n",
        "  if(satisfies_demand(scenarios[i], output, hydro, pumping_hydro)):\n",
        "    satisfied+=1\n",
        "\n",
        "print(f\"The percentage of scenarios satisfied by the determinist solution is {(satisfied/len(scenarios))*100}%\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "The percentage of scenarios satisfied by the determinist solution is 97.33333333333334%\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K83HLpGqibk4"
      },
      "source": [
        "We get similar results for the 95% reliability and the 90% reliability."
      ]
    }
  ]
}