{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "bDLJBWVorXGZ" }, "source": [ "```{index} single: application; energy systems\n", "```\n", "```{index} single: solver; highs\n", "```\n", "```{index} pandas dataframe\n", "```\n", "```{index} network optimization\n", "```\n", "```{index} stochastic optimization\n", "```\n", "```{index} SAA\n", "```\n", "\n", "# Extra: Two-stage energy dispatch optimization with wind curtailment\n", "\n", "This notebook illustrates a two-stage stochastic optimization problem in the context of energy systems where there are recourse actions that depends on the realization of the uncertain parameters. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# install dependencies and select solver\n", "%pip install -q amplpy matplotlib networkx numpy pandas\n", "\n", "SOLVER = \"highs\"\n", "\n", "from amplpy import AMPL, ampl_notebook\n", "\n", "ampl = ampl_notebook(\n", " modules=[\"highs\"], # modules to install\n", " license_uuid=\"default\", # license to use\n", ") # instantiate AMPL object and register magics" ] }, { "cell_type": "markdown", "metadata": { "id": "-lIzixzsd_do" }, "source": [ "# Preliminaries\n", "\n", "The equilibrium point for the European power network, which operates on alternating current, is at a frequency of 50 Hertz with a tolerance threshold of +- 0.05 Hertz. Power network operators are responsible for maintaining a constant frequency on the network. The way that operators do this is mainly by perfectly balancing power supply and demand.\n", "\n", "But what happens if power supply does not meet demand? When _supply exceeds demand_, then the electrical frequency increases. If _demand exceeds supply_, then the frequency drops. Since power plants are designed to operate within a certain frequency range, there is a risk that they will disconnect from the grid after a period of time. If the frequency deviates too much from the 50 Hertz, then there is a risk that the power plants switch off one after another,potentially leading to a complete power blackout.\n", "\n", "Using conventional power generators such as coal and gas, power supply and demand are often matched by increasing/decreasing the production rate of their power plants and taking generating units on or off line. With the advent of renewable energy sources, it has become increasingly more difficult to match supply and demand. Besides not being controllable power sources, network operators rely on forecasting models to predict the power generated by renewable sources. In practice, this prediction is fairly accurate for solar energy, but wind energy is particularly difficult to predict correctly. \n", "\n", "\n", "\n", "The goal of this notebook is to ensure that power demand meets supply while taking into account wind fluctuations. We will introduce two optimization problems and solve them as stochastic optimization problems. **Read first the [energy dispatch problem](../04/power-network.ipynb) from the Chapter 4 for the preliminaries on power networks and the Optimal Power Flow (OPF) problem**. An important difference from the setting there, is that here we *will not assume that the wind generation is a decision variable*. Instead, wind generation is a random variable, as will be explained later on. We do assume that solar and hydro power are decision variables, since the former is accurately predicated whereas the latter is controllable." ] }, { "cell_type": "markdown", "metadata": { "id": "LeaVlnscEkFx" }, "source": [ "\n", "## Unit Commitment" ] }, { "cell_type": "markdown", "metadata": { "id": "3vJPW_Cyx2-k" }, "source": [ "We now consider another optimization problem relevant for energy systems named the _Unit Commitment (UC)_ problem. UC is an extended version of the OPF problem in which an additional binary decision variable $x_i$ is introduced to decide whether generator $i$ should be activated or not. In this formulation we include a fixed cost $\\kappa_i$ that is incurred for the activation of a generator $i$, which is also added as a new parameter `c_fixed` to the network instance. \n", "\n", "In practice, the UC problem is often considered as a two-stage problem. Since it takes time to activate gas and coal generators before they can produce energy, this decision must be made in advance, i.e., as here-and-now decision in the first stage. In the second stage, we then decide the power generation levels of the (activated) coal and gas generators. Note that, in particular, *we cannot produce from generators we did not already turn on!* Lastly, the UC still includes the same physical constraints as the OPF, i.e., power generation limits and line capacity constraints.\n", "\n", "All solar and hydro generators are activated by default. Solar power is treated as deterministic in the sense that in the instance that you are given it holds that $p_{\\min}=p_{\\max}$ for solar generators. Hydro generators are still controllable within their nontrivial generation limits.\n", "\n", "The uncertainty will be described in terms of the _wind speed_ $v$ instead of the wind power. To determine the wind power resulting from a given wind speed, you need to use a so-called _power curve_ $g_i(\\cdot)$ for a wind generator $i$, which is an explicit function that maps the wind speed $v$ to a wind power $g_i(v)$. In the network instance that you are given, there are two wind generators, one in node 64 which is an on-shore wind park, and one in node 65 which is a off-shore wind park. Being structurally different, they have different power curves. See below a plot of the power curve function for the on-shore wind generator. \n", "\n", "\n", "
\n", "\n", "\n", "An analytical description of the power curves is given as follows\n", "\n", "$$\n", "\\begin{align}\n", "g_{64}(v) = \n", "\\begin{cases}\n", "0 & \\text{if } v \\leq 3 \\\\ \n", "0.16563 \\cdot v^3 - 4.4718 & \\text{if } 3 \\leq v \\leq 14 \\\\\n", "450 & \\text{if } v \\geq 14\n", "\\end{cases}\n", "\\end{align}\n", "$$\n", "\n", "$$\n", "\\begin{align}\n", "g_{65}(v) = \n", "\\begin{cases}\n", "0 & \\text{if } v \\leq 3.5 \\\\ \n", "0.18007 \\cdot v^3 - 7.72049 & \\text{if } 3.5 \\leq v \\leq 15 \\\\\n", "600 & \\text{if } v \\geq 15\n", "\\end{cases}\n", "\\end{align}\n", "$$\n", "\n", "We now present the two-stage UC problem formulation. In the first stage, we decide which coal and gas generators to activate (all wind, solar and hydro generators are active by default). In the second stage, the exact power output of the wind parks is calculated using the power curves knowing the realization of the two wind speeds and the power outputs of the already activated generators needs to be adjusted to match demand exactly.\n", "\n", "$$\n", "\\begin{align}\n", "\\begin{array}{llll}\n", "\\min & \\sum_{i \\in \\mathcal{G}^{\\text{coal}} \\cup \\mathcal{G}^{\\text{gas}}} \\kappa_i x_i + \\mathbb{E}_v Q(x, v) \\\\\n", "\\text{s.t.} & x_i \\in \\{0,1\\} & \\forall i \\in \\mathcal{G}^{\\text{coal}} \\cup \\mathcal{G}^{\\text{gas}}, \n", "\\end{array}\n", "\\end{align}\n", "$$\n", "\n", "where\n", "\n", "$$\n", "\\begin{align}\n", "\\begin{array}{lllll}\n", "Q(x,v) := &\\min & \\sum_{i \\in \\mathcal{G}^{\\text{coal}} \\cup \\mathcal{G}^{\\text{gas}}} c_i(p_i) \\\\\n", "&\\text{s.t.}\n", "& p_i = g_i(v_i) & \\forall i \\in \\mathcal{G}^{\\text{wind}}\\\\\n", "&& x_i p_{i}^{\\min } \\leq p_{i} \\leq x_i p_{i}^{\\max } & \\forall i \\in \\mathcal{G}^{\\text{coal}} \\cup \\mathcal{G}^{\\text{gas}} \\\\\n", "&& p_{i}^{\\min } \\leq p_{i} \\leq p_{i}^{\\max } & \\forall i \\in V \\setminus (\\mathcal{G}^{\\text{coal}} \\cup \\mathcal{G}^{\\text{gas}} \\cup \\mathcal{G}^{\\text{wind}}) \\\\\n", "&& \\sum_{j: (i, j) \\in E} f_{ij} - \\sum_{j: (j, i) \\in E} f_{ji} = p_i - d_i & \\forall \\, i \\in V\\\\\n", "&& f_{ij} = b_{ij}(\\theta_i - \\theta_j), & \\forall \\, (i, j) \\in E \\\\\n", "&& -f_{ij}^{max} \\leq f_{ij} \\leq f_{ij}^{\\max} & \\forall (i, j) \\in E\\\\\n", "&& \\theta_i \\in \\mathbb{R} & \\forall i \\in V \\\\\n", "&& f_{ij} \\in \\mathbb{R} & \\forall (i, j) \\in E\\\\\n", "&& p_{i} \\geq 0 & \\forall i \\in V\n", "\\end{array}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "id": "BOJKv8sSuuV1" }, "source": [ "## Package and data import" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "id": "bxwjVztvKQpk" }, "outputs": [], "source": [ "# Load packages\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import networkx as nx\n", "import time\n", "\n", "# Download the data\n", "base_url = (\n", " \"https://raw.githubusercontent.com/ampl/mo-book.ampl.com/dev/notebooks/10/\"\n", ")\n", "nodes_df = pd.read_csv(base_url + \"nodes.csv\").set_index(\"node_id\")\n", "edges_df = pd.read_csv(base_url + \"edges.csv\").set_index([\"node_id1\", \"node_id2\"])\n", "\n", "# Replace 'na' by an empty string\n", "nodes_df.fillna(\"\", inplace=True)\n", "\n", "network = {\"nodes\": nodes_df, \"edges\": edges_df}" ] }, { "cell_type": "markdown", "metadata": { "id": "DYkkFfFaKMRG" }, "source": [ "### Network data\n", "\n", "This notebook uses the same network as in [energy dispatch problem](../04/power-network.ipynb) and hence the same data structure. In particular, the nodes and edges of the network are stored in the `nodes` and `edges` dataframes, respectively. The `edges` dataframe contains the following columns: \n", "* `node_id1`: first node of a given edge;\n", "* `node_id2`: second node of a given edge;\n", "* `f_max`: describing the maximum power flow capacity of the edge; and\n", "* `b`: the susceptance of the edge." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "id": "wKwPFK-PKi7-" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
bf_max
node_id1node_id2
0110.0100270.705509
223.5849415.734756
34125.3133265.273978
249.2593400.159230
4518.5185217.852748
............
646528.90591200.000000
676828.90591200.000000
807928.90591200.000000
86854.8216602.814908
67115246.91361200.000000
\n", "

179 rows × 2 columns

\n", "
" ], "text/plain": [ " b f_max\n", "node_id1 node_id2 \n", "0 1 10.0100 270.705509\n", " 2 23.5849 415.734756\n", "3 4 125.3133 265.273978\n", "2 4 9.2593 400.159230\n", "4 5 18.5185 217.852748\n", "... ... ...\n", "64 65 28.9059 1200.000000\n", "67 68 28.9059 1200.000000\n", "80 79 28.9059 1200.000000\n", "86 85 4.8216 602.814908\n", "67 115 246.9136 1200.000000\n", "\n", "[179 rows x 2 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "edges_df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The network includes 18 generators of different type, which is described in the `energy_type` field. We distinguish between conventional generators (coal, gas) and renewable generators (hydro, solar, wind). Every conventional generator node has two parameters:\n", "\n", "- `c_fixed`, describing the activation cost of the conventional generators;\n", "- `c_var`, describing the unit cost of producing energy for each conventional generator\n", "\n", "Renewable generators are assumed to have zero marginal cost and zero activation cost. Nodes `64` and `65` correspond to the two wind generators that this notebook focuses on." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
dp_minp_maxc_varis_generatorenergy_typec_fixed
node_id
044.2303440.00.00False0
117.6903990.00.00False0
235.5142090.00.00False0
335.2489770.00.00False0
40.0000000.00.00False0
........................
1136.7848560.00.00False0
11419.6640450.00.00False0
115163.0092060.00.00False0
11616.1756330.00.00False0
11731.5044590.00.00False0
\n", "

118 rows × 7 columns

\n", "
" ], "text/plain": [ " d p_min p_max c_var is_generator energy_type c_fixed\n", "node_id \n", "0 44.230344 0.0 0.0 0 False 0\n", "1 17.690399 0.0 0.0 0 False 0\n", "2 35.514209 0.0 0.0 0 False 0\n", "3 35.248977 0.0 0.0 0 False 0\n", "4 0.000000 0.0 0.0 0 False 0\n", "... ... ... ... ... ... ... ...\n", "113 6.784856 0.0 0.0 0 False 0\n", "114 19.664045 0.0 0.0 0 False 0\n", "115 163.009206 0.0 0.0 0 False 0\n", "116 16.175633 0.0 0.0 0 False 0\n", "117 31.504459 0.0 0.0 0 False 0\n", "\n", "[118 rows x 7 columns]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nodes_df" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "id": "lvjm6z_RqIWM" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
dp_minp_maxc_varis_generatorenergy_typec_fixed
node_id
90.00.000000400.0000000Truehydro0
110.00.000000200.0000000Truehydro0
240.00.000000422.08643128Truecoal1689
250.00.000000227.38437018Truecoal1057
300.00.000000235.30623919Truecoal1837
450.00.000000371.34967519Truecoal1456
480.0227.262510227.2625100Truesolar0
530.097.52601297.5260120Truesolar0
580.0284.753966284.7539660Truesolar0
600.098.69380898.6938080Truesolar0
640.016.05035216.0503520Truewind0
650.077.74725777.7472570Truewind0
790.00.000000360.55486738Truegas1504
860.00.000000445.39709938Truegas1751
880.00.000000298.39068333Truegas2450
990.00.000000440.53405733Truegas2184
1020.00.000000454.66058135Truegas1617
1100.00.000000451.13388334Truegas1627
\n", "
" ], "text/plain": [ " d p_min p_max c_var is_generator energy_type c_fixed\n", "node_id \n", "9 0.0 0.000000 400.000000 0 True hydro 0\n", "11 0.0 0.000000 200.000000 0 True hydro 0\n", "24 0.0 0.000000 422.086431 28 True coal 1689\n", "25 0.0 0.000000 227.384370 18 True coal 1057\n", "30 0.0 0.000000 235.306239 19 True coal 1837\n", "45 0.0 0.000000 371.349675 19 True coal 1456\n", "48 0.0 227.262510 227.262510 0 True solar 0\n", "53 0.0 97.526012 97.526012 0 True solar 0\n", "58 0.0 284.753966 284.753966 0 True solar 0\n", "60 0.0 98.693808 98.693808 0 True solar 0\n", "64 0.0 16.050352 16.050352 0 True wind 0\n", "65 0.0 77.747257 77.747257 0 True wind 0\n", "79 0.0 0.000000 360.554867 38 True gas 1504\n", "86 0.0 0.000000 445.397099 38 True gas 1751\n", "88 0.0 0.000000 298.390683 33 True gas 2450\n", "99 0.0 0.000000 440.534057 33 True gas 2184\n", "102 0.0 0.000000 454.660581 35 True gas 1617\n", "110 0.0 0.000000 451.133883 34 True gas 1627" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nodes_df[nodes_df.is_generator]" ] }, { "cell_type": "markdown", "metadata": { "id": "t971BeuikraX" }, "source": [ "Wind turbines risk to break if the wind speed is too high, so wind turbines need to be _curtailed_ (i.e., deactivated) once the wind speed exceeds a specified maximum speed $v_{\\max}$. An additional second stage decision is therefore needed to decide whether to curtail the wind turbines or not. We include this additional requirement in the second-stage formulation, assuming that we have $T$ samples available.\n", "\n", "Assume that \n", " - $v_{64} \\sim \\text{Weibull}$ with scale parameter 15 and shape parameter 2.6.\n", " - $v_{65} \\sim \\text{Weibull}$ with scale parameter 18 and shape parameter 3.0. \n", "\n", "We sample $T=100$ data points from these two distributions. Using these data points, compute the average wind speeds. Then compute the objective value. Analyze the differences in the objective value and explain the differences." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "id": "_uE_NmVgKdlU" }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Power Output')" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABerklEQVR4nO3dd3gU1f7H8fduOqkESAESDJ1A6AgRpBnpCoIFLiIoNm5QAbFgQ/FesVyv7afgVQG9tmtXUEBEmhA6SJUmEFoSWirpO78/1iyuoSSQZJLN5/U8++zsnDOz38nssl/OnDnHYhiGgYiIiIiLspodgIiIiEh5UrIjIiIiLk3JjoiIiLg0JTsiIiLi0pTsiIiIiEtTsiMiIiIuTcmOiIiIuDR3swOoDGw2G0ePHsXf3x+LxWJ2OCIiIlIChmGQkZFB3bp1sVrP336jZAc4evQoERERZochIiIil+DQoUPUr1//vOVKdgB/f3/A/scKCAgwORoREREpifT0dCIiIhy/4+ejZAccl64CAgKU7IiIiFQxF+uCog7KIiIi4tKU7IiIiIhLU7IjIiIiLk19dkrIZrORl5dndhim8fDwwM3NzewwRERESk3JTgnk5eWxf/9+bDab2aGYKigoiLCwMI1FJCIiVYqSnYswDINjx47h5uZGRETEBQctclWGYXDmzBlSUlIACA8PNzkiERGRklOycxEFBQWcOXOGunXrUqNGDbPDMY2Pjw8AKSkphISE6JKWiIhUGdWvmaKUCgsLAfD09DQ5EvMVJXv5+fkmRyIiIlJySnZKSP1U9DcQEZGqScmOiIiIuDTTk50jR45w6623UqtWLXx8fIiJiWH9+vWOcsMweOqppwgPD8fHx4e4uDj27NnjtI9Tp04xcuRIAgICCAoKYuzYsWRmZlb0oYiIiEglZGqyc/r0abp27YqHhwfz589nx44dvPzyy9SsWdNR58UXX+T1119n5syZrFmzBl9fX/r27UtOTo6jzsiRI9m+fTuLFi1i3rx5LF++nLvvvtuMQxIREZFKxmIYhmHWmz/66KOsXLmSFStWnLPcMAzq1q3Lgw8+yOTJkwFIS0sjNDSUOXPmMHz4cHbu3El0dDTr1q2jY8eOACxYsIABAwZw+PBh6tate9E40tPTCQwMJC0trdhEoDk5Oezfv5+oqCi8vb0v84gr1pEjR3jkkUeYP38+Z86coXHjxsyePdvxd/qze++9l7fffptXXnmFCRMmnHN/VflvISJVWE6a/SFVm384uHmU6S4v9Pv9Z6beev7dd9/Rt29fbrrpJpYtW0a9evX4+9//zl133QXA/v37SUpKIi4uzrFNYGAgnTt3JiEhgeHDh5OQkEBQUJDTD3hcXBxWq5U1a9Zwww03FHvf3NxccnNzHa/T09PL8SjNUdRq1qtXL+bPn0+dOnXYs2ePU6tZka+//prVq1eXKDEUEalQmz6CufeDrcDsSORyjd8AtRub8tamJju///47M2bMYNKkSTz22GOsW7eO+++/H09PT0aPHk1SUhIAoaGhTtuFhoY6ypKSkggJCXEqd3d3Jzg42FHnr6ZPn84zzzxzSTEbhkF2fuElbXu5fDzcSnxH1AsvvEBERASzZ892rIuKiipW78iRI9x3330sXLiQgQMHllmsIiKXzVYIPz19NtFxV4tylWbiHb2mJjs2m42OHTvy3HPPAdCuXTu2bdvGzJkzGT16dLm975QpU5g0aZLjdXp6OhERESXaNju/kOinFpZXaBe0Y1pfaniW7JRdrNUM7H//UaNG8dBDD9GyZcvyCltE5NIcXAVZKeAdCA/tK/NLIFJ9mNpBOTw8nOjoaKd1LVq0IDExEYCwsDAAkpOTneokJyc7ysLCwhzTGBQpKCjg1KlTjjp/5eXlRUBAgNPD1RS1mjVp0oSFCxcybtw47r//ft5//31HnRdeeAF3d3fuv/9+EyMVETmPHd/Yn1tcp0RHLoupLTtdu3Zl165dTut2795NgwYNAPtll7CwMBYvXkzbtm0BeyvMmjVrGDduHACxsbGkpqayYcMGOnToAMDPP/+MzWajc+fOZR6zj4cbO6b1LfP9lvS9S+pirWYbNmzgtddeY+PGjRosUEQqH1sh7PjOvhxdvO+lSGmYmuxMnDiRq666iueee46bb76ZtWvX8p///If//Oc/gH3E3gkTJvCPf/yDJk2aEBUVxZNPPkndunUZMmQIYG8J6tevH3fddRczZ84kPz+f8ePHM3z48HLpcGuxWEp8KclM52s1+/LLLwFYsWIFKSkpREZGOsoLCwt58MEHefXVVzlw4EBFhisi4iwx4Y9LWEHQsIfZ0UgVZ+qvdqdOnfj666+ZMmUK06ZNIyoqildffZWRI0c66jz88MNkZWVx9913k5qaSrdu3ViwYIHTrc8fffQR48eP55prrsFqtTJs2DBef/11Mw6p0rhYq9moUaOc7nID6Nu3L6NGjeL222+vsDhFRM5p+9f25+aDdAlLLpvpTRSDBg1i0KBB5y23WCxMmzaNadOmnbdOcHAwH3/8cXmEV2VdrNWsVq1a1KpVy2kbDw8PwsLCaNasmRkhi4jYFRbAjm/tyy11CUsun+nTRUj5KGo1++STT2jVqhXPPvtssVYzEZFK6cAKyDoOPsG6hCVlwvSWHSk/F2s1+yv10xGRSmH7V/bn6Ot1CUvKhFp2RESk8ijIg51z7cuthpkbi7gMJTsiIlJ5/L4Usk+DXyg06Gp2NOIilOyIiEjl4biENQSsJR9bTORClOyIiEjlkJ8DO+fZl1sNNTcWcSlKdkREpHLY+xPkZUBAfah/pdnRiAtRsiMiIpVD0SWslkPAqp8nKTv6NImIiPnysmDXfPuyLmFJGVOyIyIi5tu9EPLPQM0roG57s6MRF6NkR0REzLfNPkkxLYeCxWJuLOJylOy4sCNHjnDrrbdSq1YtfHx8iImJYf369Y7yMWPGYLFYnB79+vUzMWIRqZZy0mHPIvuyBhKUcqDpIlzU6dOn6dq1K7169WL+/PnUqVOHPXv2ULNmTad6/fr1Y/bs2Y7XXl5eFR2qiFR3u+ZDYS7UbgqhLc2ORlyQkh0X9cILLxAREeGUyERFRRWr5+XlRVhYWEWGJiLiTJewpJzpMlZpGYb9rgEzHoZR4jC/++47OnbsyE033URISAjt2rXjnXfeKVZv6dKlhISE0KxZM8aNG8fJkyfL8q8lInJhZ07Bvp/ty7oLS8qJWnZKK/8MPFfXnPd+7Ch4+pao6u+//86MGTOYNGkSjz32GOvWreP+++/H09OT0aNHA/ZLWEOHDiUqKop9+/bx2GOP0b9/fxISEnBz0zDtIlIBfpsHtnwIbQV1mpkdjbgoJTsuymaz0bFjR5577jkA2rVrx7Zt25g5c6Yj2Rk+fLijfkxMDK1bt6ZRo0YsXbqUa665xpS4RaSa2fbHQIJq1ZFypGSntDxq2FtYzHrvEgoPDyc6OtppXYsWLfjyyy/Pu03Dhg2pXbs2e/fuVbIjIuUvIxn2L7Mvt7zB3FjEpSnZKS2LpcSXkszUtWtXdu3a5bRu9+7dNGjQ4LzbHD58mJMnTxIeHl7e4YmI2KeHMGxQvxMENzQ7GnFh6qDsoiZOnMjq1at57rnn2Lt3Lx9//DH/+c9/iI+PByAzM5OHHnqI1atXc+DAARYvXszgwYNp3Lgxffv2NTl6EakWtnxmf4652dw4xOUp2XFRnTp14uuvv+aTTz6hVatWPPvss7z66quMHDkSADc3N7Zs2cL1119P06ZNGTt2LB06dGDFihUaa0dEyt/JfXB0I1jcdAlLyp0uY7mwQYMGMWjQoHOW+fj4sHDhwgqOSETkD1s/tz836gV+dcyNRVyeWnZERKRiGYYuYUmFUrIjIiIV6+hGOLUP3H2g+UCzo5FqQMmOiIhUrC1/XMJqPgC8/MyNRaoFJTsiIlJxCgvOzoWlS1hSQZTslJBRinmpXJX+BiJy2Q4sh6wU8AmGxhq8VCqGkp2LKJojKi8vz+RIzHfmzBkAPDw8TI5ERKqsoktYLW8AN/1bIhVDt55fhLu7OzVq1OD48eN4eHhgtVa//NAwDM6cOUNKSgpBQUGaJFRELk1+Nuyca1+OucncWKRaUbJzERaLhfDwcPbv38/BgwfNDsdUQUFBhIWFmR2GiFRVu+ZDXgYERkJEZ7OjkWpEyU4JeHp60qRJk2p9KcvDw0MtOiJyeYoGEoy5EaphK7mYR8lOCVmtVry9vc0OQ0SkajpzCvYssi+31l1YUrGUWouISPnb8S3Y8iG0FYS0MDsaqWaU7IiISPlzXMJSx2SpeEp2RESkfKUegoMrAYu9v45IBVOyIyIi5WvL/+zPDbpCYH1zY5FqScmOiIiUH8OAXz+1L7cdYW4sUm0p2RERkfJzZAOc3GOf4bzF9WZHI9WUkh0RESk/v35if25xHXgHmBuLVFtKdkREpHwU5J6d4bzNcHNjkWpNyY6IiJSP3Qsh+zT4h0PDnmZHI9WYkh0RESkfRR2TW98MVk03I+ZRsiMiImUv6wTsWWhfbqO7sMRcSnZERKTsbfsSbAUQ3lbTQ4jplOyIiEjZK7oLS606UgmYmuw8/fTTWCwWp0fz5s0d5Tk5OcTHx1OrVi38/PwYNmwYycnJTvtITExk4MCB1KhRg5CQEB566CEKCgoq+lBERKRIym9wdBNY3TU9hFQK7mYH0LJlS3766SfHa3f3syFNnDiR77//ns8//5zAwEDGjx/P0KFDWblyJQCFhYUMHDiQsLAwVq1axbFjx7jtttvw8PDgueeeq/BjERERzrbqNOkDvrXNjUWESpDsuLu7ExYWVmx9Wloa7733Hh9//DG9e/cGYPbs2bRo0YLVq1fTpUsXfvzxR3bs2MFPP/1EaGgobdu25dlnn+WRRx7h6aefxtPTs6IPR0SkerMVwpbP7MsaW0cqCdP77OzZs4e6devSsGFDRo4cSWJiIgAbNmwgPz+fuLg4R93mzZsTGRlJQkICAAkJCcTExBAaGuqo07dvX9LT09m+fft53zM3N5f09HSnh4iIlIH9yyDjKHgHQdN+ZkcjApic7HTu3Jk5c+awYMECZsyYwf79+7n66qvJyMggKSkJT09PgoKCnLYJDQ0lKSkJgKSkJKdEp6i8qOx8pk+fTmBgoOMRERFRtgcmIlJdFY2t02oYuHuZG4vIH0y9jNW/f3/HcuvWrencuTMNGjTgs88+w8fHp9zed8qUKUyaNMnxOj09XQmPiMjlys2AnXPty23/Zm4sIn9i+mWsPwsKCqJp06bs3buXsLAw8vLySE1NdaqTnJzs6OMTFhZW7O6sotfn6gdUxMvLi4CAAKeHiIhcph3fQf4ZqNUY6nUwOxoRh0qV7GRmZrJv3z7Cw8Pp0KEDHh4eLF682FG+a9cuEhMTiY2NBSA2NpatW7eSkpLiqLNo0SICAgKIjo6u8PhFRKq1zR/Zn9uMAIvF3FhE/sTUy1iTJ0/muuuuo0GDBhw9epSpU6fi5ubGiBEjCAwMZOzYsUyaNIng4GACAgK47777iI2NpUuXLgD06dOH6OhoRo0axYsvvkhSUhJPPPEE8fHxeHnpWrGISIU5uQ8OrgSLVQMJSqVjarJz+PBhRowYwcmTJ6lTpw7dunVj9erV1KlTB4BXXnkFq9XKsGHDyM3NpW/fvrz11luO7d3c3Jg3bx7jxo0jNjYWX19fRo8ezbRp08w6JBGR6qmoVafRNRBYz9xYRP7CYhiGYXYQZktPTycwMJC0tDT13xERKS1bIbzSEjKOwU3vQ8shZkck1URJf78rVZ8dERGpgvYutic6PsHQrP/F64tUMCU7IiJyeTb91/7cZrjG1pFKScmOiIhcuqwTsGu+fbndrebGInIeSnZEROTSbfkMbPlQtx2EtjQ7GpFzUrIjIiKXxjBg04f2ZbXqSCWmZEdERC7N0U2Qsh3cvaHVjWZHI3JeSnZEROTSFLXqtLgefIJMDUXkQpTsiIhI6eVnw9Yv7Mu6hCWVnJIdEREpvZ1zITcNgiLhiqvNjkbkgpTsiIhI6RWNrdP2VrDqp0QqN31CRUSkdE4fgP3LAQu01aSfUvkp2RERkdIp6pjcsKf9MpZIJadkR0RESq6w4Gyy02G0ubGIlJCSHRERKbk9P9on/axRG5oNNDsakRJRsiMiIiW3YY79ud1IcPc0NRSRklKyIyIiJZN6CPYusi+31yUsqTqU7IiISMls+hAMG0R1h1qNzI5GpMSU7IiIyMUVFpwdW0etOlLFKNkREZGL2/sTpB8Bn2BocZ3Z0YiUirvZAYiISBVQ1DG57d/A3atYcaHN4FhadsXGJFVKaIA3Hm7mtLEo2RERkQtLOwJ7FtqXO4wpVpySkcPwt1fz+4msio1LqpSfH+xBwzp+pry3kh0REbmwoo7JDbpB7SbFit9dsd+R6Hi5q3eEnJvFYjHtvZXsiIjI+dkKz3ZMPkerTnZeIf9bdwiAd27ryLXRoRUYnEjJKAUXEZHz2/czpB0Cn5rn7Jg899ejpGXnU7+mD72bh5gQoMjFKdkREZHzK+qY3OZv4OHtVGQYBu8nHADg1i4NcLOad5lC5EKU7IiIyLmlH4Nd8+3L55j0c2NiKtuPpuPlbuWWjhEVHJxIySnZERGRc9v4ARiFEBkLdZoVK/7vH60617WpS01fzZMllZeSHRERKa4wHzbMti93urNY8fGMXL7fegyA22IbVGRkIqWmZEdERIrbNR8yjoFvnXN2TP7fukTyCw3aRgTRun5QxccnUgpKdkREpLh179qf299WbMTkgkIbH61JBNSqI1WDkh0REXF2Yg/sXwZYzjm2zk87kzmWlkMtX08GxIRXeHgipaVkR0REnK2fZX9u2g+CIosVv7/qIAC3dIrA28OtIiMTuSRKdkRE5Ky8LNj0kX35HB2Tf0tKJ+H3k7hZLYzsoktYUjUo2RERkbO2fQm5aVDzCmjUu1jx7F8OANC3ZSj1gnwqNjaRS6RkR0RE7AwD1r5jX+44FqzOPxGnsvL4ZvMRAG7vGlXR0YlcMiU7IiJid2QDJG0BNy9od2ux4k/WJpJbYKNVvQA6NqhpQoAil0bJjoiI2BXdbt5qKNQIdirKL7Tx3wR7x+Tbr4rCYtE8WFJ1KNkRERHIOgnbvrIvn6Nj8vxtSSSl51Dbz4tBbXS7uVQtSnZERAQ2fwiFuRDeBup1KFY8e+V+AG7tEomXu243l6pFyY6ISHVns50dW6fTnfCXS1SbD6WyKTEVTzcrIzvrdnOpepTsiIhUd3sXwekD4BUIrW4sVlzUqjOoTTh1/L2KlYtUdkp2RESquzUz7c/tR4FnDaei5PQcvt9in938Dt1uLlWUkh0Rkers+C7Y9zNggSvvKlb84eqDFNgMOl1Rk1b1Ais+PpEyoGRHRKQ6W/sf+3OzAfZRk/8kJ7+Qj/+Y3VyDCEpVpmRHRKS6yk6FzZ/YlzvfU6z4u1+PcjIrj3pBPvSJDq3Y2ETKUKVJdp5//nksFgsTJkxwrMvJySE+Pp5atWrh5+fHsGHDSE5OdtouMTGRgQMHUqNGDUJCQnjooYcoKCio4OhFRKqgTR9CfhaERENUd6ciwzCY9Yu9Y/Ko2Aa4u1WanwuRUqsUn95169bx9ttv07p1a6f1EydOZO7cuXz++ecsW7aMo0ePMnToUEd5YWEhAwcOJC8vj1WrVvH+++8zZ84cnnrqqYo+BBGRqsVWePYSVud7it1u/sveE/yWlEENTzdGdIo0IUCRsmN6spOZmcnIkSN55513qFnz7FwraWlpvPfee/z73/+md+/edOjQgdmzZ7Nq1SpWr14NwI8//siOHTv48MMPadu2Lf379+fZZ5/lzTffJC8vz6xDEhGp/HYvhNSD4B0EMTcXK/7P8t8BuKVTBIE1PCo4OJGyZXqyEx8fz8CBA4mLi3Nav2HDBvLz853WN2/enMjISBISEgBISEggJiaG0NCz15L79u1Leno627dvP+975ubmkp6e7vQQEalWim437zC62O3mO4+ls2LPCawW3W4ursHdzDf/9NNP2bhxI+vWrStWlpSUhKenJ0FBQU7rQ0NDSUpKctT5c6JTVF5Udj7Tp0/nmWeeuczoRUSqqJSdsH8ZWKznnAfr3RX2vjr9W4UTEVyjWLlIVWNay86hQ4d44IEH+Oijj/D29q7Q954yZQppaWmOx6FDhyr0/UVETLXmbftz84EQ5NwfJzk9h+9+PQLAnVerVUdcg2nJzoYNG0hJSaF9+/a4u7vj7u7OsmXLeP3113F3dyc0NJS8vDxSU1OdtktOTiYsLAyAsLCwYndnFb0uqnMuXl5eBAQEOD1ERKqF7NPw66f25c73Fiues+oA+YX2QQTbRdYsVi5SFZmW7FxzzTVs3bqVzZs3Ox4dO3Zk5MiRjmUPDw8WL17s2GbXrl0kJiYSGxsLQGxsLFu3biUlJcVRZ9GiRQQEBBAdHV3hxyQiUult/AAKsiG0FTTo6lSUlVvAR6sPAnDX1Q3NiE6kXJjWZ8ff359WrVo5rfP19aVWrVqO9WPHjmXSpEkEBwcTEBDAfffdR2xsLF26dAGgT58+REdHM2rUKF588UWSkpJ44okniI+Px8tLk9WJiDgpLIC179qXz3G7+WfrD5GeU0BUbV/iWmgQQXEdpnZQvphXXnkFq9XKsGHDyM3NpW/fvrz11luOcjc3N+bNm8e4ceOIjY3F19eX0aNHM23aNBOjFhGppH6bC2mJ4BMMMTc5FRUU2pj1x+zmd3SLwmq1nGsPIlWSxTAMw+wgzJaenk5gYCBpaWnqvyMirskw4N04OLIeuj8MvR93Kv5+yzHiP95IzRoerHr0Gnw83UwKVKTkSvr7bfo4OyIiUgEOrbUnOm6exWY3NwyD/6ywDyI4qksDJTricpTsiIhUBwlv2J9b3wx+IU5F6w+e5tdDqXi6WxkVe0XFxyZSzpTsiIi4ulO/w8559uXY8cWK315mb9UZ2q4edfx1c4e4nlInO8uXLz/nrOIFBQUsX768TIISEZEytHomYEDjOAhp4VS0OzmDn3YmY7HAnbrdXFxUqZOdXr16cerUqWLr09LS6NWrV5kEJSIiZST7NGz60L4cG1+seOayfQD0jQ6jcYhfRUYmUmFKnewYhoHFUvyWxJMnT+Lr61smQYmISBnZMAfysyCkJTR0/g/p4dNn+G7zUQDu7dnIhOBEKkaJx9kZOnQoABaLhTFjxjgN2ldYWMiWLVu46qqryj5CERG5NAV5Z+fBio0vNojguyv2U2AzuKpRLdpGBFV8fCIVpMTJTmBgIGBv2fH398fHx8dR5unpSZcuXbjrrrvOt7mIiFS07V9DxjHwC4WYG52KTmXl8em6RAD+3rOxGdGJVJgSJzuzZ88G4IorrmDy5Mm6ZCUiUpkZxtnbza+8G9yd77Kas+oAOfk2YuoF0rVxLRMCFKk4pZ4uYurUqeURh4iIlKUDKyBpK7j7QMc7nIqycgt4f9UBAMb1bHTOfpgirqTUyU5UVNQFvxi///77ZQUkIiJlYNX/2Z/bjYQawU5Fn6xNJC07n6javvRtGWZCcCIVq9TJzoQJE5xe5+fns2nTJhYsWMBDDz1UVnGJiMilSt4BexYCFujyd6eivAIb766wT/h5T/eGuGnCT6kGSp3sPPDAA+dc/+abb7J+/frLDkhERC7Tytfsz9HXQy3nW8q/2XyEpPQcQvy9uKF9PROCE6l4ZTZdRP/+/fnyyy/LanciInIpUg/Bti/sy10nOBXZbIZjEME7r47Cy10Tfkr1UGbJzhdffEFwcPDFK4qISPlJeBNsBRDVHeq1dyr6cUcSvx/PIsDbnRFXRpoUoEjFK/VlrHbt2jl1UDYMg6SkJI4fP85bb71VpsGJiEgpnDkFG9+3L3eb6FRkGAZv/LwXgNtir8Df26OioxMxTamTnSFDhji9tlqt1KlTh549e9K8efOyiktEREpr7X8g/wyEtS42NcSSXSlsP5pODU837ugWZVKAIubQODsiIq4gL+vs1BDdJjhNDWEYBq8vtrfqjOrSgGBfTxMCFDFPqZMdsM+F9fXXX7Nz504AoqOjGTx4MO7ul7Q7ERG5XBv/C9mnoGYUtBjsVPTL3hNsPpSKt4eVO69uaFKAIuYpdXayfft2rrvuOpKTk2nWrBkAL7zwAnXq1GHu3Lm0atWqzIMUEZELKMyHhD8GEbzqPnBz/qf9jT9adUZcGUkdf6+/bi3i8kp9N9add95Jq1atOHz4MBs3bmTjxo0cOnSI1q1bc/fdd5dHjCIiciHbvoK0Q+AbAm1HOhWt/v0kaw+cwtPNyj3dG51nByKurdQtO5s3b2b9+vXUrFnTsa5mzZr885//pFOnTmUanIiIXIRhwMpX7ctd7gUPb6fi1xfvAeDmTvUJC/RGpDoqdctO06ZNSU5OLrY+JSWFxo0bl0lQIiJSQnt+hJQd4OkPHcc6FW04eIpV+07ibrVwbw+16kj1VepkZ/r06dx///188cUXHD58mMOHD/PFF18wYcIEXnjhBdLT0x0PEREpR4YBK/5tX+54O/gEORUX3YE1rH196tesUcHBiVQeFsMwjNJsYLWezY+KBhcs2sWfX1ssFgoLC8sqznKVnp5OYGAgaWlpBAQEmB2OiEjJ7F8B7w8CNy944FcICHcU/XoolcFvrsTNauHnB3vQoJaviYGKlI+S/n6Xus/OkiVLLiswEREpI8tfsj+3H+WU6ACO0ZIHt6mrREeqvVInO1FRUURERDhNGQH21pxDhw4RGan5VkREyt2htbB/GVjdi034ueNoOj/tTMZigb/3Ul9KkVL32YmKiuL48ePF1p86dYqoKA1BLiJSIYpaddqMgKAIp6JXf9oNwMCYcBqH+FV0ZCKVTqmTnaL+OH+VmZmJt7duaxQRKXdHN9nvwrJY4epJTkVbD6fx4w57q86EuCYmBShSuZT4MtakSfYvlMVi4cknn6RGjbM9+wsLC1mzZg1t27Yt8wBFROQvlv/L/hxzEwQ7T/9Q1KozuE1dGof4V3RkIpVSiZOdTZs2AfaWna1bt+LpeXYiOU9PT9q0acPkyZPLPkIRETkreTv8Ng+wwNUPOhVtSjzN4t9SsFrg/mvUqiNSpMTJTtFdWLfffjuvvfaabtEWETHDipftz9GDoU4zp6JXfrKPljy0fX0a1lFfHZEipb4ba/bs2eURh4iIXMyJPfZ5sAC6O7ekrz9wiuW7j+NmtXB/b7XqiPxZqZOd3r17X7D8559/vuRgRETkAlb8GzCg2QAIi3EqeuWPvjo3dahPZC2NlizyZ6VOdtq0aeP0Oj8/n82bN7Nt2zZGjx5dZoGJiMifnNoPW/5nX77auVVn9e8nWbn3JB5uFuI1ro5IMaVOdl555ZVzrn/66afJzMy87IBEROQcfnkFjEJo1Bvqd3CsNgyDfy+yt+rc3DGCiGC16oj8VanH2TmfW2+9lVmzZpXV7kREpMip/bD5I/tyj0ecilbtO8na/afwdLOqVUfkPMos2UlISNCggiIi5WH5v8BWYG/VieziWP3nVp2/dY6kbpCPWRGKVGqlvow1dOhQp9eGYXDs2DHWr1/Pk08+WWaBiYgIcHIf/PqJfbnnY05FS3cfZ8PB03i5WxnXs5EJwYlUDaVOdgIDA51eW61WmjVrxrRp0+jTp0+ZBSYiItjnwDIKoUkfiOjkWG2zGby4YBcAo7o0IDRALesi56NxdkREKqsTe87egdVzilPR3C1H2XksHX8vd/XVEbmIUic72dnZLFq0iN277deJmzVrRlxcHD4+ulYsIlKmlr0Ahg2a9od67R2r8wttjr46d3dvSE1fz/PtQUQoZbLz3Xffceedd3LixAmn9bVr1+a9997juuuuK9PgRESqrZTfYOsX9uVezq06/1t3iIMnz1Dbz5M7ukWZEJxI1VLiu7FWrVrFjTfeSPfu3Vm5ciWnTp3i1KlT/PLLL1x99dXceOONrF69ujxjFRGpPpY9DxjQfBCEnx3MNTuvkNcX2+fAGt+rMb5epW6gF6l2LIZhGCWpOGDAACIiInj77bfPWX7PPfdw6NAhfvjhhzINsCKkp6cTGBhIWlqaJjgVEfMlb4cZV9mX710JYa0cRTOW7uOFBb9Rv6YPix/sgZe7m0lBipivpL/fJW7ZWb16NePHjz9veXx8PAkJCaUKcsaMGbRu3ZqAgAACAgKIjY1l/vz5jvKcnBzi4+OpVasWfn5+DBs2jOTkZKd9JCYmMnDgQGrUqEFISAgPPfQQBQUFpYpDRKRSWTrd/hw9xCnRSTuTz4ylewGYdG1TJToiJVTiZCc7O/uCWVNgYCA5OTmlevP69evz/PPPs2HDBtavX0/v3r0ZPHgw27dvB2DixInMnTuXzz//nGXLlnH06FGncX4KCwsZOHAgeXl5rFq1ivfff585c+bw1FNPlSoOEZFK49gW2DkXsEDPR52KZi7fR3pOAU1D/Rjctp458YlURUYJxcTEGLNmzTpv+XvvvWfExMSUdHfnVbNmTePdd981UlNTDQ8PD+Pzzz93lO3cudMAjISEBMMwDOOHH34wrFarkZSU5KgzY8YMIyAgwMjNzS3xe6alpRmAkZaWdtnxi4hclg9vMoypAYbx+e1Oq5PTso1mT/xgNHhknvHj9qTzbCxSvZT097vELTu33347kydPPmefnO+//56HH36YMWPGXHLSVVhYyKeffkpWVhaxsbFs2LCB/Px84uLiHHWaN29OZGSk43JZQkICMTExhIaGOur07duX9PR0R+vQueTm5pKenu70EBEx3cEE2LMQLG7FRkt+/ec95OTbaB8ZRFyLEJMCFKmaStyN/4EHHmDVqlUMGjSIZs2a0aJFCwzDYOfOnezZs4chQ4YwYcKEUgewdetWYmNjycnJwc/Pj6+//pro6Gg2b96Mp6cnQUFBTvVDQ0NJSkoCICkpySnRKSovKjuf6dOn88wzz5Q6VhGRcmMY8NPT9uV2t0LtswMF7j+RxadrDwHwSL/mWCwWEwIUqbpK3LJjtVr5/PPP+eSTT2jWrBm//fYbu3btonnz5nz00Ud8+eWXWK2ln1e0WbNmbN68mTVr1jBu3DhGjx7Njh07Sr2f0pgyZQppaWmOx6FDh8r1/URELmrPj3BoNbh7F5vZ/IX5v1FgM+jVrA6dG9YyKUCRqqvUAzTccsst3HLLLWUWgKenJ40b2/8H06FDB9atW8drr73GLbfcQl5eHqmpqU6tO8nJyYSFhQEQFhbG2rVrnfZXdLdWUZ1z8fLywsvLq8yOQUTksthssHiaffnKuyHwbOfjdQdOsWB7ElYLTBnQwqQARaq20jfFlDObzUZubi4dOnTAw8ODxYsXO8p27dpFYmIisbGxAMTGxrJ161ZSUlIcdRYtWkRAQADR0dEVHruIyCXZ9iUkbwOvAOg20bHaMAz++f1OAG7pFEnTUH+zIhSp0kwdenPKlCn079+fyMhIMjIy+Pjjj1m6dCkLFy4kMDCQsWPHMmnSJIKDgwkICOC+++4jNjaWLl26ANCnTx+io6MZNWoUL774IklJSTzxxBPEx8er5UZEqoaCPFjyD/ty1/uhRrCj6Putx9h8KJUanm5MvLaJSQGKVH2mJjspKSncdtttHDt2jMDAQFq3bs3ChQu59tprAXjllVewWq0MGzaM3Nxc+vbty1tvveXY3s3NjXnz5jFu3DhiY2Px9fVl9OjRTJs2zaxDEhEpnU0fwOkD4FsHOo9zrM4tKOSFBb8BcE/3RoT4e5sUoEjVV+LpIlyZposQEVPknYHX20JmMvR/CTrf7Sh6d8Xv/OP7nYT4e7H0oZ7U8NQcWCJ/VebTRQDk5+fj7u7Otm3bLjtAEZFqb81Me6ITFAkdxjhWp57J442f7dNCPNinqRIdkctUqmTHw8ODyMhICgsLyyseEZHqIfs0rHzVvtzrcXD3dBT93897ScvOp1moPzd2iDAnPhEXUuq7sR5//HEee+wxTp06VR7xiIhUDyv+DTlpUKcFxNzkWJ148gzvJxwAYMqA5rhZNYCgyOUqddvo//3f/7F3717q1q1LgwYN8PX1dSrfuHFjmQUnIuKSTh+0X8ICuHYaWM/OXv7iwt/ILzS4ukltejStY1KAIq6l1MnOkCFDyiEMEZFqZPE0KMyDqB7Q5FrH6nUHTjFvyzEsFpjSv4WmhRApI6VOdqZOnVoecYiIVA9HNsC2LwAL9HkW/khobDaDZ+baJzAe3imC6Lq6M1SkrFzSCMqpqam8++67TJkyxdF3Z+PGjRw5cqRMgxMRcSmGAT8+aV9uMxzC2ziKvthwmG1H0vH3cufBPs1MClDENZW6ZWfLli3ExcURGBjIgQMHuOuuuwgODuarr74iMTGRDz74oDziFBGp+nb9AAdX2if77P2EY3VGTj4vLrQPIPhAXBNq+2kEeJGyVOqWnUmTJjFmzBj27NmDt/fZET0HDBjA8uXLyzQ4ERGXUZgPi56yL8fGQ2B9R9H//byXE5l5NKzty22xV5gTn4gLK3Wys27dOu65555i6+vVq0dSUlKZBCUi4nI2zIGTe6FGbeg6wbF6/4ksZq3cD8ATg1rg6V7p5mcWqfJK/a3y8vIiPT292Prdu3dTp45ukxQRKSYnDZZOty/3fBS8z3Y+/uf3O8gvNOjRtA69moWYFKCIayt1snP99dczbdo08vPzAbBYLCQmJvLII48wbNiwMg9QRKTK++VVOHMSajVxmhZi+e7j/LQzBXerhScH6VZzkfJS6mTn5ZdfJjMzk5CQELKzs+nRoweNGzfG39+ff/7zn+URo4hI1ZV6CFa/ZV++dhq4eQCQX2jj2Xk7ALgt9goah/ibFaGIyyv13ViBgYEsWrSIX375hS1btpCZmUn79u2Ji4srj/hERKq2RU9BQQ406AbN+jtWf7T6IHtSMgn29eSBa5qYGKCI6yt1spOTk4O3tzfdunWjW7du5RGTiIhrOLgKtn8FFiv0m+4YQPBEZi4vL9oNwKRrmxJYw8PMKEVcXqmTnaCgIK688kp69OhBr169iI2NxcfHpzxiExGpumyFMP9h+3L70RDe2lE0/YffyMgpoFW9AEZcGWlSgCLVR6n77Pz000/069ePNWvWcP3111OzZk26devG448/zqJFi8ojRhGRqmfTfyFpK3gFOg0guO7AKb7ceBiAZwe30qzmIhXAYhiGcakbFxQUsG7dOt5++20++ugjbDYbhYWFZRlfhUhPTycwMJC0tDQCAjQfjYhcpuxUeKMDnDkB/Z6HLuMAKCi0MeiNX/gtKYPhnSJ4fljrC+9HRC6opL/fpb6MBfYxdZYuXep45ObmMmjQIHr27Hmp8YqIuI5lL9oTndpNodOdjtXvJxzkt6QMgmp48HC/5iYGKFK9lDrZqVevHtnZ2fTs2ZOePXvyyCOP0Lp1a40PISICcHw3rH3bvtxvuuNW85T0HF75o1PyI/2aE+zraVaEItVOqfvs1KlThzNnzpCUlERSUhLJyclkZ2eXR2wiIlWLYcDCKWArgKb9ofHZITn++cNOMnMLaBMRxC0dI0wMUqT6KXWys3nzZpKSknj00UfJzc3lscceo3bt2lx11VU8/vjj5RGjiEjVsOdH2PsTWD2g79lBVhP2neTbzUexWOAfg1thVadkkQp1WR2UT548ydKlS/n222/55JNP1EFZRKqvglx4KxZO7YOr7oc+zwL2kZIHvLaCPSmZjOrSgGeHtDI5UBHXUW4dlL/66itHx+QdO3YQHBxMt27dePnll+nRo8dlBS0iUmWtesOe6PiGQPeHHKtn/bKfPSmZ1PL1ZHKfZiYGKFJ9lTrZuffee+nevTt33303PXr0ICYmpjziEhGpOk4fhOX/si/3/adjVvNDp87w6k97AHi0f3ONlCxiklInOykpKeURh4hI1bXgUSjIhiuuhpibADAMgye+2UZ2fiFdGgZzY4f6JgcpUn1d0jg7hYWFfPPNN+zcuROA6OhoBg8ejJubW5kGJyJS6e1aALt+AKs7DPiXY/6r7349yrLdx/F0t/LcDTEankPERKVOdvbu3cuAAQM4cuQIzZrZrz9Pnz6diIgIvv/+exo1alTmQYqIVEp5Z2D+H/1zYuMhxD5QYOqZPKbN3QHA/b0b07COn1kRigiXcOv5/fffT6NGjTh06BAbN25k48aNJCYmEhUVxf33318eMYqIVE6//BtSEyGgHnR/2LH6uR92cjIrj6ahftzdXf8BFDFbqVt2li1bxurVqwkODnasq1WrFs8//zxdu3Yt0+BERCqtE3th5Wv25X7TwcveerNq3wk+W2+f6HP60Bg83Uv9f0oRKWOl/hZ6eXmRkZFRbH1mZiaenhr+XESqAcOAHyZDYZ59lOQW1wOQk1/I419vA+DWLpF0aBB8ob2ISAUpdbIzaNAg7r77btasWYNhGBiGwerVq7n33nu5/vrryyNGEZHKZcc38PsScPOC/i86OiW/uWQv+09kEeLvpYk+RSqRUic7r7/+Oo0aNSI2NhZvb2+8vb3p2rUrjRs35rXXXiuPGEVEKo+cdFjwmH252wSoZe+Tszs5gxlL9wEwbXBLArw1po5IZVHqPjtBQUF8++237N2713HreYsWLWjcuHGZByciUuksfgYyjkLNKOg2EYCCQhsPfbGFAptBXItQ+rYMMzlIEfmzEic7NpuNl156ie+++468vDyuueYapk6dio+PT3nGJyJSeSSugXXv2ZevexU87P/+vffLfn49lIq/tzv/GNJKY+qIVDIlvoz1z3/+k8ceeww/Pz/q1avHa6+9Rnx8fHnGJiJSeRTkwdwHAAPajoSGPQHYm5LJy4t2A/DkwGjCAr3Ni1FEzqnEyc4HH3zAW2+9xcKFC/nmm2+YO3cuH330ETabrTzjExGpHFa+Csd3Qo3a0OcfABTaDB7+4lfyCmx0b1qHmzpqSgiRyqjEyU5iYiIDBgxwvI6Li8NisXD06NFyCUxEpNI4sQeWv2Rf7vc81LDfUj575X42Jqbi5+XO80M1JYRIZVXiZKegoABvb+fmWQ8PD/Lz88s8KBGRSsNms1++KhpTJ+ZGAPafyOKlhbsAeHxgC+oGqf+iSGVV4g7KhmEwZswYvLy8HOtycnK499578fX1daz76quvyjZCEREzbfoADq4Ejxow8N9gsWD74/JVboGNbo1rM7xThNlRisgFlDjZGT16dLF1t956a5kGIyJSqWQkwY9P2Zd7PQ41GwDwfsIB1h04ja+nG88P0+UrkcquxMnO7NmzyzMOEZHKxTBg3iTITYPwttD5XgAOnszihQW/ATBlQAvq16xhYpAiUhKaoU5E5Fy2fQm7vgerOwx+E9zcKSi0MfF/m8nJtxHbsBZ/uzLS7ChFpASU7IiI/FVmin2iT4DuD0NYKwDeXv674+6rl25qjdWqy1ciVYGSHRGRPzMMmDcRsk9DWAxcPQmAbUfSeOWPwQOfub6lLl+JVCGmJjvTp0+nU6dO+Pv7ExISwpAhQ9i1a5dTnZycHOLj46lVqxZ+fn4MGzaM5ORkpzqJiYkMHDiQGjVqEBISwkMPPURBQUFFHoqIuIrtX8Fv8/64fPUWuHmQk1/IxP9tpsBm0L9VGEPb1zM7ShEpBVOTnWXLlhEfH8/q1atZtGgR+fn59OnTh6ysLEediRMnMnfuXD7//HOWLVvG0aNHGTp0qKO8sLCQgQMHkpeXx6pVq3j//feZM2cOTz31lBmHJCJVWeZx+P6Py1dXT4bw1gC8uGAXe1Iyqe3nxT9v0N1XIlWNxTAMw+wgihw/fpyQkBCWLVtG9+7dSUtLo06dOnz88cfceKN9IK/ffvuNFi1akJCQQJcuXZg/fz6DBg3i6NGjhIaGAjBz5kweeeQRjh8/jqen50XfNz09ncDAQNLS0ggICCjXYxSRSuyz22DHtxAaA3f9DO6erNx7gpHvrgFg9phO9GoeYnKQIlKkpL/flarPTlpaGgDBwfah2Dds2EB+fj5xcXGOOs2bNycyMpKEhAQAEhISiImJcSQ6AH379iU9PZ3t27ef831yc3NJT093eohINbftK3uiY3WHIW+Buydp2flM/vxXAEZ2jlSiI1JFVZpkx2azMWHCBLp27UqrVvY7H5KSkvD09CQoKMipbmhoKElJSY46f050isqLys5l+vTpBAYGOh4RERr9VKRay0iC7x+0L1/9oOPy1dRvt3EsLYcratXg8YEtTAxQRC5HpUl24uPj2bZtG59++mm5v9eUKVNIS0tzPA4dOlTu7ykilZRhwLfjIfsUhLW299UBvt18hG82H8XNauGVW9pSw7PEY7CKSCVTKb6948ePZ968eSxfvpz69es71oeFhZGXl0dqaqpT605ycjJhYWGOOmvXrnXaX9HdWkV1/srLy8tpji8RqcbWz4K9i8DNC4a+A+6eHDyZxeNfbwNgfK/GtIusaXKQInI5TG3ZMQyD8ePH8/XXX/Pzzz8TFRXlVN6hQwc8PDxYvHixY92uXbtITEwkNjYWgNjYWLZu3UpKSoqjzqJFiwgICCA6OrpiDkREqqYTe+HHJ+zL1z4DIc3JK7Bx/yebyMwt4Morgrmvd2NzYxSRy2Zqy058fDwff/wx3377Lf7+/o4+NoGBgfj4+BAYGMjYsWOZNGkSwcHBBAQEcN999xEbG0uXLl0A6NOnD9HR0YwaNYoXX3yRpKQknnjiCeLj49V6IyLnV1gAX98N+WcgqgdceQ8ALy/axa+H0wj08eDV4W1xd6s0V/tF5BKZmuzMmDEDgJ49ezqtnz17NmPGjAHglVdewWq1MmzYMHJzc+nbty9vvfWWo66bmxvz5s1j3LhxxMbG4uvry+jRo5k2bVpFHYaIVEUrXoYjG8A7EIbMAKuV5buP8/ay3wF48cbW1A3yMTlIESkLlWqcHbNonB2RaubIBnj3WjAKYdh7EHMjxzNy6f/aCk5k5jKqSwOeHdLK7ChF5CKq5Dg7IiLlLi8Lvrrbnui0GgYxN2KzGUz6bDMnMnNpHuav28xFXIySHRGpXuY/Aif3gn9dGPAvAN795XdW7DmBt4eVN0a0w9vDzeQgRaQsKdkRkepj25ew6b+ABYa+DTWC2Zh4mhcX2CcgnnpdS5qE+psbo4iUOSU7IlI9nD4AcyfYl7tPhqjunMrKI/6jjRTYDAa2Dmd4J42mLuKKlOyIiOsrzIcvxkJuOkR0hh6PUmgzeODTTRxLy6FhHV9eGNZas5mLuCglOyLi+pY8B0fW228zH/YuuLnzxs97HP10ZozsgJ9XpRhQXkTKgZIdEXFtvy+FX16xL1/3OgRFsnz3cV5bvAeA526IoVmY+umIuDIlOyLiurJO2G8zx4AOY6DlEI6mZvPAp5swDBhxZSRD29e/2F5EpIpTsiMirslWaE90MpOhTgvoO528AhvxH2/k9Jl8WtULYOp1mj9PpDpQsiMirmn5v2DfYnD3gRtngWcNps/fyabEVAK83ZkxsoPG0xGpJpTsiIjr2fczLJ1uXx70CoRG89XGw8xeeQCAl29uS0RwDfPiE5EKpWRHRFxL2hH48k7AgPajoe0IthxO5dGvtgIwvldjro0ONTdGEalQSnZExHUU5sMXt8OZkxAWA/1f5HhGLvf8dwN5BTauaR7CpGubmh2liFQwJTsi4jp+ehoOrQGvQLj5A/Isnvz9ow2OgQNfGd4Wq1UDB4pUN0p2RMQ17PgOEv7PvjzkLQhuyLR521l34DT+Xu68c1tHArw9zI1RREyhZEdEqr7ju+Gbv9uXr7oPWgzik7WJfLg6EYsFXhvRlkZ1/MyNUURMo2RHRKq2nDT49G+QlwGRV8E1U9lw8BRPfbsNgAevbUrv5uqQLFKdKdkRkarLZrMPHHhyDwTUg5vf53B6Pvf8dyP5hQYDYsKI79XY7ChFxGRKdkSk6lo6HXYvADcvuOVDMtxrMnbOek5k5tIiPICXbmyjmcxFRMmOiFRRO76D5S/al697jYKwttz3ySZ2JWcQ4u/Fe6M74quZzEUEJTsiUhWl7IRvxtmXO4+DtiN4dt4Olu46jreHlXdHd6RukI+5MYpIpaFkR0SqluzTf3RIzoQrroY+z/L+qgO8n3AQgFdvaUvr+kHmxigilYqSHRGpOgrz4bPRcOp3CIyEm+awZO9pnpm7HYBH+jWnX6twk4MUkcpGyY6IVA2GAT9Mhv3LwMMXRnzMznRP7vt4EzYDbu5Yn3t7NDQ7ShGphJTsiEjVsHoGbJgDWODG9zjs1YjRs9aSmVtAl4bB/GNIjO68EpFzUrIjIpXfrgWw8DH7cp9/kBpxDWNmryMlI5dmof68Paojnu7650xEzk3/OohI5Za0Fb64AzCg/WhyOt7Lne+vZ29KJuGB3sy5oxOBPprzSkTOT8mOiFReGUnw8XDIz4Ko7hT2/xf3f7qZ9QdPE+Dtzvt3XEl4oG4xF5ELU7IjIpVTbiZ8MhzSD0OtJhg3fcDT3+/mxx3JeLpbeee2jjQN9Tc7ShGpApTsiEjlU5gPn4+Bo5vAJxj+9j/eXH2C/64+iMViH0unc8NaZkcpIlWEkh0RqVwMA+ZNgL2LwN0H/vYZH+x2418/7gZg6qBoBsRoLB0RKTklOyJSuSx9HjZ9CBYr3DSbr46H89S39kED77+mCWO6RpkcoIhUNUp2RKTy2DAHlj1vXx74MgsL2vHQF1sAuL3rFUyMa2JebCJSZSnZEZHKYdcCmDfRvtz9YVYGXc99H2+i0GZwY4f6PDkwWoMGisglUbIjIuZLXGPvkGzYoO2tbGg4jrs+WE9eoY1+LcN4fmgMVqsSHRG5NEp2RMRcx7bARzdBQTY0vpZtHZ7h9jnrOJNXyNVNavPaiLa4u+mfKhG5dPoXRETMc2Iv/PcGyE2DyFi2X/1/jJy1kfScAjo2qMnbozrg5e5mdpQiUsUp2RERc6Qegg8Gw5kTENaa33q/y8j3t5CWnU/7yCBm396JGp7uZkcpIi5AyY6IVLzMFPjvEPvoyLWbsrvPB4z4YAepZ/JpGxHEnDuuxN9b812JSNlQsiMiFSv7NPx3KJzcC4GR7Ov3IcM/2svpM/m0qR/IB2OvJECJjoiUISU7IlJxctLgwxsheSv4hrB/4Efc/Ekip7LyiKkXyAdjOyvREZEyp2RHRCpGTjp8OAyOrAefmuzp9yFDP03iZFYereoF8OHYzgT6KNERkbKn3n8iUv5yM+yJzuF14B3Ejms/5JYvU8nIKaBN/UDev+NKAmso0RGR8qFkR0TKV26G/dLV4bXgHcSW3h8w/JtMzuQV0umKmswa00mdkUWkXCnZEZHyk5tpHzDw0GrwDmR9j9mM/O4MuQU2rm5Sm7dHddDt5SJS7kzts7N8+XKuu+466tati8Vi4ZtvvnEqNwyDp556ivDwcHx8fIiLi2PPnj1OdU6dOsXIkSMJCAggKCiIsWPHkpmZWYFHISLnlJsBH98MiQngFUhC13cZMS+H3AIbcS1CeOe2jkp0RKRCmJrsZGVl0aZNG958881zlr/44ou8/vrrzJw5kzVr1uDr60vfvn3Jyclx1Bk5ciTbt29n0aJFzJs3j+XLl3P33XdX1CGIyLlkn4YPhsDBleAVwKKOMxj5Qz75hQYDW4cz49YOeHtoZGQRqRgWwzAMs4MAsFgsfP311wwZMgSwt+rUrVuXBx98kMmTJwOQlpZGaGgoc+bMYfjw4ezcuZPo6GjWrVtHx44dAViwYAEDBgzg8OHD1K1bt0TvnZ6eTmBgIGlpaQQEBJTL8YlUG5nH7VNAJG/F8KnJ5y1e5+FV9sTm5o71mT60NW6a1FNEykBJf78r7a3n+/fvJykpibi4OMe6wMBAOnfuTEJCAgAJCQkEBQU5Eh2AuLg4rFYra9asOe++c3NzSU9Pd3qISBlIPwpzBtgTHd8QZlxxNtGJ79WIF4Yp0RGRildpk52kpCQAQkNDndaHhoY6ypKSkggJCXEqd3d3Jzg42FHnXKZPn05gYKDjERERUcbRi1RDpw/ArH5wYjdGQD2m1fkXL26yJzpTr4vmob7NsViU6IhIxau0yU55mjJlCmlpaY7HoUOHzA5JpGo7vhtm9YfUg9iCophQYzqzf3PHw83C6yPacXvXKLMjFJFqrNLeChEWFgZAcnIy4eHhjvXJycm0bdvWUSclJcVpu4KCAk6dOuXY/ly8vLzw8vIq+6BFqqNDa+13XWWfJj+4KbflP0bCAXd8Pd14e1RHujWpbXaEIlLNVdqWnaioKMLCwli8eLFjXXp6OmvWrCE2NhaA2NhYUlNT2bBhg6POzz//jM1mo3PnzhUes0i1s2sBvH89ZJ8mq05b+qVNIeG4JyH+Xnx6d6wSHRGpFExt2cnMzGTv3r2O1/v372fz5s0EBwcTGRnJhAkT+Mc//kGTJk2IioriySefpG7duo47tlq0aEG/fv246667mDlzJvn5+YwfP57hw4eX+E4sEblEG/8Lcx8Ao5CUsB70PXI7p/M9aR7mz6wxnagb5GN2hCIigMnJzvr16+nVq5fj9aRJkwAYPXo0c+bM4eGHHyYrK4u7776b1NRUunXrxoIFC/D29nZs89FHHzF+/HiuueYarFYrw4YN4/XXX6/wYxGpNgwDVvwLfv4HAL+FXcegAzdRgDu9mtXhjb+1x8+r0l4hF5FqqNKMs2MmjbMjUkKFBTD/YVj/HgA/1xnFHYf6ARbGXHUFTwxsgbtbpb06LiIupqS/3/rvl4iUTE46fHE77P0JAwuz/O/h2UPdsVrgqUHRjNEdVyJSSSnZEZGLO30QPr4Fju/E5ubNo5b7+ex4W/y93Hl9RDt6NQ+5+D5EREyiZEdELuzQOvh0BGQdJ9urDn/LnMimwitoHOLHf0Z1oGEdP7MjFBG5ICU7InJ+W7+Ab/4Ohbkc82nCDafvJ4la9IkO5d+3tFVHZBGpEvQvlYgUZyuEJc/Z77oC1nl1YfTpu8m2eDMprinjezXGqjmuRKSKULIjIs6yT8OXd8HeRQB8YLmOp9NuwdfLk3eHt+WaFqEX2YGISOWiZEdEzkreAZ/+DU7vJ9/qxeScsXxr60Z0eABvjmxPVG1fsyMUESk1JTsiYrf9a/gmHvKzSHEL4fYzE9huXMGtXSJ5YmA03h5uZkcoInJJlOyIVHeF+bB4Gqyyjzy+hhjuzYon3yuYN4bGcF0bTb0iIlWbkh2R6iz1EHxxBxxeC8DMgkG8VHALzcJr6rKViLgMJTsi1dWuBfDNvZB9mkx8mZx3FwtsVzKqSwMeH9hCl61ExGUo2RGpbgrzYfEzsOoNALbYGhKffx/ZvpHMvrG1RkMWEZejZEekOjm1H766Cw6vA2BWQT+eLxhBz+j6TB8aQy0/L5MDFBEpe0p2RKoDw4DNH2HMfwRLXibpRg0eyr+HFe5deHZYNDd3jMBi0SCBIuKalOyIuLqskzD3fvhtHhZgja05D+aPIzSyKfNvbkODWuqELCKuTcmOiCvb8xPGt3/HkplMnuHGvwtu4r/WwUwe1ILbYq/ATVM+iEg1oGRHxBXlpMOip2DDbCzAHls9JuTHE9SoA/NvaE1krRpmRygiUmGU7Ii4mt0/Ypv7ANaMowDMLujLm26jmHxDG27ppL45IlL9KNkRcRVnTmEseATLls+wAgdsoTxacBcBzXsxd3BLwgN9zI5QRMQUSnZEXMH2byic9yBu2ScoNCzMKuzPZ/6jeWxwe42bIyLVnpIdkars1H4Kvn8Y930/4gbsttXjMdu9XNWjH3N7NtIoyCIiKNkRqZrycyj85TWMFS/jbssl33BjRuF1bI66i5cGt9OcViIif6JkR6Sq2fsT2d9OwifjIAArC1vyH/9x3DqoD/e1CFEHZBGRv1CyI1JVnD5A+ndTCNj/Az5AshHEvy1jaN53NO90uQJPd6vZEYqIVEpKdkQqu5w0Un98Hr9N7xBg5FNgWPnA1o/jHScx5do2BNXwNDtCEZFKTcmOSGVVmE/qL//BffkLBBWmAfBLYUuWNZzI364boH45IiIlpGRHpLIxDE5v/o7CBU9QOzcRsI+A/H3YOK4dPIrH6wWZG5+ISBWjZEekEkn59UeyFz5DgzPbADhp+PNN0GjaDnmACVEaL0dE5FIo2RGpBA5t/pncH5+h8ZnNAGQbniz0G0LogMe4I/oK3WElInIZlOyImMQwDHZuWIpt8XO0yl4LQK7hzlL/QYQMeIwh0c1MjlBExDUo2RGpYHn5haxd8i2+616nXf4mAPINNxIC+1FnwBP0bR5tcoQiIq5FyY5IBTmdmUPC/A+J3DGTbsYeAAoMK7/W7EPwgCfo3jTG5AhFRFyTkh2RcmQYBpv2J7N78Rw6HP4vAyyHAcjBk93hg4kY9Agd6jUxOUoREdemZEekHKSdyWfB6k0UrHmXPjnzaW9JBwtkWWpwuPFIogZOpnVQmNlhiohUC0p2RMqIzWaw7sApElb8SMN9H3KDJQFPSyFYINW9DpltbqdeXDzNfILMDlVEpFpRsiNymfamZDBv/W5yNn5Ov7yFTLD+Dn9MU5US1BbfHvcR1HowQW4e5gYqIlJNKdkRuQQpGTnM3XyUHet+ptOpudzlloCvJResUGDxIK3RdQT3up+Qeu3MDlVEpNpTsiNSQsfSslmwLYnVv24n4ugPDLMuZ6z1kONblOkXhVfnMXi0G0ktvzrmBisiIg5KdkQuIPHkGeZvO8ayLfuol7SIwdaVjLbuwOpuAFBg9aKw+WC8Ot+OX2QsaKRjEZFKR8mOyJ/kFdhYf+AUS3alkPDbIeqdXMX1bquYbd2El0e+o15u3SvxancL7q2G4e5T08SIRUTkYpTsSLV3LC2bZbuOs2RXClv3HKBLwTr6ua1jknUrPp55jnoFtZrh3uZmiLkJr5oNTIxYRERKQ8mOVDspGTkk7DvJ6t9Psvr3U+Sf3E8v62ZGWdfRxboTd0+bo64tsAHWltdDzM24h8XoMpWISBWkZEdcmmEYHEnNZmNiKmv3nyRh30kOHz9NrHUHPay/cqd1C428jjlvExKNpcX10GIQ1tBWSnBERKo4JTviUrLzCtl6JI2NiafZlHiajYmpnMo4Q4xlP52tO3nKup3OXjvxtpztf2NY3bFEdIamfaH5ICy1Gpl4BCIiUtaU7EiVlZGTz85jGew4msaOY+lsP5rOrqQMsOXT2vI7na2/8TfrDjp47cbPkuO8cUB9aBIHjeOwRHUH70BzDkJERMqdyyQ7b775Ji+99BJJSUm0adOGN954gyuvvNLssKQM5BXYOHAyi30pmexJyWTnsXR2HEvn4MkzgEFdTtLGuo/rrXtp676P1tb9+JDrvBPvIGjQFa7oCo2ugTrNdHlKRKSacIlk53//+x+TJk1i5syZdO7cmVdffZW+ffuya9cuQkJCzA5PSsBmMziemcuhU2f4/UQW+45nsi8lk33Hs0g8dYZCmwEYhHOKFtaDDLIk0sZjH+3d9lGb1OI7rFELGlwFDbrBFd0gJBqs1oo+LBERqQQshmEYZgdxuTp37kynTp34v//7PwBsNhsRERHcd999PProoxfdPj09ncDAQNLS0ggICCizuI4fPUBBfi7n+gsbtqKVxh8PHPUMx2vDsdKxC6ed/WW7P5c5tjv7PsXrOVacO0bjr/GdrWQUVfjTvoxi2xWP2WbA6axckjNySE7LISU9h+SMHFLSc8kvNLD8aTsPCmhkPUoLSyKt3A7R3JqIv5FZPFCrO4S2hHodoF5HqN8RajVRciMi4uJK+vtd5Vt28vLy2LBhA1OmTHGss1qtxMXFkZCQcM5tcnNzyc09e5kjPT29XGLLfncAkbYj5bJvl+POxT+NBvbEpnZTCG0F4W3siU14G/DwqYAgRUSkKqryyc6JEycoLCwkNDTUaX1oaCi//fbbObeZPn06zzzzTLnHlm/xJNvwdLw2OH8fkb+WGVgoWvXXRhcDi6P2ObdzLP+V5Zz1ite90D4vcAyW82/3Z1aLBavVirvVgpvVgpvb2WX7kVnOhlEzyp7YhLWyP9dpBu5e541BRETkr6p8snMppkyZwqRJkxyv09PTiYiIKPP3afTk5jLfp4iIiJROlU92ateujZubG8nJyU7rk5OTCQsLO+c2Xl5eeHmpdUBERKQ6qPI9OD09PenQoQOLFy92rLPZbCxevJjY2FgTIxMREZHKoMq37ABMmjSJ0aNH07FjR6688kpeffVVsrKyuP32280OTUREREzmEsnOLbfcwvHjx3nqqadISkqibdu2LFiwoFinZREREal+XGKcnctVXuPsiIiISPkp6e93le+zIyIiInIhSnZERETEpSnZEREREZemZEdERERcmpIdERERcWlKdkRERMSlKdkRERERl6ZkR0RERFyakh0RERFxaS4xXcTlKhpEOj093eRIREREpKSKfrcvNhmEkh0gIyMDgIiICJMjERERkdLKyMggMDDwvOWaGwuw2WwcPXoUf39/LBZLme03PT2diIgIDh065LJzbrn6Mer4qj5XP0YdX9Xn6sdYnsdnGAYZGRnUrVsXq/X8PXPUsgNYrVbq169fbvsPCAhwyQ/wn7n6Mer4qj5XP0YdX9Xn6sdYXsd3oRadIuqgLCIiIi5NyY6IiIi4NCU75cjLy4upU6fi5eVldijlxtWPUcdX9bn6Mer4qj5XP8bKcHzqoCwiIiIuTS07IiIi4tKU7IiIiIhLU7IjIiIiLk3JjoiIiLg0JTuX6c033+SKK67A29ubzp07s3bt2gvW//zzz2nevDne3t7ExMTwww8/VFCkpTd9+nQ6deqEv78/ISEhDBkyhF27dl1wmzlz5mCxWJwe3t7eFRRx6Tz99NPFYm3evPkFt6lK5++KK64odnwWi4X4+Phz1q8K52758uVcd9111K1bF4vFwjfffONUbhgGTz31FOHh4fj4+BAXF8eePXsuut/Sfo/Ly4WOLz8/n0ceeYSYmBh8fX2pW7cut912G0ePHr3gPi/lc16eLnYOx4wZUyzefv36XXS/VeEcAuf8TlosFl566aXz7rMyncOS/C7k5OQQHx9PrVq18PPzY9iwYSQnJ19wv5f63S0pJTuX4X//+x+TJk1i6tSpbNy4kTZt2tC3b19SUlLOWX/VqlWMGDGCsWPHsmnTJoYMGcKQIUPYtm1bBUdeMsuWLSM+Pp7Vq1ezaNEi8vPz6dOnD1lZWRfcLiAggGPHjjkeBw8erKCIS69ly5ZOsf7yyy/nrVvVzt+6deucjm3RokUA3HTTTefdprKfu6ysLNq0acObb755zvIXX3yR119/nZkzZ7JmzRp8fX3p27cvOTk5591nab/H5elCx3fmzBk2btzIk08+ycaNG/nqq6/YtWsX119//UX3W5rPeXm72DkE6Nevn1O8n3zyyQX3WVXOIeB0XMeOHWPWrFlYLBaGDRt2wf1WlnNYkt+FiRMnMnfuXD7//HOWLVvG0aNHGTp06AX3eynf3VIx5JJdeeWVRnx8vON1YWGhUbduXWP69OnnrH/zzTcbAwcOdFrXuXNn45577inXOMtKSkqKARjLli07b53Zs2cbgYGBFRfUZZg6darRpk2bEtev6ufvgQceMBo1amTYbLZzllelc2cYhgEYX3/9teO1zWYzwsLCjJdeesmxLjU11fDy8jI++eST8+6ntN/jivLX4zuXtWvXGoBx8ODB89Yp7ee8Ip3rGEePHm0MHjy4VPupyudw8ODBRu/evS9YpzKfw7/+LqSmphoeHh7G559/7qizc+dOAzASEhLOuY9L/e6Whlp2LlFeXh4bNmwgLi7Osc5qtRIXF0dCQsI5t0lISHCqD9C3b9/z1q9s0tLSAAgODr5gvczMTBo0aEBERASDBw9m+/btFRHeJdmzZw9169alYcOGjBw5ksTExPPWrcrnLy8vjw8//JA77rjjgpPdVqVz91f79+8nKSnJ6RwFBgbSuXPn856jS/keVyZpaWlYLBaCgoIuWK80n/PKYOnSpYSEhNCsWTPGjRvHyZMnz1u3Kp/D5ORkvv/+e8aOHXvRupX1HP71d2HDhg3k5+c7nY/mzZsTGRl53vNxKd/d0lKyc4lOnDhBYWEhoaGhTutDQ0NJSko65zZJSUmlql+Z2Gw2JkyYQNeuXWnVqtV56zVr1oxZs2bx7bff8uGHH2Kz2bjqqqs4fPhwBUZbMp07d2bOnDksWLCAGTNmsH//fq6++moyMjLOWb8qn79vvvmG1NRUxowZc946VencnUvReSjNObqU73FlkZOTwyOPPMKIESMuOLliaT/nZuvXrx8ffPABixcv5oUXXmDZsmX079+fwsLCc9avyufw/fffx9/f/6KXeCrrOTzX70JSUhKenp7FEvCL/TYW1SnpNqWlWc+lROLj49m2bdtFrxPHxsYSGxvreH3VVVfRokUL3n77bZ599tnyDrNU+vfv71hu3bo1nTt3pkGDBnz22Wcl+p9WVfLee+/Rv39/6tate946VencVXf5+fncfPPNGIbBjBkzLli3qn3Ohw8f7liOiYmhdevWNGrUiKVLl3LNNdeYGFnZmzVrFiNHjrzojQCV9RyW9HehMlDLziWqXbs2bm5uxXqYJycnExYWds5twsLCSlW/shg/fjzz5s1jyZIl1K9fv1Tbenh40K5dO/bu3VtO0ZWdoKAgmjZtet5Yq+r5O3jwID/99BN33nlnqbarSucOcJyH0pyjS/kem60o0Tl48CCLFi26YKvOuVzsc17ZNGzYkNq1a5833qp4DgFWrFjBrl27Sv29hMpxDs/3uxAWFkZeXh6pqalO9S/221hUp6TblJaSnUvk6elJhw4dWLx4sWOdzWZj8eLFTv87/rPY2Fin+gCLFi06b32zGYbB+PHj+frrr/n555+Jiooq9T4KCwvZunUr4eHh5RBh2crMzGTfvn3njbWqnb8is2fPJiQkhIEDB5Zqu6p07gCioqIICwtzOkfp6emsWbPmvOfoUr7HZipKdPbs2cNPP/1ErVq1Sr2Pi33OK5vDhw9z8uTJ88Zb1c5hkffee48OHTrQpk2bUm9r5jm82O9Chw4d8PDwcDofu3btIjEx8bzn41K+u5cSuFyiTz/91PDy8jLmzJlj7Nixw7j77ruNoKAgIykpyTAMwxg1apTx6KOPOuqvXLnScHd3N/71r38ZO3fuNKZOnWp4eHgYW7duNesQLmjcuHFGYGCgsXTpUuPYsWOOx5kzZxx1/nqMzzzzjLFw4UJj3759xoYNG4zhw4cb3t7exvbt2804hAt68MEHjaVLlxr79+83Vq5cacTFxRm1a9c2UlJSDMOo+ufPMOx3pURGRhqPPPJIsbKqeO4yMjKMTZs2GZs2bTIA49///rexadMmx91Izz//vBEUFGR8++23xpYtW4zBgwcbUVFRRnZ2tmMfvXv3Nt544w3H64t9jyvL8eXl5RnXX3+9Ub9+fWPz5s1O38nc3NzzHt/FPucV7ULHmJGRYUyePNlISEgw9u/fb/z0009G+/btjSZNmhg5OTmOfVTVc1gkLS3NqFGjhjFjxoxz7qMyn8OS/C7ce++9RmRkpPHzzz8b69evN2JjY43Y2Fin/TRr1sz46quvHK9L8t29HEp2LtMbb7xhREZGGp6ensaVV15prF692lHWo0cPY/To0U71P/vsM6Np06aGp6en0bJlS+P777+v4IhLDjjnY/bs2Y46fz3GCRMmOP4eoaGhxoABA4yNGzdWfPAlcMsttxjh4eGGp6enUa9ePeOWW24x9u7d6yiv6ufPMAxj4cKFBmDs2rWrWFlVPHdLliw552ey6DhsNpvx5JNPGqGhoYaXl5dxzTXXFDv2Bg0aGFOnTnVad6HvcUW60PHt37//vN/JJUuWOPbx1+O72Oe8ol3oGM+cOWP06dPHqFOnjuHh4WE0aNDAuOuuu4olLVX1HBZ5++23DR8fHyM1NfWc+6jM57AkvwvZ2dnG3//+d6NmzZpGjRo1jBtuuME4duxYsf38eZuSfHcvh+WPNxURERFxSeqzIyIiIi5NyY6IiIi4NCU7IiIi4tKU7IiIiIhLU7IjIiIiLk3JjoiIiLg0JTsiIiLi0pTsiIiIiEtTsiMiFW7p0qVYLJZikwWW1pgxYxgyZEiZxGSGnj17MmHCBLPDEHF5SnZE5JLNnDkTf39/CgoKHOsyMzPx8PCgZ8+eTnWLEpx9+/Zx1VVXcezYMQIDA8s9xnfeeYc2bdrg5+dHUFAQ7dq1Y/r06eX+viJSebibHYCIVF29evUiMzOT9evX06VLFwBWrFhBWFgYa9asIScnB29vbwCWLFlCZGQkjRo1AiAsLKzc45s1axYTJkzg9ddfp0ePHuTm5rJlyxa2bdtW7u8tIpWHWnZE5JI1a9aM8PBwli5d6li3dOlSBg8eTFRUFKtXr3Za36tXL8fyny9jzZkzh6CgIBYuXEiLFi3w8/OjX79+HDt2zLF9YWEhkyZNIigoiFq1avHwww9zsan9vvvuO26++WbGjh1L48aNadmyJSNGjOCf//yno07RpbBnnnmGOnXqEBAQwL333kteXp6jjs1mY/r06URFReHj40ObNm344osvnN5r27Zt9O/fHz8/P0JDQxk1ahQnTpxwlGdlZXHbbbfh5+dHeHg4L7/8csn/0CJyWZTsiMhl6dWrF0uWLHG8XrJkCT179qRHjx6O9dnZ2axZs8aR7JzLmTNn+Ne//sV///tfli9fTmJiIpMnT3aUv/zyy8yZM4dZs2bxyy+/cOrUKb7++usLxhYWFsbq1as5ePDgBestXryYnTt3snTpUj755BO++uornnnmGUf59OnT+eCDD5g5cybbt29n4sSJ3HrrrSxbtgyA1NRUevfuTbt27Vi/fj0LFiwgOTmZm2++2bGPhx56iGXLlvHtt9/y448/snTpUjZu3HjBuESkjJTZ/OkiUi298847hq+vr5Gfn2+kp6cb7u7uRkpKivHxxx8b3bt3NwzDMBYvXmwAxsGDBw3DMIwlS5YYgHH69GnDMAxj9uzZBmDs3bvXsd8333zTCA0NdbwODw83XnzxRcfr/Px8o379+sbgwYPPG9vRo0eNLl26GIDRtGlTY/To0cb//vc/o7Cw0FFn9OjRRnBwsJGVleVYN2PGDMPPz88oLCw0cnJyjBo1ahirVq1y2vfYsWONESNGGIZhGM8++6zRp08fp/JDhw4ZgLFr1y4jIyPD8PT0ND777DNH+cmTJw0fHx/jgQceuNCfV0TKgPrsiMhl6dmzJ1lZWaxbt47Tp0/TtGlT6tSpQ48ePbj99tvJyclh6dKlNGzYkMjIyPPup0aNGo7+PADh4eGkpKQAkJaWxrFjx+jcubOj3N3dnY4dO17wUlZ4eDgJCQls27aN5cuXs2rVKkaPHs27777LggULsFrtjdtt2rShRo0aju1iY2PJzMzk0KFDZGZmcubMGa699lqnfefl5dGuXTsAfv31V5YsWYKfn1+xGPbt20d2djZ5eXlO8QcHB9OsWbPzxi4iZUfJjohclsaNG1O/fn2WLFnC6dOn6dGjBwB169YlIiKCVatWsWTJEnr37n3B/Xh4eDi9tlgsF+2TU1KtWrWiVatW/P3vf+fee+/l6quvZtmyZRe8rFYkMzMTgO+//5569eo5lXl5eTnqXHfddbzwwgvFtg8PD2fv3r1lcBQicqnUZ0dELluvXr1YunQpS5cudbrlvHv37syfP5+1a9eWKLE4n8DAQMLDw1mzZo1jXUFBARs2bCj1vqKjowF7h+Eiv/76K9nZ2Y7Xq1evxs/Pj4iICKKjo/Hy8iIxMZHGjRs7PSIiIgBo374927dv54orrihWx9fXl0aNGuHh4eEU/+nTp9m9e3ep4xeR0lPLjohctl69ehEfH09+fr6jZQegR48ejB8/nry8vMtKdgAeeOABnn/+eZo0aULz5s3597//fdFBCceNG0fdunXp3bs39evX59ixY/zjH/+gTp06xMbGOurl5eUxduxYnnjiCQ4cOMDUqVMZP348VqsVf39/Jk+ezMSJE7HZbHTr1o20tDRWrlxJQEAAo0ePJj4+nnfeeYcRI0bw8MMPExwczN69e/n0009599138fPzY+zYsTz00EPUqlWLkJAQHn/8ccdlNBEpX0p2ROSy9erVi+zsbJo3b05oaKhjfY8ePcjIyHDcon45HnzwQY4dO8bo0aOxWq3ccccd3HDDDaSlpZ13m7i4OGbNmsWMGTM4efIktWvXJjY2lsWLF1OrVi1HvWuuuYYmTZrQvXt3cnNzGTFiBE8//bSj/Nlnn6VOnTpMnz6d33//naCgINq3b89jjz0G2C/ZrVy5kkceeYQ+ffqQm5tLgwYN6NevnyOheemllxyXu/z9/XnwwQcvGLuIlB2LUVYXxUVEqqAxY8aQmprKN998Y3YoIlJO1IYqIiIiLk3JjoiIiLg0XcYSERERl6aWHREREXFpSnZERETEpSnZEREREZemZEdERERcmpIdERERcWlKdkRERMSlKdkRERERl6ZkR0RERFza/wOOh7fxmnZngQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def g_offshore(v):\n", " return np.piecewise(\n", " v, [v <= 3.5, v > 15], [0, 600, lambda v: 0.18007 * v**3 - 7.72049]\n", " )\n", "\n", "\n", "def g_onshore(v):\n", " return np.piecewise(\n", " v, [v <= 3, v > 14], [0, 450, lambda v: 0.16563 * v**3 - 4.4718]\n", " )\n", "\n", "\n", "# assign wind generators\n", "g_func = {64: g_onshore, 65: g_offshore}\n", "\n", "# show plots\n", "v = np.linspace(0, 20, 1000)\n", "for i in g_func.keys():\n", " plt.plot(v, g_func[i](v), label=i)\n", "\n", "plt.legend()\n", "plt.xlabel(\"Wind Speed\")\n", "plt.ylabel(\"Power Output\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting uc_windcurtailment.mod\n" ] } ], "source": [ "%%writefile uc_windcurtailment.mod\n", "\n", "set T;\n", "set V;\n", "set E within {V, V};\n", "set W;\n", "set NW;\n", "\n", "param d{V};\n", "param p_min{V};\n", "param p_max{V};\n", "param c_var{V};\n", "param is_generator{V};\n", "param energy_type{V} symbolic;\n", "param c_fixed{V};\n", "\n", "param b{E};\n", "param f_max{E};\n", "\n", "param vmax{W};\n", "\n", "param samples{W, T};\n", "param g{W, T};\n", "\n", "param M;\n", "\n", "# Declare decision variables\n", "# Binary variable indicating whether a generator is on or off\n", "var x{V} binary;\n", "# Power generation of each generator\n", "var p{V, T} >= 0;\n", "# Voltage angle of each node\n", "var theta{V, T};\n", "# Power flow on each edge\n", "var f{(i,j) in E, T} <= f_max[i,j], >= -f_max[i,j];\n", "# Binary variable indicating whether a wind turbine is curtailed (0) or not (1)\n", "var y{W, T} binary;\n", "\n", "# Objective function, with the goal of minimizing the total cost of generation\n", "var term1 =\n", " sum{i in V: energy_type[i] in {\"coal\", \"gas\"}} M * x[i];\n", "\n", "var term2{t in T} =\n", " sum{i in V: is_generator[i] == 1} c_var[i] * p[i, t];\n", "\n", "minimize objective: \n", " term1 + sum{t in T} term2[t] / card(T);\n", "\n", "# Declare constraints\n", "# Wind power output is equal to the wind speed times the wind turbine capacity, if the wind turbine is not curtailed\n", "s.t. wind_speed_to_power{i in W, t in T}:\n", " p[i, t] == (1 - y[i, t]) * g[i, t];\n", " \n", "s.t. wind_curtailment{i in W, t in T}:\n", " samples[i, t] <= vmax[i] + y[i, t] * M;\n", "\n", "s.t. generation_upper_bound{i in NW, t in T}:\n", " p[i, t] <= x[i] * p_max[i];\n", " \n", "s.t. generation_lower_bound{i in NW, t in T}:\n", " p[i, t] >= x[i] * p_min[i];\n", "\n", "# Net flow in every node must be equal to the power generation minus the demand in that node\n", "var incoming_flow{i in V, t in T} = sum{j in V: (j, i) in E} f[j, i, t];\n", "var outgoing_flow{i in V, t in T} = sum{j in V: (i, j) in E} f[i, j, t];\n", "s.t. flow_conservation{i in V, t in T}:\n", " incoming_flow[i, t] - outgoing_flow[i, t] == p[i, t] - d[i];\n", "\n", "s.t. susceptance{(i,j) in E, t in T}:\n", " f[i, j, t] == b[i, j] * (theta[i, t] - theta[j, t]);" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "id": "lnGh_DTyqZbY" }, "outputs": [], "source": [ "def UC_windcurtailment(network, samples, g, vmax):\n", " df_nodes = network[\"nodes\"]\n", " df_edges = network[\"edges\"]\n", "\n", " # Define a model\n", " model = AMPL()\n", " model.read(\"uc_windcurtailment.mod\")\n", "\n", " model.set_data(df_nodes, \"V\")\n", " model.set_data(df_edges, \"E\")\n", "\n", " model.set[\"T\"] = range(samples.shape[1])\n", " model.set[\"W\"] = df_nodes.index[df_nodes[\"energy_type\"] == \"wind\"]\n", " model.set[\"NW\"] = df_nodes.index[df_nodes[\"energy_type\"] != \"wind\"]\n", "\n", " model.param[\"M\"] = 1000\n", " model.param[\"samples\"] = samples\n", " model.param[\"g\"] = g\n", " model.param[\"vmax\"] = vmax\n", "\n", " model.option[\"solver\"] = SOLVER\n", "\n", " return model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We sample from two different continuous Weibull distributions with different scale and shape parameters and then solve the problem." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "id": "RfIVAHLtJxs3" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective value: 29662.84\n", "Number of scenarios wind tubine 64 is active: 18 out of 100\n", "Number of scenarios wind tubine 65 is active: 16 out of 100\n" ] } ], "source": [ "seed = 0\n", "rng = np.random.default_rng(seed)\n", "\n", "scale64 = 15\n", "shape64 = 2.6\n", "scale65 = 18\n", "shape65 = 3\n", "\n", "nsamples = 100\n", "samples = [\n", " {64: scale64 * rng.weibull(shape64), 65: scale65 * rng.weibull(shape65)}\n", " for i in range(nsamples)\n", "]\n", "\n", "g = [{k: g_func[k](s[k]) for k in s} for s in samples]\n", "\n", "df_samples = pd.DataFrame(samples).T\n", "df_g = pd.DataFrame(g).T\n", "\n", "# Maximum wind speeds for each wind generator\n", "vmax = {64: 20.0, 65: 22.0}\n", "\n", "model = UC_windcurtailment(network, df_samples, df_g, vmax)\n", "model.get_output(\"solve;\")\n", "\n", "print(f\"Objective value: {model.obj['objective'].value():.2f}\")\n", "\n", "y = model.var[\"y\"].to_pandas()\n", "y.index.names = [\"W\", \"T\"]\n", "y = y.rename(columns={y.columns[0]: \"y\"})\n", "\n", "W = [64, 65]\n", "for w in W:\n", " q = f\"W == {w}\"\n", " v = y.query(q)[\"y\"].sum()\n", " print(f\"Number of scenarios wind tubine {w} is active: {v:2.0f} out of {nsamples}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Assessing the quality of the solution\n", "\n", "How can we assess the quality of the solution? For example, we could investigate how it performs compared to some baseline solution: we could have assumed 'average' values for the wind speed and solve the problem like a single-stage problem where all the data is known. That would give us some initial here-and-now commitment decisions, for which we could check in how many cases is it not even possible to schedule the production $p_i$ to meet the demand.\n", "\n", "To do so, we need a variant of our model building function in which we can fix the initial decisions $x_i$ and seek the best possible decisions for the other variables given an observed wind speed, for each of the $T$ scenarios, and see how often is the problem infeasible. The following is a modified version of our function." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def UC_windcurtailment_fixed_x(network, samples, g, vmax, fixed_x=None):\n", " model = UC_windcurtailment(network, samples, g, vmax)\n", "\n", " df_nodes = network[\"nodes\"]\n", " G = df_nodes.query('is_generator and energy_type != \"wind\"').index\n", "\n", " if fixed_x is not None:\n", " for i in G:\n", " model.var[\"x\"][i].fix(fixed_x[i])\n", "\n", " return model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we construct our average sample, solve for it, and then optimize for the second stage decisions for each of the simulated scenarios separately." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Nominal objective = 15101.25\n" ] } ], "source": [ "T = len(samples)\n", "\n", "mean_sample = [\n", " {i: np.array([samples[t][i] for t in range(T)]).mean() for i in samples[0].keys()}\n", "]\n", "g = [{k: g_func[k](s[k]) for k in s} for s in mean_sample]\n", "\n", "df_samples = pd.DataFrame(mean_sample).T\n", "df_g = pd.DataFrame(g).T\n", "\n", "m_nominal = UC_windcurtailment_fixed_x(network, df_samples, df_g, vmax)\n", "m_nominal.get_output(\"solve;\")\n", "\n", "print(f\"Nominal objective = {m_nominal.obj['objective'].value():0.2f}\")\n", "fixed_x = m_nominal.var[\"x\"].to_dict()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, for each of our scenarios we will solve the problem by fixing the initial decisions to be the one for the 'average' scenario and see in how many of them we cannot obtain a feasible second-stage decision, which we extract using the try-except mechanism of Python catching an error in case we are trying to call a solution that does not exist." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The second stage problem is infeasible in 65 out of 100 instances\n" ] } ], "source": [ "# Variable to count the number of infeasible second stage problems\n", "n_infeasible = 0\n", "\n", "for t in range(T):\n", " sample = [samples[t]]\n", " g = [{k: g_func[k](s[k]) for k in s} for s in sample]\n", "\n", " df_samples = pd.DataFrame(sample).T\n", " df_g = pd.DataFrame(g).T\n", "\n", " m_single = UC_windcurtailment_fixed_x(network, df_samples, df_g, vmax, fixed_x)\n", " m_single.get_output(\"solve;\")\n", " solve_result = m_single.get_value(\"solve_result\")\n", " n_infeasible += 1 if solve_result == \"infeasible\" else 0\n", "\n", "print(f\"The second stage problem is infeasible in {n_infeasible} out of {T} instances\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As it turns out, the 'nominal' commitment solution would lead to infeasibility, thus a blackout, in 65% of the cases. Since a blackout is the very last thing a network operator would like to happen, this clearly demonstrates the value of the stochastic solution. It is important to note that a `very robust` solution that avoids blackouts could be explained by simply assuming very poor performance of the wind turbines. Such a solution, however, could be unnecessarily conservative and expensive, and it is therefore better to prepare for a set of `realistic` scenarios." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" }, "vscode": { "interpreter": { "hash": "96c80cd78af1433053b8f057723d5c1af154f81d1718cfecf4b67cabd3ea98cf" } } }, "nbformat": 4, "nbformat_minor": 4 }