diff --git a/docs/source/examples/notebooks/creating_models/1-an-ode-model.ipynb b/docs/source/examples/notebooks/creating_models/1-an-ode-model.ipynb index 535aec383a..9130f35137 100644 --- a/docs/source/examples/notebooks/creating_models/1-an-ode-model.ipynb +++ b/docs/source/examples/notebooks/creating_models/1-an-ode-model.ipynb @@ -5,7 +5,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Creating a simple ODE model\n", + "# Creating a simple ODE model: the reservoir model\n", "\n", "In the [Getting Started](../getting_started/tutorial-1-how-to-run-a-model.ipynb) series, we learnt how to use in-built PyBaMM models. In this series, our focus will be on how to build PyBaMM models from scratch, from the very simple reservoir model all the way up to the Doyle-Fuller-Newman model. We strongly advise to familiarise yourself with the Getting Started notebooks, if you have not done so, before following this series.\n", "\n", @@ -16,14 +16,14 @@ "In this notebook we create and solve the following simple ODE model:\n", "$$\n", "\\begin{align*}\n", - "\\frac{dx_n}{dt} &= -\\frac{I(t)}{Q_n}, \\\\\n", - "\\frac{dx_p}{dt} &= \\frac{I(t)}{Q_p}, \\\\\n", - "V(t) &= U_p(x_p) - U_n(x_n) - I(t)R, \\\\\n", - "x_n(0) &= x_{n0}, \\\\\n", - "x_p(0) &= x_{p0}, \\\\\n", + "\\frac{\\mathrm{d} x_\\mathrm{n}}{\\mathrm{d} t} &= -\\frac{I(t)}{Q_\\mathrm{n}}, \\\\\n", + "\\frac{\\mathrm{d} x_\\mathrm{p}}{\\mathrm{d} t} &= \\frac{I(t)}{Q_\\mathrm{p}}, \\\\\n", + "V(t) &= U_\\mathrm{p}(x_\\mathrm{p}) - U_\\mathrm{n}(x_\\mathrm{n}) - I(t)R, \\\\\n", + "x_\\mathrm{n}(0) &= x_{\\mathrm{n}0}, \\\\\n", + "x_\\mathrm{p}(0) &= x_{\\mathrm{p}0}, \\\\\n", "\\end{align*}\n", "$$\n", - "where $x_n$ and $x_p$ are the dimensionless stochiometries of the negative and positive electrodes, $I(t)$ is the current, $Q_n$ and $Q_p$ are the capacities of the negative and positive electrodes, $U_p(x_p)$ and $U_n(x_n)$ are the open circuit potentials of the positive and negative electrodes, and $R$ is the internal resistance of the battery.\n", + "where $x_n$ and $x_p$ are the dimensionless stochiometries of the negative and positive electrodes, $I(t)$ is the current, $Q_\\mathrm{n}$ and $Q_\\mathrm{p}$ are the capacities of the negative and positive electrodes, $U_\\mathrm{p}(x_\\mathrm{p})$ and $U_\\mathrm{n}(x_\\mathrm{n})$ are the open circuit potentials of the positive and negative electrodes, and $R$ is the internal resistance of the battery.\n", "\n", "We begin by importing the PyBaMM library into this notebook, along with NumPy and Matplotlib, which we use for plotting:\n" ] @@ -37,8 +37,18 @@ "name": "stdout", "output_type": "stream", "text": [ + "\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.3.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m24.0\u001b[0m\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n", "Note: you may need to restart the kernel to use updated packages.\n" ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n" + ] } ], "source": [ @@ -70,12 +80,12 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "x_n = pybamm.Variable(\"Negative electrode stochiometry\")\n", - "x_p = pybamm.Variable(\"Positive electrode stochiometry\")" + "x_n = pybamm.Variable(\"Negative electrode stoichiometry\")\n", + "x_p = pybamm.Variable(\"Positive electrode stoichiometry\")" ] }, { @@ -97,7 +107,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -124,7 +134,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -152,230 +162,311 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We first initialise the model using the `BaseModel` class. This sets up the required structure for our model. " + "Now that we have defined the parameters and variables for the reservoir model, we can define the model itself. In PyBaMM a model can be defined using the [`pybamm.BaseModel`](https://docs.pybamm.org/en/stable/source/api/models/base_models/base_model.html) class. PyBaMM models can be composed by differential and/or algebraic equations. Differential equations have the form $\\frac{\\mathrm{d} x}{\\mathrm{d} t} = f(x, y, t)$, where $x$ and $y$ are state variables, $t$ is time, and $f(x, y, t)$ is an arbitrary function that will refer as the right-hand-side (rhs) of the differential equation. Algebraic equations have the form $g(x, y, t) = 0$. For example, let's consider the following arbitrary system:\n", + "\n", + "$$\n", + "\\begin{align*}\n", + "\\frac{\\mathrm{d} x}{\\mathrm{d} t} &= f(x, y, t), \\\\\n", + "g(x, y, t) &= 0, \\\\\n", + "x(0) &= x_{0},\n", + "\\end{align*}\n", + "$$\n", + "with an output variable (i.e. a variable that is computed directly from the state and indepedent variables) $z = h(x, y, t)$. \n", + "\n", + "\n", + "We will use this example to illustrate the main attributes of the `BaseModel` class. These four main attributes are:\n", + "1. `rhs` - a Python dictionary of the right-hand-side term for the differential equations, where the key is the state variable we are solving for (e.g. $x$) and the value is the rhs function (e.g. $f(x, y, t)$). For our example, we would write `{x: f(x, y, t)}`.\n", + "2. `algebraic` - a Python dictionary of the algebraic equations. The key is still a state variable (just for indexing purposes, e.g $y$) and the value is the function that is equal to zero (e.g. $g(x, y, t)$). For our example, we would write `{y: g(x, y, t)}`. It is worth reiterating that this imposes $g(x, y, t) = 0$, not $y = g(x, y, t)$.\n", + "3. `initial_conditions` - a Python dictionary of the initial conditions. Again, the key is the variable we are imposing an initial condition for and the value is the initial value. Even though mathematically initial conditions are only needed for differential equations, PyBaMM also requires them for algebraic equations and the value is used as an initial guess for the solver. For our example, we would write `{x: x_0, y: y_0}` (where `y_0` would be a reasonable value for y at the initial state).\n", + "4. `variables` - a Python dictionary of the output variables. The key is a string with the variable name (as a string) and the value is the expression for that variable. For our example, we would write `{\"z\": h(x, y, t)}`.\n", + "\n", + "Now we have all the pieces we need to define the reservoir model. We start by initialising the model:" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [] - }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ - "model = pybamm.BaseModel()" + "model = pybamm.BaseModel(\"reservoir model\")" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Next, we define the variables in the model using the `Variable` class. In more complicated models we can give the variables more informative string names, but here we simply name the variables \"x\" and \"y\"" + "Note that we can pass the name of the model as a string when initialising it. This will be used by PyBaMM's plotting functionality and comes handy when plotting various models together. Next we define the rhs and initial conditions for our model:" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ - "x = pybamm.Variable(\"x\")\n", - "y = pybamm.Variable(\"y\")" + "model.rhs[x_n] = -i / (Q_n * 3600)\n", + "model.initial_conditions[x_n] = x_n_0\n", + "model.rhs[x_p] = i / (Q_p * 3600)\n", + "model.initial_conditions[x_p] = x_p_0" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "We can now use the symbols we have created for our variables to write out our governing equations. Note that the governing equations must be provied in the explicit form `d/dt = rhs` since pybamm only stores the right hand side (rhs) and assumes that the left hand side is the time derivative. " + "Note the factor of 3600 in the equation to convert the Ah to C.\n", + "\n", + "Finally, we define the output variables for our model. Note that these can be any variables we might need to use (e.g. to analyse, plot...) after solving the model, including the state variables. For the reservoir model we define the stoichiometry of each electrode and the voltage:" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ - "dxdt = 4 * x - 2 * y\n", - "dydt = 3 * x - y" + "model.variables = {\n", + " \"Negative electrode stoichiometry\": x_n,\n", + " \"Positive electrode stoichiometry\": x_p,\n", + " \"Voltage [V]\": U_p - U_n - i * R,\n", + "}" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "The governing equations must then be added to the dictionary `model.rhs`. The dictionary stores key and item pairs, where the key is the variable which is governed by the equation stored in the corresponding item. Note that the keys are the symbols that represent the variables and are not the variable names (e.g. the key is `x`, not the string \"x\")." + "The naming convention of the variables is up to the user, but we suggest sticking to the PyBaMM convention of the name followed by the units in square brackets. The only constraint is that state variables need to be given the same name as when they were defined." ] }, { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "model.rhs = {x: dxdt, y: dydt}" - ] - }, - { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "The initial conditions are also stored in a dictionary, `model.initial_conditions`, which again uses the variable as the key" + "### PyBaMM expressions\n", + "\n", + "It is worth pausing here and discussing the concept of an \"expression\" in PyBaMM. Notice that in the model definition we have used expressions like `-i / Q_n` and `U_p - U_n - i * R`. These expressions do not evaluate to a single value like similar expressions involving Python variables would. Instead, they are symbolic expressions that represent the mathematical equation.\n", + "\n", + "The fundamental building blocks of a PyBaMM model are the expressions, either those for the RHS rate equations, the algebraic equations of the expression for the output variables such as $V(t)$. PyBaMM encodes each of these expressions using an \"expression tree\", which is a tree of operations (e.g. addition, multiplication, etc.), variables (e.g. $x_\\mathrm{n}$, $x_\\mathrm{p}$, $I(t)$, etc.), and functions (e.g. $\\exp$, $\\sin$, etc.), which can be evaluated at any point in time.\n", + "\n", + "As an example, lets consider the expression for $\\frac{\\mathrm{d} x_\\mathrm{n}}{\\mathrm{d} t}$, given by `-i / Q_\\mathrm{n}`. We can visualise the expression tree for this expression using the `render` method:" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 8, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/\n", + "├── -\n", + "│ └── Current function [A]\n", + "│ └── time\n", + "└── *\n", + " ├── 3600.0\n", + " └── Negative electrode capacity [A.h]\n" + ] + } + ], "source": [ - "model.initial_conditions = {x: pybamm.Scalar(1), y: pybamm.Scalar(2)}" + "model.rhs[x_n].render()" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Finally, we can add any variables of interest to our model. Note that these can be things other than the variables that are solved for. For example, we may want to store the variable defined by $z=x+4y$ as a model output. Variables are added to the model using the `model.variables` dictionary as follows:" + "or print it to a .png file (which renders it more nicely) by typing: `model.rhs[x_n].visualise(\"x_n_rhs.png\")`. You can also print the expression as a string using the `print` method:" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 9, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-Current function [A] / (3600.0 * Negative electrode capacity [A.h])\n" + ] + } + ], "source": [ - "model.variables = {\"x\": x, \"y\": y, \"z\": x + 4 * y}" + "print(model.rhs[x_n])" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Note that the keys of this dictionary are strings (i.e. the names of the variables). The string names can be different from the variable symbol, and should in general be something informative. The model is now completely defined and is ready to be discretised and solved!" + "In any case, you can see that the expression tree is a symbolic representation of the mathematical equation, which can then be later on used by the PyBaMM solvers to solve the model equations over time." ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "## Using the model" + "### Events\n", + "\n", + "An event is a condition that can be used to stop the solver, for example when the voltage reaches a certain value. Even though they are not a fundamental part of the model (and they are definitely optional), events are very useful to prevent the model reaching unfeasible regions (e.g. negative concentrations). In PyBaMM an event can be defined using the [`pybamm.Event`](https://docs.pybamm.org/en/stable/source/api/events/event.html) class. The class takes two required arguments: a string with the name of the event (for identification purposes) and an expression that is the condition to trigger the event. For example, let's consider the following event for the maximum stoichiometry in the negative electrode:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "max_stoichiometry = pybamm.Event(\"Maximum negative stochiometry\", 1 - x_n)" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "We first discretise the model using the `pybamm.Discretisation` class. Calling the method `process_model` turns the model variables into a `pybamm.StateVector` object that can be passed to a solver. Since the model is a system of ODEs we do not need to provide a mesh or worry about any spatial dependence, so we can use the default discretisation. Details on how to provide a mesh will be covered in the following notebook." + "The expression `1 - x_n` is the condition that the event is looking for, and the solver will stop when this expression becomes zero. Note that the\n", + "expression must evaluate to a non-negative number for the duration of the simulation, and then become zero when the event is triggered.\n", + "\n", + "PyBaMM models have an `events` attribute, which is of all the events of the model. When one of these events is triggered the solver will stop. For example, in the reservoir model we want the stoichiometries to be between 0 and 1, so we define:" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ - "disc = pybamm.Discretisation() # use the default discretisation\n", - "disc.process_model(model);" + "model.events = [\n", + " pybamm.Event(\"Minimum negative stochiometry\", x_n),\n", + " pybamm.Event(\"Maximum negative stochiometry\", 1 - x_n),\n", + " pybamm.Event(\"Minimum positive stochiometry\", x_p),\n", + " pybamm.Event(\"Maximum positive stochiometry\", 1 - x_p),\n", + "]" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "Now that the model has been discretised it is ready to be solved. Here we choose the ODE solver `pybamm.ScipySolver` and solve, returning the solution at 20 time points in the interval $t \\in [0, 1]$" + "## Using the model\n", + "\n", + "Now that we have defined the model, we can proceed to solve it. However, before actually solving it, we need to define the parameter values we will use. \n", + "\n", + "### Parameter values\n", + "You should already be familiar with the [`pybamm.ParameterValues`](https://docs.pybamm.org/en/stable/source/api/parameters/parameter_values.html)\n", + "class, which is used to define the values of the parameters in the model (see [Tutorial 4 - Setting parameter values](../getting_started/tutorial-4-setting-parameter-values.ipynb) for a full recap). Recall that the parameter values are passed to the `pybamm.ParameterValues` class as a dictionary, where the key is the parameter name (i.e. the string we used when defining the parameter) and the value is the parameter value. For parameter functions, we can define them either as a normal or a lambda function. For the reservoir model, we define the following parameter values:" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ - "solver = pybamm.ScipySolver()\n", - "t = np.linspace(0, 1, 20)\n", - "solution = solver.solve(model, t)" + "def graphite_LGM50_ocp_Chen2020(sto):\n", + " u_eq = (\n", + " 1.9793 * np.exp(-39.3631 * sto)\n", + " + 0.2482\n", + " - 0.0909 * np.tanh(29.8538 * (sto - 0.1234))\n", + " - 0.04478 * np.tanh(14.9159 * (sto - 0.2769))\n", + " - 0.0205 * np.tanh(30.4444 * (sto - 0.6103))\n", + " )\n", + "\n", + " return u_eq\n", + "\n", + "def nmc_LGM50_ocp_Chen2020(sto):\n", + " u_eq = (\n", + " -0.8090 * sto\n", + " + 4.4875\n", + " - 0.0428 * np.tanh(18.5138 * (sto - 0.5542))\n", + " - 17.7326 * np.tanh(15.7890 * (sto - 0.3117))\n", + " + 17.5842 * np.tanh(15.9308 * (sto - 0.3120))\n", + " )\n", + " \n", + " return u_eq\n", + "\n", + "parameter_values = pybamm.ParameterValues({\n", + " \"Current function [A]\": lambda t: 1 + 0.5 * pybamm.sin(t / 100),\n", + " \"Initial negative electrode stochiometry\": 0.9,\n", + " \"Initial positive electrode stochiometry\": 0.3,\n", + " \"Negative electrode capacity [A.h]\": 1.2,\n", + " \"Positive electrode capacity [A.h]\": 1,\n", + " \"Electrode resistance [Ohm]\": 0.1,\n", + " \"Positive electrode OCV\": nmc_LGM50_ocp_Chen2020,\n", + " \"Negative electrode OCV\": graphite_LGM50_ocp_Chen2020,\n", + "})" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "After solving, we can extract the variables from the solution. These are automatically post-processed so that the solutions can be called at any time $t$ using interpolation. The times at which the model was solved at are stored in `solution.t` and the statevectors at those times are stored in `solution.y`" + "### Solving the model\n", + "\n", + "Finally we can proceed to solve the model. For this simple case, we can use the `pybamm.Simulation` class directly. We define the simulation by passing the model and the parameter values and then solve it for a given time interval (in this case 1 second):" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ - "t_sol, y_sol = solution.t, solution.y # get solution times and states\n", - "x = solution[\"x\"] # extract and process x from the solution\n", - "y = solution[\"y\"] # extract and process y from the solution" + "sim = pybamm.Simulation(model, parameter_values=parameter_values)\n", + "sol = sim.solve([0, 3600])" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "We then plot the numerical solution against the exact solution." + "We can plot the solution by using the PyBaMM's standard plotting capabilities. In this case, as it is a custom model, we need to specify which variables we want to plot:" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 14, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABYyUlEQVR4nO3dd3hUVf7H8fdJMmmkQQotQOg9SEdFpFiwoVhYbFgQrGvb/a11ERFdd11FEXtDxQJ27CAWFGkB6UhvoYYWAqTn/P6YwFImkDKZO5l8Xs8zz8zce3PvZ644Z763nGOstYiIiIiIiIhUtiCnA4iIiIiIiEj1oAJUREREREREfEIFqIiIiIiIiPiEClARERERERHxCRWgIiIiIiIi4hMhTgc4kYSEBJuSkuJ0DBERkaPMmzdvp7U20ekcaidFRMRfldRW+nUBmpKSQlpamtMxREREjmKM2eB0BlA7KSIi/quktlKX4IqIiIiIiIhPqAAVERERERERn1ABKiIiIiIiIj7h1/eAioj4g/z8fNLT08nJyXE6ivhYeHg4ycnJuFwup6OIiPg1tZXVV1nbShWgIiInkZ6eTnR0NCkpKRhjnI4jPmKtZdeuXaSnp9O4cWOn44iI+DW1ldVTedpKXYIrInISOTk5xMfHq0GtZowxxMfH62i+iEgpqK2snsrTVnqlADXGvGmM2WGMWVLCfGOMGWuMWW2MWWSM6eSN7ZbKokkwph2MjHM/L5rkX+sTkSpBDWr1pP/uIj6k31hVnr4zq6ey/nf31hnQ8UD/E8w/D2he/BgOvOSl7Z7Yoknw5Z2QuQmw7ucv7yz/F5q313fkevWFKyIiIlWJN3+/VNZvLBHxO165B9RaO90Yk3KCRS4G3rHWWmCWMSbOGFPXWrvVG9sv0bRRkJ999LT8bHZPfph/rWhJkDEEBbmr9iCD+70xGAOG4mlBBlewwRUcxI1z/0mMh/Ud/PYRpgf1IjTEvdyhR4QrmIjQYPezK5jw0CBCg4OOPkpw6Av30HoPfeECpA4q3+deNMn92TPTITYZ+o0o/7pExHGbNm1iyJAhbN++HWMMw4cP56677jo8f/z48fTu3ZtGjRp5/ejzs88+S61atRgyZAjjx4/nnHPOoV69egAMHjyYxx57jObNm3t1m6VlraVfv358/vnnZGZmlmkf+dtnEal0XvxtUFBYRM78D4n8/l6CCv73+6Xgi78yc1UGfyaeR25BIbkFRe5Hvvt1QZE9vA5rj17nQ6seppaH31j7vxnBVzk9iIlwERPuIiYihNji19HhIYQEn+Bcin4PVSs5OTn06tWL3NxcCgoKuPzyy3n00Uex1mKMYeTIkYwcOfLwe2+6++67ufTSS+nVqxfPPvssw4cPJzIyEoCzzjqLjz76iJo1a3p1m6WVnZ1N//79+fHHH1m8eDG33nor+/btIzg4mIceeoi//OUvJe6jyvosxh77DVDeFbkL0K+ste08zPsKeNJa+1vx+2nAfdbaNA/LDsd9lpSGDRt23rBhQ/lDjYwDjv98RRh6hn1CkQWLdT9b93ORtRQVWeyh1xYKiorIL7SsDbuKIA//XousoUnue6WKFGQoLkxDiAgN4uPs4dS2Gccttze0Dm93+4ro8BCiwkOICQ8hOtxFVFgI0cWvo8NDCAs5SUEL4IqAi8bqS1eknJYvX07r1q0d2/7WrVvZunUrnTp1Iisri86dO/P5558TGxvLI488QqNGjWjcuDG//vorr7zyite2W1BQQKdOnZg/fz4hISH07t2b//73v3Tp0gWAX375hQkTJvDaa695bZtl8fXXX/PDDz8wZsyYMu+jsnwWT//9jTHzrLVdfPJBT6BLly42Le24plTkaKX4bZBbUMi2zBy2Zuawc38uO7Ny2bk/j537c8nIynVP25/H3oN5HMgr5LfQO0kO2nncptKLEuiZNxYAYyA8JJgwl/sAvOsExeKvOZcS5Ok320l+Y0WHh5AUHUZSdDhJMWGHX3fMnErHBY8QXKjfQ77idFtpreXAgQNERUWRn59Pz549ee655wgLC+Ott94CoG/fvsyZM4cnnnjCa9vdtWsXF1xwAbNmzQIgJSWFtLQ0EhISAHj77bdJT0/noYce8to2y+KFF16goKCAu+66i5UrV2KMoXnz5mzZsoXOnTuzfPly1q1b53EfleWzlKWt9LtecK21rwKvgrthrdDKYpOLL+U4WlBsMr/f06+sueDZZPdRtGMURtfjm5vPIL+wiPzCIvIKi8grKCInv4ic/EIO5hWSnV9ITn4h2cWvs/MLyckrJGn58V/eADG52xnzw8qT5nIFG2IjQqkZ6aJmZCgv7XyI+ILjjyDmfj+STXXOp1aNUOIiXAR5qqRLoiOIIo6qW7cudevWBSA6OprWrVuzefNm2rRpw+OPP0737t1p164dkydPBmDNmjXcfvvtZGRkEBkZyWuvvUazZs049dRTeeqpp+jduzcPPPAAQUFBPP7446SkpDBo0CC+/fZbIiIieP/992nWrBk//vgjnTp1IiQkhI8//pi0tDSuvvpqIiIimDlzJmeccQbXX389BQUFhISUvjmZMGECY8eOJS8vj+7du/Piiy8yf/58hg4dypw5cygsLKRbt25MnDiRnTt3MmLECKKjo1m9ejV9+vThxRdfJCgoiPfee4/hw4eXeR9587OIVAklXBG2a/LDDP21Ppv3ZpORlXvcnwUHGWrVCCUhKozE6DCaJkVRMzKUmHAX9X/b5XFT9YN2sXjkOYSFBOMKNqU/0zTG8282Yuvz+4192ZeTz77sAvZl55OZnc++HPfzngN57MjKZUdWLvM37mHHvlxyC4r4LfQ/BAcd/5n3fvkwn+7rSvPaUTRPiqZ2TJjuWwwQxhiioqIA95Aw+fn5GGPo2LEjERERnHrqqeTn5/PSS+47AefNm8e9997L/v37SUhIYPz48URGRtKtWzcmT55My5YtufLKK+nbty/Dhg0jKiqKYcOGMWXKFOrUqcOHH35IYmIin3zyCf37u+9EHDt2LFu2bKFPnz4kJCTw008/MWDAAM4444wyF6BPPfUUkyZNIjc3l4EDB/Loo4/y2WefMW7cOH744Qe2bdvGmWeeyfTp0/nuu+/47LPPyMzMZPPmzVxzzTU88sgjALz33nu8//77ALRo0eLw+uvVq0dSUhIZGRke95E3P8uxfNXKbgYaHPE+uXha5eo3wvMRv34jyrwqYwz0e8Tj+lznjKRNvZjyZSzhCzcoLplVd57HgdwCsnIOPfLJyilgf27x6+J5ew+6v4D3HMyjZsHxZ1MBXPu3cNYzvwDuBiUhKpTE6DASixuV/70O/9/76DBqrPgU4+1LhEWqsEe/XMqyLfu8us429WJ45KK2pVp2/fr1/PHHH3Tv3p0tW7bwyCOPcOONN9K4cWNuv/12XnrpJYYPH87LL79M8+bNmT17Nrfddhs//vgj48eP5/LLL+f555/nu+++Y/bs2YfXGxsby+LFi3nnnXe4++67+eqrr5gxYwadO3cG4PLLL2fcuHFHnTUEaNasGQsXLjy83MksX76ciRMnMmPGDFwuF7fddhvvvfceQ4YMYcCAATz88MNkZ2dzzTXX0K5dO37++WfmzJnDsmXLaNSoEf379+fTTz/l8ssvZ8aMGR7P+JZmH3njs4j4m7yCIjbuPsj6nQdYt/MA63YdYP3OA0zITPfY6UfNgh1EhYXQp2Ui9eMiqRcXTt3YCBKiQ0mMCqNmZGjJB6wXe/79YmKTiQ4vx7i5JfxmCzrrEerFRVCPiFKtxlpLVm4B0U96LpBj8nYw6qtlh99Hh4fQonY0zZOiaFE7mnb1Y2lbL4YaYR5+IuuAfKk51VYWFhbSuXNnVq9eze2330737t1ZsGABb731Ftdeey19+/bl4Ycf5pFHHuGvf/0rX3zxBYmJiUycOJGHHnqIN998k3HjxnH99ddz1113sWfPHoYNGwbAgQMH6NKlC2PGjGHUqFE8+uijjBs3jhkzZnD55ZcDcOedd/LMM8/w008/HT5rWLNmTXJzc9m1axfx8fGl+qxTpkxh1apVzJkzB2stAwYMYPr06QwcOJBPPvmEF154ge+++45HH32UOnXqADBnzhyWLFlCZGQkXbt25YILLiA1NZW1a9eSkpJy3DbmzJlDXl4eTZs29biPRo8e7ZXP4omvCtDJwB3GmA+B7kBmpd//Cf/7UvDWl4W31wcnLJJdwUHERYYSFxla+vWVUNDm1ajLcwNPYfeBPHbtzyMjK5eM4stqlm/NYuf+3KPuzThkRtgD1DfHH0HM+e4R0uucR93YCM9f0iLidfv37+eyyy7j2WefJSYmhpiYGF577TXGjx/PGWecwTXXXMP+/fv5/fffueKKKw7/XW6u+8xG27Ztufbaa7nwwguZOXMmoaH/+2658sorDz/fc889gPvS35NdTpWUlHT4Mp7SmDZtGvPmzaNr166A+96UpKQkAEaMGEHXrl0JDw9n7Nixh/+mW7duNGnS5HC+3377jcsvv5zdu3cTHR1d5n3krc8iUqlOUOwUFlk27DrAyu1ZrNi2n5Xbs/hz2z7W7zpI4RFteVyki8YJNdgbkkStgu3HbSIoNpkJN3UvXz4vHuQHvPYbyxhDTLirxKvgTFwyacPOYuX2LFbvcO+7ldv38/3SbXw4d1PxOqBZYhTtk2NpXz+W1ORY2u+eQug3d+uAvJ8LDg5mwYIF7N27l4EDB7JkyRI6dOjAc889x8iRI7nkkku4+OKLWbp0KUuWLOHss88G3IXroatozj77bD766CNuv/12Fi5ceHjdQUFB/OUvfwHgmmuu4dJLLwXcbWViYuIJcx1qX8pSgE6ZMoWOHTsC7rZt1apV9OrVi+eff5527drRo0ePw233odyH1n/ppZfy22+/Ua9ePeLi4o5b/9atW7n22mt5++23CQoK8riPvPVZPPFK5WCM+QDoDSQYY9KBRwAXgLX2ZeAb4HxgNXAQuMEb2y2V1EHe/WKojPWB94raEhqE8P6PcnFq/RL/rKjIsjc7nx1ZOe7itPhR7yfPRxBDD2zlrGemAxATHkK9uAjqxobTsFYkDWpF0rBWJI3ia9CgVgSRocf8M9MRRKnCSnum0tvy8/O57LLLuPrqqw83eodcf/31h18XFRURFxfHggULPK5n8eLFxMXFsWPHjqOmH3kJ2qHXERERJx3XKycnh4iIo89MzJ49m5tvvhmAUaNGMWDAgMPzrLVcd911/Otf/zpuXbt27WL//v3k5+eTk5NDjRo1jst25PuQkBCKiooICnKf2yntPirLZ/EWY0wD4B2gNu7OCV611j53zDJXA/cBBsgCbrXWLjx2XVINeOigMP/zvzJx1gY+zO3Bqu37yS0oAtzFUqNakbSoHc157erSJLEGjRPcj8MHsBeN9m6xCJVzUN6bv7FK+D1k+o0gISqMhKgwTmuacHiWtZaMrFyWbMlkUXomi9Mz+XXVTj6d775gb0bYQx4PyDNtlH7DeOBUW3lIXFwcffr04bvvvqNdO3f3NCNHjgTcbYi1lrZt2zJz5szj/raoqIjly5cTGRnJnj17SE5O9riNiraVn332GY8++igAr7/++lFX5FhreeCBBw63pUdKT08nKCiI7du3H9UGemorPWXbt28fF1xwAY8//jg9evQ46m+P3Edl+Sxl5a1ecK88yXwL3O6NbQUkb37hlrNBCCq+z6NWjVBa1TlixnzPRxDzo+rx3MBT2LI3h62Z2WzZm8OWvdmkrd9DVm7BUcsmRofRsLgoPbvgF85Z+wQhhcX/M+gIoshJWWsZOnQorVu35t577z3hsjExMTRu3JiPPvqIK664AmstixYtokOHDnz66afs3r2b6dOnc+GFFzJnzpzDR0YnTpzI/fffz8SJEzn11FMBaN26NatXrz687ujoaLKyso7a3sqVKw837occuuTJk379+nHxxRdzzz33kJSUxO7du8nKyqJRo0bcfPPNPPbYY6xbt4777ruPcePGAe7LhNatW0ejRo2YOHHi4fs+W7Zsydq1a2nWrFmZ9lFZPosXFQB/s9bON8ZEA/OMMVOttcuOWGYdcKa1do8x5jzc/SGU8/SUVEV7D+axZPM+Onw9guhj7tl0FeXQb8srTGl4JkNOjadF7Wha1ommWVLU8Qd6j1UZxeKh9fpr213Gz2yMISkmnL4x4fRtVRtwf/du35fLwvS91PvI8wF5m5nOb6sy6JpSi3BXcKV8FCmdjIwMXC4XcXFxZGdnM3XqVO677z6Py7Zs2ZKMjAxmzpx5+L7HlStX0rZtW8aMGUPr1q154oknuOGGG5g5cyYul4uioiI+/vhjBg8ezPvvv0/Pnj2B/7WVvXv3Bv7Xvhy6bNVay7Zt2467DHbgwIEMHDjQY75zzz2Xf/7zn1x99dVERUWxefNmXC4XtWrV4sYbb+SDDz7g7bff5plnnuHvf/87AFOnTmX37t1ERETw+eef8+abb1KzZk0KCwvJyckhPDycvLw8Bg4cyJAhQw5fNnwipf0sZaVrJwORD44ghp070uMZVWstew/ms3H3QTbsPsim3QfZuOsgG3YfYM663fwteywh5pijRPnZ7PnyYT7a25mmiVE0SYyiQc2IE3etLlKNzJgxg3fffZf27dtzyimnAPDEE09w/vnne1z+vffe49Zbb2X06NHk5+czePBg6tevz/3338+0adNo0KABd9xxB3fddRdvv/02AHv27CE1NZWwsDA++OADAM477zyuvfbaw+u9/vrrueWWWw533LNv3z4iIiIO339SGm3atGH06NGcc845FBUV4XK5eOGFF/jll19wuVxcddVVFBYWctppp/Hjjz8SFBRE165dueOOOw53QnSowb7gggv4+eefadasWZn3kTc+S1kU33aytfh1ljFmOVAfWHbEMr8f8SezcPeXIAEqJ7+QxZszmb9hDws27WXx5kzS97jb2rVhW93nwY9Rl528c2O38m3Qn4vFylLBz2yMoU5sOHVi68AUzwfkt9h4rn1jDuGuIHo0iadX80TObJlIk4Qa6tzIx7Zu3cp1111HYWEhRUVFDBo0iAsvvNDjsqGhoXz88cfceeedZGZmUlBQwN13301ISAivv/46c+bMITo6ml69ejF69GgeffRRatSowZw5cxg9ejRJSUlMnDgRcLdFr7zyCjfddBMAw4cPp3///tSrV4+ffvqJefPm0aNHjzJ1cHfOOeewfPnywweEo6KimDBhAi+//DJnnHEGPXv2pEOHDofv9QT37SqXXXYZ6enpXHPNNYfPqJ5zzjn89ttvnHXWWUyaNInp06eza9cuxo8fD7iHKTvUbh7LG5/FE68Nw1IZ1L28n/DiJbN2ZBymhKFxmuT8r5t1V7AhJb4GTROjaFE7ipZ1YmhVN5qU+BoEe+oQQZf1SiVyumv5ynZsN+tHGjhwIP/5z388jo85ZswYYmJiGDp0aKVl+/nnn/nvf//LV199ddy8rVu3MmTIEKZOnVrh7Zzos3h7GJbiYcumA+2stR576TDG/B1oZa29ycM87w1XJt5zgnbIWsuWzBzmb9jD/I17mL9xL8u2ZJJf6G4PG9SKILV+HO3qu+85PPXLMwned3yv+8Q2gHuW+PJTySElDGWTe96zzIjsw/SVO5m+MoO1Ow8A0LBWJOe2rc25bevQqWHNso0+UEUFelsZFRXF/v37Pc7r2bMnX331lcf7Le+66y4GDBhAv35lG4GjLMaPH09aWtrhK4eONH/+fMaMGcO7775b4e2c6LNU6WFYxA958aipOcHQOAv/cQ5rdu5nzY79rMk4wJoMd+cAU5Zt41CfCmEhQYcvO2pVJ5rWdWNov3sKMVPvVccAIpXgySefZOvWrR4L0Li4uKPOkPpa3bp1GTZsGPv27SMmppw9kRfz1WcxxkQBnwB3n6D47AMMBXp6mu/V4crEOzzcs1n0xZ38/OcOPs4/lXkb9rB9n7szsHBXEKnJcQzt2YRODePo2LAmidFhR6/vLM+97lfonk2pmBIu6Q1LHURfOHzZ7qbdB/l5ZQbTlm9n/O/ree3XdSREhXFOcTF6apN4QkN0hVegefrpp9m4caPHArRdu3aVWnyeTKdOnejTpw+FhYUEB1fsMnFvfRadARXfKsVg2MfKyS9k9Y79/Lktiz+37mPF9qzDPfcCJQ6GbWMbYHSkWLwg0I/qyol56wyoMcYFfAV8b619poRlUoHPgPOstScdDFrtpH+wY9piPIwTnl6UwOAar9GpYU06NYyjc6NatKobjas0t5joyp4qb19OPj/9uYMpS7fz04odHMwrpGakiwtS63LJKfXp3KgmZvFHAfPfWW1l9aYzoOK/ytERQrgrmHb1Y2lXP/ao6Tv357JiWxb1J5TcMcB1b84htfhv2yfHUi82XPdkiIjPGfcXzxvA8hMUnw2BT4FrS1N8inPyC4tYlJ7JrLW7mL1uN+P3bsZT01I/aBe/3de3fBupjvdsBpiYcBcXn1Kfi0+pT05+Ib+u2snkhVv4eF46E2Zt5IboOTxY+DKuInXMKNWLClDxPS81qglRYSQ0CytxrK9MVxI79uXw0uqdh8dFi68RSseGcXRqVJNODWvSITmOiFAPlyPoyLOIeNfpwLXAYmPMguJpDwIN4fCQZSOAeODF4gNlBeW9z1S8y1rLyu37+XVVBr+t3smcdbs5mFcIQPOkKDJDk6iZf/w4myZW/UiJW7grmLPb1ObsNrXZn1vA90u2ceY3d/+v+DxEQ7tINaACVKq+EnrqrXnRaL5L7UVOfiHLt+5j8eZMFm7K5I+Ne/hhuXsMxJAgQ+u6MXQ6oihNTv8Kc8y9PDoiKSIVYa39DY/9mh61zE3AcZ0OiTN27Mvht9U7+W3VTn5bvZMdWe7bPpok1OCyTsmc1jSebo1rER8VVjnjbErAigoL4bLOyfBlhsf5NjOdFdv20apOxe5tF/FXKkCl6jvJZb3hrmA6NqxJx4Y1wd2bNbsP5PHHxuLeCDfsZVJaOm/PdPckOTP8AeqiwaZFRAKWh6tcsltdxqx1u9wF56qdrNjuHie2ZqSL05slcEbzBHo2T6R+nIcB2CtrnE0JbCVcwbXFxtP/2V/p2DCOG05vzHnt6pTuvmGRKkIFqASGMl7WW6tGKP1a16Zfa3evdQWFRfy5LYv5G/dQ5/uS7yldsW0fLZKiq0V36iIiAclDj7W5n93BQwUL+TT/NEJDguiaUpNLOrbijOYJtKkbU7rvfN2zKWVVwhVcsec8xj9z2zBh1gbu/OAP6saGc91pKVzZtSGxkS7n8op4iQpQESAkOOh/HR3N8nxEcnOR+4hkzUgX3RrXokeTeM5onkjTRA02LcfQPcQifimvoIii7x8hPP/oq1zCbC4jIj7mkmvvplvjWoS7KjZUgUiplHDmPCp1EEOBG05L4acVO3jjt3U8+e2fPPfDKq7okswNpzemcUINR6N7hdrKaksFqMixSjgiGdH3Uf4b2oFZa3cxa+0uvl/q7nCiXmw4ZzRP5IwWCfRslkBcZKhDwcUveDi7UtF7iOfOncvQoUOZM2cOhYWFdOvWjYkTJ9KuXTsvhRYJXDv25fDTih389Ke7A6FFbPF4N25c3nZ6tUj0fUCp3k5w5jwoyBy+WmvZln28OWMdH87ZxIRZG7ggtR6392lade8TrYS2csSIEdSqVYu7774bgIceeoikpCTuuusuLwQWb9I4oCKelOKo3KbdB/l11U6mr8xgxpqdZOUUYAykJsfRq3kCvVokckqDOPd9GzrKV6WVaWyzMe08nkEntgFUYFzahx9+mJycHLKzs0lOTuaBBx4o97qkbLw1DmhlUDt5PGstK7ZnMWXpdqYu287izZkA1I0Np0+rJP65ahARB7cc/4cV/H9UxBd2ZOXw1oz1vPP7eg7kFXJOm9rc0bcZqclxTkdzvK1cv349l156KfPnz6eoqIjmzZszZ84c4uPjy7U+KRuNAypSUaW4l6dBrUiu6t6Qq7o3pKCwiIXpmfy6KoPpKzN44afVPP/jaqLDQri3zkKu3fk0IYUa56ta8DAY/Qmnl9KIESPo2rUr4eHhjB07tkLrEgk0BYVFzNuwhynLtjNl2TY27c7GGOjYII5/9G9J31ZJtKwd7b5dYtGj6rFWqqyk6HDu69+Km3s14a0Z63lrxjqmLNvOmS0SuefsFpyyZ0rVOOBdCW1lSkoK8fHx/PHHH2zfvp2OHTuq+PRTKkBFvCAkOIjOjWrSuVFN7j6rBZnZ+cxcs5Of/syg/5LbCOH4cb7stFEYf2wUpGJK6NWQCo4HuGvXLvbv309+fj45OTnUqBEA9/+IlFYJvdb+uiqDKcu2M235dvYczCc0OIjTm8VzW+9m9GudRFJ0+PHrUo+1EgDiIkO55+wW3HRGY96dtYHXf13Hmy/9m6fC3iDMuocM8usD3pXUVt50002MHz+ebdu2ceONN1ZoXVJ5dAmuSCWzI+MwHP//WRGG0Z1/56w2SXRLqUWIulj3W2W6rOjY+1rAfXblorEV+gEwYMAABg8ezLp169i6dSvjxo0r97qkbHQJrsM8/D+Va8J4sHAYn+SdRkx4CH1bJXFO2zr0apFIVJiOrUv1sz+3gKJn2hKTu+34mT66vNwf2sq8vDzat29Pfn4+q1atIjhYHYr5ii7BFfEjpoSjfLuDE5kwewNvzlhHfI1Qzmlbh/Pb16FHk3iN91WVVcLZlXfeeQeXy8VVV11FYWEhp512Gj/++CN9+/b1UmgR/1X4w6MEe+i19p9hH3HpkHvo1riWvjOl2osKC4Hc7R7n2cx0T/1uOauSrkQIDQ2lT58+xMXFqfj0YypARSpbCb3qJlz0OAtanc0vKzL4Zsk2Ji/YzAdzNhIX6eKcNrU5r31dTm+aQGiIflhVOV4eD3DIkCEMGTIEgODgYGbPnu21dYv4o537c/luyTa+WbyVCZnpnnutzd/B6c0SfB9OxF+VcMB7K/FMm7WBK7s28K+rrSph7NyioiJmzZrFRx995NX1inepABWpbCc4yhcJnNe+Lue1r0tOfiHTV2bw7ZJtfLt4G5PS0okOD+Hs1rW5qEM9ejZP0FF+EQlYO7Jy+H7JNr5evJU563ZTZKFJQg2ywuoQm+fpssKK3SsmEnA8HPAuCongk6gbefrzJUyYuYF/XtiGns0D88DNsmXLuPDCCxk4cCDNmzd3Oo6cgApQEV8oxVG+cFcw57Stwzlt65BbUMiM1Tv5dvE2pizbzqd/bCa+RigXpNblko716dggDrP4I3WiISJV2o59OXxbfKZzzvrdWAtNE2twR59mnJ9a191z7eLH1GutSGl4OOAd1G8Ed7S/guZLt/H4N8u55o3ZnNW6No9c1IYGtSKdzetlbdq0Ye3atU7HkFJQASrih8JCgunbqjZ9W9Xm8YIiflmZwecLNjNx7ibembmBG2Pm8kDBS7iKNLSLr1hr3UM4SLXizx31VVV7DuTxzZKtTF6w5XDR2aJ2FHf2bc4FqXVpUTv66D9Qr7UipefhgLcB+rerS++WSbw1Yz3P/7iKs8f8wp39mjPsjCZevbpKbWX1VNa2Ur3gilQhWTn5fLdkG72/7Uti4Y7jF9BA6pVi3bp1REdHEx8fr4a1GrHWsmvXLrKysmjcuPFR89QL7kkcM2xKzpkP8505g8kLtzB9ZQYFRZamiTW4qEM9LkytS7Ok6JOvU0S8YsvebEZOXsqUZdtpUTuKJwa2p0tKrQqvV21l9VSetlIFqEhVNDIOPAztYjH8fOVKejVPJDhIX/7ekp+fT3p6Ojk5OSdfWAJKeHg4ycnJuFyuo6arAD0BD8MrZNtQ7su/ibTos7jolHoM6FCPNnVj9CNVxEFTl23nkS+WsCUzh8FdG3D/ea2Iiwwt9/rUVlZfZW0rdQmuSFV0gp7ubnhrLnViwrm8czJXdEmmUXwNBwIGFpfLddxRPRE5XmGRJf/7Rwg/ZtiUCJPHUzW/wPW3JwjSwTERv3B2m9qc1jSe56at4o3f1jHtzx08MbA9Z7epXa71qa2U0lKXmiJVUb8R7k44juSKIOmSJ3jp6k60qhvNiz+v5synfmbwqzP57I90cvILnckqIgHNWssfG/fw6JdL6fGvaYTu3+JxubADW1R8iviZGmEhPHh+a764/XTia4Qy7J007pm4gL0H85yOJgFMZ0BFqqISOuUISR3EebiHdtmamc0n89KZlJbOPRMX8thXyxnUpQFXd28YcD3fiYjvbdx1kM/+2Mxnf6SzftdBQoOD6NMqkZwt9YjM9lCEatgUEb/Vrn4sk+/oybgfV/HCz2uYsXonTwxsz1nlPBsqciK6B1QkwBUVWX5fs4sJszYwdfl2iqylT8skru3RiDNbJOqMhEg5VNd7QDOz8/lm8VY+nZ/O3PV7ADi1STwDO9bn3HZ1iI1webwHFFcEXDRWPdeKVAFLNmfy948W8ue2LEY3Wc5V+98iaN9m9UAtZaZOiESErZnZfDB7Ix/M3URGVi4NakVwTfdGXNGlAbXWfK5hDkRKqToVoPmFRUxfmcGn8zczdfl28gqKaJpYg0s7JXNJx/rUj4s4/o+O6QVX3yciVUteQRFTJz5Pn5WjiTRHXI6rg0lSBpVagBpj+gPPAcHA69baJ4+Zfz3wFLC5eNI4a+3rJ1uvClCRypFXUMSUZdt4d+YGZq/bzaWu33ky5DVCbe7/FlIjI1KishagxpgGwDtAbdxdWL9qrX3umGUM7rb0fOAgcL21dv6J1uuVdtJDsWjbX8GSzfv4ZH46Xy7cwq4DedSqEcqADvUY2LE+qcmx6sFWJNCNaeexw0MN+SalVWm94BpjgoEXgLOBdGCuMWaytXbZMYtOtNbeUdHtiUjFhYYEcWFqPS5MrcfK7Vkkvn4Pofm5Ry+Un+3+UaoCVMQbCoC/WWvnG2OigXnGmKnHtJXnAc2LH92Bl4qfK8+xl8tmbiL/87/y1DfLeXVvF0KDgzirTRIDOyZzZotEQkPUd6FItZGZ7nGyzUxHh5+kIrzRCVE3YLW1di2AMeZD4GLg2AJURPxQi9rRkL/D4zybmU5eQSFhIcE+TiUSWKy1W4Gtxa+zjDHLgfoc3VZeDLxj3ZcmzTLGxBlj6hb/beWYNuroezUBV1EON+VNIGXgDVzQvi6xka4S/lhEAtoJhnxbumx7uYdrEfHGocz6wJH/OtOLpx3rMmPMImPMx8WXInlkjBlujEkzxqRlZGR4IZ6InFQJvVNuLorn9Cd/YtyPq8g8mO/jUCKByRiTAnQEZh8zq1TtqVfbyRLOcCQVZXBV94YqPkWqMw9DvhWFRPBejesY9k4aT377JwWFRQ6Fk6rMV9fSfAmkWGtTganA2yUtaK191VrbxVrbJTEx0UfxRKo5D42MdUVw4IwHaVc/hv9OWcnp//6RJ7/9k4ys3BJWIiInY4yJAj4B7rbW7ivPOrzaTpY0NIqGTBGR1EHuviBiGwAGYhsQNGAsf737Ia7q3pCXf1nDNW/M1u8CKTNvXIK7GTjyjGYy/+tsCABr7a4j3r4O/McL2xURb/EwrqjpN4KWqYMYDyzbso+XflnDq9PX8NaMdfylawOG92pCck2NJypSWsYYF+7i8z1r7aceFjlpe+p1/UZ4HjKl34hK3ayIVBGpg47rCyIceGJgezo1rMlDny3mgrG/8sLVneiaUsuZjFLlVLgXXGNMCLAS6Ie7oZwLXGWtXXrEMofvYTHGDATus9b2ONm61QuuiH9Zt/MAr/yyhk/mp2MtXHxKfW7t3YRmSdFORxPxqXL0gmtwX/2z21p7dwnLXADcgbsX3O7AWGtttxOtt7J6wVXnYyJSGsu37uPWCfPYtCebB89vzY2np6iHbDmssodhOR94FvcwLG9aax83xowC0qy1k40x/wIG4O4FcDdwq7X2z5OtVwWoiH/ampnNa9PX8f6cDeQWFHF+u7rcfVZzmtdWISrVQzkK0J7Ar8Bi4NBNUw8CDQGstS8XF6njgP64h2G5wVp7wkZQ7aSIOG1fTj5/n7SQKcu2c3nnZB4f2E6dFwpQyQVoZVHDKuLfdu3P5c0Z6xg/Yz0H8wsZ0KEed/ZrTtPEKKejiVSqshaglUXtpIj4g6Iiy7PTVjF22io6N6rJK9d2JiEqzOlY4rCS2koN6CUi5RYfFcb/nduKX+/ry829mjJl6XbOfuYX/jZpIRm/v+sexHpknPt50SSn44qIiEglCAoy3Ht2C164qhNLt2Ry8bgZLNtSrn7WpBpQASoiFVarRij3n9eKX+/rw42nN4bFk6jx/b3F44dZ9/OXd6oIFRERCWAXpNbl41tOo7DIctlLv/Pdkm1ORxI/pAJURLwmISqMhy9sw3/iPifS5B09Mz/b3dGJiIiIBKx29WOZfMfptKwTzS0T5vHKL2vw51v+xPdUgIqI1wVneR45wmamczCvwMdpRERExJeSYsL5cHgPLkyty7++/ZORk5dSWKQiVNxUgIqI95UwiP3monh6P/Uz78/eSEFhkcdlREREpOoLdwUzdnBHhvdqwtszN3DrhHnk5Bc6HUv8gApQEfG+fiPcg9kfyRVB3pkP0aBWJA9+tphzn53O90u36bIcERGRABUUZHjw/NY8clEbpi7fzpWvzWL/3PfVSWE1F+J0ABEJQIcGsT9mcPsmqYP4uK9lyrLt/Pu7P7n53Xl0S6nFiIva0K5+rLOZRUREpFLccHpj6saGM3XiOEK+fg3Idc841Ekh/O+3gwQ8jQMqIo4oKCxiYtomnpmykt0H8/hLlwb87ZyWJEZr3DDxfxoHVESk7HKfakPYAQ/9RMQ2gHuW+D6QVCqNAyoifiUkOIiruzfix7/3Zujpjfl4Xjp9//szr05fQ16B7g8VEREJNGEHtniekZnu2yDiKBWgIuKo2AgXD1/Yhu/v6UXXxrV44ps/OffZ6Uxbvl33h4qIiASSEjopLHG6BCQVoCLiF5omRvHm9V1564auGAND307jurfmsm7nAaejiYiIiDd46KTwoA1lcau7HAokTlABKiJ+pU/LJL6/uxcPX9CaPzbs4dxnpzNm6kp31+2LJqnnPBERkaoqdRBcNNZ9zyeGwphkXoy+k0t+rc8XCzyPIS6BR73giojfcQUHcdMZTRjQoR6Pfb2c56atYv/c93mw8CWCC3PcC6nnPBERkaonddDhdjsYuDknn7lvp3H3xAUUFFou66zLcQOdzoCKiN9Kignn+Ss78u7QbgzLn/C/4vOQ/Gz3UC8iIiJSJUWHu3j7xm6c3jSBv3+8kE/mqUOiQKcCVET83hnNE6ltd3qeqZ7zREREqrRwVzCvDemiIrSaUAEqIlWCKaGHvPyoej5OIiIiIt4WEeouQk9rGs/fP17Ip/NVhAYqFaAiUjV46Dkvm1D+sfcSnpm6UmOHioiIVHERocG8PqQrpzWN528fqQgNVCpARaRqOKbnPGIbUHjBc9B+EGOnreLC539lwaa9TqcUERGRCji2CP38D/WOG2jUC66IVB1H9JwHEAWM6QoXdajLQ58t4dIXZ3Dj6Y352zktiQgNdi6niIiIlNuhIvSG8XP420cLqREWwtltajsdS7xEZ0BFpMrr26o2U+7pxVXdG/L6b+s499npzF67y+lYIiIiUk4RocG8fl1X2tWP5fb35zNjdQmdEUqVowJURAJCdLiL0Ze058PhPTAGBr82i8e/XkZOfqHT0UQwxrxpjNlhjFlSwvxYY8yXxpiFxpilxpgbfJ1RRMTfRIWF8PYNXWkcX4Nh76Qxf+MepyOJF6gAFZGA0qNJPN/ceQZXd2/Ia7+u46Lnf2PJ5kynY4mMB/qfYP7twDJrbQegN/C0MSbUB7lERPxaXGQo7w7tRmJ0GNe/OYflW/c5HUkqSAWoiAScGmEhjL6kPW/f2I19Oflc8sIMnp+2ioJC9ZQrzrDWTgd2n2gRINoYY3Df3rwbKPBFNhERf5cUE86Eod2pERbCtW/MYd3OA05HkgpQASoiAevMFolMuftMLkity9NTV3LZyzNZk7Hf6VginowDWgNbgMXAXdZaj0dMjDHDjTFpxpi0jIwMX2YUEXFMg1qRvDu0O0XWcs3rs9mxL8fpSFJOKkBFJKDFRrp4bnBHxl3VkQ27DnDB2F+ZMGsD1lqno4kc6VxgAVAPOAUYZ4yJ8bSgtfZVa20Xa22XxMRE3yUUEXFYs6Qo3r6hG3sO5nH9W3PJysl3OpKUgwpQEakWLkytx5S7e9E1pRYPf76EWyfM50Da+zCmHYyMcz8vmuR0TKm+bgA+tW6rgXVAK4cziYj4nfbJsbx4dSdWbM/itvfmk6/ba6ocFaAiUm0kxYTz9g3deOj81kSs+ISgr+6GzE2AdT9/eaeKUHHKRqAfgDGmNtASWOtoIhERP9W7ZRL/urQ9v67aycQ3nsbqYHKV4pUC1BjT3xizwhiz2hhzv4f5YcaYicXzZxtjUryxXRGRsgoKMgzr1YR/x31OBLlHz8zPhmmjnAkmAc0Y8wEwE2hpjEk3xgw1xtxijLmleJHHgNOMMYuBacB91loNeiciUoJBXRrwUuoaLt38H4wOJlcpIRVdgTEmGHgBOBtIB+YaYyZba5cdsdhQYI+1tpkxZjDwb+AvFd22iEh5he7f4nlGZrpvg0i1YK298iTztwDn+CiOiEhA6L/9VYzJO3rioYPJqYOcCSUn5Y0zoN2A1dbatdbaPOBD4OJjlrkYeLv49cdAv+Ku5kVEnBGb7HFydmRdHwcRERGR8jAlHTTWwWS/5o0CtD6w6Yj36cXTPC5jrS0AMoF4L2xbRKR8+o0AV8RRk3II4769l/Cvb5drzFARERF/V8LB5BKni1/wu06INL6ZiPhE6iC4aCzENgAMxDYg+JKxxHS7ild+WctVGmNMRETEv3k4mJxNGPtOf9ChQFIaFb4HFNgMNDjifXLxNE/LpBtjQoBYYJenlVlrXwVeBejSpYsG6hORypM66Kh7RFzA6FOgS6NaPPDpYs4f+xtjrzyF05omOBZRRERESnCoDZ82CjLTyatRj4f3DWTDvBTe61RIWEiws/nEI2+cAZ0LNDfGNDbGhAKDgcnHLDMZuK749eXAj1ajwIuIn7qkY30m33E6sREhXPP6bF74aTVFRfrKEhER8Tupg+CeJTByL6H/t4y+V9xB2oY9PPDJYlRu+KcKF6DF93TeAXwPLAcmWWuXGmNGGWMGFC/2BhBvjFkN3AscN1SLiIg/aV47msl39OSC1Ho89f0KbnonjcyD+U7HEhERkRO4ILUu957dgk//2MxLv6xxOo544I1LcLHWfgN8c8y0EUe8zgGu8Ma2RER8pUZYCGMHn0LXlJo89tUyLnlxBq9e25nmtaOdjiYiIiIl+GvfZqzesZ//fLeC5knRnN2mttOR5Ah+1wmRiIg/McYw5NQUPhjWg6ycAga++DtTlm5zOpaIiIiUwBjDfy5PJTU5lnsmLmD1jv1OR5IjqAAVESmFLim1+PKvp9MksQbD353Hcz+s0n2hIiIifircFczL13QmLCSIm99NIytHt9H4CxWgIiKlVDc2gkk3n8qlHesz5oeV3PrePPbnFjgdS0RERDyoFxfBuKs6sX7XQf42aaEOHPsJFaAiImUQ7grm6UEd+OeFbfhh+Q4ufXEG63cecDqWiIiIeHBq03gePL81U5Zt58WfVzsdR1ABKiJSZsYYhvZszDs3dmNHVi4XvzCD31fvdDqWiIiIeHDj6Slccko9np66kp9W7HA6TrWnAlREpJxOb5bA5Nt7khQdxpA35zDz85dgTDsYGed+XjTJ6YgiIiLVnjGGf12aSus6Mdz1wR9s2KUrl5ykAlREpAIaxkfyyW2ncW+dBXT4YwRkbgKs+/nLO1WEioiI+IGI0GBeubYzxhhue28+OfmFTkeqtlSAiohUUEy4i1sL3yPS5B09Iz8bpo1yJpSIiIgcpUGtSJ4Z1IGlW/Yx+utlTseptlSAioh4gcnc7HlGZrpvg4iIiEiJ+rWuzc29mjBh1kYmL9zidJxqSQWoiIg3xCZ7nJwXVc/HQURERORE/n5uS7o0qskDnyxiTcZ+p+NUOypARUS8od8IcEUcNSmbMB7OulQ95IqIiPgRV3AQz1/VkTBXMLfrflCfUwEqIuINqYPgorEQ2wAwENuAnP7PsCDubK57a44u8xEREfEjdWMjeGZQB/7clsXIyUudjlOthDgdQEQkYKQOcj+K1QQ+6pDPsHfSuPODP9ixL4ebzmjiXD4RERE5rHfLJG7v05QXflpDjybxXNKxvtORqgWdARURqUSxES7eubEb57Wrw+ivlzP6q2UUFVmnY4mPGWPeNMbsMMYsOcEyvY0xC4wxS40xv/gyn4hIdXXPWS3omlKThz9fwsZdB52OUy2oABURqWThrmDGXdWJ605txOu/reOuiQvILdD9JtXMeKB/STONMXHAi8AAa21b4ArfxBIRqd5CgoMY85dTMAbumvgHBYVFTkcKeCpARUR8IDjIMHJAW+7r34ovF27h+jfnkpWT73Qs8RFr7XRg9wkWuQr41Fq7sXj5HT4JJiIiJNeM5ImB7flj417GTlvldJyApwJURMRHjDHc2rspzwzqwNz1u7nqtdnsPpDndCzxDy2AmsaYn40x84wxQ0pa0Bgz3BiTZoxJy8jI8GFEEZHAdVGHelzeOZmNv4wn96nWMDIOxrSDRZOcjhZwVICKiPjYpZ2SeeXazqzcnsVfXpnJtswcpyOJ80KAzsAFwLnAP40xLTwtaK191VrbxVrbJTEx0ZcZRUQC2uimy3nS9QZhB7YAFjI3wZd3qgj1MhWgIiIO6Ne6NuNv6MaWvdlc8crv6vhA0oHvrbUHrLU7gelAB4cziYhUK+G/jCac3KMn5mfDtFHOBApQKkBFRBxyatN43h/Wg6ycAi5/+XdWbMtyOpI45wugpzEmxBgTCXQHljucSUSkeslML9t0KRcVoCIiDurQII5JN58KwF9encnCTXudDSSVwhjzATATaGmMSTfGDDXG3GKMuQXAWrsc+A5YBMwBXrfWljhki4iIVILY5LJNl3JRASoi4rAWtaP5+JbTiA4P4arXZjFzzS6nI4mXWWuvtNbWtda6rLXJ1to3rLUvW2tfPmKZp6y1bay17ay1zzoYV0Skeuo3AlwRR03KNWEU9R3hUKDApAJURMQPNIyP5ONbTqNeXATXvzWH6SvVu6mIiIhPpQ6Ci8ZCbAPAcCC8Lv+XO5R39nd1OllAUQEqIuInaseEM/HmU2mcUIOb3knj5xUaClJERMSnUgfBPUtg5F4i71vOvuYDefK7P1mbsd/pZAFDBaiIiB+pVSOUD4b1oFliFMPfmcfS7153j0Om8chERER8yhjDvy9LJSwkmL99tJCCwiKnIwUEFaAiIn6mZo1Q3h/WneFxaTSe+YB7HDKNRyYiIuJztWPCGXVxW/7YuJdXpq91Ok5AUAEqIuKH4iJDuTf4QyJN3tEzNB6ZiIiITw3oUI8L2tfl2R9WsnzrPqfjVHkqQEVE/FTQvs2eZ2g8MhEREZ8xxjD6knbERrj4x8eLdCluBVWoADXG1DLGTDXGrCp+rlnCcoXGmAXFj8kV2aaISLWh8chERET8Qs0aoYy6uB2LN2fy2q/rnI5TpVX0DOj9wDRrbXNgWvF7T7KttacUPwZUcJsiItWDh/HIDtpQFre6y6FAIiIi1df57evSv20dxvywUr3iVkBFC9CLgbeLX78NXFLB9YmIyCHHjEdWFJPMSzF3cdlvyfyicUJFRER8btTFbQkPCeK+TxZRVGSdjlMlVbQArW2t3Vr8ehtQu4Tlwo0xacaYWcaYS060QmPM8OJl0zIy9ANLRKq5I8YjC7p3KTfdeh/NkqIY/k4av6/Z6XQ6ERGRaiUpJpx/XtiGuev3MGH2BqfjVEknLUCNMT8YY5Z4eFx85HLWWguUdBigkbW2C3AV8KwxpmlJ27PWvmqt7WKt7ZKYmFiWzyIiEvBiI128O7QbjeIjGTo+jbnrdzsdSUREpFq5vHMyZzRP4N/f/kn6noNOx6lyTlqAWmvPsta28/D4AthujKkLUPy8o4R1bC5+Xgv8DHT02icQEalm4qPCmHBTd+rGhnPDW3NZsGmv05FERESqDWMMTwxsjwUe/GwJ7vNwUloVvQR3MnBd8evrgC+OXcAYU9MYE1b8OgE4HVhWwe2KiFRrSdHhvD+sB7VqhDLkjdks3ZLpdCQREZFqo0GtSP5xbkumr8zg0/klDJsmHlW0AH0SONsYswo4q/g9xpguxpjXi5dpDaQZYxYCPwFPWmtVgIqIVFCd2HDeH9adqLAQrn1jDqt3qEc+ERERXxlyagodG8bx+DfL2Xswz+k4VUaFClBr7S5rbT9rbfPiS3V3F09Ps9beVPz6d2tte2tth+LnN7wRXEREILlmJO8P60GQMQx5Yzab92Y7HUlERKRaCAoyPH5JezKz8/n3d386HafKqOgZUBERcVhKQg3evrErWTkFXPvGbHbtz3U6koiISLXQpl4MN56ewgdzNpGmjgFLRQWoiEgAaFsvltev68LmPdncMH4u+3MLnI4kIiJSLdx9VgvqxYbz0GdLyC8scjqO31MBKiISILo3ieeFqzqxdMs+hr+TRm5BodORREREAl6NsBBGDmjLiu1ZvPHbOqfj+D0VoCIiAeSsNrX5z2Wp/L5mF3d9sIDCInUNLyIiUtnOaVuHs9vU5tkfVrJpt8YGPREVoCIiAeayzsn888I2fLd0Gw99tljjk/kBY8ybxpgdxpglJ1muqzGmwBhzua+yiYiId4wc0JYgYxg5eana3hNQASoiEoCG9mzMHX2a8eHcTfx3ygqn4wiMB/qfaAFjTDDwb2CKLwKJiIh31Y+L4J6zWjDtzx18v3S703H8lgpQEZEA9bdzWnBltwa88NMa3pu9wek41Zq1djpwsu4R/wp8Auyo/EQiIlIZrj89hVZ1onn0y6UczFOHgJ6oABURCVDGGB67uB19Wibyz8+XMG25jsb6K2NMfWAg8JLTWUREpPxcwUE8dkk7tmbm8NLPa5yO45dUgIqIBLCQ4CDGXdWJtvViueP9P1i4aa/TkcSzZ4H7rLUn7b/fGDPcGJNmjEnLyMio/GQiIlImXVNqcckp9Xhl+lo27lKHRMdSASoiEuBqhIXwxvVdiI8KZdKbT1PwdFsYGQdj2sGiSU7HE7cuwIfGmPXA5cCLxphLPC1orX3VWtvFWtslMTHRhxFFRKS07j+vNSFBhse+XuZ0FL+jAlREpBpIig7n49PTedi+TEhWOmAhcxN8eaeKUD9grW1srU2x1qYAHwO3WWs/dzaViIiUV53YcP7atzlTl23nl5W6WuVIKkBFRKqJOnP/QwR5R0/Mz4Zpo5wJVI0YYz4AZgItjTHpxpihxphbjDG3OJ1NREQqx409U2icUINfP3kRO0ZXHx0S4nQAERHxkcz0sk0Xr7HWXlmGZa+vxCgiIuIjYSHBjGu3msYzx2Fyiw8AH7r6CCB1kHPhHKQzoCIi1UVsctmmi4iISIW0Xf4skUZXHx1JBaiISHXRbwS4Io6adNCGMqvx7Q4FEhERCXC6+ug4KkBFRKqL1EFw0ViIbQAYbGwyb9W6h2vnNmLOut1OpxMREQk8uvroOCpARUSqk9RBcM8SGLkXc89Srhn2fzSoGcnN76axYdcBp9OJiIgEFg9XH1lXhHt6NaUCVESkGouNdPHG9V0psjD07TT25eQ7HUlERCRwHHH1kcWQXpTAgg6jqm0HRKACVESk2mucUIOXr+nM+p0HuP29+RQUFjkdSUREJHAUX31UNGIPN9V6izuXNSO3oNDpVI5RASoiIpzaNJ7Rl7Tj11U7eeyrZU7HERERCTjBQYYHz2/Npt3ZvPP7BqfjOEYFqIiIADC4W0OGndGYt2duYMKs6tswioiIVJZeLRLp3TKRsT+uYveBvJP/QQBSASoiIofdf15rerdM5NEvlzJ3vXrGFRER8bYHz2/NgdwCxk5b5XQUR6gAFRGRw4KDDM8N7khyzUhunTCfrZnZTkcSEREJKC1qRzO4W0MmzNrA2oz9TsfxORWgIiJylNgIF69e25nsvAJueXceOfnVt6MEERGRynDPWS0ICwniX9/+6XQUn1MBKiIix2leO5pn/nIKC9MzeeizJVhrnY4kIiISMBKjw7itTzOmLtvOzDW7nI7jUypARUTEo3Pb1uGufs35ZH46b/++3uk4IiIiAWVoz8bUiw3niW+WV6sDvSpARUSkRHf1a85ZrWvz2NfLq90RWhERkcoU7grmnrNbsHhzJt8u2eZ0HJ+pUAFqjLnCGLPUGFNkjOlyguX6G2NWGGNWG2Pur8g2RUTEd4KCDGP+0oGU+Ehuf38+6XsOOh1JREQkYFzaKZnmSVH89/sVFBQWOR3HJyp6BnQJcCkwvaQFjDHBwAvAeUAb4EpjTJsKbldERHwkOtzFa0O6kF9QxG3vzSe3QJ0SiYiIeENwkOHv57Zk7c4DfDwv3ek4PlGhAtRau9xau+Iki3UDVltr11pr84APgYsrsl0REfGtJolR/HdQBxalZzLqy2VOxxEREQkY57SpTceGcTz7w6pq0fO8L+4BrQ9sOuJ9evE0j4wxw40xacaYtIyMjEoPJyIipXNu2zrc3KsJ783eyGd/VI+jtCIiIpXNGMN9/VuxbV9Otej076QFqDHmB2PMEg+PSjmLaa191VrbxVrbJTExsTI2ISIi5fR/57akW+NaPPDpYv7cts/pOCIiIgGhR5N4zmyRyIs/ryEzO9/pOJXqpAWotfYsa207D48vSrmNzUCDI94nF08TEZEqJiQ4iHFXdiQ63MWtE+aTlRPYjaSIiIiv/N+5LcnMzufV6WucjlKpfHEJ7lyguTGmsTEmFBgMTPbBdkVEpBIkxYQz7sqObNx9kH98vKhajV1WXsaYN40xO4wxS0qYf7UxZpExZrEx5ndjTAdfZxQREWe1qx/LgA71ePO39ezYl+N0nEpT0WFYBhpj0oFTga+NMd8XT69njPkGwFpbANwBfA8sByZZa5dWLLaIiDipe5N4/nFuS1zLPmb/v1vDyDgY0w4WTXI6mr8aD/Q/wfx1wJnW2vbAY8CrvgglIiL+5d6zW5BfWMTzP652OkqlCanIH1trPwM+8zB9C3D+Ee+/Ab6pyLZERMS/DK85j7ywNwjLyXVPyNwEX97pfp06yLlgfshaO90Yk3KC+b8f8XYW7ttVRESkmklJqMHgbg34YM5Ghp3RhIbxkU5H8jpfXIIrIiIByEwbRZjNPXpifjZMG+VMoMAxFPi2pJnqLV5EJLD9tW9zgoIMz/+4yukolUIFqIiIlE9mCUOxlDRdTsoY0wd3AXpfScuot3gRkcBWOyaca7o34tM/NrNu5wGn43idClARESmf2BKuEi1pupyQMSYVeB242Fq7y+k8IiLinFt6N8EVbHh+WuCdBVUBKiIi5dNvBLgijpp00IaysePfHQpUdRljGgKfAtdaa1c6nUdERJyVFB3OkFNT+HzBZtZk7Hc6jlepABURkfJJHQQXjYXYBoChMCaZf4fcxrVzG2l80GMYYz4AZgItjTHpxpihxphbjDG3FC8yAogHXjTGLDDGpDkWVkRE/MLNvZoQ7gpmbICdBa1QL7giIlLNpQ463ONtMHDR+t1MeHUWD362hLGDT8EY42w+P2GtvfIk828CbvJRHBERqQLio8K47rQUXv5lDXf0aUbz2tFOR/IKnQEVERGv6ZJSi3vPbsGXC7cwce4mp+OIiIhUacPPaEKkK5hnA+gsqApQERHxqlvPbErPZgmM/HIpK7dnOR1HRESkyqpZI5Qbezbm60Vb+XPbPqfjeIUKUBER8aqgIMMzf+lAVFgId7w/n5z8QqcjiYiIVFk39WxCdFgIz04NjLOgKkBFRMTrkqLDeWbQKazcvp8nv/3T6TgiIiJVVmykixt7Nua7pdtYsjnT6TgVpgJUREQqRa8Widxwegrjf1/PLysznI4jIiJSZd3YszEx4SGM+3G101EqTAWoiIhUmvv6t6JF7Sj+/tFCdu3PdTqOiIhIlRQb4eL601L4bum2Kt+/ggpQERGpNOGuYJ4b3JHMg/nc/+lirLVORxIREamSbji9MZGhwbzwU9U+C6oCVEREKlXrujH8o39Lpi7bzocamkVERKRcatYI5doejfhy4RbW7zzgdJxyUwEqIiKV7sbTG9OzWQKjvlzG2oz9TscRERGpkoae0ZiBIb8T+0pHGBkHY9rBoklOxyoTFaAiIlLpgoIM/72iA2GuIO6euID8wiKnI4mIiFQ5Sesm86+Q16iZvx2wkLkJvryzShWhKkBFRMQn6sSG86+B7VmUnsmzP6x0Oo6IiEjVM20UofaYTv3ys2HaKGfylIMKUBER8Znz2tflis7JvPTzGuZt2ON0HBERkaolM71s0/2QClAREfGpERe1oW5sBH//aCHZeYVOxxEREak6YpPLNt0PqQAVERGfig538dQVqazbeYB/f/en03FERESqjn4jwBVx1CQbEuGeXkWoABUREZ87rWkC15+Wwvjf1/P7mp1OxxEREakaUgfBRWMhtgEWQ3pRAt80fsA9vYpQASoiIo74R/+WpMRH8o+PF7E/t8DpOCIiIlVD6iC4Zwlm5F5GN5/E/ataVal2VAWoiIg4IjI0hKcHdWDL3mwe/3q503FERESqnFt7NyUrp4AP52x0OkqpqQAVERHHdG5Ui2FnNOGDORv5ecUOp+OIiIhUKR0axNGjSS3e+G0deQVVY4xtFaAiIuKoe85uQfOkKO7/ZDGZB/OdjiMiIlKl3HxmU7Zm5vDlwi1ORykVFaAiIuKocFcwTw/qQMb+XB79cqnTcURERKqU3i0SaVk7mlemr8Fa63Sck1IBKiIijktNjuP23k359I/NTF223ek4XmeMedMYs8MYs6SE+cYYM9YYs9oYs8gY08nXGUVEpGoyxnDzmU1YuX0/P6/IcDrOSVWoADXGXGGMWWqMKTLGdDnBcuuNMYuNMQuMMWkV2aaIiASmO/o2p1WdaKZ/8gJFz7SFkXEwph0smuR0NG8YD/Q/wfzzgObFj+HASz7IJCIiAeKiDvWoFxvOy7+scTrKSVX0DOgS4FJgeimW7WOtPcVaW2KhKiIi1VdoSBCvdFjLAwUvEbQvHbCQuQm+vLPKF6HW2unA7hMscjHwjnWbBcQZY+r6Jp2IiFR1ruAgbuzZmNnrdvPHxj1OxzmhChWg1trl1toV3gojIiLVW6MF/yXS5B09MT8bpo1yJpDv1Ac2HfE+vXjacYwxw40xacaYtIwM/7/USkREfGNwt4bEhIfw6vS1Tkc5IV/dA2qBKcaYecaY4SdaUA2riEg1lpletunVkLX2VWttF2ttl8TERKfjiIiIn4gKC+HaUxvx3dJtrNt5wOk4JTppAWqM+cEYs8TD4+IybKentbYT7ntcbjfG9CppQTWsIiLVWGxy2aYHjs1AgyPeJxdPExERKbXrTkvBFRzEa7/671nQkxag1tqzrLXtPDy+KO1GrLWbi593AJ8B3cofWUREAla/EeCKOGpSUXCEe3pgmwwMKe4NtweQaa3d6nQoERGpWpKiw7msUzIfz0snIyvX6TgeVfoluMaYGsaY6EOvgXNwd14kIiJytNRBcNFYiG2AxbCVRP4deit5bS53OlmFGGM+AGYCLY0x6caYocaYW4wxtxQv8g2wFlgNvAbc5lBUERGp4oad0Zj8wiLenbXB6SgehVTkj40xA4HngUTga2PMAmvtucaYesDr1trzgdrAZ8aYQ9t731r7XQVzi4hIoEodBKmDMMCy5dt55e00In5ezd1ntXA6WblZa688yXwL3O6jOCIiEsCaJEbRr1US783awG29mxLuCnY60lEq2gvuZ9baZGttmLW2trX23OLpW4qLT6y1a621HYofba21j3sjuIiIBL5+rWtz8Sn1eOGn1azYluV0HBERkSrhxtMbs+tAHl8s8L/uBHzVC66IiEi5jLiwDVFhITzw6SKKiqzTcURERPzeqU3jaVUnmjd/W4/7Ihv/oQJURET8WnxUGA+e35r5G/fy4dxNJ/8DERGRas4Yw9CejVmxPYsZq3c5HecoKkBFRMTvXd45me6Na/Hkt8v9tlc/ERERf3JRh3okRIXyxm/+NSSLClAREfF7xhgeH9ie7PxCRn+9zOk4IiIifi/cFcw1PRrx04oMVu/Y73Scw1SAiohIldAsKYpbezfjiwVb+HVVhtNxRERE/N41PRoRGhLE+N/XOR3lMBWgIiJSZdzWuymNE2rw8OdLyMkvdDqOiIiIX0uICuOSU+rxybzN7D2Y53QcQAWoiIhUIeGuYEZf0o4Nuw7ywk+rnY4jIiLi927s2Zjs/ELen7PR6SiAClAREaliTm+WwMCO9Xn5lzWs3qGxQUVERE6kVZ0YTm8Wzzu/byC/sMjpOCpARUSk6nnogtZEhobw4GdL/G58MxEREX8ztGdjtu3L4ZvFW52OogJURESqnoSoMB44rxVz1u3mo3npTscRERHxa71bJNEkoQZvzljvdBQVoCIiUjUN6tKALo1q8uS3f5J5MN/pOCIiIn4rKMgw5NRGLNy0l0Xpe53N4ujWRUREyikoyDDq4nbsPZjH01NXOB1HRETEr13aOZnI0GDembnB0RwqQEVEpMpqUy+Ga3s0YsKsDSzdkul0HBEREb8VE+7i0k71mbxwC7sPODckiwpQERGp0u49pyU1I0MZ8cVSiorUIZGIiEhJhpyaQl5BEZPSNjmWQQWoiIhUabERLu47rxXzNuzh0z82Ox1HRETEb7WoHU2PJrV4d+YGCh06aKsCVEREqrzLOyXTsWEcT367nMxsdUgkIiJSkutOTWHz3mx++nOHI9tXASoiIlVeUJBh1IB27DqQx7M/rHQ6joiIiN86u01t6sSE8/bM9Y5sXwWoiIgEhPbJsVzVrSHvzNzAn9v2OR1HRETEL4UEB3F194b8umonazP2+3z7KkBFRCRg/N+5LYkJD2HE50uxVh0SiYiIeDK4W0NcwYZ3Z/l+SBYVoCIiEjDiIkP5R/9WzFm/my8WbHE6zlGMMf2NMSuMMauNMfd7mN/QGPOTMeYPY8wiY8z5TuQUEZHAlxgdxvnt6/JxWjoHcgt8um0VoCIiElD+0qUBHZJjeeKb5T5vVEtijAkGXgDOA9oAVxpj2hyz2MPAJGttR2Aw8KJvU4qISHUy5NQUsnIL+HyBb3uQVwEqIiIBJSjI8MiAtuzIyuXFn1c7HeeQbsBqa+1aa20e8CFw8THLWCCm+HUs4F+ncEVEJKB0ahhH23oxvPP7Bp/etqICVEREAk6nhjW5tGN9Xvt1HRt3HXQ6DkB94MhRv9OLpx1pJHCNMSYd+Ab4q6cVGWOGG2PSjDFpGRkZlZFVRESqAWMM152awortWcxdv8dn21UBKiIiAekf/VsREmR44pvlTkcprSuB8dbaZOB84F1jzHHttLX2VWttF2ttl8TERJ+HFBGRwHFRh3pEh4fw3mzfdUakAlRERAJSndhwbuvdlNDlH5PzVGsYGQdj2sGiSU7E2Qw0OOJ9cvG0Iw0FJgFYa2cC4UCCT9KJiEi1FBEazMhGS/nH8suxPmonQyp17SIiIg4aXnMehaFvEH4g1z0hcxN8eaf7deogX0aZCzQ3xjTGXXgOBq46ZpmNQD9gvDGmNe4CVNfYiohI5Vk0iYGb/0OQyXa/90E7qTOgIiISsEJ/Hk0EuUdPzM+GaaN8msNaWwDcAXwPLMfd2+1SY8woY8yA4sX+BgwzxiwEPgCutxrMVEREKtO0UQQVZB89rZLbyQqdATXGPAVcBOQBa4AbrLV7PSzXH3gOCAZet9Y+WZHtioiIlEpmetmmVyJr7Te4Oxc6ctqII14vA073dS4REanGHGgnK3oGdCrQzlqbCqwEHjh2gVKOfSYiIuJ9scllmy4iIlKdONBOVqgAtdZOKb6sCGAW7k4VjlWasc9ERES8r98IcEUcPc0V4Z4uIiJS3TnQTnrzHtAbgW89TC/N2GciIiLelzoILhoLsQ0A436+aKyvOyASERHxTw60kye9B9QY8wNQx8Osh6y1XxQv8xBQALxX0UDGmOHAcICGDRtWdHUiIlLdpQ5SwSkiIlISH7eTJy1ArbVnnWi+MeZ64EKgXwm99ZVm7LMjt/cq8CpAly5d1PufiIiIiIhIgKjQJbjFvdv+AxhgrT1YwmKHxz4zxoTiHvtsckW2KyIiIiIiIlVPRe8BHQdEA1ONMQuMMS8DGGPqGWO+gZLHPqvgdkVERERERKSKqdA4oNbaZiVM3wKcf8T748Y+ExERERERkerFm73gioiIiIiIiJRIBaiIiIiIiIj4hPHcca1/MMZkABu8tLoEYKeX1lVdaJ+VnfZZ2WmflY/2W9l5c581stYmemld5aZ20nHaZ+Wj/VZ22mdlp31Wdt7eZx7bSr8uQL3JGJNmre3idI6qRPus7LTPyk77rHy038pO++zEtH/KTvusfLTfyk77rOy0z8rOV/tMl+CKiIiIiIiIT6gAFREREREREZ+oTgXoq04HqIK0z8pO+6zstM/KR/ut7LTPTkz7p+y0z8pH+63stM/KTvus7Hyyz6rNPaAiIiIiIiLirOp0BlREREREREQcpAJUREREREREfCLgClBjTH9jzApjzGpjzP0e5ocZYyYWz59tjElxIKZfKcU+u9cYs8wYs8gYM80Y08iJnP7kZPvsiOUuM8ZYY0y17wa8NPvMGDOo+N/aUmPM+77O6G9K8f9mQ2PMT8aYP4r//zzfiZz+xBjzpjFmhzFmSQnzjTFmbPE+XWSM6eTrjE5TO1l2aifLTu1k+aitLDu1lWXjF+2ktTZgHkAwsAZoAoQCC4E2xyxzG/By8evBwESnc1eBfdYHiCx+fav22cn3WfFy0cB0YBbQxenc/r7PgObAH0DN4vdJTueuAvvsVeDW4tdtgPVO53b6AfQCOgFLSph/PvAtYIAewGynM/t4/6idrJx9pnayjPuseDm1k2Xcb2ory7XP1FYevT8cbycD7QxoN2C1tXattTYP+BC4+JhlLgbeLn79MdDPGGN8mNHfnHSfWWt/stYeLH47C0j2cUZ/U5p/ZwCPAf8GcnwZzk+VZp8NA16w1u4BsNbu8HFGf1OafWaBmOLXscAWH+bzS9ba6cDuEyxyMfCOdZsFxBlj6vomnV9QO1l2aifLTu1k+aitLDu1lWXkD+1koBWg9YFNR7xPL57mcRlrbQGQCcT7JJ1/Ks0+O9JQ3EdFqrOT7rPiyxUaWGu/9mUwP1aaf2ctgBbGmBnGmFnGmP4+S+efSrPPRgLXGGPSgW+Av/omWpVW1u+8QKN2suzUTpad2snyUVtZdmorva/S28kQb65MApsx5hqgC3Cm01n8mTEmCHgGuN7hKFVNCO5Li3rjPnsw3RjT3lq718lQfu5KYLy19mljzKnAu8aYdtbaIqeDiVRHaidLR+1khaitLDu1lX4m0M6AbgYaHPE+uXiax2WMMSG4T8Xv8kk6/1SafYYx5izgIWCAtTbXR9n81cn2WTTQDvjZGLMe9/Xzk6t5Bwul+XeWDky21uZba9cBK3E3stVVafbZUGASgLV2JhAOJPgkXdVVqu+8AKZ2suzUTpad2snyUVtZdmorva/S28lAK0DnAs2NMY2NMaG4O0+YfMwyk4Hril9fDvxoi++4raZOus+MMR2BV3A3qtX9XgM4yT6z1mZaaxOstSnW2hTc9wMNsNamORPXL5Tm/83PcR/RxRiTgPsyo7U+zOhvSrPPNgL9AIwxrXE3qhk+TVn1TAaGFPfy1wPItNZudTqUD6mdLDu1k2WndrJ81FaWndpK76v0djKgLsG11hYYY+4AvsfdK9ab1tqlxphRQJq1djLwBu5T76tx34A72LnEzivlPnsKiAI+Ku6HYqO1doBjoR1Wyn0mRyjlPvseOMcYswwoBP7PWlttz7qUcp/9DXjNGHMP7k4Wrq/mhQLGmA9w/zhLKL7f5xHABWCtfRn3/T/nA6uBg8ANziR1htrJslM7WXZqJ8tHbWXZqa0sO39oJ0013v8iIiIiIiLiQ4F2Ca6IiIiIiIj4KRWgIiIiIiIi4hMqQEVERERERMQnVICKiIiIiIiIT6gAFREREREREZ9QASoS4IwxccaY25zOISIi4q/UVor4jgpQkcAXB6hRFRERKVkcaitFfEIFqEjgexJoaoxZYIx5yukwIiIifkhtpYiPGGut0xlEpBIZY1KAr6y17ZzOIiIi4o/UVor4js6AioiIiIiIiE+oABURERERERGfUAEqEviygGinQ4iIiPgxtZUiPqICVCTAWWt3ATOMMUvUsYKIiMjx1FaK+I46IRIRERERERGf0BlQERERERER8QkVoCIiIiIiIuITKkBFRERERETEJ1SAioiIiIiIiE+oABURERERERGfUAEqIiIiIiIiPqECVERERERERHzi/wGMJFmGva32iwAAAABJRU5ErkJggg==", + "application/vnd.jupyter.widget-view+json": { + "model_id": "6010f8d88e4544dd973089a9bc96fb60", + "version_major": 2, + "version_minor": 0 + }, "text/plain": [ - "
" + "interactive(children=(FloatSlider(value=0.0, description='t', max=2519.87471894461, step=25.1987471894461), Ou…" ] }, - "metadata": { - "needs_background": "light" - }, + "metadata": {}, "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "t_fine = np.linspace(0, t[-1], 1000)\n", - "\n", - "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))\n", - "ax1.plot(t_fine, 2 * np.exp(t_fine) - np.exp(2 * t_fine), t_sol, x(t_sol), \"o\")\n", - "ax1.set_xlabel(\"t\")\n", - "ax1.legend([\"2*exp(t) - exp(2*t)\", \"x\"], loc=\"best\")\n", - "\n", - "ax2.plot(t_fine, 3 * np.exp(t_fine) - np.exp(2 * t_fine), t_sol, y(t_sol), \"o\")\n", - "ax2.set_xlabel(\"t\")\n", - "ax2.legend([\"3*exp(t) - exp(2*t)\", \"y\"], loc=\"best\")\n", - "\n", - "plt.tight_layout()\n", - "plt.show()" + "sol.plot([\"Voltage [V]\", \"Negative electrode stoichiometry\", \"Positive electrode stoichiometry\"])" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "In the [next notebook](./2-a-pde-model.ipynb) we show how to create, discretise and solve a PDE model in pybamm." + "In this tutorial we have introduced the main concepts to define a PyBaMM model and used them to build the reservoir model. In the [next notebook](./2-a-pde-model.ipynb) we will build on what we have learn to define a simple partial differential equation (PDE) model: the Single Particle Model for a half cell." ] }, { @@ -390,7 +481,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 15, "metadata": {}, "outputs": [ { @@ -399,8 +490,7 @@ "text": [ "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n", "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n", - "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.\n", - "[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n", + "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n", "\n" ] } @@ -426,7 +516,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.11.6" } }, "nbformat": 4,