![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\" "\\"Open"){kind=icon}
![\"Download\"](\"https://img.shields.io/badge/Github-Download-blue.svg\" "\\"Download"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[2.6 Dynamic Optimization with Pyomo.DAE](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6-Dynamic-Optimization-with-Pyomo.DAE)",
"section": "2.6 Dynamic Optimization with Pyomo.DAE"
}
},
"source": [
"# 2.6 Dynamic Optimization with Pyomo.DAE\n",
"\n",
"Adapted from https://github.com/Pyomo/pyomo/blob/master/examples/dae/car_example.py"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpages": {
"level": 1,
"link": "[2.6 Dynamic Optimization with Pyomo.DAE](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6-Dynamic-Optimization-with-Pyomo.DAE)",
"section": "2.6 Dynamic Optimization with Pyomo.DAE"
}
},
"outputs": [],
"source": [
"import sys\n",
"if \"google.colab\" in sys.modules:\n",
" !wget \"https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py\"\n",
" import helper\n",
" #!pip install casadi\n",
" helper.install_idaes()\n",
" helper.install_ipopt()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[2.6.1 Car Example](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1-Car-Example)",
"section": "2.6.1 Car Example"
}
},
"source": [
"## 2.6.1 Car Example\n",
"\n",
"You are a race car driver with a simple goal. Drive distance $L$ in the minimal amount of time **but** come to a complete stop at the finish line.\n",
"\n",
"### Optimal Control Problem Formulation\n",
"\n",
"Mathematically, you want to solve the following optimal control problem:\n",
"\n",
"$$\\begin{align*}\n",
"\\min_{u} \\quad & t_f \\\\\n",
"\\mathrm{s.t.} \\quad & \\frac{dx}{dt} = v \\\\\n",
"& \\frac{dv}{dt} = u - R v^2 \\\\\n",
"& x(t=0) = 0, ~~ x(t=t_f) = L \\\\\n",
"& v(t=0) = 0, ~~ v(t=t_f) = 0 \\\\\n",
"& -3 \\leq u \\leq 1\n",
"\\end{align*}$$\n",
"\n",
"where $a$ is the acceleration/braking (your control variable) and $R$ is the drag coefficient (parameter).\n",
"\n",
"### Scale Time Scale Time\n",
"\n",
"Let $t = \\tau \\cdot t_f$ where $\\tau \\in [0,1]$. Thus $dt = t_f d\\tau$. The optimal control problem becomes:\n",
"\n",
"$$\\begin{align*}\n",
"\\min_{u} \\quad & t_f \\\\\n",
"\\mathrm{s.t.} \\quad & \\frac{dx}{d\\tau} = t_f v \\\\\n",
"& \\frac{dv}{d\\tau} = t_f (u - R v^2) \\\\\n",
"& x(\\tau = 0) = 0, ~~ x(\\tau = 1) = L \\\\\n",
"& v(\\tau = 0) = 0, ~~ v(\\tau = 1) = 0 \\\\\n",
"& -3 \\leq u \\leq 1\n",
"\\end{align*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.1-Orthogonal-Collocation-on-Finite-Elements:-Manual-Approach)",
"section": "2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach"
}
},
"source": [
"### 2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach\n",
"\n",
"**Indices and Sets**\n",
"* Finite elements: $i \\in \\mathcal{I} = \\{1,2,...,N_{FE}\\}$\n",
"* Finite elements (except first): $\\mathcal{I} \\setminus 1$\n",
"* Internal collocation points: $j,k \\in \\mathcal{J} = \\{1,2,...,N_{C}\\}$\n",
"\n",
"\n",
"**Parameters**\n",
"* Coefficients in collocation/Runge-Kutta formula: $a_{j,j}$\n",
"* Drag coefficient: $R$\n",
"* Race length: $L$\n",
"* Scaled time for each finite element: $h_i = \\frac{1}{N_{FE}}$\n",
"\n",
"**Variables**\n",
"* Position (internal collocation points): $x_{i,j}$\n",
"* Position (beginning of each finite element): $\\bar{x}_{i}$\n",
"* Position derivative (internal collocation points): $\\dot{x}_{i,j}$\n",
"* Velocity (internal collocation points): $v_{i,j}$\n",
"* Velocity (beginning of each finite element): $\\bar{v}_{i}$\n",
"* Velocity derivative (internal collocation points): $\\dot{v}_{i,j}$\n",
"* Time (internal collocation points): $t_{i,j}$\n",
"* Time (beginning of each finite element): $\\bar{t}_{i}$\n",
"* Acceleration (control degrees of freedom): $u_i$\n",
"\n",
"**Objective and Constraints**\n",
"\n",
"$$\\begin{align*}\n",
"\\min \\quad & t_f \\\\\n",
"\\mathrm{s.t.} \\quad & \\mathrm{Differential~equation~for~x:} \\\\ \n",
"& x_{i,j} = \\bar{x}_{i} + h_i \\sum_{k \\in \\mathcal{J}} a_{k,j} \\dot{x}_{i,k},\\quad \\forall i \\in \\mathcal{I},~j \\in \\mathcal{J} \\\\\n",
"& \\bar{x}_{i} = \\bar{x}_{i-1} + h_{i-1} \\sum_{k \\in \\mathcal{J}} a_{k,N_C} \\dot{x}_{i-1,k},\\quad \\forall i \\in \\mathcal{I} \\setminus 1\\\\\n",
"& \\dot{x}_{i,j} = t_f v_{i,j}, \\quad \\forall i \\in \\mathcal{I} \\setminus 1,~j \\in \\mathcal{J} \\\\ \\\\\n",
"& \\mathrm{Differential~equation~for~v:} \\\\ \n",
"& v_{i,j} = \\bar{v}_{i} + h_i \\sum_{k \\in \\mathcal{J}} a_{k,j} \\dot{v}_{i,k},\\quad \\forall i \\in \\mathcal{I},~j \\in \\mathcal{J} \\\\\n",
"& \\bar{v}_{i} = \\bar{v}_{i-1} + h_{i-1} \\sum_{k \\in \\mathcal{J}} a_{k,N_C} \\dot{v}_{i-1,k},\\quad \\forall i \\in \\mathcal{I} \\setminus 1\\\\\n",
"& \\dot{v}_{i,j} = t_f (u_{i} - R v_{i,j}^2), \\quad \\forall i \\in \\mathcal{I} \\setminus 1,~j \\in \\mathcal{J} \\\\ \\\\\n",
"& \\mathrm{Differential~equation~for~time:} \\\\\n",
"& t_{i,j} = \\bar{t}_{i} + h_i \\sum_{k \\in \\mathcal{J}} a_{k,j} t_f,\\quad \\forall i \\in \\mathcal{I},~j \\in \\mathcal{J} \\\\\n",
"& \\bar{t}_{i} = \\bar{t}_{i-1} + h_{i-1} \\sum_{k \\in \\mathcal{J}} a_{k,N_C} t_f,\\quad \\forall i \\in \\mathcal{I} \\setminus 1\\\\ \\\\\n",
"& \\mathrm{Initial~conditions:} \\\\ \n",
"& \\bar{x}_{1} = 0, \\quad \\bar{v}_{1} = 0, \\quad \\bar{t}_1 = 0\\\\ \\\\\n",
"& \\mathrm{Boundary~conditions:} \\\\ \n",
"& L = \\bar{x}_{N_{FE}} + h_{N_{FE}} \\sum_{k \\in \\mathcal{J}} a_{k,N_C} \\dot{x}_{N_{FE},k}, \\\\\n",
"& 0 = \\bar{v}_{N_{FE}} + h_{N_{FE}} \\sum_{k \\in \\mathcal{J}} a_{k,N_C} \\dot{v}_{N_{FE},k}, \\\\ \\\\\n",
"& \\mathrm{Bounds~on~acceleration:} \\\\\n",
"& -3 \\leq a_i \\leq 1, \\quad \\forall i \\in \\mathcal{I}\n",
"\\end{align*}$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.1-Orthogonal-Collocation-on-Finite-Elements:-Manual-Approach)",
"section": "2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach"
}
},
"outputs": [],
"source": [
"from pyomo.environ import *\n",
"import matplotlib.pyplot as plt\n",
"from pyomo.dae import *\n",
"\n",
"m2 = ConcreteModel()\n",
"\n",
"# Define model parameters\n",
"m2.R = Param(initialize=0.001) # Friction factor\n",
"m2.L = Param(initialize=100.0) # Final position\n",
"\n",
"# Define finite elements and collocation points\n",
"NFE = 15\n",
"NC = 3\n",
"m2.I = Set(initialize=RangeSet(1,NFE)) # Finite elements\n",
"m2.J = Set(initialize=RangeSet(1,NC)) # Internal collocation points\n",
"\n",
"# Define first order derivative collocation matrix\n",
"A = {}\n",
"A[1,1] = 0.19681547722366\n",
"A[1,2] = 0.39442431473909\n",
"A[1,3] = 0.37640306270047\n",
"A[2,1] = -0.06553542585020\n",
"A[2,2] = 0.29207341166523\n",
"A[2,3] = 0.51248582618842\n",
"A[3,1] = 0.02377097434822\n",
"A[3,2] = -0.04154875212600\n",
"A[3,3] = 0.11111111111111\n",
"\n",
"m2.a = Param(m2.J, m2.J, initialize=A)\n",
"\n",
"# Define step for each finite element\n",
"m2.h = Param(m2.I,initialize=1/NFE)\n",
"\n",
"# Define objective (final time)\n",
"m2.tf = Var(domain=NonNegativeReals)\n",
"\n",
"# Variables for x\n",
"m2.x0 = Var(m2.I)\n",
"m2.x = Var(m2.I,m2.J)\n",
"m2.xdot = Var(m2.I,m2.J)\n",
"\n",
"# Variables for v\n",
"m2.v0 = Var(m2.I)\n",
"m2.v = Var(m2.I,m2.J)\n",
"m2.vdot = Var(m2.I,m2.J)\n",
"\n",
"# Variables for t\n",
"m2.t0 = Var(m2.I)\n",
"m2.t = Var(m2.I,m2.J)\n",
"\n",
"# Acceleration\n",
"m2.u = Var(m2.I, bounds=(-3,1))\n",
"\n",
"### Finite Element Collocation Equations\n",
"# position\n",
"def FECOLx_(m2,i,j):\n",
" return m2.x[i,j] == m2.x0[i] + m2.h[i]*sum(m2.a[k,j]*m2.xdot[i,k] for k in m2.J)\n",
"m2.FECOLx = Constraint(m2.I,m2.J,rule=FECOLx_)\n",
"\n",
"# velocity\n",
"def FECOLv_(m2,i,j):\n",
" return m2.v[i,j] == m2.v0[i] + m2.h[i]*sum(m2.a[k,j]*m2.vdot[i,k] for k in m2.J)\n",
"m2.FECOLv = Constraint(m2.I,m2.J,rule=FECOLv_)\n",
"\n",
"# time\n",
"def FECOLt_(m2,i,j):\n",
" return m2.t[i,j] == m2.t0[i] + m2.h[i]*sum(m2.a[k,j]*m2.tf for k in m2.J)\n",
"m2.FECOLt = Constraint(m2.I,m2.J,rule=FECOLt_)\n",
"\n",
"\n",
"### Continuity Equations\n",
"# position\n",
"def CONx_(m2,i):\n",
" if i == 1:\n",
" return Constraint.Skip\n",
" else:\n",
" return m2.x0[i] == m2.x0[i-1] + m2.h[i-1]*sum(m2.a[k,NC]*m2.xdot[i,k] for k in m2.J)\n",
"m2.CONx = Constraint(m2.I, rule=CONx_)\n",
"\n",
"# velocity\n",
"def CONv_(m2,i):\n",
" if i == 1:\n",
" return Constraint.Skip\n",
" else:\n",
" return m2.v0[i] == m2.v0[i-1] + m2.h[i-1]*sum(m2.a[k,NC]*m2.vdot[i,k] for k in m2.J)\n",
"m2.CONv = Constraint(m2.I, rule=CONv_)\n",
"\n",
"# time\n",
"def CONt_(m2,i):\n",
" if i == 1:\n",
" return Constraint.Skip\n",
" else:\n",
" return m2.t0[i] == m2.t0[i-1] + m2.h[i-1]*sum(m2.a[k,NC]*m2.tf for k in m2.J)\n",
"m2.CONt = Constraint(m2.I, rule=CONt_)\n",
"\n",
"### Differential equations\n",
"# position\n",
"def ODEx_(m2,i,j):\n",
" return m2.xdot[i,j] == m2.tf*m2.v[i,j]\n",
"m2.ODEx = Constraint(m2.I,m2.J,rule=ODEx_)\n",
"\n",
"# velocity\n",
"def ODEv_(m2,i,j):\n",
" return m2.vdot[i,j] == m2.tf*(m2.u[i] - m2.R*m2.v[i,j]**2)\n",
"m2.ODEv = Constraint(m2.I,m2.J,rule=ODEv_)\n",
"\n",
"### Initial conditions\n",
"m2.xIC = Constraint(expr=m2.x0[1] == 0)\n",
"m2.vIC = Constraint(expr=m2.v0[1] == 0)\n",
"m2.tIC = Constraint(expr=m2.t0[1] == 0)\n",
"\n",
"### Boundary conditions\n",
"m2.xBC = Constraint(expr=m2.L == m2.x0[NFE] + m2.h[NFE]*sum(m2.a[k,NC]*m2.xdot[NFE,k] for k in m2.J))\n",
"m2.vBC = Constraint(expr=0 == m2.v0[NFE] + m2.h[NFE]*sum(m2.a[k,NC]*m2.vdot[NFE,k] for k in m2.J))\n",
"\n",
"### Set objective\n",
"m2.obj = Objective(expr=m2.tf)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.1-Orthogonal-Collocation-on-Finite-Elements:-Manual-Approach)",
"section": "2.6.1.1 Orthogonal Collocation on Finite Elements: Manual Approach"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ipopt 3.13.2: \n",
"\n",
"******************************************************************************\n",
"This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
" Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
" For more information visit http://projects.coin-or.org/Ipopt\n",
"******************************************************************************\n",
"\n",
"This is Ipopt version 3.13.2, running with linear solver ma27.\n",
"\n",
"Number of nonzeros in equality constraint Jacobian...: 1093\n",
"Number of nonzeros in inequality constraint Jacobian.: 0\n",
"Number of nonzeros in Lagrangian Hessian.............: 105\n",
"\n",
"Total number of variables............................: 286\n",
" variables with only lower bounds: 1\n",
" variables with lower and upper bounds: 15\n",
" variables with only upper bounds: 0\n",
"Total number of equality constraints.................: 272\n",
"Total number of inequality constraints...............: 0\n",
" inequality constraints with only lower bounds: 0\n",
" inequality constraints with lower and upper bounds: 0\n",
" inequality constraints with only upper bounds: 0\n",
"\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 0 9.9999900e-03 1.00e+02 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n",
" 1r 9.9999900e-03 1.00e+02 9.99e+02 2.0 0.00e+00 - 0.00e+00 1.05e-07R 2\n",
" 2r 3.0578540e+00 9.99e+01 9.31e+02 2.0 4.49e+01 - 1.35e-01 6.80e-02f 2\n",
" 3r 2.2146807e+00 9.91e+01 5.68e+02 1.3 2.16e+00 - 1.00e+00 3.91e-01f 1\n",
" 4r 2.0626784e+00 9.87e+01 4.21e+02 0.6 2.94e+00 - 6.36e-01 2.59e-01f 1\n",
" 5r 2.7501672e+00 9.77e+01 2.42e+02 0.6 8.74e+00 - 5.56e-01 4.24e-01f 1\n",
" 6r 5.8696920e+00 9.43e+01 1.91e+02 -0.1 8.62e+00 - 4.13e-01 7.57e-01f 1\n",
" 7r 8.6536528e+00 8.80e+01 1.65e+02 -0.1 1.29e+01 - 2.98e-01 7.27e-01f 1\n",
" 8r 9.0137628e+00 8.55e+01 1.09e+02 -0.1 3.77e+00 - 5.42e-01 1.00e+00f 1\n",
" 9 9.4134225e+00 8.42e+01 9.84e-01 -1.0 1.13e+02 - 3.63e-01 1.56e-02h 7\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 10 9.7719760e+00 8.29e+01 9.61e-01 -1.0 1.13e+02 - 6.94e-01 1.56e-02h 7\n",
" 11 1.0069154e+01 8.16e+01 1.10e+00 -1.0 1.20e+02 - 6.26e-01 1.52e-02h 7\n",
" 12 1.0350049e+01 8.04e+01 1.27e+00 -1.0 1.24e+02 - 1.00e+00 1.56e-02h 7\n",
" 13 1.0613260e+01 7.91e+01 1.29e+00 -1.0 1.25e+02 - 6.03e-01 1.56e-02h 7\n",
" 14 1.0858870e+01 7.79e+01 1.31e+00 -1.0 1.25e+02 - 1.00e+00 1.56e-02h 7\n",
" 15 1.1089714e+01 7.66e+01 1.31e+00 -1.0 1.23e+02 - 7.24e-01 1.56e-02h 7\n",
" 16 1.1306240e+01 7.55e+01 1.31e+00 -1.0 1.24e+02 - 1.00e+00 1.56e-02h 7\n",
" 17 1.1714847e+01 7.31e+01 1.28e+00 -1.0 1.22e+02 - 8.51e-01 3.12e-02h 6\n",
" 18 1.2077191e+01 7.08e+01 1.24e+00 -1.0 1.23e+02 - 1.00e+00 3.12e-02h 6\n",
" 19 2.2531170e+01 7.10e+01 8.96e-01 -1.0 1.19e+02 - 1.00e+00 1.00e+00w 1\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 20 1.7706786e+01 4.49e+00 1.75e-01 -1.0 2.28e+01 - 1.00e+00 1.00e+00w 1\n",
" 21 1.7349863e+01 1.01e-01 7.40e-04 -1.0 8.28e+00 - 1.00e+00 1.00e+00h 1\n",
" 22 1.6478946e+01 1.06e+00 3.16e-02 -2.5 1.68e+01 - 8.21e-01 1.00e+00h 1\n",
" 23 1.6333670e+01 1.11e-01 9.02e-04 -2.5 1.27e+01 - 9.80e-01 1.00e+00h 1\n",
" 24 1.6324255e+01 5.10e-04 4.40e-06 -2.5 6.64e-01 - 1.00e+00 1.00e+00h 1\n",
" 25 1.6289923e+01 3.05e-03 1.05e-05 -3.8 1.03e+00 - 1.00e+00 1.00e+00h 1\n",
" 26 1.6289753e+01 1.16e-07 1.47e-09 -3.8 5.97e-03 - 1.00e+00 1.00e+00h 1\n",
" 27 1.6287822e+01 9.84e-06 3.34e-08 -5.7 5.95e-02 - 1.00e+00 1.00e+00f 1\n",
" 28 1.6287798e+01 1.59e-09 5.42e-12 -8.6 7.53e-04 - 1.00e+00 1.00e+00h 1\n",
"\n",
"Number of Iterations....: 28\n",
"\n",
" (scaled) (unscaled)\n",
"Objective...............: 1.6287797665173553e+01 1.6287797665173553e+01\n",
"Dual infeasibility......: 5.4235099871471539e-12 5.4235099871471539e-12\n",
"Constraint violation....: 1.5949126463965513e-09 1.5949126463965513e-09\n",
"Complementarity.........: 2.5440456242856559e-09 2.5440456242856559e-09\n",
"Overall NLP error.......: 2.5440456242856559e-09 2.5440456242856559e-09\n",
"\n",
"\n",
"Number of objective function evaluations = 106\n",
"Number of objective gradient evaluations = 23\n",
"Number of equality constraint evaluations = 106\n",
"Number of inequality constraint evaluations = 0\n",
"Number of equality constraint Jacobian evaluations = 30\n",
"Number of inequality constraint Jacobian evaluations = 0\n",
"Number of Lagrangian Hessian evaluations = 28\n",
"Total CPU secs in IPOPT (w/o function evaluations) = 0.014\n",
"Total CPU secs in NLP function evaluations = 0.001\n",
"\n",
"EXIT: Optimal Solution Found.\n",
"final time = 16.29\n"
]
}
],
"source": [
"solver = SolverFactory('ipopt')\n",
"solver.solve(m2,tee=True)\n",
"\n",
"print(\"final time = %6.2f\" %(value(m2.tf)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.6.1.2 Orthogonal Collocation on Finite Elements: Pyomo.dae](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.2-Orthogonal-Collocation-on-Finite-Elements:-Pyomo.dae)",
"section": "2.6.1.2 Orthogonal Collocation on Finite Elements: Pyomo.dae"
}
},
"source": [
"### 2.6.1.2 Orthogonal Collocation on Finite Elements: Pyomo.dae\n",
"\n",
"#### Declare Model"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.6.1.2 Orthogonal Collocation on Finite Elements: Pyomo.dae](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.2-Orthogonal-Collocation-on-Finite-Elements:-Pyomo.dae)",
"section": "2.6.1.2 Orthogonal Collocation on Finite Elements: Pyomo.dae"
}
},
"outputs": [],
"source": [
"m = ConcreteModel()\n",
"\n",
"m.R = Param(initialize=0.001) # Friction factor\n",
"m.L = Param(initialize=100.0) # Final position\n",
"\n",
"m.tau = ContinuousSet(bounds=(0,1)) # Unscaled time\n",
"m.time = Var(m.tau) # Scaled time\n",
"m.tf = Var()\n",
"m.x = Var(m.tau,bounds=(0,m.L+50))\n",
"m.v = Var(m.tau,bounds=(0,None))\n",
"m.u = Var(m.tau, bounds=(-3.0,1.0),initialize=0)\n",
"\n",
"m.dtime = DerivativeVar(m.time)\n",
"m.dx = DerivativeVar(m.x)\n",
"m.dv = DerivativeVar(m.v)\n",
"\n",
"m.obj = Objective(expr=m.tf)\n",
"\n",
"def _ode1(m,i):\n",
" if i == 0 :\n",
" return Constraint.Skip\n",
" return m.dx[i] == m.tf * m.v[i]\n",
"m.ode1 = Constraint(m.tau, rule=_ode1)\n",
"\n",
"def _ode2(m,i):\n",
" if i == 0 :\n",
" return Constraint.Skip\n",
" return m.dv[i] == m.tf*(m.u[i] - m.R*m.v[i]**2)\n",
"m.ode2 = Constraint(m.tau, rule=_ode2)\n",
"\n",
"def _ode3(m,i):\n",
" if i == 0:\n",
" return Constraint.Skip\n",
" return m.dtime[i] == m.tf\n",
"m.ode3 = Constraint(m.tau, rule=_ode3)\n",
"\n",
"def _init(m):\n",
" yield m.x[0] == 0\n",
" yield m.x[1] == m.L\n",
" yield m.v[0] == 0\n",
" yield m.v[1] == 0\n",
" yield m.time[0] == 0\n",
"m.initcon = ConstraintList(rule=_init)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 4,
"link": "[2.6.1.2.1 Discretize/Transcribe and Solve](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.2.1-Discretize/Transcribe-and-Solve)",
"section": "2.6.1.2.1 Discretize/Transcribe and Solve"
}
},
"source": [
"#### 2.6.1.2.1 Discretize/Transcribe and Solve"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"nbpages": {
"level": 4,
"link": "[2.6.1.2.1 Discretize/Transcribe and Solve](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.2.1-Discretize/Transcribe-and-Solve)",
"section": "2.6.1.2.1 Discretize/Transcribe and Solve"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ipopt 3.13.2: \n",
"\n",
"******************************************************************************\n",
"This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
" Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
" For more information visit http://projects.coin-or.org/Ipopt\n",
"******************************************************************************\n",
"\n",
"This is Ipopt version 3.13.2, running with linear solver ma27.\n",
"\n",
"Number of nonzeros in equality constraint Jacobian...: 1145\n",
"Number of nonzeros in inequality constraint Jacobian.: 0\n",
"Number of nonzeros in Lagrangian Hessian.............: 135\n",
"\n",
"Total number of variables............................: 319\n",
" variables with only lower bounds: 46\n",
" variables with lower and upper bounds: 91\n",
" variables with only upper bounds: 0\n",
"Total number of equality constraints.................: 305\n",
"Total number of inequality constraints...............: 0\n",
" inequality constraints with only lower bounds: 0\n",
" inequality constraints with lower and upper bounds: 0\n",
" inequality constraints with only upper bounds: 0\n",
"\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 0 0.0000000e+00 1.00e+02 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n",
" 1 3.6882930e+01 9.96e+01 6.45e+02 -1.0 1.00e+04 - 9.91e-05 3.69e-03h 9\n",
" 2 3.9470634e+01 1.07e+01 1.19e+03 -1.0 1.53e+02 0.0 2.42e-03 9.90e-01f 1\n",
" 3 3.8532247e+01 2.51e-01 2.03e+02 -1.0 2.54e+00 -0.5 9.65e-01 9.90e-01h 1\n",
" 4 3.8410228e+01 7.80e-04 3.81e-01 -1.0 1.22e-01 -1.0 1.00e+00 1.00e+00h 1\n",
" 5 1.1753203e+01 1.06e+02 6.26e+03 -2.5 2.54e+03 - 6.71e-02 1.02e-01f 1\n",
" 6 1.4082983e+01 7.02e+01 4.81e+04 -2.5 1.08e+02 - 7.78e-02 3.38e-01h 1\n",
" 7 1.5801574e+01 3.64e+01 8.78e+04 -2.5 8.55e+01 - 8.81e-02 4.94e-01h 1\n",
" 8 1.6404903e+01 1.44e+01 8.17e+04 -2.5 7.55e+01 - 2.84e-01 6.29e-01h 1\n",
" 9 1.6434448e+01 6.13e+00 2.33e+04 -2.5 4.53e+01 - 6.87e-01 5.78e-01h 1\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 10 1.6664583e+01 1.13e-01 2.73e+03 -2.5 7.97e+00 - 9.14e-01 1.00e+00h 1\n",
" 11 1.6635032e+01 8.41e-03 1.02e-05 -2.5 4.85e+00 - 1.00e+00 1.00e+00h 1\n",
" 12 1.6528482e+01 1.77e-02 3.21e-04 -3.8 1.54e+00 - 1.00e+00 1.00e+00h 1\n",
" 13 1.6526502e+01 7.15e-06 5.43e-08 -3.8 2.76e-02 - 1.00e+00 1.00e+00h 1\n",
" 14 1.6520275e+01 6.32e-05 1.05e-06 -5.7 9.58e-02 - 1.00e+00 1.00e+00h 1\n",
" 15 1.6520269e+01 7.57e-11 1.88e-11 -5.7 1.33e-04 - 1.00e+00 1.00e+00h 1\n",
" 16 1.6520192e+01 9.76e-09 1.62e-10 -8.6 1.19e-03 - 1.00e+00 1.00e+00f 1\n",
"\n",
"Number of Iterations....: 16\n",
"\n",
" (scaled) (unscaled)\n",
"Objective...............: 1.6520191521080921e+01 1.6520191521080921e+01\n",
"Dual infeasibility......: 1.6218795328535750e-10 1.6218795328535750e-10\n",
"Constraint violation....: 9.7603845006233314e-09 9.7603845006233314e-09\n",
"Complementarity.........: 2.5852537039546602e-09 2.5852537039546602e-09\n",
"Overall NLP error.......: 9.7603845006233314e-09 9.7603845006233314e-09\n",
"\n",
"\n",
"Number of objective function evaluations = 27\n",
"Number of objective gradient evaluations = 17\n",
"Number of equality constraint evaluations = 27\n",
"Number of inequality constraint evaluations = 0\n",
"Number of equality constraint Jacobian evaluations = 17\n",
"Number of inequality constraint Jacobian evaluations = 0\n",
"Number of Lagrangian Hessian evaluations = 16\n",
"Total CPU secs in IPOPT (w/o function evaluations) = 0.009\n",
"Total CPU secs in NLP function evaluations = 0.001\n",
"\n",
"EXIT: Optimal Solution Found.\n",
"final time = 16.52\n"
]
}
],
"source": [
"#discretizer = TransformationFactory('dae.finite_difference')\n",
"#discretizer.apply_to(m,nfe=15,scheme='BACKWARD')\n",
"\n",
"discretizer = TransformationFactory('dae.collocation')\n",
"discretizer.apply_to(m,nfe=15,scheme='LAGRANGE-RADAU',ncp=3)\n",
"#discretizer.apply_to(m,nfe=15,scheme='LAGRANGE-LEGENDRE',ncp=3)\n",
"\n",
"# force piecewise constant controls (acceleration) over each finite element\n",
"m = discretizer.reduce_collocation_points(m,var=m.u,ncp=1,contset=m.tau)\n",
"\n",
"solver = SolverFactory('ipopt')\n",
"solver.solve(m,tee=True)\n",
"\n",
"print(\"final time = %6.2f\" %(value(m.tf)))"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 4,
"link": "[2.6.1.2.2 Plot Results](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.2.2-Plot-Results)",
"section": "2.6.1.2.2 Plot Results"
}
},
"source": [
"#### 2.6.1.2.2 Plot Results"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"nbpages": {
"level": 4,
"link": "[2.6.1.2.2 Plot Results](https://ndcbe.github.io/CBE60499/02.06-Pyomo-DAE.html#2.6.1.2.2-Plot-Results)",
"section": "2.6.1.2.2 Plot Results"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+SklEQVR4nO2dd5xc5XX3v2e2qu2q7q567wJRFhCIajUQRfAax8YNl4SQgOMUYxPbKZ83Vozj5I0LSWxi48hxwAGbIkAFSfQqJBDSruqqS9u1KiuttszMef+4d1aj1ay0ZWbunZnz/XzmM7fNvWfmmXvmmd9znnNEVTEMwzDSi4DXBhiGYRjxx5y7YRhGGmLO3TAMIw0x524YhpGGmHM3DMNIQ8y5G4ZhpCEZ59xFZJ+IzE/i9b4tIr9I1vWM7iEiXxKRt3p5jutEZEe8bDLig4ioiExK4vVOisiEZF3vQmScc08kInKjiByK3qaq/6iqf+iVTUbiUdU3VXVqZD3ZHQgj+YjIayJy1n2tqv1VdY9XNnXEnLthGEYUIpLttQ3xIGOdu4jkiciPRKTSffxIRPKi9i8RkU0ickJEdovIze72L4vINhFpFJE9IvLH7vZ+wEpghPv37KSIjBCRvxeR30Sd9w4RKReRY+6v//SofftE5BsisllEjovI/4pIfvI+ldRFRB4Wkd912PZjEfmJiBSKyC9FpEpEDovI90Qkq5PzXCMiH7if/wcick3UvsEi8iv3+3JURJ5zt7f/YxOR/wbGAC+434FvishLIvK1DtfZLCJ3xvdT8D9uO+1275+tInJX1L4/irq3torIZe720SLyjIjUicgREXk06jVfcV9zVERWi8jYTq6bJyL/LCIHRKRGRH4mIn3cfTeKyCER+ZaIVAO/EpFBIvKie82j7vIo9/ilwHXAo24bP+pub5eB3O/cr93X7xeR74pIwN33JRF5y7XnqIjsFZFb4v5hq2pGPYB9wHzg/wLvAUXAMOAd4B/cY64EjgMLcH4ARwLT3H23AhMBAW4AmoDL3H03Aoc6XO/vgd+4y1OAU+55c4BvAhVAbpRt64ERwGBgG3C/159ZKjyAsW5bFLjrWUAVMAd4Dvg50M9t7/XAH7vHfQl4y10eDBwFvgBkA/e460Pc/S8B/wsMctvvhljtHvmORa3/AfB+1Pps4Eik3TPpAXzK/X4HgE+798Nwd/th4Ar33prktmkW8DHwr2775QPXuue6071/prvt9V3gnahrKTDJXf4RsNxt4wHAC8D3o9ovCPwAyAP6AEOATwJ93eOfBp6LOvdrwB92eG/R1/s18Lz72nHATuCrUd+5NuCP3Pf3J0AlIHH9rL1ubA++XPtwnPtuYHHU9kXAPnf558C/dvF8zwFfj/qSnM+5/w3wVNS+gPuFvjHKts9H7f8n4Gdef2ap8gDeAr7oLi9w27gYaAH6RB13D/Cqu/wlzjj3LwDrO5zzXfeY4UAYGBTjume1O+c69zygAZjsrv8z8O9ef15+eACbgCXA6sh91GH/1UAdkB1j38qIw3TXAzg/8GPddcX5kRCcH5GJHc67N6r9WoH889h5CXA0av01OnHursNuAWZE7ftj4LWo71xF1L6+7mtL4vnZZqwsg9N72B+1vt/dBjAaxzGcg4jcIiLviUiDiBwDFgNDe3JNVQ0DB3H+GUSojlpuAvp38dwGPIHjuAE+666PxellV7lS2DGcH++iGK/v+J3AXR+J851oUNWj3TVKVVuAp4DPu3/N7wH+u7vnSQdE5Iuu3Blpi1k4909n99xoYL+qBmPsGwv8OOpcDTiOfGSH44bhONCNUceucrdHqFPV5ig7+4rIz11J5QTwBjCwMzmvA0OBXM71LzHvc1Vtchfjeq9nsnOvxPlyRBjjbgPH4U7s+AJxNPnf4/S8ilV1ILAC5wsFzq9vl68pIoLz5T3cffONGDwN3Ohqo3fhOPeDOL2ooao60H0UqOrMGK/v+J0A53tx2D3PYBEZ2AU7Yn0PlgGfA+YBTar6blfeUDrh6uH/CTyII3UNBMpw7p+Y95y7fYzEHuQ8iCOvDYx69FHVdzocVw+cBmZGHVeoqtHOtGOb/RUwFbhKVQuA6yNvo5PjO16vjXP9S1Lv80x27k8C3xWRYSIyFPhbIDLw+UvgyyIyT0QCIjJSRKbh/Brn4fxNDLqDIAujzlkDDBGRwk6u+RRwq3veHJwvUAuO3m/0ElWtw/m7/Cucv9zbVLUKeBn4FxEpcNtzoojcEOMUK4ApIvJZEckWkU8DM4AX3fOsBP7dHWzLEZHrY5wDnO/BWfHOrjMPA/9ChvbacTRzxbl/EJEv4/TcAX4BfENELheHSe6PwXqcsZNHRKSfiOSLyFz3NT8D/lpEZrrnKxSRT3W8qPsP+T+BfxWRIvfYkSKy6Dy2DsD5QTgmIoOBv+uw/5w2jrpeCOdeXyoiA9z38Zec8S9JIZOd+/eADcBmYAvwobsNVV0PfBlnEOc48DqOjtcI/BlOwx3F+eu/PHJCVd2O86Oxx/37N4IoVHUH8Hngpzi/7rcDt6tqa+LeZsbxBM6YyhNR276I88O8FafdfoejoZ+Fqh4BbsP50T2CM+B9m6rWu4d8AadHth2oBf68Exu+j9NxOCYi34ja/mvgIpJ8k/sFVd2K8+P2Lo5zvAh42933NLAUp90accayBruO8nYcLfsAcAhnIBZVfRZnEPS3rnRSBnQWdfItnMHX99xj1+L0zDvjRzgDq/U4gRerOuz/MXC3G+3ykxiv/xqOzr8HZyzoCeDx81wv7ogr6BuGkWBE5IvAfap6rde2GOlPJvfcDSNpiEhf4E+Bx7y2xcgMzLkbRoJxtd06HCniiQscbhhxwWQZwzCMNMR67sY5iMjjIlIrImVR234oItvFmTb/bBdDAg3D8Ahf9NyHDh2q48aN89oMA9i4cWM9zrTrk8CvVXUWgIgsBF5R1aCI/ABAVb91vnNZu/qHjRs31qvqsAsf2TWsbf3B+drVF9nPxo0bx4YNG7w2wwBEZL+qviEi46K3q+rLUavvAXdf6FzWrv5BRDrOvO0V1rb+4HztarKM0RO+gjOhxzAMn2LO3egWIvIdnAx6/9PJ/vtEZIOIbKirq0uucYZhtGPO3egyInIvzgzOz2kngzWq+piqlqpq6bBhcZN4DcPoJr7Q3A3/I06xkm/h5DBvutDxhmF4i/XcjXMQkSdx8n9MFadCzVeBR3GSKa1xU7b+zFMjDcM4Lxd07p3EPA8WkTUisst9HhS1769FpEJEdlwg65rhU1T1HlUdrqo5qjpKVX+pqpNUdbSqXuI+7vfaTuPCfOUrX6GoqAggVopj3AyMP3Hv2c3ilrYzUp+u9Nz/C7i5w7aHgXWqOhlY564jIjOAz+B8kW7GSY/aleT2hmEkgC996UusWtUxoeFZ3AJMdh/3Af+RDLuMxHNBzT1WzDNOWawb3eVlODm0v+Vu/61beWaviFTg1CPNuMIEfkVVOdDQxKaDx9haeYKHb5mGUzPESAbhsPL0xoMsvmg4A/JzEn6966+/nn379p3vkCU4k9UUJx3uQBEZ7uav7zKbDx1j7daaXlia+ogId18+itGD+3ptCtDzAdXiSOOralUkAT5OGan3oo47xLklrwAnZA6np8CYMWN6aIbRHWobm3no6c28vtMJUeyTk8WX546npDDfY8syh9d31fGt32/hZEuIr1473mtzwLk/D0atR+7Zc5z7+e7Zn6yrYO22GjK5n6AKwXCYhxZN89oUIP7RMrGattOQOdz0p6Wlpd7nQEhzXt1Ry0NPf0xjc5Bv3jyVG6YMY2rxALKzbEw9mTz7oVNp7cP9R/3i3ONyz7YEQ1wyeiDPPTA31kszgml/s5Jg2D+urKfOvSby101EhuNUpQHnV3901HGjOFOX1PCIX7y5h++9tI1pJQN44o/mMKV4gNcmZSQnW4K8vNWpi/zhgW7X2U4UcblngyElO5DB3XYgIIIPUnW109Nu23LgXnf5XuD5qO2fEZE8ERmPM0izvncmGr1hVVkV33tpG7fMKuG5B+aaY/eQlVuqaG4Lc8fsEVQdb6by2GmvTQLnnv2iGzUzBzjeXb0dHDkiO8uce9hHPfeuhELGinl+BFggIruABe46qlqOU190K07NwQfcGoiGB1Qfb+Zbv9/C7FGF/Ogzl5CfY4FLXvLsR4cZN6Qvf3idI8ds2J/43vs999zD1VdfDZAXuX9F5H4RiYSyrsCp81mBU0T6T3tynWBYyclwiU8EQj7qunclWuaeTnbN6+T4pTiFbg0PCYeVbzz9Ma3BMD/6zKXkZZtj95Kq46d5d88Rvj5vMtOHF9AnJ4sP9x/ljtkjLvziXvDkk08CICIfqmppx/1ulMwDvb2OyTLpI8sYPud/1h/grYp6/ua2GYwf2s9rczKe5zdVogp3XTqSnKwAs0cX+kl37zVtoXDGD84HBMI+8u6Z3RppSlNrkB+v3clV4wdzz5WjL/wCI+E899FhLhszkLFDnB/ay8cOorzyBE2tQY8tiw/BsPXcAyLm3I3Esuyd/dSfbOWbN0+1CUo+YGvlCbZXN3LXpWemfFw+dhChsLL50HEPLYsfobBmfM9dRPDReKo593TjRHMbP3t9NzdNHcblYwd7bY4BPLfpMNkB4daLz+jrl4520jFtTMKgajJoC4XJyfieuzMD3C+Yc08zfvnmXo6fbuMvF0z12hQDp0f7/KbD3Di1iMH9ctu3D+qXy8Rh/fgwTZx7MKQWCilCOOy1FWcw555GHD3Vyi/f2svNM0u4aFSh1+YYwHt7jlBzouUsSSbC5WMHsfHAUV/19npKMBwmK5DZ7iQrYJq7kSB+9c4+TrUG+YsFU7w2xXB55sPDDMjLZt70onP2lY4dzLGmNvbUn/LAsvjSFlJyMrznLoJp7kb8aQmGeOL9/XxiahFTS2wWqh843RpiVVkVt1xUEnMC2WVjXd19X+pLM6Gwkp3hPXcnzt0/3j2zWyONWLmlmvqTrXxp7jivTTFc1myr4VRriDsviZkYlQlD+zGwb05aDKq2hcIZ33MP+GyGqjn3NOGpDQcZM7gv104a6rUphsvzHx1meGE+cyYMibk/EBAuG+Po7qlOMGwDqgELhTTizcGGJt7ZfYS7Lx9lce0+oeFUK6/vrOOO2SMInCdE8PKxg6ioPcmxptYkWhdfVJVQWDN+QFVshqoRb37/4SFE4JOXj/LaFMPlpc2VBMPKkk4kmQiXjXF091RORdAWchyaxbmb5m7EkXBY+d3GQ8ydOJSRA/t4bY7h8tymSqYU92f68PMPbl8yeiDZAUlp3T3kahGZPkPV4tyNuPLe3iMcOnqaT5Var90vHGxoYuP+oyy5ZOQFZbI+uVnMHFHAhhSOmGlzPVqmD6iaLGPEleWbKumXm8WimSVem2K4PL/JKaW35JKupfO9bOwgPj50jLaQj7p93SDoyjKWOMwGVI04EQora7bWcNO0IivE4RNUlec2VXLFuEGMGtS3S68pHTuY5rYw5ZUnEmxdYgi6P0pZGS7LZAVMczfixMb9RzlyqtV67T5ia9UJKmpPXnAgNZrScYlLIrZq1SqmTp0KMEtEHu64X0RuFJHjIrLJffxtd68RKQptA6omyxhxYnV5NblZAW6cOsxrUwyX5ZsqnQyQFw3v8muKC/IZObAPG/c3xNWWUCjEAw88wMqVKwHKgXtEZEaMQ99U1Uvcx//t7nXaZZkM77lbyl8jLqgqq8urmTtpCAPyc+J6bhF5XERqRaQsattgEVkjIrvc50FxvWgaEA4ryz+u5IYpwxgUlQGyK5SOG8SGffFNIrZ+/XomTZrEhAkTABT4LbAkbhdwsQFVB+u5G3Fha9UJDh09nShJ5r+AmztsexhYp6qTgXXuuhHFB/saqDrezB1dHEiNpnTsIGobWzh09HTc7Dl8+DCjR59ViesQEEsvulpEPhaRlSIys7Pzich9IrJBRDbU1dW1bz8zoJrZ7sQqMRlxYXV5DQGB+TOK435uVX0D6KgRLAGWucvLgDvjfuEU5/mPK+mTk8WCHrTJpe5kpnhWZurkX0DHjR8CY1V1NvBT4LnznO8xVS1V1dJhw85IgZEon6yM19wtzt2IAy+XV1M6djBD++cl65LFqloF4D6fm8OWznt36U5rMMyKLVUsnFlM39zsbr9+1CBnAlrV8fj13EeNGsXBgwfP2gRURm9Q1ROqetJdXgHkiEi3EhRFJjFluixjce5Gr9l/5BTbqxtZODP+vfbe0lnvLt15Y2cdx5raOs0AeSEK++SQmx2gtrElbjZdccUV7Nq1i7179wII8BlgefQxIlIi7kwrEbkSxycc6c51gm53NdMHVJ30A15bcYbudzEMz1mztQYg2SGQNSIyXFWrRGQ4UJvMi/ud5z+uZFDfHK6d3LOsnCJCSUE+NSea42ZTdnY2jz76KIsWLQKYCfyDqpaLyP0Aqvoz4G7gT0QkCJwGPqPdHNW13DIOgcCZz8IPZPZPbYry+s46JhX1Z/Tgrk2SiRPLgXvd5XuB55N5cT9zqiXImq3V3HrxcHJ60XstLsiLq3MHWLx4MTt37gQoU9Wl4Dh117Gjqo+q6kxVna2qc1T1ne5ew0IhHWxA1egVp1tDvL+3gRumJE7yEJEngXeBqSJySES+CjwCLBCRXcACd93A+SfV3Bbmjtk9k2QiFBXkU3sifrJMsojIMjag6q84d5NlUoz39h6hNRjm+gQ6d1W9p5Nd8xJ20RRm+ceVjCjMp3Rs70L/iwfk8+qJWlQ1pfLyR3rumT6gGpBOI5Q8wXruKcbbu+rJzQ5w1fjBXptiAEdPtfLGzjpuv0BRjq5QUphHU2uIky3BOFmXHNoHVC3O3crsGT3nrYp6rhg3yBKF+YQVZVUEw8rts7s/cakjxQX5ANSkmDTTZj13wE0/YHHuRk+obWxme3Ujc61Oqm9YvqmSCcP6MXNEQa/PVTQg4tzjO6iaaCwU0sHSDxg95t3dTvixFcH2B9XHm1m/r4E7Zo+Ii0ZeXOBMSEs552753AH/xbn3yrmLyF+ISLmIlInIkyKSbwmmEsdbu+op7JPDzBGFXptiAC9urkQV7oiDJAOpK8sE28vsZbhzD6RJz11ERgJ/BpSq6iwgC2cGnCWYSgCqytsV9VwzcUjGh5z5heUfV3LRyEImDOsfl/P1y8tmQF52CvbcbUAVIil/08C5u2QDfUQkG+iLk7fCEkwlgL31p6g83tzjGZBGfNlXf4rNh47HrdceoSgBE5kSjQ2oOqSNLKOqh4F/Bg4AVcBxVX0ZSzCVEN6uqAdg7kRz7n7ghY+d/Fu3Xtz1ohxdoaQwvikIkoENqDpkpcuAqqulLwHGAyOAfiLy+a6+PlMTTPWUd/ccYeTAPowdktSUA0YnvLC5kivHDWbEwD5xPW/xgPzU1dwzXC702wzV3vzUzgf2qmqdqrYBzwDX4CaYArAEU/FBVXlvTwNXTRicUjMX05Ud1Y3srDnJ7bPj22sHNwVBYzNhP3mJC2DRMg7ppLkfAOaISF83Zeg8YBuWYCru7Ko9ScOpVuZMGOK1KQaOJBMQuKUbdVK7SklBHm0h5WhTa9zPnSiCVqwDcOPcffSj3OPcMqr6voj8DqeSSxD4CHgM6A885SabOgB8Kh6GZjLv7XHi26825+45qk6d1LmThiakUEp0OOSQ5BVi6RVtYSUnSzL+X6XfZJleJQ5T1b8D/q7D5hYswVRceXe3o7dHqvUY3rH50HEONDTx4CcmJeT8RQVnZqnOiMOs12QQCmvG99ohjeLcjeQQDivv721gzoQhGd8z8gMvfFxJTpYkrFBKSWHqpSBoC4XJyfAYd4ho7l5bcQZrEZ9zRm+3LJBeEw4rL26u4oYpRRT2yUnINYa5Ukx1nJz7qlWrmDp1KsAsETlnQqE4/EREKkRks4hc1t1rBEOa8bNTwVL+Gt0korfbYKr3bNh/lOoTzQmJkomQmx1gaP/cuIRDhkIhHnjgAVauXAlQDtwjIjM6HHYLMNl93Af8R3evEwyHMz7GHawSk9FN3nPj25NcUs+IwQsfV5KfE2D+9MQWJi8akE9tHHru69evZ9KkSUyYMAFAgd/izE2JZgnwa3V4DxgYCWXuKm0hzfj6qZBmA6pGYono7TdNjTnJ10giwVCYlWVVzJteTL+8xN42xQV5cZFlDh8+zOjRo6M3HQKu6nDYSOBgh2NG4sw6PwsRuQ+nd8+YMWPat/+fy0ZyzUT7Z+m3nrs5dx9jert/eG9PA/UnW7k9zukGYlFSmM+Wwyd6fZ5O9N+OG2N1uWO/UPUxnHBnSktL24+5xlJiABHN3WsrzmCyjI8xvd0/vLi5kn65WdyYhH9RRQPyOXKqhbZQ78r6jBo1ioMHD561CSe5XzSHgNEXOMboAoGAv3ru5tx9zPp9DYwozDe93WNag2FWllWzYEZxUsobFhfkowp1jb0bVL3iiivYtWsXe/fuBaeH/hmcGeTRLAe+6EbNzMFJAHiOJGNcGBEn5t8vmHP3KarK+r0NXGmFsD3n7Yp6jp9ui0ud1K5QUhifikzZ2dk8+uijLFq0CGAm8JSqlovI/SJyv3vYCmAPUAH8J/CnvbpoBuO3lL+mufuU/UeaqGts4Qpz7p7zwseVFORnc93k5GQvPVNLtffhkIsXL2bx4sWISJmqLgVQ1Z9F9qsjzD/Q6wsZVkPV6Brr9zYAcOU4c+5e0twW4uWtNSyaWUJudnJul+KC1JulavgvWsacu09Zv6+Bwf1ymVQUnxJuRs94fWcdJ1uC3JYkSQZgSL9csgNizj3FsPQDRpdYv7eBK8YN8l0+mVhF0b22KZG8tLmKQX1zkhrHHQgIRQPyUq5oR6YTmcfllxQE5tx9SPXxZg40NHGFzySZ8xRFT0tOt4ZYu62Gm2eVkJPk6fVFBalXbi/TyXI7Yn7pvZtz9yHr9zl6+1XjfRnfHqsoelry2o5amlpD3HZx8iSZCCXm3FOOQCDi3P3h3c25+5D1e4/QLzeL6cMHeG3KWZynKHo76VT4/MXNVQztn8tVHkQsFRfkmXNPMSIKqjl3o1M+2HuUy8cN9l2mva4URU+XwudNrUHWbXckGS/aoaggnxPNQU63hpJ+baNnBCKyTO8mFscNf3kPg6OnWtlR0+hJb7ELdFYUPe1Yt62W5rawJ5IMOLIMxC+vu5F4AtZzN87Hhv1HAXw3mOrSWVH0tOOlzVUMG5DnWTtYrHvq0d5zN+duxOL9PUfIzQ5w8ahCr005B1V9H4gURd+C8/15zFOjEsDJliCv7qjl1ouGe1YbNF4pCIzkIT6LlrH0Az7jvb1HuGzMwKQkqOoJnRRFTyvWbauhJRjm1iSk9+2MIuu5pxwW5250yvHTbZRXnrAUvx7z4uYqSgryuXzMIM9sGJCXTd/cLJvIlEIEfNZzN+fuIz7Y24Cq5W/3ksbmNl7fWcctF5W0xy17gYhQXJBvA6ophA2oGp3ynqu3XzJ6oNemZCxrt9XQGvQuSiaa4oK8uNRSNZKDTWIyOsXvensm8NLmKkYU5nOpD35greeeWkRkGZ/4dhtQ9QsRvf3r8yZ7bUrGcqK5jTd21vPFq8d6KslEcFIQtKCq3U4g19DQwKc//Wn27dvHuHHjwMkDdA4isg9oBEJAUFVLe2d15mKyjBET09u9Z+3WGlpDYRZ7GCUTTVFBPq3BMMea2rr92kceeYR58+axa9cu5s2bB1BynsNvUtVLzLH3jsgPsF9K7Zlz9wmmt3vPS5urGDmwjy8kGejdLNXnn3+ee++9FyDy7F3oT4bgN1nGnLtPML3dW46fbuONXXUsvqjENzn0ezORqaamhuHDnX8g7nNnEqwCL4vIRhG573znTKekcInAb7KMae4+wPR271m7tYa2kLL4In9IMhBdSzW2c58/fz7V1dXnbF+6dGl3LjNXVStFpAhYIyLbVfWNWAeq6mO4M5JLS0v94cF8hN/i3Hvl3EVkIPALYBZOD+ArwA7gf4FxwD7gD1T1aG+uk+5E9Haf5m/PCF7a4kgyfpLFigoiPffYE5nWrl3b6WuLi4upqqpi+PDhVFVVAQRjHaeqle5zrYg8C1wJxHTuxvlJt5S/PwZWqeo0YDZOEqmHgXWqOhlY564b5+HdPUfIyw5w6ZiBXpuSkRw/3cabPpNkAPKysxjcL7dHmvsdd9zBsmXLACLPxzoeIyL9RGRAZBlYCJT1wuSM5ozmnuLOXUQKgOuBXwKoaquqHsPJ973MPWwZcGfvTEx/3tl9hMvHDjK93SP8KMlEKC7I79FEpocffpg1a9YwefJk1qxZA05xFURkhIisiJweeEtEPgbWAy+p6qo4mZ5xpJMsMwGoA34lIrOBjcDXgWJVrQJQ1SpXyzsHd/DmPoAxY8b0wozUpuFUK9uqTvCNhVO8NiVjWeFDSSZCcUFej3ruQ4YMYd26de3rIhKCdhlmsbu8B+cftxEHIjVd0kGWyQYuA/5DVS8FTtENCSZdKvb0lnd3HwHg6olDPbYkMznR3Mabu+q5ZZa/JJkIkYlMhv+RNKrEdAg45Ob4BifP92VAjYgMB3Cfa3tnYnrz9u56+udlM9uH+dszAb9NXOpIUUE+9SdbaAv5xGMYnZI2xTpUtRo4KCJT3U3zgK3AcuBed9u9wPO9sjDNeaeinqvG+69eaqawYot/csnEoqQgH1Woa7Teu99Jtzj3rwH/IyK5wB7gyzg/GE+JyFdxyrJ9qpfXSFsOHW1i35Emvnj1OK9NyUgiuWS+cPVYX0oy4Gju4MS6jxjYx2NrjPORTgOqqOomIFY+inm9OW+m8E6Fo7fPnWR6uxes2+ZKMj6MkolgtVRTh3SLczd6wVsV9Qztn8eU4v5em5KRvLS5muE+lmQASgojzt1kGb+TNnHuRu9QVd7ZfYS5k4b4VhJIZxqbnVwyt8wa7ov0vp0xuG8uOVlied1TAL/JMubcPWJnzUnqT7Yw10IgPeGV7bW0BsPcevH5MuF6TyAgFA3Ip+a4OXe/0z6g6hPvbs7dI96qqAfgmkmWT8YLXnKLYF862v+ZcIsL8qhpNOfud86U2fPYEBdz7h7x5q46Jgzrx6hBfb02JeM42RLktZ113DzL2yLYXaW4IJ9q67n7HtPcDZrbQry35wjXT87cmbleEpFk/BwlE02xzVJNCc7EuXtrRwRz7h6wcf9RmtvCXDfZ9HYvWLmliqIBeZSO9b8kA07EzMmWICdbYmbtNXxCe5k967lnLm/sqiMnS1KyXqqIDBSR34nIdhHZJiJXe21Td2hqDfLqjtqUkWTg7IlMhn/x2wxVc+4e8ObOei4fO4h+eSlZCCtWDv+U4dXtdTS3hbl5lr+jZKJpn8hkuruvMc09w6lrbGFr1QmuS0G9/Tw5/FOGFWVVDOmXm1JVr3pTKNtIHoE0ygpp9IC3KpzCwik6mBqdw/8jEfmFW8GnHT8XUT7dGuLV7bUsmlVCVopIMnCm527O3d9Y+oEM582d9Qzul8vMEQVem9ITLpjD3895+l/fWUdTa4hbUyRKJkK/vGwG5GVT242ImaeffpqZM2cSCATYsGFDp8eJyM0iskNEKkTESmL2ApuhmsGoKm/squfaSUNTZjCvA53l8E8JVpZVMahvDleNH+y1Kd2muLB7se6zZs3imWee4frrr+/0GBHJAv4NuAWYAdwjIjN6a2umEnC9qV8095Qc0UtVtlc3Un+yJWVDIFW1WkQOishUVd3BmRz+vqclGGLdtlpuu3h4SubOLynI75YsM3369K4cdiVQ4ZbbQ0R+i1MDOSXa1G9kWc89c3l9p6NBp+JgahSRHP6bgUuAf/TWnK7x1q56TrYEuSXFJJkIzkSmuGvuI4GDUeuH3G0x8fN4ih8Qn1Visp57Enl9Rx3TSga0p3FNRc6Tw9/XrNhSTWGfHK6ZmDpRMtGUFOZR29hCKKztg8Hz58+nurr6nGOXLl3KkiVLunLaWNpgp55JVR8DHgMoLS31hwfzEX6LczfnniROtgTZsL+Br1w73mtTMo7WYJg1W6tZOLOEnBSUZMDpuYfCypFTLRQNcDoHa9eu7e1pDwGjo9ZHAZW9PWmmkjY1VI3u8XZFPW0h5YYpKS3JpCTv7K7nRHOQW1Jo4lJHzkxkimuOmQ+AySIy3i2V+RmcGshGD7A49wzltR219M/L5opxqRepkeqs3FJN/7xsrk3RgWzo/kSmZ599llGjRvHuu+9y6623AkwGEJERIrICQFWDwIPAapyZxk+pann8rc8M/BbnbrJMElBVXt1ex3WTh6asLJCqBENhXt5azbzpReRlZ3ltTo+JjNN01bnfdddd3HXXXe3rIrILQFUrgcWR7aq6AlgRR1Mzlkh4s098u/Xck8G2qkaqTzRz07Qir03JON7f28DRpjZumZWaUTIRhvbPIyCWX8bP+G1A1Zx7Enh1Ry0AN041vT3ZrCyrok9OVsqPdWQFhGED8iwzpI+xGaoZyKvba7loZGF7lIORHEJhZVVZDZ+YVkSf3NSVZCJ0dyKTkVz8prmbc08wx5pa+fDAUW6yXnvS2bj/KPUnW1Iqve/5SNBEJiNOZFnK38zi9Z11hBVuNL096azYUkVediBtxjpKuplfxkguJstkGGu31TK0fy6XjBrotSkZRTisrC6v5vopw+ifmkVRzqG4IJ8TzUFOt4a8NsWIQcS5h3zi3c25J5DWYJjXdtQyb1pxqmaBTFk2HTpG1fHmlJ641JH2iUwmzfgScb3p6bZQe83bky1Bmlq9qX2bHl0an/LBvgYam4PMn1HstSkZx6qyanKyhHnT0+ezj57ING5ovwscbSSbHDfn7w9X7+CHq3ecte+Hd1/Mp0pHx3pZwjDnnkDWbK0hLzvAtZNSd2ZkKqKqrCyr4pqJQynsk+O1OXGjpNAKZfuZPrlZ/PvnLuPw0dNnbV+6YhsHO2xLBubcE4SqsnZbDddNHpoWYXipRHnlCQ42nOaBGyd5bUpcaS+3Z4OqvmVxjJTS31+5jbAHOrxp7gliR00jh46eZn4ayQKpwuryagICC9JMDhuQn0O/3CyLdU8xsgJCyIPwyF47dxHJcoslv+iuDxaRNSKyy30e1HszU4+1W2sA+MT09AjDSyVWllVz1fghDOmf57Upcae4IL9btVQN7wmIeDKxKR4996/jZJSL8DCwTlUnA+voUEA5U1izrZZLRg+0WalJpqK2kYrak9xyUfpEyURTbLNUU46ASOrJMiIyCrgV+EXU5iXAMnd5GXBnb66RitSeaObjg8fSThZIBVZucSoTLZqZns7dJjKlHlkBIeRBjvfe9tx/BHwTiDa9WFWrANznmLpEOtdjXLPNkWTmmSSTdFaWVXP52EHtg4/pRnFBPrWNzZ70BI2eERBv8s302LmLyG1Arapu7MnrVfUxVS1V1dJhw9Ir78rq8hrGDenL1OIBXpuSURw40sTWqhNpNXGpIyUFebSFlIamVq9NMbqI03NPIecOzAXuEJF9wG+BT4jIb4AaERkO4D7X9trKFOJEcxvv7q5n0cyS9mroRnJYWVYFpK8kA10Ph3z66aeZOXMmgUCADRs2dHqciOwTkS0isklEOj/Q6DFZgRQbUFXVv1bVUao6Dqf24iuq+nmcGoz3uofdCzzfaytTiFe21dIWUhamsYPxKyvLqrloZCGjB/f12pSEUexWZKptPL9znzVrFs888wzXX399V057k6peoqqlvbfQ6Ih4FC2TiElMjwBPichXgQPApxJwDd+yYksVJQX5XDp6oNemZBRVx0+z6eAxHlo01WtTEkp7CoILFMqePn16MswxukCWeCPLxMW5q+prwGvu8hFgXjzOm2qcagny+s467rlyjCUKSzKry5womXTJ3d4ZwwbkIdL1WqpdQIGXRUSBn6vqY50dKCL3AfcBjBkzJl7XT3u8ipax9ANx5NUdtbQEw2k9oOdXVpVXM7moPxOH9ffalISSkxVgaP88ao43M3/+fKqrq885ZunSpSxZsqSrp5yrqpUiUgSsEZHtqvpGrANdx/8YQGlpqYXrdJFAwJtoGXPucWRlWTVD++dSOm6w16YkFBHJAjYAh1X1Nq/tOXKyhfV7G3jwpvTKJdMZkXJ7a9eu7fW5VLXSfa4VkWeBK4GYzt3oGak8Q9UAmttCvLq9lkUzS8hKf0mm46xkT3l5aw1hhUUZ8o8pXuX2RKSfiAyILAMLgbJen9g4C680d3PuceK1HXU0tYa4Zda5WeHSiU5mJXvKqrJqxgzuy4zhBV6bkhSKC/Iu6NyfffZZRo0axbvvvsutt94KMBlAREaIyIrIqYC3RORjYD3wkqquSqDpGUnAo1BIk2XixIubKxncL5c5E9JbkuHMrOSYM7SSPeh2/HQb7+yu58tzx2fMvIKSgnyONrXR3BYiPyd2Oum77rqLu+66q31dRHZBuwyz2F3eA8xOvMWZjfXcU5im1iDrttWy+KISsrPS9yPtyqzkZM88fmV7DW0hTfsomWjaY90tO2RK4PTcPbhu8i+ZfqzZWsPpthC3XzzCa1MSTWezkj1jVVk1xQV5GVWAPLrcnuF/AkLqZYU0HF74uIrigjyuSPMomfPMSvaEplZnXsGimSUZNa+gpNCceyqRssU6Mp3jTW28vrOW2y4ekVEOxg+8sbOO5rZwRkkycCa/TI2l/k0JAqk8QzWTWb21mraQcvvstJdkziJ6VrJXrCqrZlDfHK5M839MHSnIz6ZPjpXbSxVSLnGY4fDCx5WMHtyH2aMKvTYlo2gJhli3rZYFM4rTehA7FiJCcUGeOfcUwdHcPbhu8i+ZPhw52cI7u49w+8UjMiYMzy+8s/sIjS3BjJNkIji1VM25pwIBMc095XhpSxWhcOZJMn5gdVk1/fOymTtpqNemeEJJodVSTRWyAilYQzXTefajw0wrGcD0DJkZ6RdCYWXN1hpumlZEXnbsSTzpTklBPjUnWlAPeoRG97BomRRjX/0pPjpwjDsvHem1KRnHB/saOHKqNaOzbxYX5NMaDHO0qc1rU4wL4CQO8+C6yb9kevDMh4cQgSWXmCSTbFaVVZOXHeCGKelVe7c7tMe6Wzik77FJTClEOKz8/sPDXDtpKMML+3htTkahqrxcXs31U4bRLy9zI3nbY91Nd/c9qVggO2N5f28Dh4+d5u7LR3ltSsax+dBxKo83c3OG16i1Waqpg+VzTyF+/+Eh+udls3BGZjsYL1hVXk12QJg3vchrUzylKFJuz2QZ32M99xShqTXIyi1V3HrRcPrkZmakhleoKqvKqrl64hAG9s312hxPyckKMKTfhfO6G95jPfcUYcWWak61hvikSTJJZ1ftSfbWn2JRhksyEc43S/Whhx5i2rRpXHzxxZG87jF7IiJys4jsEJEKEXk4geZmLJbyN0V4cv0BJgzrxxXjBnltSsaxqqwaEVg4o9hrU3xBSUF+p7LMggULKCsrY/PmzUyZMgXgnF9EtxbuvwG3ADOAe0RkRgJNzkiyBJNl/M7OmkY27j/KZ64YbekGPGBVWTWXjRlEkRspkukUF3ZeS3XhwoVkZzvRRHPmzAGIpWNdCVSo6h5VbcXJ0b8kMdZmLgHT3P3Pb9cfJCdL+ORlJskkm4MNTWytOsGimdZrjxApt9cSDJ33uMcffxzgeIxdI4GDUeuH3G0xEZH7RGSDiGyoq6vrgcWZSZZHmnvmBgp3k+a2EM98dIiFM0sY0j/Pa3MyjtXl1QDcPDO9C5B3h//86y9TufcgFz/fj5yozJhLly5lyZIl7ctuD74hxili/f3s1Aup6mPAYwClpaWW96CLeDWgas69i6wur+ZYUxv3XJH4os/Guawqq2b68ALGDOnrtSm+4dHfPMu9j6/nN/dfHbMK2LJly3jxxRdZt24d/fr1i3WKQ8DoqPVRQGVirM1cHFnGg+sm/5KpyZPrDzB6cB+umTjEa1MyjtrGZjYeOJrxE5c60l5LNcag6qpVq/jBD37A8uXL6du30x/ED4DJIjJeRHJxSicuT5C5GUtWAOu5+5U9dSd5b08DDy2aaqX0PGDN1hpUydjc7Z1Rcp4UBA8++CAtLS0sWLAgsmkMgIiMAH6hqotVNSgiDwKrcUIlH1fV8mTYnklkWZk9//I/7x8gOyB8ymLbPWFVWTXjhvRlSnF/r03xFQV9ssnPCcTsuVdUVJy1LiIHAFS1Elgc2a6qK4AVibU0swlYmT1/cro1xNMbDnLzrBILwfOA401tvLv7CItmlVj4aQdExIl1t1mqviYgKVasQ0RGi8irIrJNRMpF5Ovu9sEiskZEdrnPKT3b55mPDnGiOcgXrx7ntSkZySs7agiG1WaldkJxQeex7oY/SMViHUHgr1R1OjAHeMCd3fYwsE5VJwPr3PWURFX5r7f3MXNEgc1I9YhVZdUUF+RxyaiBXpviS6zcnv9xeu4eXLenL1TVKlX90F1uBLbhTIBYAixzD1sG3NlLGz3j7Yoj7Ko9yZfnjjdJwANOt4Z4fWcdC2eU2EB2J1i5Pf+TFSDleu7tiMg44FLgfaBYVavA+QEAYuZmTYXZbr96ey9D++dy+2ybOOMFr++so7ktbFEy58HK7fmflM0KKSL9gd8Df66qJ7r6OlV9TFVLVbV02DD/lUvbW3+KV3bU8tmrxmZsEWavWV1ezcC+OVw5/twJOoaDldvzPwERVEn6v6teOXcRycFx7P+jqs+4m2tEZLi7fzhQ2zsTveE/39xDTlaAz8+xGale0BoMs25bDfOmFZ81td44Gyu353+yXEkx2bHuvYmWEeCXwDZV/X9Ru5YD97rL9wLP99w8b6htbOZ3Gw/xyctGUTTAwh8jdBYhlQje23OEE81Bk2QugJXb8z/tzj3JPffeTGKaC3wB2CIim9xt3wYeAZ4Ska8CB4BP9cpCD/jV2/sIhsL88fUTvDbFb0QipD4UkQHARhFZo6pb432h1eXV9M3N4rrJQ+N96rTCyu35n4AbjJFs2b3Hzl1V3yJ2VjmAeT09r9c0Nrfxm/f2c8us4YwbGjPZUsbiDpBHBssbRSQSIRVX5x4OKy9vreHGqcPIz7HxjvNh5fb8TyTQK2VkmXTlifcP0Ngc5P4bJnptiq/pECEVvb3XUVAfHTxKXWOLTVzqIiWFnZfbM7zHK1nGnHsULcEQv3xrL9dOGspFowq9Nse3nC9CKh5RUKvLa8jJEm6aFjOK1ujA+crtGd4TkWWSnYLAnHsUT204RG1jC39yo/XaO6OTCKm4oaqsKqtm7qShFOTnxPv0aYmlIPA3KRctk260BsP8x6sVXD52kOVs74TzREjFje3VjRxoaDJJphtEyu01t52/3J7hDRHNPdm5w8y5u/z2gwNUHm/m6/MmW6qBzolESH1CRDa5j8UXelF3WFVWjQjMn261UrtKsRsOWXuixWNLjFhEUmcke5aq5XPHyWHy01cquHL8YAu9Ow8XiJCKC6vLqykdO4hhA6xObVcZHhXrbmUI/UeWeCPLmHMHlr27j7rGFv79c5dZr91DDhxpYnt1I9+9dbrXpqQU7eX2onT3hx56iBdeeIHc3FwmTpwITqWlcxCRfUAjEAKCqlqaaHszjYBp7t5wormNn72+mxunDotZZNhIHqvLqwFMb+8mEVmmJipiZsGCBZSVlbF582amTJkCcL4P9SZVvcQce2LISrVJTOnCv71awbGmNr6xcKrXpmQ8q8qrmTG8gNGDTVroDgPysumbm3VWz33hwoXty3PmzAHITb5lBkDA7UJbnHsS2X/kFL96ax93Xz6KWSMtrt1Lahub+fDAUeu194ALldt7/PHHAY538nIFXhaRjSJy3wWu4/s03X4kYJp78vn+iu1kZwkPLbJeu9es2VqDKiyaZVEyXWX+/PlUVztS1qGjp9moyuv/2JelS5eyZMkSAJYuXUp2djZAQyenmauqlSJSBKwRke2q+kasA1X1MeAxgNLSUqsO0kWyLFomuby7+wiryqv5xsIp7WlTDe9YXV7D2CF9mVo8wGtTUoa1a9e2L//F/27ig30NvPWtT7RvW7ZsGS+++CLr1q2jX7/YeZJUtdJ9rhWRZ4ErgZjO3egZXkXLZKQsEwor//DiVkYO7MMfXmeZH73m+Ok23qmoZ9HMEotW6iHFBfnURpXbW7VqFT/4wQ9Yvnw5ffvGHsMQkX5udk9EpB+wEChLls2ZgsW5J5En3t/P1qoT/PSeSy3roA94dXstwbCa3t4LSgryaA2FaTjVypD+eTz44IO0tLSwYMGCyCFjAERkBPALVV0MFAPPuj+o2cATqrrKC/vTmTO5ZZJ73Yxz7jUnmvmn1TuYO2kIt11stVH9wOryaooG5HHp6IFem5KyRBftGNI/j4qKirP2i8gBaJdhFrvLe4DZybU088iyaJnEo6p859kyWoNhvnfnRSYB+IDmthCv7ahjwYzi9r+vRvexcnv+xatomYxy7i9srmLtthr+auEUxlshDl/w5q56TreFrJxeLzlTKNvyy/gNr6JlMsa51zW28PfLy5k9qpCvzB3vtTmGy+ryagrys5kzwTJx9oZh/fMIiNVS9SOWzz2BhMPKXz61iZMtQf7p7tlkZ2XE2/Y9wVCYtdtqmDe9mBxrk16RnRVgaP+8s1IQGP6gXZaxnnv8eezNPby5q56/vW0GU0ssjtovrN/bwLGmNhbNtIlL8aCksPNZqoZ3tMsySY6WSXvn/uGBo/zz6h0svqiEz101xmtzjChWl1eTlx3g+ik9K8dnnI1VZPInFi2TAI41tfK1Jz6ipDCf7/+fiy06xkeoKi9vreH6KcPom5txEbkJ4Xz5ZQzvaNfczbnHh7ZQmAef+Ii6xhZ+es+lFPaxepx+YvOh41Qdb7aJS3GkpDCfY1Zuz3fYgGocceLZt/BWRT3fu3MWl44Z5LVJRgdWl1eTFRDmTy/y2pS0wWLd/YkVyI4jP163i6c2HOLPPjGJP7hitNfmGDFYXV7NnAmDGdjX0ozHi/aKTBYx4ytMlokTP3t9Nz9au4tPXjaKv1gwxWtzjBhU1J5kd90pk2TiTEmhU3fWdHd/cabnntzrps1Ilqryo7W7+PG6Xdx28XAe+aSlF/ArkXJ6C2eYc48nJsv4k0hWDcsK2QNUlUdWbufnb+zh7stH8YNPXtz+a2n4j9Xl1cwePbB9yrwRHwbk59AvN8tSEPgMr1L+prwsc7IlyNee/Iifv7GHz88Zwz+ZY/c1lcdOs/nQcZu4lCCKCy3W3W94VawjpXvuWw4d52tPfsiBhia+dfM07r9hgkkxPudlV5IxvT0xlBTkU3X8tNdmGFF4FS2Tks79WFMrP163i/9+dz/DBuTx5B/N4SpLPJUSrC6vYVJRfyYO6++1KWlJcUE+6/d2Vi7V8IKILJNkVSZxzl1EbgZ+DGThVH55pLfnPNkS5LfrD/DTVypobG7j01eM5puLpjGon4XTJYvetOvRU62s39fA/TdYacNEEUlB8N3vfpfly5cTCAQoKioCiDmLLxH3qXE2EZU42ekHEuLcRSQL+DdgAXAI+EBElqvq1u6cpy0UZv+RU2ytauTl8mpe2V5LU2uI6yYP5Tu3TmdaSUEizDc6obftunZbDaGwcvNMq4CVKEoK8giGla888HW+973vAfCTn/yENWvWnPOhx+s+Nc5PumnuVwIVbhkvROS3wBKgS1+aOx59i2NNbVQeO03Q/UCG9MvlzktHcvflo7jMZpx6Ra/adXV5DSMH9mHWSPtRThSRCKQmzvybPXXqVGeH96o9ja4RkWV+vG4Xy97Z16tz/fn8KdzaxfKgiXLuI4GDUeuHgKuiDxCR+4D7AMaMOTtb46Rh/QmrctvFw5lU1J/JRQOYPnyA5WH3nl6165Ti/lw6ZqANeieQCcP6s/iiEnKzA3znO9/h17/+NYWFhQCVMQ6/YHtGc762NTpnSL9cvnrt+LgMdBf06brLTpRzj3X3nvWfRFUfAx4DKC0tPWvf//v0JQkyy+glvWrXb948LXGWZSDz58+nurr6nO1Lly5lSvHlLF26lKVLl/L973+fb3/727GS+FywPc/acZ62NTpHRPib22Yk/bqJcu6HgOikLqOI3XMwUgtrVx+xdu3aLh332c9+lm9/+9uxtExrzzQmUTrHB8BkERkvIrnAZ4DlCbqWkTysXVOEXbt2tS8vX74cIJYmYO2ZxiSk566qQRF5EFiNE2L1uKqWJ+JaRvKwdk0dHn74YXbs2EEgEGDs2LHgausiMgIn5HGxtWd6k7A4d1VdAaxI1PkNb7B2TQ1+//vfn7UuIm0AqloJLI5st/ZMXyz8xDAMIw0x524YhpGGmHM3DMNIQ8y5G4ZhpCGiyU5VFssIkTpgf4fNQ4F6D8xJNH5/X2NVdVg8TmTt6ivi1q4Qs239/v57it/fV6ft6gvnHgsR2aCqpV7bEW/S9X11lXR9/+n6vrpKur7/VH5fJssYhmGkIebcDcMw0hA/O/fHvDYgQaTr++oq6fr+0/V9dZV0ff8p+758q7kbhmEYPcfPPXfDMAyjh5hzNwzDSEN859xF5GYR2SEiFSLysNf29AYReVxEakWkLGrbYBFZIyK73OeMqBlo7Zq+pEvbplu7+sq5RxXsvQWYAdwjIskvYRI//gu4ucO2h4F1qjoZWOeupzXWrulLmrXtf5FG7eor505UwV5VbQUiBXtTElV9A2josHkJsMxdXgbcmUybPMLaNX1Jm7ZNt3b1m3OPVbB3pEe2JIpiVa0CcJ9j1bZMN6xd05d0b9uUbVe/OfduFew1UgZr1/TF2tan+M25Z0LB3hoRGQ7gPtd6bE8ysHZNX9K9bVO2Xf3m3DOhYO9y4F53+V7geQ9tSRbWrulLurdt6rarqvrqgVPfcSewG/iO1/b08r08CVQBbTg9nK8CQ3BG3Xe5z4O9ttPa1drV2jb92tXSDxiGYaQhfpNlDMMwjDhgzt0wDCMNMeduGIaRhphzNwzDSEPMuRuGYaQh5twvgIgMFJE/dZdHiMjvvLbJ6D3WrumJtesZLBTyAojIOOBFVZ3ltS1G/LB2TU+sXc+Q7bUBKcAjwEQR2YQzkWG6qs4SkS/hZIjLAmYB/wLkAl8AWoDFqtogIhNxUqIOA5qAP1LV7cl+E8Y5WLumJ9auEbyeReX3BzAOKIux/CWgAhiA80U4Dtzv7vtX4M/d5XXAZHf5KuAVr9+TPaxd0/Vh7XrmYT333vGqqjYCjSJyHHjB3b4FuFhE+gPXAE+LtCfPy0u+mUY3sXZNTzKqXc25946WqOVw1HoY57MNAMdU9ZIk22X0DmvX9CSj2tWiZS5MI85fuW6jqieAvSLyKQBxmB1P44weY+2anli7uphzvwCqegR42y2a+8MenOJzwFdF5GOgnBQtQZZuWLumJ9auZ7BQSMMwjDTEeu6GYRhpiDl3wzCMNMScu2EYRhpizt0wDCMNMeduGIaRhphzNwzDSEPMuRuGYaQh/x+/2THlMBtctQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
}
],
"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.8.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}