tutorial01-2.ipynb 31.6 KB
Newer Older
sp2668's avatar
sp2668 committed
1 2 3 4 5 6
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
sp2668's avatar
sp2668 committed
7
    "# Energy System Modelling - Tutorial I.2"
sp2668's avatar
sp2668 committed
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use approximations to seasonal variations of wind and solar power generation $W(t)$\n",
    "and $S(t)$ and load $L(t)$:\n",
    "$$W(t) = 1 + A_W \\cos \\omega t$$\n",
    "$$S(t) = 1 - A_S \\cos \\omega t$$\n",
    "$$L(t) = 1 + A_L \\cos \\omega t$$\n",
    "\n",
    "The time series are normalized to $\\langle{W}\\rangle = \\langle{S}\\rangle = \\langle{L}\\rangle := \\frac{1}{T} \\int_0^T L(t)\n",
    "d t = 1$, and the constants have the values\n",
    "\n",
    "$$\\omega = \\frac{2\\pi}{T} $$\n",
    "$$T = 1 \\text{ year}$$\n",
    "$$ A_W = 0.4 $$ $$ A_S = 0.75 $$ $$ A_L = 0.1 $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
sp2668's avatar
sp2668 committed
32
    "## Imports"
sp2668's avatar
sp2668 committed
33 34 35 36
   ]
  },
  {
   "cell_type": "code",
sp2668's avatar
sp2668 committed
37
   "execution_count": 2,
sp2668's avatar
sp2668 committed
38 39 40
   "metadata": {},
   "outputs": [],
   "source": [
sp2668's avatar
sp2668 committed
41 42 43 44 45
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate as integrate\n",
    "from scipy.optimize import minimize\n",
    "%matplotlib inline"
sp2668's avatar
sp2668 committed
46 47 48 49 50 51
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
sp2668's avatar
sp2668 committed
52
    "## Parameters"
sp2668's avatar
sp2668 committed
53 54 55 56
   ]
  },
  {
   "cell_type": "code",
sp2668's avatar
sp2668 committed
57
   "execution_count": 30,
sp2668's avatar
sp2668 committed
58 59 60
   "metadata": {},
   "outputs": [],
   "source": [
sp2668's avatar
sp2668 committed
61 62 63 64 65
    "A_w = 0.4\n",
    "A_s = 0.75\n",
    "A_l = 0.1\n",
    "T = 1\n",
    "phi = 0"
sp2668's avatar
sp2668 committed
66 67 68 69 70 71 72 73 74
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
sp2668's avatar
sp2668 committed
75 76 77 78 79 80 81
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Remark:** The `lambda` operator is a way to create small anonymous functions, i.e. functions without a name, in Python. "
   ]
  },
sp2668's avatar
sp2668 committed
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def w(t, phi=0):\n",
    "    return 1 + A_w * np.cos(2*t*np.pi/T-phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def s(t):\n",
    "    return 1 - A_s * np.cos(2*t*np.pi/T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def l(t):\n",
    "    return 1 + A_l * np.cos(2*t*np.pi/T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def c(gamma):\n",
    "    return 1 - gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mismatch(???):\n",
    "    return ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(alpha):\n",
    "    return 1/T * integrate.quad(lambda t: mismatch(???)**2, 0, T)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0,1,0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Availability time series as shown on the worksheet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VMX6wPHvpPdeIYSEEmoghBAIUsUCiCg2mghYEEXs2K7XguWqV3/XgoCAgChS7YgVkSIJJAFC7yQQSO+97M7vjwMRNZAQNjnZzXyeJ4+GnT3nDey+mZ3yjpBSoiiKolgWK70DUBRFUUxPJXdFURQLpJK7oiiKBVLJXVEUxQKp5K4oimKBVHJXFEWxQCq5K4qiWCCV3BVFUSyQSu6KoigWyKauBkKIxcAoIFNK2b2Wx92Bz4Dgc9d7W0q5pK7r+vj4yJCQkMsOWFEUpSVLTEzMllL61tWuzuQOLAXmAMsu8vgM4ICU8kYhhC9wWAixXEpZeamLhoSEkJCQUI/bK4qiKOcJIVLq067OYRkp5WYg91JNAFchhABczrWtrs/NFUVRlMZRn557XeYA3wJnAVdgrJTSaILrKoqiKA1kignV64HdQCsgApgjhHCrraEQYpoQIkEIkZCVlWWCWyuKoii1MUXPfSrwhtRqBx8TQpwEOgM7/t5QSrkAWAAQFRWlag0rinJZqqqqSE1Npby8XO9QGp2DgwNBQUHY2to26PmmSO6ngGHAFiGEP9AJOGGC6yqKovxFamoqrq6uhISEoE3zWSYpJTk5OaSmphIaGtqga9RnKeQKYAjgI4RIBV4EbM8FMB94BVgqhNgLCOBpKWV2g6JRFEW5hPLycotP7ABCCLy9vbmS4es6k7uUcnwdj58FrmtwBIqiKJfB0hP7eVf6c5piWEZRzF9lCWQcgPwUKC/QvowGcHDXvlz8ICAcnH30jlRR6kUld6VlqiqHE7/DoXVwejtkH0XbslEH10Bo3RvChkOnESrZK38xcuRIPv/8czw8POrVPjk5mVGjRrFv3z6Tx6KSu9KynNkJ2+fDwXVQVQL2bhAyALrfqvXMvdqDo4fWW7ey+bMXX5AK6XshfQ+kxGq/FIQVtL0K+twLnUeBtXo7tXTr16/XO4Qa6tWoWD4p4fB6+ON9OB0Hdq7Q43bociOEDAIbu4s/19lH+/JuD+0G/3m99D3aL4g9q2DNZPAIhr7TIepusHVsmp9LaXJvvfUWDg4OPPzwwzz22GMkJSXx22+/sWHDBpYsWcLWrVtJSEiguLiYESNGMGDAALZt20br1q355ptvcHR0JDExkbvvvhsnJycGDBjQaLGq5K5YtjM74ad/walt4NEWhr8BERPBodZ9dvUjBAT21L6GPKP94oidCz89p/33mpcg/DatndJoXv5uPwfOFpr0ml1bufHijd0u+vigQYN45513ePjhh0lISKCiooKqqiq2bt3KwIED2bp1a03bo0ePsmLFChYuXMgdd9zBF198wZ133snUqVP54IMPGDx4MLNmzTJp/BdSJX8Vy1SaC189AAuHQs5RGPUuzNwJ/R64ssT+d1bW2ieAu3+AyevAyQu+vBcWXQPpph9HVfTVu3dvEhMTKSoqwt7enpiYGBISEtiyZQsDBw78S9vQ0FAiIiJqnpecnExBQQH5+fkMHqx9Cpw0aVKjxap67orlOfoLfPMQlGbDgMdgwOOmTegXEzoQpm2CPSvhlxdgwRAY+hxc9Yj2S0AxqUv1sBuLra0tISEhLFmyhP79+9OjRw82btzI8ePH6dKly1/a2tvb1/y/tbU1ZWVlSCmbbCmn6rkrlqOqDL57FJbfpvWg7/tNGyJpisR+npUVREyAB7dD55Gw4WVYPBzy6lWlVTEDgwYN4u2332bQoEEMHDiQ+fPnExERUa+k7eHhgbu7e83wzfLlyxstTpXcFctQkKol0cQl0H8m3LdRGxPXi7M33P4J3LIIsg5rvfgTm/SLRzGZgQMHkpaWRkxMDP7+/jg4OPxjSOZSlixZwowZM4iJicHRsfEm34VW76vpRUVFSXVYh2ISyX/A6rugugJuWaD1mJuTnOOwYjzkHIPrXtXG/dVka4McPHjwH8Mflqy2n1cIkSiljKrruarnrpi3PWtg2Whw9NSGYZpbYgdtGeV9G7RNTz89C98/ru1+VZRGpJK7Yr7i5msrU4JjtOTpG6Z3RBdn7wp3fKpN8CYshrVTtU8aitJI1GoZxfxICb+9Clve1naG3vox2DroHVXdrKy0CV5nX21NfFkejPtcS/yKYmKq566YFym1ZYZb3oZek7RJS3NI7BeKmQE3z9fmCj67DSqK9Y5IsUAquSvmZeNrsO19iLoHRn9gvvVcIsbDbYshNR5WjIPKUr0jUiyMSu6K+dj0X9j8X4i8C0a+bf4rTrrdDGM+guStsGqiVqlSUUxEJXfFPOxYCBtfhR7jtFICVhby0u1xO9z0IRz/Db64R62isSBDhgxBz+Xedb5DhBCLhRCZQoiLFsoQQgwRQuwWQuwXQqidGoppHVwH62dBp5FaIrS0rfy9JsLwN7Uywj8+q80rKC2OwWDaX+z16f4sBYZf7EEhhAcwFxgtpewG3G6a0BQFOL1D69G27q2tijHXMfa69JsOMQ/Bjo8gdo7e0SgXUVJSwg033EDPnj3p3r07q1atYsOGDfTq1Yvw8HDuvvtuKir+ucT1gQceICoqim7duvHiiy/W/HlISAizZ89mwIABrFmzxqSx1ucM1c1CiJBLNJkAfCmlPHWufaZpQlNavJzj8PlYcGsFE1aBnZPeETWua1+BwjPw8/Pg1hq636J3RM3bD89oB6iYUkA4jHjjog//+OOPtGrViu+//x6AgoICunfvzoYNGwgLC+Ouu+5i3rx5PProo3953muvvYaXlxcGg4Fhw4axZ88eevToAYCDg8NfSgWbiikGLsMATyHE70KIRCHEXSa4ptLSVRTBygna/09c2zKOs7Oy0pZIBsfA1w9CWpLeESl/Ex4ezq+//srTTz/Nli1bSE5OJjQ0lLAwbQPd5MmT2bx58z+et3r1aiIjI+nVqxf79+/nwIEDNY+NHTu2UWI1xWdcG6A3MAxwBGKFEHFSyiN/byiEmAZMAwgODjbBrRWLZDTCV9O1c00nfaVt328pbB20nawLBsPKiTDt95bxi60hLtHDbixhYWEkJiayfv16nn32Wa677ro6n3Py5Enefvtt4uPj8fT0ZMqUKZSX/7kyytnZuVFiNUXPPRX4UUpZIqXMBjYDtZbjk1IukFJGSSmjfH19TXBrxSJtfkubXLz+tT+PtmtJXHxh3HIoyYLVk8FQpXdEyjlnz57FycmJO++8kyeffJJt27aRnJzMsWPHAPj0009rDuI4r7CwEGdnZ9zd3cnIyOCHH35oklhN0XP/BpgjhLAB7IC+wP9McF2lJTr8A/z+H+g5QTuTtKVq1UvbpPXlfdoY/Ig39Y5IAfbu3cusWbOwsrLC1taWefPmUVBQwO233051dTV9+vRh+vS/vm579uxJr1696NatG+3ateOqq65qkljrLPkrhFgBDAF8gAzgRcAWQEo5/1ybWcBUwAgsklK+W9eNVclf5R/yUuCjgeAZAnf/bH5lBRrDj89C3Fy4Yxl0vUnvaHSnSv7Wv+RvfVbLjK9Hm/8C/62rnaJcVHUlrL1bW+N9+1KV2M+75mU4vV07NjCgB3iF6h2RYiYsZJufYvY2vAxnErShCK92ekfTfNjYwW1LtFILa6aoMsFKvankrujv8A/axp0+92n1VpS/8mwLN82FtN1aRUxFqQeV3BV9FWdqQw7+4doRdErtuozSJpi3z4djv+odjWIGVHJX9COlltgriuDWhWqcvS7XvAS+neHrGVCaq3c0SjOnkruin8QlcPQnuHY2+LWcFRANZusItyyE0hz47hFVYEy5JJXcFX1kH4Of/gXthkL0NL2jMR+BPeDq5+Hgt5C0Qu9oWiQXFxeTXOell17i7bffNsm1aqOSu9L0jAb4+gGwtoOb51pObfam0n8mtB0A65+CglS9o1GaKfWuUppe3DxI3QEj/6tVfFQuj5U13DQHpEENz+hISsmsWbPo3r074eHhrFq1CoDi4mKGDRtGZGQk4eHhfPPNNzXPee211+jUqRPXXHMNhw8fbtT4LLQ4ttJs5RyH316BsBEQrkr/N5hXqDbB+sNTsPtz7cCPFubNHW9yKPeQSa/Z2aszT0c/Xa+2X375Jbt37yYpKYns7Gz69OnDoEGD8PX15auvvsLNzY3s7Gz69evH6NGj2blzJytXrmTXrl1UV1cTGRlJ7969TRr/hVTPXWk6RqO2OsbGHkb9z/zPQNVbn/sguD/89CwUpukdTYuzdetWxo8fj7W1Nf7+/gwePJj4+HiklDz33HP06NGDa665hjNnzpCRkcGWLVsYM2YMTk5OuLm5MXr06EaNT/XclaYTvwhObdM25LgF6h2N+bOy0oZn5vWHdY/B+BUt6hdmfXvYjeVidbmWL19OVlYWiYmJ2NraEhISUlPiVzThv4/quStNo+CMVmKg/TCImKB3NJbDu722eubID9oKGqXJDBo0iFWrVmEwGMjKymLz5s1ER0dTUFCAn58ftra2bNy4kZSUlJr2X331FWVlZRQVFfHdd981anyq5640jR+e0lbJjPq/FtW7bBJ9H4A9q7XVM+2GgIO73hG1CGPGjCE2NpaePXsihOCtt94iICCAiRMncuONNxIVFUVERASdO3cGIDIykrFjxxIREUHbtm0ZOHBgo8ZXZ8nfxqJK/rYgB9fBqolahcMBj9bdXrl8Z3bComEQdQ/c0Hhrp/WmSv7Wv+SvGpZRGldFkdZr9+sGMTP0jsZytY7UNoPFL4JU1WlSVHJXGtvG16HwLNz4Hljb6h2NZRv6L3AN1Na+G6r1jkbRmUruSuNJ36dVMYyaCm366B2N5XNw047jy9gH8Qv1jqbR6DWU3NSu9OdUyV1pHFLC+ifBwQOu/rfe0bQcXW7UViRtfB2KMvSOxuQcHBzIycmx+AQvpSQnJwcHh4ZXSq1ztYwQYjEwCsiUUna/RLs+QBwwVkq5tsERKZZhz2o4FQs3vg9OXnpH03IIASPegrn94NcXYcx8vSMyqaCgIFJTU8nKytI7lEbn4OBAUFBQg59fn6WQS4E5wLKLNRBCWANvAj81OBLFcpQXwM/PQ+ve0GuS3tG0PD4d4KqHYcs7EDkZ2sboHZHJ2NraEhqqzpGtj/ockL1ZCBFSR7OZwBdAow+sns4tJfZ4Dva2VtjbWONkZ42Piz3+bvZ4OtlhZaXWUOvu9zehJAsmrFIVH/Uy8AlIWqUNjU3bBNZqS4uepJQUV1STUVhBZmE5fm4OdPAzTengi7nif3EhRGtgDHA1dSR3IcQ0YBpAcHBwg+63+3Q+T32xp9bH7KytCPVxpqO/C538XYkI9qBXsCcu9uqF3WSyDsOOj6D3ZG15nqIPO2cY/jqsvgt2fgJ97tE7ohajymDkwNlCElPyOJxexJHMIo5lFFNU8ecKpvsHtePZkY27Xr9em5jO9dzX1TbmLoRYA7wjpYwTQiw9167OMfeGbmIqrzKQU1JJeZWB8ioDpZUGsosqyCgsJ62gnGOZxRzJLOJ0bhkAVgK6tnJjQAdfru3qR0QbT6xV777xfHYrnI6Hh3eCs4/e0bRsUsLSUZB5QPv3cPTUOyKLdSKrmA0HM/n9SCY7U/IpqzIA4O1sR0d/Fzr6uRLk6UiAuwN+rg6E+jgT4N6wydL6bmIyRZc2Clh5riCODzBSCFEtpfzaBNf+Bwdba1p7ONbZrqi8il2n8klIziXuZC4Lt5xg/qbjeDnbMbx7ALf0ak3vtp5NWsjH4h35WTu8+brXVGJvDoSA4f+BBYNh01va/ysmk5JTwle7zvBt0llOZJUAEObvwh1RQfQJ9SKqrVeDE7gpXHHP/W/tltLIPfeGKiirYtORLH45kMGvBzIoqzIQ7OXE2D5tGNenDd4u9k0Wi0WqrtSqEyLhgViwsdM7IuW87x6BXZ9p/y6+YXpHY9Yqqg2s35vG8rhTJKTkIQT0DfVieLcAhnXxp42XU6PHUN+ee53JXQixAhiC1ivPAF4EbAGklPP/1nYpzTS5X6i4opof96XzRWIqsSdysLOxYnTPVtwzIJQugW66xGT2Yj+En56DCash7Hq9o1EuVJwFH0RCm75wp1ql3BDZxRUsi03h8+0pZBdX0s7Hmduigrg5ojWt6jGSYEomS+6NpbkUDjuWWcQn21L4YmcqpZUGru3qz8NXdyQ8SFXWq7fSXHg/AlpHwZ1fqKqPzdG2OfDzv2DiF9DxGr2jMRsZheV8tOkEn+9IoaLayNWd/JhyVQgDOvjoNqSrkvtlKiitYsm2kyzeepLC8mqGdfbj6RGdCfN31Tu05u+HZ7QVMtP/AP+uekej1Ka6Ej6MBltHmL5VO4dVuai8kkrmbDzGp3EpGIySmyJaMWNoB9r7Nu7yxfpQyb2BisqrWBabwvxNxympqOaOqDY8dm0Y/m76TYw0aznHtaQRMRFGv693NMql7P8a1kyG0R9A5F16R9MslVcZ+GRbMnM2HqOkoppbI4OYeXVHgr0bfyy9vlRyv0J5JZV88NsxPo1Lxtbaikev6cjUq0KxtVabcv5i1Z1w7Dd4eBe4+usdjXIpUsLH10F+CszcCfb690Kbk42HMnnx2/2cyi1laCdfnh7Rmc4BzW8OTtVzv0Kezna8cGNXfn18MDHtvHl9/SFueH8L20/k6B1a85ESCwe/0w7gUIm9+RMCrn8NijNg2wd6R9NsnMkv4/5PE5i6NB4ba8Fn9/RlydToZpnYL4fqudfTLwcyeOnb/ZzJL2NSv7Y8PaJzy975KiUsugYKz8DMRG1HpGIeVk+Goz9rvfcWfFC50Sj5fMcp/rP+IEYJM4d14N4B7bCzad59XtVzN7Fru/rzy+ODuPuqUD7bnsL1/9vM1qPZeoelnwPfwJkE7YAIldjNyzUvgqEKNr2hdyS6OZ1bysRF23n+631EBHvw82ODeHBIh2af2C+H5fwkTcDJzoYXbuzKmvtjsLex4s6PtzP7uwOUn9tq3GIYqmDDbPDtAhET9I5GuVxe7SDqbtj5KWQd0TuaJiWl5IvEVIa/u5m9Zwr4zy3hfHZP3ybZfNTUVHJvgKgQL9Y/MpDJMW1Z/MdJbv7wD45kFOkdVtPZuQxyj2s9QLWkzjwNmqUti/xttt6RNJmCsioeXrmbJ9Yk0a2VOz8+OpDx0cEWW4JEJfcGcrC15uWburNkSh+yiyu48YOtrI4/rXdYja+iGH5/A4JjIGy43tEoDeXiC/0f1ibET8frHU2jSzqdzw3vb2H93jRmXd+JFdP6EeRpeb31C6nkfoWGdvbjh0cG0SfEi6e+2MOTa5Ioq7TgYZq4eVCSCdfOVjtRzV3MDHD2g19e0CbILZCUkmWxydw+PxYpYc30GGYM7dAiKsOq5G4Cvq72fHJ3NA8P68gXO1MZM/cPkrNL9A7L9Epy4I/3oPMoaBOtdzTKlbJ3gSFPw6lt2uoZC1NaWc0jK3fzwjf7uaqDN+tmDiAyuOWUPVbJ3USsrQSPXxvG0qnRpBeWM3rOVjYfsbBzHrf+H1SVwLAX9I5EMZXIyeAZChteAaNR72hM5nRuKbfOi+W7PWeZdX0nPp7cB0/nllWpVCV3Exsc5su3MwbQysORKUt2sGDzccs4qb3gDOxYCD3Hg28nvaNRTMXaVlvOmrEX9n+pdzQmEXs8h9FztpKaV8riKX2YMbRDizx+UyX3RhDs7cQXD/RnePcAXl9/iFlr91BZbea9ok1vgjTCkGf0jkQxte63gn932PiatszVjK3YcYpJH2/H28Web2ZcxdBOfnqHpBuV3BuJs70NH06I5JFhHVmbmMqkj7eTX1qpd1gNk31MO+wh6m7waNjZt0ozZmUFV/8bck9o/85myGiU/Gf9QZ79ci/9O/jw5YP9adcMKjjqSSX3RiSE4LFrw3h3bAS7TuUzZu42UnLMcKJ142tg4wCDntQ7EqWxhF2vHeax6U2oKtM7mstSXmXggeWJfLT5BJP6tWXx5CjcHGz1Dkt3Krk3gZt7tebz+/qSX1rJLXO3sSc1X++Q6i/93Fhsv+ng0nI/4lo8IbSJ8qI0iF+kdzT1lldSyYSFcfx8IIMXRnVl9k3dsFGVW4F6JHchxGIhRKYQYt9FHp8ohNhz7mubEKKn6cM0f1EhXqx9oD+OdtaMWxDHxsOZeodUPxtfBwd3bcOLYtlCBkD7q2Hr/7TNas3c6dxSbp2/jX1nC5k7IZK7B4Ra7G7ThqjPr7ilwKW2Ip4EBkspewCvAAtMEJdFau/rwpcP9ifUx5l7P0ngi8RUvUO6tNREOLwe+s8ERw+9o1GawtDnoTQHts+vu62ODqYVcsu8bWQXVfDZPX0ZEd5yq1teTJ3JXUq5Gci9xOPbpJR5576NA4JMFJtF8nN1YNX9MfRr58UTa5JYvPWk3iFd3MZXwckb+k7XOxKlqQT1hk4jYdv7UNY8hw8TU3IZ+1Es1kKw9oH+RId66R1Ss2Tqwal7gB9MfE2L42Jvw+IpfRjeLYDZ6w7wfz8fbn5r4ZP/gOO/wYDHwF6dI9uiDH0Oygsg9kO9I/mH3w9nMnGRttRx7QMx6ozjSzBZchdCDEVL7k9fos00IUSCECIhK8vCdm9eJnsba+ZM6MXYqDa8/9sxXv7uAEZjM0nwUsJvr4JLAETdo3c0SlMLCIeuN0PcXK3kRDPxw9407luWQDsfF1bfH2Pxhb+ulEmSuxCiB7AIuElKedFXg5RygZQySkoZ5evra4pbmzUbayveuDWceweEsnRbMs9+uRdDc0jwJzZq9UYGPgF26g3UIg19DqpK4Y//6R0JAF/tSmXG5zvpGeTByvv74etqr3dIzd4VJ3chRDDwJTBJStmyKv+bgBCCf93QhYev7sCqhNM8vno31QYdd7NKqa2QcQuC3pP1i0PRl28nCL8ddiyCYn1Xdn2+/RSPr06iXztvlt0Trdaw11N9lkKuAGKBTkKIVCHEPUKI6UKI87NsLwDewFwhxG4hhPkcjNpMCCF4/LpOPDW8E9/sPsvMFbuo0ivBH9sAqfEw6AmwUb2jFm3w02Co1CqB6uSTbck899VehoT5snhKH5zsWvC5xZdJHZDdzHy89SSvrDvA9d38+WB8ZNOe6SglLLwaSrK1Q69tWlYVPaUWXz8I+76AR5LANaBJb71460lmrzvAdV39mTOhid8LzZg6INtM3TMglBdv7MpP+zOY8fnOpi04duQnOLsTBs9SiV3RDJqlFRPb2rRj74u2nGD2ugMM7xbAhxNVYm8I9TfWDE29KpSXR3fjlwMZPLi8iRK8lFoNGc8QrayvogB4hUKviZCwRCv73AQWbTnBq98fZGR4AB9M6IWtKifQIOpvrZma3D+EV27qxq8HM5i5Ymfjj8EfXg/pe2DQU1qNb0U5b9Asrdzz1v9r9Fst3nqyJrG/N04l9iuh/uaasUkxIbx0bojm4cacZJUSfv8PeLWDHmMb5x6K+fIIhl53ws5lUNB4JTOWxSYz+9x8k0rsV0797TVzU64K5d+juvLDvnQeW7W7cdbBH/peq/446CmwVqsRlFoMfFzrBDTS2Pvy7Sm88M1+ru2qLSRQif3Kqb9BM3DPgFCeG9mZdXvSmLU2ybQ7WY1G+P0N8GqvrWtWlNo0Yu99bWIq//pqH0M7+TJnQi81eWoi6m/RTEwb1J7Hrw3jy51neP6bfaarRXP4e+38zMGq167UYeATWu99i+nG3r9LOstTa5MY0MGHeXf2xt7G2mTXbunM7t28JXULb8a/iZWwwgorrKyssLWyxdbKFjtrOxysHXC0ccTBxgEXWxdc7FxwtXXF3d4dD3sPPBw88HLwwsfRBycbJ7Oq/zzz6g5UVBv4cONx7G2seGFU1yuL32iE39/Ueu3dbzNdoIpl8mgDkZO03vvAx8H9ygrA/rQ/nUdX7SYqxIuFd0XhYGteib3SUElOWQ455TnkleeRX5FPQUUBRZVFFFUVUVJVQmlVKeXV5ZRVl1FprKTSUEmVsYqbO9zMpK6TGjU+s0vubvZudPXqihEjRmnEYDRQLaupMlRRaawktzyXsuoyyqrLKK4qpriyGEntvVwHawd8HH3wd/bHz8mPAOcAWju3JtAlkCCXIFq7tsbeuvns0hRC8OR1nSivMvLx1pM429nw5PWdGn7B8732MQtUr12pnwGPw85PYcs7MKrh4++bj2Qx8/Nd9AhyZ/GUPjjaNa/EbjAaSC9N50zRGc4UnyGtJI20kjQySjLILM0ksyyTosqiiz7f0cYRV1tXnGydajqbdtZ2ONs6Y2dlh4d945+PYHbv6J6+Pek5uP6HPRmlkdKqUgoqC8ivyCe/PJ/c8lxyynLIKssiqyyLjJIM9mTt4ZeUX6g2Vtc8VyDwc/Ij2C2YELcQQtxCCHUPpb1HewKcA7ASTT+qJYTg+Ru6UFpZzZyNx3Cyt+bBIR0u/0J/6bXfavpAFcvk0ebc2Pun2jBNA3rvO07mMu3TBDr4ubB0SjQu9vqlobzyPI7nH+dEwQmSC5NJKUwhpTCFM8Vn/pELfBx9CHAOIMQ9hOjAaHwcffBx9MHLwQsPew88HTxxt3PHxc4FGyv9U6v+ETQyK2GFi502PNPapfUl2xqlkeyybM4WnyW1OJXTRadJLUolpTCFn1N+pqCioKato40j7d3b09GzIx09OxLmGUZnr86427s39o+EEIJXbw6ntNLAWz8extnOhsn9Qy7vIofXn+u1f6R67crlGfg47PpMWzlzwzuX9dSk0/ncvTSe1h6OLLsnGnenptlTUWGo4GjeUY7kHeFI3hGO5h3lWP4xcsv/PIfIwdqBYLdgwjzDGBY8jDaubQhyDaK1c2sCnAOwNbP9H6q2zGXIK8/jRMEJ7Sv/BEfzj3I07+hfXiCtnFvR2aszXb270s2nG129u+Ll0DgnxVQZjDy4fCe/HMjgndt7cmvvevaipISPBkJlCcyIV8lduXzfPQK7P4eHd4P7pTtN5x3JKOKOj2Jxsbdh7fT+BLg7NEoLaIchAAAgAElEQVRoZdVlHM49zP6c/RzIOcCBnAOcLDiJQRoArWPW0aMjHTw70N69Pe092hPqHqrbp/HLVd/aMiq5m0B2WTZHco9wKO8Qh3IOcTD3IMmFyTWPt3ZpTbhPOOE+4fTw7UFX767YWZumdkt5lYF7Pokn9ngOcyf2Znj3ehR3OvQ9rJwAN8+DiAkmiUNpYfJS4INI6D0Vbni7zuanckq5bf42ANZO70+wt2nOCTBKI8kFySRlJbEnew/7svdxNO9oTSL3cfShi1cXOnt1pot3Fzp5diLINcgskvjFqOSus6LKIg7lHmJf9j72Zu9lX/Y+0krSALC1sqWrd1cifCPo5deLCL8IvB29G3yvkopqJn28nX1nClk0OYpBYZc4CEVK+GgQVBTBQwmq16403LcPQ9IKrWKkW6uLNksvKOf2j7ZRXF7Nqvuv7Gi80qpS9mXvY1fmLnZn7SYpK6lmYtPV1pXuPt3p7tOdcJ9wuvl0w8/Jr8H3aq5Ucm+Gskqz2JO1h6SsJHZn7WZ/9n4qjZUAhLiF0Nu/d81XK5eLv1lqU1BaxbiFcSRnl/DZvdH0bnuRoaBD62HleLhprlYQSlEa6nzvPepuGPnfWpvkllRyx0expBeU8/l9fekRdHmrRAorC9mVsYvEjEQSMxI5kHOAaqlNdHbw6KAtsPDtSU+/noS4hZh1j7y+VHI3A5WGSg7kHGBn5k52ZuxkZ+bOml5IK+dWRAVEEeUfRXRgdJ2TwQBZRRXc8VEsOcUVrLo/hi6Bbn9tICUsGKwdfvxQgioQply5bx6CPavP9d4D//JQcUU1ExbGcTi9iE/ujqZfu7o/nRZUFJCQkUBCegIJGQkczj2MRGJrZUt3n+5E+kUS6R9JT9+eTbJ4oTlSyd0MGaWRo3lHSchIIDEjkYT0BPIq8gBt3D46IJrowGiiA6Iv+nEzNa+U2+bFUm2UrJkeQ6iP858PHv4RVoyF0XO0zSiKcqVyT8IHvSH6PhjxZs0fl1cZmLoknh3JuSyY1JthXfxrfXppVSmJGYnsSN/B9rTtHMo9hERib21PT9+eNR2ccJ9wHGwaZwLW3JgsuQshFgOjgEwpZfdaHhfAe8BIoBSYIqXcWdeNVXKvm1EaOZ5/nB3pO4hPjyc+PZ7CykIAQt1D6RfYj76BfekT0Ac3uz976ccyi7h9fixOdjZ88cC5VQnnT1kqzdFOWVK9dsVUvp4B+9bWnNZUbTAy/bOd/Howg3fHRnBzrz8/dVYZqkjKSmJ7+na2p21nb9ZeqmU1tla29PTtSXRANH0C+tDDt4fJFh1YGlMm90FAMbDsIsl9JDATLbn3Bd6TUvat68YquV8+g9HAkbwjbE/bTlx6HDszdlJWXYaVsKK7d3f6BvalX2A/IvwiOJRWyvgFcbTycGT1/TF4nt0Ey2+DG99XB18rppV7Aj6Igr7TMV73Gk+uTeLLnWd4eXQ37oppy5G8I8SlxRGXFkdiRmLNa7abdzf6Bvalb2BfInwjVM+8nkw6LCOECAHWXSS5fwT8LqVcce77w8AQKWXapa6pkvuVu7AXFHc2jr3ZezFIA442jkT6RxJo24PPNznQyasjX9m9iHVJJszcqY7QU0zvq+nI/V/zdpc1zN11mhuii3D2OMH2tO3klOcAl/60qdRffZO7KdbBtQZOX/B96rk/u2RyV66crbWtNiYZEMWMiBkUVxYTnx5PXFocsWmx/FHwB3Zt4XS1A/8qzyG62xgGVObjZ2N5y8MU/ZRUlZDQaQjbTv/IDzmP4dKxgk154FXmRb/AfsS0iqFfYD8CnJv2gO2WzhTJvbayhLV+HBBCTAOmAQQHB5vg1sqFXOxcGBo8lKHBQwFIL0kn9mwsm35/la2OTnyf9jOs+Zn27u1r3nBRAVE42zrXcWVF+VO1sZp92fuITYsl7mwce7L2UC2rsXF1I6qikP7dH6Z/6LV09OzYIpYmNlemSO6pQJsLvg8CztbWUEq5AFgA2rCMCe6tXEKAcwBjrDwYc+YYWzs9y12HvInslImXUwprjqzhs4OfYSNs6OHbg36B/ejXqh/dfbpja6UmW5U/SSk5WXiSuLPauHl8ejzFVcUIBF28uxDjcws/JbhxR6tAXs+4H5GbCb2voFqpYhKmSO7fAg8JIVaiTagW1DXerjQRKWHTm+DaigG3P8bMjSm8t+Eo9w4Yw/vj2pGUlUTs2Vji0uKYlzSPuUlzcbJxIiogir4B2kSX6n21TOkl6TXLE+POxpFZlglAkEsQw0OHaxOhAX3ZebKSaZ8mEh3ixYuT+yC+/Rl2LIL+j4Bzw3ddK1euzuQuhFgBDAF8hBCpwIuALYCUcj6wHm2lzDG0pZBTGytY5TKd3AynYmHEf8HGnkev6UhBWRWLtp7Ew8mWh67WEjhom0fOv5m3p21nc+pmADztPekT0EdbohbYh1C3ULM64ESpn+yybBLSE9iRvoMd6TtIKUwBwMvBi+iA6JpVLW1c//yQHncihweX76RroBsL7uqtHbYxaBbsXQuxH8A1L+nzwyiA2sRk2ZaM1JapPbwbbLVlZkaj5Mk1SXy56wyv3NSNSTEhtT71wp7b9rTtZJRmAODr6EuUf1TN5pJQd5XszdH5ZJ6QkUB8ejwnCk4A4GzrrO2KPpfQL/bJbW9qAeMXxhHg7sDq+2Pwcr5gBdaaqXD0Z3h0Lzg1TkXUlqwpV8sozVHyVkj5A4a/WZPYAaysBG/e1oPC8mpe+HY/bo623BTxz9IGAc4BjG4/mtHtRyOl5HTR6ZpeXWJ6Ij8k/wBoPbtIv0h6+/eml38vOnl2ahYHFSh/klKSWpTKrqw/a7Sc75k72TjRy68Xo9uPJjogmi7eXer89zuWWczkJTtwd7Tl03ui/5rYQTuPd/9XEDsHhr3QWD+WUgfVc7dUS0dB9hFt16Ct4z8eLq8yMGXJDhKS85h/Z2+u6Vr79vDanE/258skJGYkcqb4DKDVyu7h24MI3wh6+vakh2+PFlsDRC/naxYlZSWRlJXErsxdZJdlA+Bm51ZTnyXKP6peyfxCqXml3D4/liqDZO30GEJ8LrLSavVkOLYBHt2jeu8mpmrLtGTJf8DSkXD9fyDmwYs2K66oZuLCOA6mF7F0ah/6t/dp8C3TS9LZlblLK8WauZvDeYcxSiOgbV45X8++u093wjzD1NZyEzFKIymFKezP2c/eLK209MHcg1QZqwCtJlEvv141paU7eHRo8AR5ZlE5d8yPJbeksvbCdBfK2A/z+mtj8Fc/36D7KbVTyb0l++RGyDyk9Zpq6bVfKK+kkrELYjmTV8Zn9/alV7CnSUI4X3c7KSuJPVl72JO9p+bEKhsrGzp6dKSrd1c6e3Wms1dnwjzDcLI1zQEOlqraWM3JgpMcyj3EodxDHMg5wMHcg5RUlQDap6au3l3p4dOjpgyuj2PDf2FfKL+0krEfxXE6r5TP7u1LZH1eJ6smwYnftU+PqvduMiq5t1Qp22DJCLj+dYiZUa+nZBaWc9v8WArKqlh1fz86B5h+W7iUkrSSNPbn7Gd/9n725+znYO7BmnNpBYJgt2A6emhn0rb3aE979/a0dWtrdmdXXimjNJJeks6x/GN/Oc7xeP7xmvr/9tb2hHmGacc5emvHObb3aN8o8x3FFdXcuWg7B84WsnhKHwZ0rOcvjPR9MP8qGPQUXP0vk8fVUqnk3lJ9MhoyD2q9Jbv694RP52pjqdVGyer7+9HO16URg9RIKUkvSedg7kEO5x6uSWIphSnIc5ucrYU1rV1aE+IeQohbCG1c29R8BToHmm3iN0ojWaVZNQexny46TUphCskFyZwqOkVZdVlNW19HX8I8w/5yEHuoe2iTTFyfn5uJT85j7sRIru92mSUEVt0JJzZpnyIdTfOpsKVTyb0lSomFJcPhuteg/0OX/fRjmcWM/SgWexsrVk+PIchTn2GS8upykguTOZ5/nOP5x0kuTCa5MJlThaeoMFTUtBMIfJ18aeXcCn9nf/yd/PFz8sPX0RdvR2+8HbzxdPDE3d69yVbwGKWRosoi8ivyySnLIac8h+yybDJLM8kszSSjJIOzJWdJL0mvGRcHsBJW2i8xtxBC3EMIdQ+lg0cH2rm3021CurLayLRPE9h0JIt3x0bUuqqqTul7Yf4AGPwMDH3W9EG2QCq5t0TLbtImsh7Zc1m99gsdOFvIuAWxeDrbsfr+GPzdmk8Z1vO93fM93bSSNM4Wn+VsydmaxFluKK/1ua52rrjZueFi64KLnQvOts442jjiaOOIvbU9tla22FrbYiNssBJWWAkrBAKDNGCURgzSQJWxiipDFVXGKkqrSymrLqOsqoyiqiKKK4spriqmoKKg5nDmC9kIG3ycfPBz8qOVcysCXQJp5dzqz08hLoHNquxDtcHIzBW7+GFfOv+5JZzx0VdQC2rlRDi55Vzv/fKO2VP+SSX3lqam1/4q9J95RZfadSqPOxdtJ9DDkZXT+uHjYm+iIBuXlJLCykJyynO0XnNZDnkVeeSX55NXkUdRZRHFVcUUVxZTUlVCuaGcsuoyKg2VVBoqqTJWUW2sxiiNNcNCoA0NWQkr7KztsLWyxc7KDkdb7ReDg7UDLnYuuNq64mrniru9Ox72Hng6eOLt4K19gnD0xsvBy2zKOBiMkidW7+br3Wd5/oYu3Duw3ZVdMG0PfDRQ9d5NRCX3lqam154Edlde5XH7iRwmL9lBiLczK6f1w8OpZS1dlFIikWaTkE3FaJQ8++VeViWcZtb1nZgxtINpLrxyolYOQ429X7H6JveW9cq1VCmx2pKzqx41SWIH6NvOm4V3RXEiu4RJH++goKyq7idZECFEi0vsUkpe+m4/qxJO8/DVHUyX2AGGPAMVhRA3z3TXVC6pZb16LdXv/wFnP4i626SXHdjRl3kTIzmUXsjkxTsoKm9ZCb4lkVIye90BlsWmMG1QOx67Nsy0NwgIhy43asm9LM+011ZqpZK7uUvZBic3wYBHGzyJeinDuvgzZ0Ik+84UMGVJPMUV1Sa/h6IvKSWvfX+QJX8kM/WqEJ4d0blxisENPtd7j51r+msr/6CSu7k732vv3XiVlq/vFsAH43ux+3Q+U5fsoEQleIshpeSNHw+xaOtJJse05YVRXRuvymdAd+gyWuu9l+Y2zj2UGiq5m7Pkrdok1YDHGqXXfqER4YG8P64XO0/lM2XJDtWDtwDnE/tHm05wZ79gXhrdrfHLNw95BiqLIPbDxr2PopK72ZISNr4OLgEQ1TTno9zQI5D3xkVoCV6NwZs1KSWvrz9Yk9hnj+7eNHX5/btBtzGwfT6U5DT+/VowldzN1cnNWr32gU/UWRzMlEb1aFUzRDN58Q4KVYI3O1JKXv3+IAu3nOSumLa8clN3rKya8MCVwc9AZQlse7/p7tkCqeRujs732t1aQ+RdTX77keGBzJnQiz2pBUxatJ380somj0FpGKNR8sI3+/l460mm9A/h5aYYivk7v84QfhvsWADFWU177xakXsldCDFcCHFYCHFMCPFMLY8HCyE2CiF2CSH2CCFGmj5Upcbx3+B03Lleuz7lAYZ3D2T+nb05mFbEuAVxZBdX1P0kRVcGo+TpL/bwaZy23PHFGxtx8rQug5+G6nL441197t8C1JnchRDWwIfACKArMF4I0fVvzZ4HVkspewHjALXWqbFICRtfA/c20GuSrqFc09Wfj6dEkZxTwtiPYskorL2ui6K/KoORR1ftZk1iKo8M69h4yx3ry6cj9BgL8R9DUbp+cViw+vTco4FjUsoTUspKYCVw09/aSOB8EXB34KzpQlT+4siPcCZRO+HGRv+SAAM7+vLJ1GjSC8q5bf42UnJK9A5J+ZvyKgPTP03ku6SzPDOiM49dG9Y8DjUf/BQYKmHL/+kdiUWqT3JvDZy+4PvUc392oZeAO4UQqcB6oNbKVUKIaUKIBCFEQlaWGmu7bEaj1mv3agcRE/SOpkbfdt6smNaP4vJqbpsfy8G0Qr1DUs4pLK/irsU7+O1wJq/c3J3pg9vrHdKfvNpBrzshcQnkn667vXJZ6pPca/sV//dqY+OBpVLKIGAk8KkQ/yzMIaVcIKWMklJG+fr6Xn60Ld3Bb7X62EOehWZ2SEWPIA/WTI/BWgju+CiWhGS1SUVv2cUVjF8Qx86UPN4b14tJ/drqHdI/DZql/Xfzf/WNwwLVJ7mnAm0u+D6Ifw673AOsBpBSxgIOgGkOb1Q0RoO2Qsa3M3S/Ve9oatXBz5W1D8Tg42LPxEXb+Wm/GkvVS3J2CbfO28bxrGIWTY5idM9WeodUO4822u7qXZ9BznG9o7Eo9Unu8UBHIUSoEMIObcL027+1OQUMAxBCdEFL7mrcxZT2roXswzD0ObCy1juaiwrydGLt9Bg6B7rxwGeJfBaXondILc7u0/ncOm8bReXVrLivH0M6+ekd0qUNfBys7WDTW3pHYlHqTO5SymrgIeAn4CDaqpj9QojZQojR55o9AdwnhEgCVgBTpF6F4i2RoUqrIRMQDp1v1DuaOnm72LPivr4M7eTH81/v478/HcJoVC+HpvDboQzGL4jDyd6atdNj6BVsBrXTXQMg+j7YswoyD+kdjcVQh3WYg/iP4fvHYcIaCLtO72jqrdpg5N/f7GPFjtPc0COQd27viYNt8/3UYc6klCz5I5lXvz9At1buLJ7SB19X8zhBC9BKEbzXE9oNhnHL9Y6mWVOHdViKylLt42qbftDxWr2juSw21la8PiacZ0d0Zv3eNMYtiCOrSG12MrVqg5EXvtnP7HUHuKaLP6vu72deiR3A2Vs71P3QOm2pr3LFVHJv7nYsgOJ0uOZFaA5rky+TEIL7B7dn3sTeHEov5OYP/2DfmQK9w7IY+aWVTF0az6dxKdw/uB3z7+yNk52N3mE1TMwMcPKGDbP1jsQiqOTenJUXwNb/QYdroW1/vaO5IsO7B7Dm/v4YpeS2+dv4ZvcZvUMye4fSCxk95w+2n8jlrVt78OyILk1bAMzU7F21khonfocTm/SOxuyp5N6cbfsAyvNh2L/1jsQkwoPc+fahAYS3dueRlbv5z/qDVBuMeodlltbvTeOWudsorzKw8v5+3NGnTd1PMgdR92gF8TbM1kptKA2mkntzVZShHUfWbQwE9tQ7GpPxdbVn+b39uLNfMB9tPsGERdvJVDVp6q2y2sjL3+3nweU76RTgynczBxBpDiti6svWQSsqdiZBG39XGkwl9+Zq81tgqICrLaPXfiE7GytevTmc/43tyd7UAka+v4Vtx7L1DqvZO5tfxtgFsSz5I5kp/UNYNS0Gfzd9qoI2qoiJ4BOm9d4N6sSvhlLJvTnKOQ6JS6H3FPBuRrVATGxMryC+eegqPJzsmPjxdt768RBVapimVuv3pjHivS0czSjmwwmRvDS6G3Y2Fvr2tbaBYS9C9hHYrZZFNpSFvjrM3IbZYG0Pg57SO5JGF+bvyjczrmJsVBvm/n6cW+dt42S2qix5XklFNbPWJPHg8p2E+DizbuYAbugRqHdYja/zDRAUrW3eqyzVOxqzpJJ7c5OaCAe+1tb8uvrrHU2TcLa34Y1bezBvYiQpOaWMfG8LS/842eJ3tcYez2HEe1v4YmcqM6/uwNrpMYT4OOsdVtMQAq59GYrSYPs8vaMxSyq5NydSwq8vgpMPxDykdzRNbkR4ID89Ooi+7bx46bsDjFsQ1yJ78cUV1Tz/9V7GL4zDSsCq+2N44rpO2Fq3sLdr2/4QNgK2vqsO026AFvZqaeaO/gzJW7RDDBzc6m5vgQLcHVgypQ//va0HB9MLGf7uZj7YcJSKaoPeoTU6KSU/7kvjuv/bxPLtp7h3QCg/PDKIPiFeeoemn2tehMpiVRK4AVRyby4M1fDzv8GrvVYCtQUTQnB7VBt+fXwww7r48c4vRxj+7hY2H7HcQqPJ2SVMWRLP9M924uZoy9rp/Xl+VFcc7Vp4LR6/LtpxkvELVUngy6SSe3Oxa5lW0vfa2c3i+LzmwN/NgbkTe/PJ3dEA3LV4B3cvjedIRpHOkZlOXkklr6w7wLX/20RiSh7/HtWVdTMH0LutBa1dv1JD/6UtMPj1Rb0jMSuqKmRzUFEE7/cC744wdb1Z1pBpbBXVBpb8kcyHG49RUlHN7b3bMHNYB4I8nfQOrUFKKqpZFpvC3N///Hkevy7MMtetm8Kmt7QjJqf+CG1j9I5GV/WtCqmSe3Pw26vamOK9v0FQb72jadbySiqZs/EYy2KTkRJujQziwaHtaettHqtIisqrWBabwqItJ8grrWJoJ1+eGdGFTgGueofWvFWWwAe9wa0V3LuhRXeAVHI3FwWp8EGUtq73to/1jsZsnM0vY/6m46yMP43BKBneLYDJ/UPoE+KJaIZv/NO5pSyLTWZl/GmKyqsZ2smXh67uqIZfLseuz+CbGXDrxxB+m97R6EYld3Ox9h6thsZD8eARrHc0ZiezsJyPt55kxY5TFJZX0zXQjfHRbRjVoxWezvrOXVRUG9h4KIu1iafZcCgTKyEY3j2A6YPaEx7krmtsZslogAVDoDQHHkoAO/MckrtSJk3uQojhwHuANbBISvlGLW3uAF4CJJAkpZxwqWuq5A6cioPF12s7Ua/+l97RmLXSymq+3nWWZbHJHEovwtZaMLSTHzf0CGRwmC8eTk2T6CuqDWw/kcvPB9JZtyeN/NIqfFzsGR/dhol92xLgrsbUr0jKNlgyAoY8C0Oe0TsaXZgsuQshrIEjwLVAKtqB2eOllAcuaNMRWA1cLaXME0L4SSkzL3XdFp/cjUZYOBSKM2FmAtiZx5hxcyel5EBaIV/tPMM3SWfJKqrA2krQu60nAzv40CfUi4g2HiY77s9glBxOLyIhJZfY4zlsPpJFSaUBB1srrusawJjI1gzs4INNS9uA1JhWT4YjP2nvG/cgvaNpcvVN7vU5siUaOCalPHHuwiuBm4ADF7S5D/hQSpkHUFdiV4CkFZC2G25ZqBK7CQkh6NbKnW6t3HluZBd2p+az4WAGGw5m8s4vRwCwtRaE+bsS5u9KR38XQryd8Xezx8/VAU9nOxxsrGqSsdEoqag2UlReRUZhBRmF5ZzJL+NIRhFHM4o5mF5IUblWuTDAzYHREa0Y1tmfqzr4qDXqjeXa2XDkR/j1Jbh1kd7RNFv1Se6tgdMXfJ8K9P1bmzAAIcQfaEM3L0kpfzRJhJaoogg2vAxBfSD8dr2jsVhWVoLIYE8igz2ZdX1n8ksrSUzJIz45jwNphcSdyOGrXbWfCGVjJbASgsqLVKl0c7ChU4ArN0W0ondbT6LaehHk6dgsJ3Mtjmdb6D9TW2HW5z4I/ns6UqB+yb22V+vfx3JsgI7AECAI2CKE6C6lzP/LhYSYBkwDCA5uwZOHv7+hDceMW9Gil3Q1NQ8nO4Z18WdYlz8LshWWV3Emr4yMwnIyCyvIL6ukvMpIRbUBgxEcbK1wsLXG2d4Gf1d7/NwcaOXugK+rvUrkehrwGOxaDuufhGm/g5X6lPR39UnuqcCFZ3gFAWdraRMnpawCTgohDqMl+/gLG0kpFwALQBtzb2jQZi3zEGyfD5GT1Jr2ZsDNwRa3QFu6BLbMWj5my84Zrn8N1k6FhMUQfZ/eETU79ZnliQc6CiFChRB2wDjg27+1+RoYCiCE8EEbpjlhykAtgpRaT8POBYa9pHc0imLeuo2B0EHw2ytQok7y+rs6k7uUshp4CPgJOAisllLuF0LMFkKMPtfsJyBHCHEA2AjMklKqGp1/t/9LrerjsH+Ds7fe0SiKeRMCRr6t7V799SW9o2l21CamplJRBHP6gIsf3LdRjREqiqn8/G/Y9j7c8yu06aN3NI2uvksh1eLbprLxdShKh5HvqMSuKKY0+ClwbQXrHgVDld7RNBsquTeFs7u0SdSou1tEz0JRmpS9K4x8CzL2QZw6ku88ldwbm6EavnsEnH1h2At6R6MolqnzKOg0UjtQOy9F72iaBZXcG9uOBZCWBCPeBEcPvaNRFMskBIz8LyC0FWk6zSU2Jyq5N6b801qt9o7XQdeb9Y5GUSybexBc/bx2FvH+r/SORncquTcWKbXhGNCWa6ndjIrS+KKnQWAErJ8FJS17NbZK7o1l93I4vgGueUmrhaEoSuOztoGb50J5AfzwlN7R6Eol98ZQeBZ+fA7aXgV97tU7GkVpWfy7waBZsG8tHPpe72h0o5K7qUkJ6x4DQyWM/gCs1F+xojS5gY+Df7j2XizN1TsaXajMY2pJK7Va08P+Dd7t9Y5GUVoma1u4+UPtSL4fntY7Gl2o5G5KeSnaRE5wf+g7Xe9oFKVlC+ypDc/sXQ37vtQ7miankrupGA3w1bmEPma+KjGgKM3BwCehdZQ2PFNQ+8Eslkold1P54z04tU3bSKFWxyhK82BtA7cs0ObAvn5AO7u4hVDJ3RTO7tYKg3W9CXqO0zsaRVEu5N0ern8dTm6C7S2n9oxK7leqvFA7DcbZF0a9qzYrKUpz1HuKVnvmlxfhzE69o2kSKrlfifO7UPNS4LbF4OSld0SKotRGCLjpQ3ANgDVToCy/zqeYO5Xcr0TCYu10pav/BW1j9I5GUZRLcfLSOmGFZ+Dbhyy+uFi9krsQYrgQ4rAQ4pgQ4plLtLtNCCGFEHWeEmL20vbAj89C+2Fw1WN6R6MoSn20iYZhL8LB77SKrRaszuQuhLAGPgRGAF2B8UKIrrW0cwUeBrabOshmpyQHVk0EJ29tJl7tQlUU8xHzEISNgJ+eg5RtekfTaOqTlaKBY1LKE1LKSmAlcFMt7V4B3gLKTRhf82OohrVToCgDxn0Gzj56R6QoyuWwstL2oniGwKpJWmluC1Sf5N4auPCnTz33ZzWEEL2ANlLKdSaMrXn6+Xk4uRlufBda99Y7GkVRGsLRA8Z9DtUV2qfwylK9IzK5+iT32tb21cxECCGsgP8BT9R5ISGmCSEShHJU5ywAAAkFSURBVBAJWVlZ9Y+yudj1mbZOtu8DEDFB72gURbkSvp3g1oXa/Nm3My1ugrU+yT0VaHPB90HA2Qu+dwW6A78LIZKBfsC3tU2qSikXSCmjpJRRvr6+DY9aD8d/05Y9thsC172qdzSKophCpxHa6U371sJvr+gdjUnZ1KNNPNBRCBEKnAHGATXdVillAVAz8CyE+B14UkqZYNpQdZS+F1bdBb6d4Y5l2pZmRVEsw8AnIP8UbHlHO6ov6m69IzKJOrOUlLJaCPEQ8BNgDSyWUu4XQswGEqSU3zZ2kLrKPw3LbwcHN5i4Bhzc9Y5IURRTEgJu+D8oSoPvnwDXVvD/7d1/kFV1Gcfx94PID5XEAsRRFjCgINqK2VFsClMJkAxohmQhAwNCceqf/mhs/MehmByasnRwEoIMjF/RD1ZLLJAfhaKsoyIwooSIW/wMxBYGBPbpj+/RdtZl79nde++595zPa+bO3LP33L3Ps+feZ7/3nO95zifGJB1Vu5kntJ+pqqrKa2tLfHD/30Pw2Feg/jBMXwtXfmgGqIikxZn68Hk/sjsM5Pp/MemImmVmL7p7znOJNEH7Qk4ehSXjw9lsU1aqsIukXefL4BurQ1fXZZNg/9akI2oXFffmnDoGSybA8TdDYVdrAZFsuKwnTF0TetA8PhHqXkw6ojZTcW+q/nAYsR/dDdW/hf4jko5IRIqpW2+Y9kToRbP0a/DWc0lH1CYq7o0d3weLR8N/9kD1chgwMumIRCQJl18Ndz4ZRvJLJ8DrTycdUaupuL/v0E5YNDrskpm6BgaqsItkWvcKmP50mAK9fDK8siLpiFpFxR3gtT+Hwm4dwqyYPtclHZGIlIJLe4QRfL8vwB/vgvVzwvWSy0C2i3tDA2z4MayYAj0GwMx10Gtw0lGJSCnp3C1MjRw2LZzotGxSWVzsI7vFvf4wLK+GTQ/AZ6bAt9aG/WwiIk117AzjHoLbHoS9G2HhTfCv0p5Jk83ivqsGHhkeNtKtP4EJj8DFXZKOSkRKXdX0sJvm3Bn41Zfhmblw/mzSUTUrW8X93QPw+5mw6puhh8Rdm+D6WbqotYjEVzEcZj8LlbfD5nmw8GZ4e1vSUX1INor7eydh4wPw8DDYtQZuvBdmrtf+dRFpm67dwwU/Jj0O9Ydg0UhYPSM0ICsR6W5veOpYuIj1Cwuh/iAMmQAj74eP9k86MhFJg8FfDW3At/wCnn04XJv1s5Nh+D2hX3yC0tc47PxZ2PcP2PkH2L4Kzp2Gj98MI76vNgIiUjgn6mDTPNi+MtSdASOhshoGjcprN9m4jcPKv7iffjecgHRwO9TVwht/hdPvQMeuUPn18B9Uu19EpFhOHoXaX0PtotBGuMPFoY1J389D70q4qjK0OGij9Bb3N/4Ga38Ap0+E2/kz/3/s0p7hv+Unbwuj9U6X5C9gEZHWaGiAum3w2hOw+6nQ1uR9N3wHRs9t06+NW9zLb5971yug99DwNadL97Dca3D4j9itt2a+iEhp6NABKq4Pt1E/CoPRgzvCld2u/FTBX778Ru4iIhmmi3WIiGRYrOJuZmPMbLeZ7TGze5t5/HtmtsvMtpvZejPrm/9QRUQkrpzF3cwuAuYDtwJDgMlm1vSacy8BVe5eCawG5uU7UBERiS/OyP06YI+773X394AVwPjGK7j7Bnc/FS1uBa7Jb5giItIacYr71cDbjZbrop9dyAzgqfYEJSIi7RNnKmRzcwubnWJjZncAVcCNF3h8FjALoKKiImaIIiLSWnFG7nVAn0bL1wD/brqSmY0E7gPGufuZpo8DuPsCd69y96qePXu2JV4REYkhTnHfBgw0s/5m1gmoBmoar2BmnwMeJRT2w/kPU0REWiPWSUxmNhb4OXARsNjd55rZHKDW3WvMbB3waeBA9JT97j4ux+88ArzVxrh7AEfb+NxylsW8s5gzZDPvLOYMrc+7r7vn3PWR2Bmq7WFmtXHO0EqbLOadxZwhm3lnMWcoXN46Q1VEJIVU3EVEUqhci/uCpANISBbzzmLOkM28s5gzFCjvstznLiIiLSvXkbuIiLSgpIt7jG6Unc1sZfT482bWr/hR5l8Wu3DmyrnRehPNzM0sFbMq4uRtZrdH23unmS0rdoz5FuP9XWFmG8zspeg9PjaJOPPJzBab2WEz23GBx83MHor+JtvNbFi7X9TdS/JGmFP/T+BaoBPwCjCkyTr3AL+M7lcDK5OOu0h53wRcEt2fXe55x8k5Wq8bsJnQnK4q6biLtK0HErquXhEt90o67iLkvACYHd0fAuxLOu485D0CGAbsuMDjYwk9uQwYDjzf3tcs5ZF7zm6U0fJvovurgVvMyv46e1nswhlnWwP8kNBO+nQxgyugOHl/G5jv7scBvPzPAI+TswMfie5fTjPtTsqNu28GjrWwynhgiQdbge5mdlV7XrOUi3ucbpQfrOPu54ATwMeKEl3hZLELZ86coxYXfdz9yWIGVmBxtvUgYJCZbTGzrWY2pmjRFUacnO8H7jCzOuAvwHeLE1qiWvu5z6mUL5Adpxtl7I6VZSRvXTjLSIs5m1kH4EHgzmIFVCRxtnVHwq6ZLxG+of3dzIa6+zsFjq1Q4uQ8GXjM3X9qZjcAS6OcGwofXmLyXstKeeQepxvlB+uYWUfCV7iWvvqUg7x14SwjuXLuBgwFNprZPsI+yZoUHFSN+x5f4+5n3f1NYDeh2JerODnPAFYBuPtzQBdC/5U0i/W5b41SLu45u1FGy9Oi+xOBZzw6OlHGstiFs8Wc3f2Eu/dw937u3o9wnGGcu9cmE27exHmP/4lwAB0z60HYTbO3qFHmV5yc9wO3AJjZYEJxP1LUKIuvBpgazZoZDpxw9wO5ntSipI8i5zjCPBZ4nXB0/b7oZ3MIH2wIG/13wB7gBeDapGMuUt7rgEPAy9GtJumYC51zk3U3koLZMjG3tQE/A3YBrwLVScdchJyHAFsIM2leBkYlHXMecl5O6Jp7ljBKnwHcDdzdaDvPj/4mr+bj/a0zVEVEUqiUd8uIiEgbqbiLiKSQiruISAqpuIuIpJCKu4hICqm4i4ikkIq7iEgKqbiLiKTQ/wBJP0YAD3kWlwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wind = plt.plot(x, w(x), label='wind')\n",
    "solar = plt.plot(x, s(x), label='solar')\n",
    "load = plt.plot(x, l(x), label='load')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
sp2668's avatar
minor  
sp2668 committed
193
    "## Problem I.2"
sp2668's avatar
sp2668 committed
194 195 196 197 198 199 200 201
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**(a) What is the seasonal optimal mix $\\alpha$, which minimizes**\n",
sp2668's avatar
sp2668 committed
202 203
    "$$\\langle\\left[ \\alpha W(\\cdot) + (1-\\alpha) S(\\cdot) - L(\\cdot) \\right]^2 \\rangle = \\frac1T \\int_0^T \\left[ \\alpha W(t) + (1-\\alpha) S(t) - L(t) \\right]^2 \\,\\mathrm d t$$\n",
    "> **Hint:** You can use the function [`scipy.optimize.minimize`](https://docs.scipy.org/doc/scipy-1.0.0/reference/generated/scipy.optimize.minimize.html) and use e.g. [`method='nelder-mead'`](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method)."
sp2668's avatar
sp2668 committed
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**(b) How does the optimal mix change if we replace $A_L \\to -A_L$?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**(c) Now assume that there is a seasonal shift in the wind signal\n",
    "$$ W(t) = 1 + A_W \\cos \\left( \\omega t - \\phi \\right).$$\n",
sp2668's avatar
sp2668 committed
235 236 237 238
    "Express the optimal mix $\\alpha$ as a function of $\\phi$.**\n",
    "> **Remark:** Note, that $\\alpha\\in [0,1]$ and you need to add this as bounds.\n",
    "\n",
    "> **Hint:** If you encounter problems (why?), try another optimisation algorithm [`method='TNC'`](https://en.wikipedia.org/wiki/Truncated_Newton_method)"
sp2668's avatar
sp2668 committed
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**(d) A constant conventional power source $C(t) = 1 - \\gamma$ is now introduced. The mismatch then becomes\n",
    "$$\\Delta(t) = \\gamma \\left[ \\alpha W(t) + (1-\\alpha) S(t) \\right] + C(t) - L(t)$$\n",
    "Analogously to (a), find the optimal mix $\\alpha$ as a function of $0 \\leq \\gamma \\leq 1$, which minimizes $\\langle{\\Delta^2}\\rangle$.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}