{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Non-linear optimization in Python\n",
"\n",
"** *\n",
"*Note: if you are visualizing this notebook directly from GitHub, some mathematical symbols might display incorrectly or not display at all. This same notebook can be rendered from nbviewer by following [this link.](http://nbviewer.jupyter.org/github/david-cortes/datascienceprojects/blob/master/optimization/nonlinearopt.ipynb)*\n",
"\n",
"This IPython notebook tries to showcase the usage of different optimization libraries, by finding the optimal solution of a real-valued toy function, subject to non-linear equality and inequality constraints of the variables, using different solver interfaces available for Python.\n",
"\n",
"The function to minimize and the constraints are non-linear, thus the most typical algorithms such as simplex or gradient descent are not applicable.\n",
"\n",
"\n",
"Some of the algorithms shown here are based on jacobians and/or hessians (matrices of first and second derivatives of the function to be minimize and its non-linear constraints), which help speed up things especially in large problems (in terms of having many variables, constraints, or functions over data points), but most of the methods are an overkill for this simple example.\n",
"\n",
"The choice of method and information passed to solver can, however, make a big impact on speed in bigger problems (such as calculating the parameters of some predictive model for a large dataset).\n",
"\n",
"\n",
"The file, unlike all the other in this page, was produced using Python 2.\n",
"** *\n",
"## Examples\n",
"\n",
"[1. Derivative-free solver with NLopt](#p1)\n",
"\n",
"[2. Gradient-based solver with SciPy](#p2)\n",
"* [Obtaining the derivatives](#p21)\n",
"* [Invoking the solver](#p22)\n",
"\n",
"[3. KKT-based solver with CVXOPT](#p3)\n",
"* [Obtaining the jacobian and hessian](#p31)\n",
"* [Invoking the solver](#p32)\n",
"\n",
"[4. Interior-point solver with Ipopt](#p4)\n",
"* [Ipopt through pyipopt](#p41)\n",
" * [Obtaining the gradient, jacobian and hessian](#p411)\n",
" * [Invoking the solver](#p412)\n",
"* [Ipopt through CasADi](#p42)\n",
"\n",
"[5. Alternative interior-point method wtih NLopt](#p5)\n",
"\n",
"[6. Benchmarking algorithms - solving a bigger problem](#p6)\n",
"* [Graphical comparison of solvability and time it took](#p61)\n",
"\n",
"** *\n",
"\n",
"## Mathematical formulation\n",
"** *\n",
"The goal is to solve the following mathematical problem:\n",
"\n",
"$$ {min \\: x_1 x_4 ( x_1 + x_2 + x_3) + x_3} $$\n",
"\n",
"$$ {s.t. \\:}$$\n",
"$$ {x_1 x_2 x_3 x_4 \\ge 25} $$\n",
"\n",
"$$ {x_1^2 + x_2^2 + x_3^2 + x_4^2 = 40} $$\n",
"\n",
"$$ {1 \\le x_1, x_2, x_3, x_4 \\le 5} $$\n",
"\n",
"(i.e. to find the values of $x_1, x_2, x_3$, and $x_4$ that minimize the function while respecting the constraints)\n",
"\n",
"the suggested starting point for algorithms is $ x_1=1, x_2=5, x_3=5, x_4=1 $\n",
"\n",
"This corresponds to problem 71 of the [Hock-Schittkowsky test suite](http://apmonitor.com/wiki/index.php/Apps/HockSchittkowski), available [here](http://apmonitor.com/wiki/uploads/Apps/hs071.apm).\n",
"** *\n",
"Note that while the function is non-convex over $\\mathbb{R}^4$, all the variables are bounded to be between 1 and 5, and in that domain the function is convex and thus can be approached with a variety of methods."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Proving that the function is not convex over R^4\n",
"\n",
"# a function is convex if, for any number a in [0,1], it satisfies this:\n",
"# f(a*x1 + (1-a)*x2) <= a*f(x1) + (1-a)*f(x2)\n",
"\n",
"import numpy as np\n",
"\n",
"f=lambda x: x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"a=np.float(0.2)\n",
"x1=np.array([1,-100,2,3],dtype='float64')\n",
"x2=np.array([-1,100,2,3],dtype='float64')\n",
"f(a*x1 + (1-a)*x2) <= a*f(x1) + (1-a)*f(x2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 1. Derivative-free solver (library NLopt)\n",
"** *\n",
"This problem can be easily solved by the algorithm COBYLA (Constrained Optimization BY Linear Approximations), available in the [NLopt library](http://ab-initio.mit.edu/wiki/index.php/NLopt) (among others), which comes with an interface to Python. This library can be easily installed from [conda-forge](https://anaconda.org/conda-forge/nlopt) (if using the Anaconda distribution for Python).\n",
"\n",
"The library includes both derivative-based and derivate-free methods, and all of its algorithms require the functions and constraints to be passed as evaluable functions with two arguments: the variables and the gradient (not used here, and passed to me modified in-place in the function when used), even when the gradients are not used.\n",
"\n",
"It can throw errors depending on the type of numbers used and outputted by functions (e.g. integers vs. floating point numbers) and the tolerances used for the objective and constraints, so oftentimes it helps to set these to values other than default. This library, as most Python-based solvers, uses NumPy's data structures.\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 1. 4.74299965 3.82114997 1.3794083 ]\n",
"Objective value: 17.0140172892\n",
"Result code: 4\n"
]
}
],
"source": [
"import nlopt, numpy as np\n",
"\n",
"############ defining the problem ##########\n",
"\n",
"# function to minimize\n",
"f0=lambda x,grad: np.float(x[0]*x[3]*(x[0]+x[1]+x[2])+x[2])\n",
"\n",
"# greater-or-equal inequality constraint\n",
"## NLopt requires inequalities to be expressed as f(x)-c<=0, thus this one first has to be flipped\n",
"f1=lambda x,grad: np.float(25.0-x[0]*x[1]*x[2]*x[3])\n",
"\n",
"# euqality constraint\n",
"f2=lambda x,grad: np.float(x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2 - 40.0)\n",
"\n",
"# upper and lower bounds for the variables\n",
"lb=np.array([1]*4, dtype='float64')\n",
"ub=np.array([5]*4, dtype='float64')\n",
"\n",
"# suggested starting point\n",
"x0=np.array([1,5,5,1], dtype='float64')\n",
"\n",
"\n",
"############ passing it to the solver ##########\n",
"\n",
"# initializing the optimization problem with 4 variables\n",
"opt = nlopt.opt(nlopt.LN_COBYLA,4)\n",
"\n",
"# setting the function to minimize\n",
"opt.set_min_objective(f0)\n",
"\n",
"# adding the constraints\n",
"opt.add_inequality_constraint(f1)\n",
"opt.add_equality_constraint(f2)\n",
"\n",
"# setting the bounds for the variables\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"# setting the numerical tolerance for the objective\n",
"opt.set_xtol_rel(1e-10)\n",
"\n",
"# running the solver\n",
"x = opt.optimize(x0)\n",
"minf = opt.last_optimum_value()\n",
"print \"Optimal x:\", x\n",
"print \"Objective value:\", minf\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 2. Gradient-based solver (library SciPy)\n",
"** *\n",
"The problem, being a smooth, differentiable function, can be solved more easily when using the partial derivatives of the function and its non-linear constraints with respect to the variables to minimize (also known as gradient or jacobian), using algorithms such as ‘Sequential Least-Squares Programming’ (SLSQP), available in many libraries such as in SciPy's optimization module. A different implementation of the same algorithm is also available in many other libraries, such as the NLopt library used above.\n",
"\n",
"\n",
"The SciPy library comes installed by default with Python’s Anaconda distribution for scientific computing."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Obtaining the derivatives\n",
"\n",
"The indefinite derivatives of the function and its constraints can be obtained manually or calculated with symbolic mathematics modules, such as Python’s SymPy (also installed by default in the Anaconda distribution), as shown here:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[x[0]*x[3] + x[3]*(x[0] + x[1] + x[2]), x[0]*x[3], x[0]*x[3] + 1, x[0]*(x[0] + x[1] + x[2])]\n"
]
}
],
"source": [
"from sympy import symbols, diff\n",
"\n",
"# This time, the functions variables will be defined as symbolic math\n",
"x=symbols('x[0] x[1] x[2] x[3]')\n",
"\n",
"# function to minimize\n",
"f0=x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"\n",
"# constraints\n",
"f1=x[0]*x[1]*x[2]*x[3]\n",
"f2=x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2\n",
"\n",
"# derivatives of the function and the constraints\n",
"d_f0_symbols=[diff(f0,x[var]) for var in range(len(x))]\n",
"d_f1_symbols=[diff(f1,x[var]) for var in range(len(x))]\n",
"d_f2_symbols=[diff(f2,x[var]) for var in range(len(x))]\n",
"\n",
"# checking some of these derivatives\n",
"print(d_f0_symbols)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These SymPy symbolic functions can be converted to evaluable functions using lambdify or passing them as a string to python's eval function:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"First derivative w.r.t. x1: x[0]*x[3] + x[3]*(x[0] + x[1] + x[2])\n",
"\n",
"Evaluating first derivative at point (1, 5, 5, 1)\n",
"Hard-coded function: 12\n",
"Evaluate string from sympy: 12\n",
"Sympy's lambdify: 12\n",
"\n",
"Evaluating first derivative at point (1, 2, 3, 4)\n",
"Hard-coded function: 28\n",
"Evaluate string from sympy: 28\n",
"Sympy's lambdify: 28\n",
"\n",
"Evaluating first derivative at point (10, 10, 10, 10)\n",
"Hard-coded function: 400\n",
"Evaluate string from sympy: 400\n",
"Sympy's lambdify: 400\n",
"\n"
]
}
],
"source": [
"from sympy.utilities.lambdify import lambdify\n",
"\n",
"d_example=lambda x: x[0]*x[3] + x[3]*(x[0] + x[1] + x[2])\n",
"d_example1=lambda x: eval(str(d_f0_symbols[0]))\n",
"d_example2=lambdify(x, d_f0_symbols[0])\n",
"\n",
"# some points to try\n",
"sample_points=[(1,5,5,1),(1,2,3,4),(10,10,10,10)]\n",
"\n",
"print \"First derivative w.r.t. x1:\",str(d_f0_symbols[0])\n",
"print \n",
"for p in sample_points:\n",
" print \"Evaluating first derivative at point\",p\n",
" print \"Hard-coded function:\",d_example(p)\n",
" print \"Evaluate string from sympy:\",d_example1(p)\n",
" print \"Sympy's lambdify:\",d_example2(*p)\n",
" print"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The examples below fill use the first method (evaluating sympy outputs as strings), as, while not the most efficient, it's the most convenient to write."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Invoking the solver"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
" fun: 17.014017245572898\n",
" jac: array([ 14.57227022, 1.37940764, 2.37940764, 9.56415073])\n",
" message: 'Optimization terminated successfully.'\n",
" nfev: 5\n",
" nit: 5\n",
" njev: 5\n",
" status: 0\n",
" success: True\n",
" x: array([ 1. , 4.74299606, 3.82115467, 1.37940764])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.optimize import minimize\n",
"import numpy as np\n",
"\n",
"############ defining the rest of the problem ##########\n",
"\n",
"# gradient of the function\n",
"grad=lambda x: np.array(eval(str(d_f0_symbols)))\n",
"\n",
"# constraints\n",
"cons=[\n",
" {'type':'ineq','fun':lambda x: eval(str(f1))-25,'jac':lambda x: np.array(eval(str(d_f1_symbols)))},\n",
" {'type':'eq','fun':lambda x: eval(str(f2))-40,'jac':lambda x: np.array(eval(str(d_f2_symbols)))}\n",
"]\n",
"\n",
"# upper and lower bounds for the variables\n",
"lb=np.array([1]*4, dtype='float64')\n",
"ub=np.array([5]*4, dtype='float64')\n",
"\n",
"# suggested starting point\n",
"x0=np.array([1,5,5,1], dtype='float64')\n",
"\n",
"############ passing everything to the solver ##########\n",
"minimize(fun=lambda x: eval(str(f0)), jac=grad, constraints=cons, bounds=[(lb[var],ub[var]) for var in range(len(x))],\n",
" x0=x0, method='SLSQP')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 3. KKT-based solver (library CVXOPT)\n",
"\n",
"As the problem is convex in the interval to be solved, its minimum can be determined to be the point that satisfies the [Karush–Kuhn–Tucker conditions (KKT)](https://en.wikipedia.org/wiki/Karush–Kuhn–Tucker_conditions) for optimality, which in turn require solving some equations with its Jacobian(matrix of partial derivatives of the function and each non-linear constraint with respect to each variable) and Hessian (matrix of second-order partial derivates of the sum of the function and each constraint multiplied by a lagrange multiplier, with respect to each pair of variables)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Obtaining the gradient, jacobian and hessian\n",
"\n",
"These will be calculated using SymPy just like in the example above:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from sympy import symbols, diff\n",
"\n",
"x=symbols('x[0] x[1] x[2] x[3]')\n",
"z=symbols('z[0] z[1] z[2] z[3]') # these will be the lagrangian multipliers\n",
"\n",
"f0=x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"f1=-(x[0]*x[1]*x[2]*x[3])\n",
"\n",
"# note that CVXOPT doesn't non-linear equality constraints directly\n",
"# the original contrainf f_2(x)==40 is equivalent to these two combined:\n",
"# f_2(x)>=40 and f_3(x)<=40\n",
"f2=x[0]**2+x[1]**2+x[2]**2+x[3]**2\n",
"f3=-(x[0]**2+x[1]**2+x[2]**2+x[3]**2)\n",
"g=[f0,f1,f2,f3]\n",
"m=z[0]*f0+z[1]*f1+z[2]*f2+z[3]*f3 # this is the sum to differentiate for the hessian\n",
"\n",
"jac_symbols=[[diff(g[con],x[var]) for var in range(len(x))] for con in range(len(g))]\n",
"hess_symbols=[[diff(diff(m,x[var1]),x[var2]) for var1 in range(len(x))] for var2 in range(len(x))]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Invoking the solver"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"{'dual infeasibility': 0.07703721666218974,\n",
" 'dual objective': 17.077745556399535,\n",
" 'dual slack': 0.6117001389000475,\n",
" 'gap': 0.3472762664951544,\n",
" 'primal infeasibility': 0.08079169139753443,\n",
" 'primal objective': 15.95031728284468,\n",
" 'primal slack': 0.015319998875008295,\n",
" 'relative gap': 0.020335018187750183,\n",
" 'sl': <0x1 matrix, tc='d'>,\n",
" 'snl': <3x1 matrix, tc='d'>,\n",
" 'status': 'unknown',\n",
" 'x': <4x1 matrix, tc='d'>,\n",
" 'y': <0x1 matrix, tc='d'>,\n",
" 'zl': <0x1 matrix, tc='d'>,\n",
" 'znl': <3x1 matrix, tc='d'>}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from cvxopt import solvers, matrix\n",
"import numpy as np\n",
"\n",
"ncons=len(g)-1\n",
"x0=np.array([1,5,5,1],dtype='float64')\n",
"\n",
"# cvxopt requires defining a single function that outputs all the required evaluations\n",
"def F(x=None,z=None): #x is the vector of variables and z are the lagrande multipliers\n",
" \n",
" # when passed with no arguemnts (i.e. F()), it should return the number of non-linear constraints and astarting point\n",
" if x is None: return ncons,matrix(x0)\n",
" \n",
" # when the variables are off-bounds, it should return None\n",
" if min(x) < 1.0: return None,None\n",
" if max(x) > 5.0: return None,None\n",
" \n",
" # otherwise, it needs to output the function and constraints evaluated at x\n",
" f=matrix(np.array([\n",
" eval(str(f0)),\n",
" 25.0+eval(str(f1)),\n",
" eval(str(f2))-40.0,\n",
" 40.0+eval(str(f3))\n",
" ],dtype='float64'))\n",
" \n",
" # and along with it, the jacobian evaluated at the same point\n",
" Df=matrix(np.array(eval(str(jac_symbols)),dtype='float64'))\n",
" if z is None: return f,Df\n",
" \n",
" # if z is provided, it also needs to output the hessian avaluated at x,z\n",
" H=matrix(np.array(eval(str(hess_symbols)),dtype='float64'))\n",
" return f,Df,H\n",
"\n",
"######### passing this function to the solver ##########\n",
"# important to set tolerances, as otherwise it can get stuch\n",
"xsol=solvers.cp(F=F, kktsolver='ldl', options={'show_progress':False,'reltol':1e-4})\n",
"xsol"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [1.0, 4.808364383142631, 3.9238828088666087, 1.302127362224565]\n",
"Objective value: 16.5965081735\n",
"Z multipliers: [0.6117001389000475, 1.2668465747904871, 1.0880138014153908]\n"
]
}
],
"source": [
"print \"Optimal x:\",[var for var in xsol['x']]\n",
"print 'Objective value:',xsol['x'][0]*xsol['x'][3]*(xsol['x'][0]+xsol['x'][1]+xsol['x'][2])+xsol['x'][2]\n",
"print 'Z multipliers:', [var for var in xsol['znl']]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Note that this solution is quite different from what the other solvers threw and is infeasible, but close to the optimum_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 4. Interior-point solver (library Ipopt)\n",
"\n",
"The problem can also be solved using interior-point methods, which are also helped by the jacobian and hessian matrices. Ipopt provides many of these, but unfortunately doesn't come with a Python interface. It is possible, however, to use third-party interfacers to it such as CasADi or pyipopt.\n",
"** *\n",
"\n",
"### Ipopt through pyipopt (linux only)\n",
"\n",
"\n",
"This library can be easily installed in 64-bit Linux systems through conda forge, first by installing [ipopt](https://anaconda.org/conda-forge/ipopt) and then [pyipopt](https://anaconda.org/imperial-college-research-computing/pyipopt) (e.g. ‘conda install -c conda-forge ipopt=3.12.7’, then ‘conda install -c imperial-college-research-computing pyipopt=0.1’). In Windows, the quickest option is to create a virtual machine under software like [VirtualBox](https://www.virtualbox.org) and run a linux environment from there).\n",
"\n",
"If the installation through conda forge doesn't succeed, then in order to install this manually, the process would look like this:\n",
"\n",
"1. Download the [source code for ipopt](http://www.coin-or.org/download/source/Ipopt)\n",
"2. Uncompress the file, go to the folder ‘Third Party’, and within each of those subfolder (except HSL), execute the files named ‘get.package’ (e.g. go the the folger Blas and run the command ./get.Blas, repeat for the other folders)\n",
"3. Configure the program for the local user (i.e. run ‘./configure –prefix=/usr/local’, NOT the usual ./configure’), then do ‘make’ and ‘make install’.\n",
"4. Download and install the pyipopt module for interfacing ipopt through python, which requires downloading its [source code](https://github.com/xuy/pyipopt), uncompressing, and installing with setup tools by running the setup script with the build and install commands (i.e. first ‘python setup.py build’, then ‘python setup.py install’, second one as superuser). In 32-bit anaconda installations, it might require the installation of the libgcc library for anaconda (‘conda install libgcc’).\n",
"\n",
"Installing pyipopt in Python3 might require additional steps."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Obtaining the gradient, jacobian and hessian\n",
"\n",
"Obtained with SymPy, just like before"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from sympy import symbols, diff\n",
"import numpy as np\n",
"\n",
"x=symbols('x[0] x[1] x[2] x[3]')\n",
"\n",
"# note that ipopt requires defining the hessian slightly differently from cvxopt\n",
"# here the objective function is multiplied by a different variable than the constraints\n",
"d=symbols('d')\n",
"z=symbols('z[0] z[1]')\n",
"\n",
"f0=x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"f1=x[0]*x[1]*x[2]*x[3]\n",
"f2=x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2 # ipopt supports equality constraints\n",
"g=[f1,f2]\n",
"m=d*f0 + z[0]*g[0] + z[1]*g[1]\n",
"\n",
"grad=[diff(f0,x[var]) for var in range(len(x))]\n",
"jac_symbols=[[diff(g[con],x[var]) for var in range(len(x))] for con in range(len(z))]\n",
"hess_symbols=[[diff(diff(m, x[v1]), x[v2]) for v1 in range(len(x))] for v2 in range(len(x))]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Invoking the solver"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 1. 4.74299964 3.82114998 1.37940829]\n",
"Objective value: 17.014017142\n",
"Bound multipliers z_L: [ 1.08787122e+00 2.42877637e-10 3.22241253e-10 2.39607541e-09]\n",
"Bound multipliers z_U: [ 2.27272724e-10 3.53731307e-09 7.71167542e-10 2.51089040e-10]\n",
"Constraint multipliers: [-0.55229366 0.16146856]\n"
]
}
],
"source": [
"import pyipopt\n",
"\n",
"nvar=len(x)\n",
"lb=np.array([1]*4, dtype='float64')\n",
"ub=np.array([5]*4, dtype='float64')\n",
"ncon=len(g)\n",
"g_lb=np.array([25.0,40.0])\n",
"g_ub=np.array([np.inf,40.0])\n",
"\n",
"eval_f=lambda x: eval(str(f0))\n",
"eval_grad_f=lambda x: np.array(eval(str(grad)))\n",
"eval_g=lambda x: np.array([eval(str(f1)),eval(str(f2))])\n",
"\n",
"## ipopt requires passing the hessian and jacobian as a flat vector\n",
"## the indices of the lower triangular part to which each entry corresponds need to be passed too, these are easy to calculate\n",
"\n",
"# structure of the jacobian to output\n",
"str_jac=list()\n",
"for v1 in range(ncon):\n",
" for v2 in range(nvar):\n",
" str_jac.append((v1,v2))\n",
"str_jac=(np.array(str_jac,dtype='int32').T[0], np.array(str_jac,dtype='int32').T[1])\n",
"\n",
"# structure of the hessian to output\n",
"str_hess=list()\n",
"for v1 in range(nvar):\n",
" for v2 in range(v1+1):\n",
" str_hess.append((v1,v2))\n",
"str_hess=(np.array(str_hess).T[0], np.array(str_hess).T[1])\n",
"\n",
"# number of entries on each\n",
"nnzj=int((len(g))*len(x))\n",
"nnzh=int(sum([i+1 for i in range(len(x))]))\n",
"\n",
"# the jacobian and hessian can be passed with the argument 'flag' (true/false)\n",
"# when they are passed with a flag, they should output the structure that was calculated above\n",
"\n",
"def eval_jac_g(x,flag):\n",
" if flag==True:\n",
" return str_jac\n",
" else:\n",
" return np.array(eval(str(jac_symbols)), dtype='float64').reshape((nnzj,1))\n",
"def eval_h(x,z,d,flag):\n",
" if flag==True:\n",
" return str_hess\n",
" else:\n",
" out=list()\n",
" for v1 in range(len(x)):\n",
" for v2 in range(v1+1):\n",
" out.append(eval(str(hess_symbols[v1][v2])))\n",
" return np.array(out)\n",
"\n",
"# these are the number of non-duplicated non-zero entries in the jacobian and hessian\n",
"nnzj=(len(g))*len(x)\n",
"nnzh=sum([i+1 for i in range(len(x))])\n",
"\n",
"# suggested starting point\n",
"x0=np.array([1,5,5,1], dtype='float64')\n",
"\n",
"# solver doesn't support keywork argument, thus arguments need to be passed in this exact order\n",
"opt = pyipopt.create(nvar, lb, ub, ncon, g_lb, g_ub, nnzj, nnzh, eval_f, eval_grad_f, eval_g, eval_jac_g,eval_h,lambda x: True)\n",
"x_opt, zl, zu, constraint_multipliers, obj, status = opt.solve(x0)\n",
"\n",
"print \"Optimal x:\",x_opt\n",
"print \"Objective value:\",obj\n",
"print \"Bound multipliers z_L:\", zl\n",
"print \"Bound multipliers z_U:\", zu\n",
"print \"Constraint multipliers:\",constraint_multipliers"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"### Ipopt through CasADi\n",
"\n",
"In comparison to all the steps done before to get Ipopt running, passing this problem to it through CasADi is easier, as it does all the work of calculating gradients, jacobians and hessians instead of needing them provided. As such, here SymPy won’t be used (since casadi itself does the same).\n",
"\n",
"Getting the library tu run is also a lot easier, as it comes with Ipopt prepackaged. In Windows, the process would look like this:\n",
"\n",
"1. Download the [CasADi binaries](https://github.com/casadi/casadi/wiki/InstallationInstructions).\n",
"2. Unzip it, and add that directory to the Python's path variable - this can be accomplished temporarily within a Python session by executing something like this 'from sys import path;path.append(\"D:\\\\Downloads\\\\casadi-py35-np1.9.1-v3.1.1-64bit\")'.\n",
"\n",
"In Linux, it can be installed through [conda forge](https://anaconda.org/conda-forge/casadi) too ('conda install -c conda-forge casadi=3.1.1')."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [1, 4.743, 3.82115, 1.37941]\n",
"Objective value: 17.014\n",
"Constraint multipliers: [-0.552294, 0.161469]\n"
]
}
],
"source": [
"from casadi import SX, vertcat, nlpsol\n",
"import numpy as np\n",
"\n",
"# casadi uses symbolic variables just like sympy - this creates an x vector with 4 values\n",
"x=SX.sym(\"x\",4)\n",
"\n",
"f0=x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"f1=x[0]*x[1]*x[2]*x[3]\n",
"f2=x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2\n",
"\n",
"# note that constraints cannot be passed as a list, they need to be converted into a vector\n",
"g=vertcat(f1,f2)\n",
"\n",
"# the problem is created as an object wit the variables, objective and constraints\n",
"# note that, by default, it prints all the outputs of each iteration, which results in large printout\n",
"solver = nlpsol(\"solver\", \"ipopt\", {'x':x,'f':f0,'g':g}, {\"print_time\":False,\"ipopt\":{\"print_level\":0}})\n",
"\n",
"# same information given before\n",
"x0=np.array([1,5,5,1], dtype='float64')\n",
"lb=np.array([1,1,1,1], dtype='float64')\n",
"ub=np.array([5,5,5,5], dtype='float64')\n",
"g_lb=np.array([25,40], dtype='float64')\n",
"g_ub=np.array([np.inf,40], dtype='float64')\n",
"\n",
"\n",
"# invoking the solver with all these parameters\n",
"res=solver(x0=x0,lbx=lb,ubx=ub,lbg=g_lb,ubg=g_ub)\n",
"\n",
"print \"Optimal x:\" , res[\"x\"]\n",
"print \"Objective value:\" , res[\"f\"]\n",
"print \"Constraint multipliers:\",res['lam_g']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 5. Alternative interior-point method (NLopt)\n",
"** *\n",
"NLopt can solve problems in a manner similar to the barrier method from ipopt, by using the AUGLAG solver helped by local, less stringent solvers (either gradient-based or derivative free). In this toy example, doing this doesn’t bring any speed benefit, but it can speed things up in some situations:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 1. 4.80254926 3.8470158 1.35314458]\n",
"Objective value: 16.9042724845\n",
"Result code: 6\n"
]
}
],
"source": [
"from sympy import symbols, diff\n",
"import nlopt, numpy as np\n",
"\n",
"# Building things in numpy, just like before\n",
"x=symbols('x[0] x[1] x[2] x[3]')\n",
"\n",
"# function to minimize\n",
"f0=x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"\n",
"# constraints\n",
"f1=25.0 - (x[0]*x[1]*x[2]*x[3]) # nlp requires constraints as <=\n",
"f2=x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2 - 40.0\n",
"\n",
"# derivatives of the function and the constraints\n",
"d_f0_symbols=[diff(f0,x[var]) for var in range(len(x))]\n",
"d_f1_symbols=[diff(f1,x[var]) for var in range(len(x))]\n",
"d_f2_symbols=[diff(f2,x[var]) for var in range(len(x))]\n",
"\n",
"# function to minimize\n",
"def f0(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(d_f0_symbols)), dtype='float64')\n",
" return 1.0*x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]\n",
"\n",
"def f1(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(d_f1_symbols)), dtype='float64')\n",
" return 25.0 - (x[0]*x[1]*x[2]*x[3])\n",
"\n",
"def f2(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(d_f2_symbols)), dtype='float64')\n",
" return x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2 - 40.0\n",
"\n",
"# upper and lower bounds for the variables\n",
"lb=np.array([1]*4, dtype='float64')\n",
"ub=np.array([5]*4, dtype='float64')\n",
"\n",
"# suggested starting point\n",
"x0=np.array([1,5,5,1], dtype='float64')\n",
"\n",
"\n",
"############ passing it to the solver ##########\n",
"\n",
"# initializing the solver with AUGLAG\n",
"# it comes in different varieties, the one ussing derivatives has preffix 'LD_'\n",
"opt = nlopt.opt(nlopt.LD_AUGLAG_EQ,4)\n",
"\n",
"# setting the function to minimize\n",
"opt.set_min_objective(f0)\n",
"\n",
"# adding the constraints\n",
"opt.add_inequality_constraint(f1)\n",
"opt.add_equality_constraint(f2)\n",
"\n",
"# setting the bounds for the variables\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"# passing a local solver to AUGLAG - MMA wouldn't be able to solve this problem alone\n",
"local_opt=nlopt.opt(nlopt.LD_MMA,4)\n",
"opt.set_local_optimizer(local_opt)\n",
"\n",
"# setting timeouts can help prevent crashes from approximation problems\n",
"opt.set_maxtime(10)\n",
"\n",
"# running the solver\n",
"x_opt = opt.optimize(x0)\n",
"minf = opt.last_optimum_value()\n",
"print \"Optimal x:\", x_opt\n",
"print \"Objective value:\", minf\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Note that the solution, while close to the optimums that the other solvers found, is not feasible, but provides a reasonably good approximation._"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## 6. Benchmarking some algorithms with a larger problem\n",
"\n",
"Now solving problem problem 103 of the [Hock-Schittkowsky test suite](http://apmonitor.com/wiki/index.php/Apps/HockSchittkowski), available [here](http://apmonitor.com/wiki/uploads/Apps/hs103.apm).\n",
"\n",
"\n",
"$$ {min \\:\\: 10\\frac{x_1x_4^2x_7^{0.5}}{x_2x_6^3}+15\\frac{x_3 \n",
" x_4}{x_1x_2^2x_5x_7^{0.5}}+20\\frac{x_2x_6}{x_1^2x_4x_5^2}\n",
" +25\\frac{x_1^2x_2^2x_5^{0.5}x_7}{x_3x_6^2}}$$\n",
"\n",
"\n",
"\n",
"$$s.t.$$\n",
"\n",
"$$ { 1-0.5\\frac{x_1^{0.5}x_7}{x_3x_6^2}-0.7\\frac{x_1^3x_2x_6x_7^{0.5}}{x_3^2 }\n",
" -0.2\\frac{x_3x_6^{\\frac{2}{3}}x_7^{0.25}}{x_2x_4^{0.5}} \\ge 0 }$$\n",
"\n",
"$$ { 1-\\frac{1.3x_2x_6}{x_1^{0.5}x_3x_5}-0.8\\frac{x_3x_6^2}{x_4x_5}\n",
" -3.1\\frac{x_2^{0.5}x_6^{\\frac{1}{3}}}{x_1x_4^2x_5} \\ge 0 }$$\n",
"\n",
"$$ { 1-2\\frac{x_1x_5x_7^{\\frac{1}{3}}}{x_3^1.5x_6}-0.1\\frac{x_2x_5} \n",
" {x_3^{0.5}x_6x_7^{0.5}} - \\frac{x_2x_3^{0.5}x_5}{x_1} - \n",
" 0.65\\frac{x_3x_5x_7}{x_2^2x_6} \\ge 0 }$$\n",
"\n",
"$$ { 1 - 0.2\\frac{x_2x_5^{0.5}x_7^{\\frac{1}{3}}}{x_1^2x_4} - 0.3\\frac{x_1^{0.5}x_2^2 \n",
" x_3x_4^{\\frac{1}{3}}x_7^{0.25}}{x_5^{\\frac{2}{3}}} - 0.4\\frac{x_3x_5x_7^{0.75}}{ \n",
" x_1^3x_2^2} - 0.5\\frac{x_4x_7^{0.5}}{x_3^2} \\ge 0 }$$\n",
"\n",
"$$ { 10\\frac{x_1x_4^2x_7^{0.5}}{x_2x_6^3} + 15\\frac{x_3 \n",
" x_4}{x_1x_2^2x_5x_7^{0.5}} + 20\\frac{x_2x_6}{x_1^2x_4x_5^2} \n",
" + 25\\frac{x_1^2x_2^2x_5^{0.5}x_7}{x_3x_6^2} \\ge 100 }$$\n",
"\n",
"$$ { 10\\frac{x_1x_4^2x_7^{0.5}}{x_2x_6^3} + 15\\frac{x_3 \n",
" x_4}{x_1x_2^2x_5x_7^{0.5}} + 20\\frac{x_2x_6}{x_1^2x_4x_5^2} \n",
" + 25\\frac{x_1^2x_2^2x_5^{0.5}x_7}{x_3x_6^2} \\le 3000 }$$\n",
"\n",
"\n",
" \n",
"$$ { 0.1 \\le x_1,x_2,x_3,x_4,x_5,x_6 \\le 10 } $$\n",
"$$ { 0.01 \\le x_7 \\le 10 } $$\n",
"\n",
"The best-known value is $543.667958$ and the suggested starting point is $x_1,x_2,x_3,x_4,x_5,x_6,x_7=6 $\n",
"\n",
"This problem, while still smallish, is more difficult to solve than the previous one and requires more iterations within the solvers. Not all of them can find the optimum for this."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# symbolic math with simpy, just like before\n",
"from sympy import symbols, diff\n",
"from scipy.optimize import minimize\n",
"from cvxopt import solvers, matrix\n",
"import numpy as np, nlopt, pyipopt\n",
"\n",
"# variables\n",
"x=symbols('x[0] x[1] x[2] x[3] x[4] x[5] x[6]')\n",
"z=symbols('z[0] z[1] z[2] z[3] z[4] z[5] z[6]')\n",
"l=symbols('l[0] l[1] l[2] l[3] l[4] l[5]')\n",
"d=symbols('d')\n",
"\n",
"# function to minimize\n",
"# BE AWARE: in Python2, x^(1/3) is not the same as x^(1.0/3.0)\n",
"f0 = 10.0*x[0]*x[3]**2*x[6]**.5/(x[1]*x[5]**3)+15.0*x[2] \\\n",
" *x[3]/(x[0]*x[1]**2*x[4]*x[6]**0.5)+20.0*x[1]*x[5]/(x[0]**2*x[3]*x[4]**2) \\\n",
" +25.0*x[0]**2*x[1]**2*x[4]**0.5*x[6]/(x[2]*x[5]**2)\n",
" \n",
"# non-linear constraints\n",
"f1 = 1.0-.5*x[0]**0.5*x[6]/(x[2]*x[5]**2)-.7*x[0]**3*x[1]*x[5]*x[6]**.5/x[2]**2 \\\n",
" -.2*x[2]*x[5]**(2.0/3.0)*x[6]**.25/(x[1]*x[3]**.5)#>=0\n",
"\n",
"f2 = 1.0-1.3*x[1]*x[5]/(x[0]**.5*x[2]*x[4])-.8*x[2]*x[5]**2/(x[3]*x[4]) \\\n",
" -3.1*x[1]**.5*x[5]**(1.0/3.0)/(x[0]*x[3]**2*x[4])#>=0\n",
"\n",
"f3 = 1.0-2.0*x[0]*x[4]*x[6]**(1.0/3.0)/(x[2]**1.5*x[5])-.1*x[1]*x[4]/ \\\n",
" (x[2]**.5*x[5]*x[6]**.5)-1.0*x[1]*x[2]**.5*x[4]/x[0]- \\\n",
" .65*x[2]*x[4]*x[6]/(x[1]**2*x[5])#>=0\n",
"\n",
"f4 = 1.0-.2*x[1]*x[4]**.5*x[6]**(1.0/3.0)/(x[0]**2*x[3])-.3*x[0]**.5*x[1]**2 \\\n",
" *x[2]*x[3]**(1.0/3.0)*x[6]**.25/x[4]**(2.0/3.0)-.4*x[2]*x[4]*x[6]**.75/ \\\n",
" (x[0]**3*x[1]**2)-.5*x[3]*x[6]**.5/x[2]**2#>=0\n",
"\n",
"f5 = 10.0*x[0]*x[3]**2*x[6]**.5/(x[1]*x[5]**3)+15.0*x[2] \\\n",
" *x[3]/(x[0]*x[1]**2*x[4]*x[6]**0.5)+20.0*x[1]*x[5]/(x[0]**2*x[3]*x[4]**2) \\\n",
" +25.0*x[0]**2*x[1]**2*x[4]**0.5*x[6]/(x[2]*x[5]**2)#>=100\n",
"\n",
"f6 = 10.0*x[0]*x[3]**2*x[6]**.5/(x[1]*x[5]**3)+15.0*x[2] \\\n",
" *x[3]/(x[0]*x[1]**2*x[4]*x[6]**0.5)+20.0*x[1]*x[5]/(x[0]**2*x[3]*x[4]**2) \\\n",
" +25.0*x[0]**2*x[1]**2*x[4]**0.5*x[6]/(x[2]*x[5]**2)#<=3000\n",
" \n",
"g_ge=[f0,f1,f2,f3,f4,f5,-f6]\n",
"g_le=[f0,-f1,-f2,-f3,-f4,-f5,f6]\n",
"g_ipopt=[f0,f1,f2,f3,f4,f5,f6]\n",
"m_cvxopt=z[0]*g_le[0]+z[1]*g_le[1]+z[2]*g_le[2]+z[3]*g_le[3]+z[4]*g_le[4]+z[5]*g_le[5]+z[6]*g_le[6]\n",
"m_ipopt=d*f0+l[0]*g_ipopt[1]+l[1]*g_ipopt[2]+l[2]*g_ipopt[3]+l[3]*g_ipopt[4]+l[4]*g_ipopt[5]+l[5]*g_ipopt[6]\n",
"\n",
"jac_symbols=[[diff(g_ge[con],x[var]) for var in range(len(x))] for con in range(len(g_ge))]\n",
"jac_nlopt=[[diff(g_le[conn],x[varr]) for varr in range(len(x))] for conn in range(len(g_le))]\n",
"jac_ipopt=[[diff(g_ipopt[conn],x[vaarr]) for vaarr in range(len(x))] for coonn in range(len(g_ipopt))]\n",
"hess_cvxopt=[[diff(diff(m_cvxopt,x[varr1]),x[varr2]) for varr1 in range(len(x))] for varr2 in range(len(x))]\n",
"hess_ipopt=[[diff(diff(m_ipopt,x[varrr1]),x[varrr2]) for varrr1 in range(len(x))] for varrr2 in range(len(x))]\n",
"\n",
"# upper and lower bounds for the variables\n",
"lb=np.array([0.1]*(len(x)-1)+[0.01],dtype='float64')\n",
"ub=np.array([10.0]*len(x),dtype='float64')\n",
"\n",
"# suggested starting point\n",
"x0=np.array([6.0]*len(x),dtype='float64')\n",
"\n",
"# upper and lower bounds for the constaints\n",
"c_ge=np.array([np.nan,0,0,0,0,100,-3000],dtype='float64')\n",
"g_lb=np.array([0,0,0,0,100,-np.inf],dtype='float64')\n",
"g_ub=np.array([np.inf]*5+[3000],dtype='float64')\n",
"\n",
"cons_scipy=[\n",
" {'type':'ineq','fun':lambda x: eval(str(g_ge[1]))-c_ge[1],'jac':lambda x: np.array(eval(str(jac_symbols[1])))},\n",
" {'type':'ineq','fun':lambda x: eval(str(g_ge[2]))-c_ge[2],'jac':lambda x: np.array(eval(str(jac_symbols[2])))},\n",
" {'type':'ineq','fun':lambda x: eval(str(g_ge[3]))-c_ge[3],'jac':lambda x: np.array(eval(str(jac_symbols[3])))},\n",
" {'type':'ineq','fun':lambda x: eval(str(g_ge[4]))-c_ge[4],'jac':lambda x: np.array(eval(str(jac_symbols[4])))},\n",
" {'type':'ineq','fun':lambda x: eval(str(g_ge[5]))-c_ge[5],'jac':lambda x: np.array(eval(str(jac_symbols[5])))},\n",
" {'type':'ineq','fun':lambda x: eval(str(g_ge[6]))-c_ge[6],'jac':lambda x: np.array(eval(str(jac_symbols[6])))}\n",
"]\n",
" \n",
"cons_ipopt=lambda x: np.array([eval(str(g_ipopt[1])),eval(str(g_ipopt[2])),eval(str(g_ipopt[3]))\n",
" ,eval(str(g_ipopt[4])),eval(str(g_ipopt[5])),eval(str(g_ipopt[6]))],dtype='float64')\n",
"\n",
"# constraints for nlopt with in-place gradient modification\n",
"def f0_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[0])), dtype='float64')\n",
" return eval(str(g_le[0]))\n",
"def f1_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[1])), dtype='float64')\n",
" return eval(str(g_le[1]))+c_ge[1]\n",
"def f2_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[2])), dtype='float64')\n",
" return eval(str(g_le[2]))+c_ge[2]\n",
"def f3_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[3])), dtype='float64')\n",
" return eval(str(g_le[3]))+c_ge[3]\n",
"def f4_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[4])), dtype='float64')\n",
" return eval(str(g_le[4]))+c_ge[4]\n",
"def f5_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[5])), dtype='float64')\n",
" return eval(str(g_le[5]))+c_ge[5]\n",
"def f6_nlopt(x,grad):\n",
" if grad.size>0:\n",
" grad[:]=np.array(eval(str(jac_nlopt[6])), dtype='float64')\n",
" return eval(str(g_le[6]))+c_ge[6]\n",
"\n",
"# # formatting other necessary info\n",
"eval_f=lambda x: eval(str(f0))\n",
"eval_grad_f=lambda x: np.array(eval(str(jac_symbols[0])))\n",
"\n",
"nvar=int(len(x))\n",
"ncon=int(len(g_ge)-1)\n",
"\n",
"str_jac_ipopt=list()\n",
"for v1 in range(ncon):\n",
" for v2 in range(nvar):\n",
" str_jac_ipopt.append((v1,v2))\n",
"str_jac_ipopt=(np.array(str_jac_ipopt,dtype='int32').T[0], np.array(str_jac_ipopt,dtype='int32').T[1])\n",
"\n",
"str_hess_ipopt=list()\n",
"for v1 in range(nvar):\n",
" for v2 in range(v1+1):\n",
" str_hess_ipopt.append((v1,v2))\n",
"str_hess_ipopt=(np.array(str_hess_ipopt,dtype='int32').T[0], np.array(str_hess_ipopt,dtype='int32').T[1])\n",
"\n",
"nnzj=int((len(g_ge))*len(x))\n",
"nnzh=int(sum([i+1 for i in range(len(x))]))\n",
"\n",
"def eval_jac_ipopt(x,flag):\n",
" if flag==True:\n",
" return str_jac_ipopt\n",
" else:\n",
" return np.array(eval(str(jac_ipopt)), dtype='float64').reshape((nnzj,1))\n",
"def eval_h_ipopt(x,z,d,flag):\n",
" if flag==True:\n",
" return str_hess_ipopt\n",
" else:\n",
" out=list()\n",
" for v1 in range(len(x)):\n",
" for v2 in range(v1+1):\n",
" out.append(eval(str(hess_ipopt[v1][v2])))\n",
" return np.array(out, dtype='float64')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 4.39410417 0.85446875 2.84323017 3.39997898 0.72292614 0.87040643\n",
" 0.02463883]\n",
"Objective value: 543.667958466\n",
"Result code: 4\n"
]
}
],
"source": [
"# COBYLA under NLopt\n",
"opt = nlopt.opt(nlopt.LN_COBYLA,7)\n",
"\n",
"# setting the function to minimize\n",
"opt.set_min_objective(f0_nlopt)\n",
"\n",
"# adding the constraints\n",
"opt.add_inequality_constraint(f1_nlopt)\n",
"opt.add_inequality_constraint(f2_nlopt)\n",
"opt.add_inequality_constraint(f3_nlopt)\n",
"opt.add_inequality_constraint(f4_nlopt)\n",
"opt.add_inequality_constraint(f5_nlopt)\n",
"opt.add_inequality_constraint(f6_nlopt)\n",
"\n",
"# setting the bounds for the variables\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"# setting the numerical tolerance for the objective\n",
"opt.set_xtol_rel(1e-10)\n",
"\n",
"# running the solver\n",
"x = opt.optimize(x0)\n",
"print \"Optimal x:\", x\n",
"print \"Objective value:\", opt.last_optimum_value()\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 loop, best of 3: 35.3 s per loop\n"
]
}
],
"source": [
"%%timeit\n",
"opt.optimize(x0)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 4.39410332 0.85446887 2.84323014 3.39997727 0.72292638 0.87040641\n",
" 0.02463884]\n",
"Objective value: 543.667363053\n",
"Result code: 4\n"
]
}
],
"source": [
"# CCSAQ under NLopt\n",
"opt = nlopt.opt(nlopt.LD_CCSAQ,7)\n",
"\n",
"# setting the function to minimize\n",
"opt.set_min_objective(f0_nlopt)\n",
"\n",
"# adding the constraints\n",
"opt.add_inequality_constraint(f1_nlopt)\n",
"opt.add_inequality_constraint(f2_nlopt)\n",
"opt.add_inequality_constraint(f3_nlopt)\n",
"opt.add_inequality_constraint(f4_nlopt)\n",
"opt.add_inequality_constraint(f5_nlopt)\n",
"opt.add_inequality_constraint(f6_nlopt)\n",
"\n",
"# setting the bounds for the variables\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"# setting the numerical tolerance for the objective\n",
"opt.set_xtol_rel(1e-10)\n",
"\n",
"# running the solver\n",
"x = opt.optimize(x0)\n",
"print \"Optimal x:\", x\n",
"print \"Objective value:\", opt.last_optimum_value()\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 loop, best of 3: 41.5 s per loop\n"
]
}
],
"source": [
"%%timeit\n",
"opt.optimize(x0)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 4.39055143 0.85458152 2.84235139 3.4044958 0.72337832 0.87131107\n",
" 0.02467146]\n",
"Objective value: 543.668347173\n",
"Result code: 3\n"
]
}
],
"source": [
"# MMA under NLopt\n",
"opt = nlopt.opt(nlopt.LD_MMA,7)\n",
"opt.set_min_objective(f0_nlopt)\n",
"opt.add_inequality_constraint(f1_nlopt)\n",
"opt.add_inequality_constraint(f2_nlopt)\n",
"opt.add_inequality_constraint(f3_nlopt)\n",
"opt.add_inequality_constraint(f4_nlopt)\n",
"opt.add_inequality_constraint(f5_nlopt)\n",
"opt.add_inequality_constraint(f6_nlopt)\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"opt.set_xtol_rel(1e-8)\n",
"opt.set_ftol_rel(1e-8)\n",
"\n",
"# running the solver\n",
"x_opt = opt.optimize(x0)\n",
"print \"Optimal x:\", x_opt\n",
"print \"Objective value:\", opt.last_optimum_value()\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 loop, best of 3: 24.6 s per loop\n"
]
}
],
"source": [
"%%timeit\n",
"opt.optimize(x0)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
" fun: 543.660277163643\n",
" jac: array([ -57.07384276, -1089.33427157, 138.87592324, 200.25885006,\n",
" -549.76482059, -495.11913883, -5071.04873243])\n",
" message: 'Positive directional derivative for linesearch'\n",
" nfev: 59\n",
" nit: 39\n",
" njev: 35\n",
" status: 8\n",
" success: False\n",
" x: array([ 4.39408088, 0.854464 , 2.84325329, 3.39998429, 0.72292814,\n",
" 0.87041065, 0.02464184])"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# SLSQP under SciPy\n",
"minimize(fun=eval_f, jac=eval_grad_f, constraints=cons_scipy, bounds=[(lb[bn],ub[bn]) for bn in range(len(x))],\n",
" x0=x0, method='SLSQP')"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 loop, best of 3: 5.4 s per loop\n"
]
}
],
"source": [
"%%timeit\n",
"minimize(fun=eval_f, jac=eval_grad_f, constraints=cons_scipy, bounds=[(lb[bn],ub[bn]) for bn in range(len(x))],\n",
" x0=x0, method='SLSQP')"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" pcost dcost gap pres dres\n",
"Couldn't solve problem\n"
]
}
],
"source": [
"# CP+LDL solver under CVXOPT\n",
"def F(x=None,z=None):\n",
" if x is None: return 6,matrix(x0)\n",
" if max(x) >= 10.0: return None,None\n",
" if min(x[0:6]) <= 0.1: return None,None\n",
" if x[6] <= 0.01: return None,None\n",
" \n",
" f=matrix(np.array([\n",
" eval(str(g_le[0]))+c_le[0],\n",
" eval(str(g_le[1]))+c_le[1],\n",
" eval(str(g_le[2]))+c_le[2],\n",
" eval(str(g_le[3]))+c_le[3],\n",
" eval(str(g_le[4]))+c_le[4],\n",
" eval(str(g_le[5]))+c_le[5],\n",
" eval(str(g_le[6]))+c_le[6]\n",
" ],dtype='float64'))\n",
" \n",
" Df=matrix(np.array(eval(str(jac_nlopt)),dtype='float64'))\n",
" if z is None: return f,Df\n",
" \n",
" # if z is provided, it also needs to output the hessian avaluated at x,z\n",
" H=matrix(np.array(eval(str(hess_cvxopt)),dtype='float64'))\n",
" return f,Df,H\n",
"\n",
"\n",
"options={'show_progress':True,'kktreg':1e-4,'reltol':1e-10} #setting refinement>=2 makes it loop forever\n",
"try:\n",
" xsol=solvers.cp(F=F, kktsolver='ldl', options=options)\n",
"except:\n",
" print \"Couldn't solve problem\""
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# pyipopt crashes when trying to solve this problem\n",
"\n",
"# solver doesn't support keywork argument, thus arguments need to be passed in this exact order\n",
"opt = pyipopt.create(nvar, lb, ub, ncon, g_lb, g_ub, nnzj, nnzh, eval_f, eval_grad_f, cons_ipopt, eval_jac_ipopt,eval_h_ipopt, lambda x: True)\n",
"opt.num_option('tol', 1e-10)\n",
"# opt.solve(x0) # --> this would invoke the solver, but it crashes"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [4.3941, 0.854469, 2.84323, 3.39998, 0.722926, 0.870406, 0.0246388]\n",
"Objective value: 543.668\n",
"Constraint multipliers: [-38.705, -834.115, -1287.32, -86.7667, -4.15984e-11, 7.49225e-12]\n"
]
}
],
"source": [
"from casadi import SX, vertcat, nlpsol\n",
"\n",
"x_casadi=SX.sym(\"x_casadi\",7)\n",
"\n",
"f0_casadi = 10.0*x_casadi[0]*x_casadi[3]**2*x_casadi[6]**.5/(x_casadi[1]*x_casadi[5]**3)+15.0*x_casadi[2] \\\n",
" *x_casadi[3]/(x_casadi[0]*x_casadi[1]**2*x_casadi[4]*x_casadi[6]**0.5)+20.0*x_casadi[1]*x_casadi[5]/(x_casadi[0]**2*x_casadi[3]*x_casadi[4]**2) \\\n",
" +25.0*x_casadi[0]**2*x_casadi[1]**2*x_casadi[4]**0.5*x_casadi[6]/(x_casadi[2]*x_casadi[5]**2)\n",
" \n",
"# non-linear constraints\n",
"f1_casadi = 1.0-.5*x_casadi[0]**0.5*x_casadi[6]/(x_casadi[2]*x_casadi[5]**2)-.7*x_casadi[0]**3*x_casadi[1]*x_casadi[5]*x_casadi[6]**.5/x_casadi[2]**2 \\\n",
" -.2*x_casadi[2]*x_casadi[5]**(2.0/3.0)*x_casadi[6]**.25/(x_casadi[1]*x_casadi[3]**.5)#>=0\n",
"\n",
"f2_casadi = 1.0-1.3*x_casadi[1]*x_casadi[5]/(x_casadi[0]**.5*x_casadi[2]*x_casadi[4])-.8*x_casadi[2]*x_casadi[5]**2/(x_casadi[3]*x_casadi[4]) \\\n",
" -3.1*x_casadi[1]**.5*x_casadi[5]**(1.0/3.0)/(x_casadi[0]*x_casadi[3]**2*x_casadi[4])#>=0\n",
"\n",
"f3_casadi = 1.0-2.0*x_casadi[0]*x_casadi[4]*x_casadi[6]**(1.0/3.0)/(x_casadi[2]**1.5*x_casadi[5])-.1*x_casadi[1]*x_casadi[4]/ \\\n",
" (x_casadi[2]**.5*x_casadi[5]*x_casadi[6]**.5)-1.0*x_casadi[1]*x_casadi[2]**.5*x_casadi[4]/x_casadi[0]- \\\n",
" .65*x_casadi[2]*x_casadi[4]*x_casadi[6]/(x_casadi[1]**2*x_casadi[5])#>=0\n",
"\n",
"f4_casadi = 1.0-.2*x_casadi[1]*x_casadi[4]**.5*x_casadi[6]**(1.0/3.0)/(x_casadi[0]**2*x_casadi[3])-.3*x_casadi[0]**.5*x_casadi[1]**2 \\\n",
" *x_casadi[2]*x_casadi[3]**(1.0/3.0)*x_casadi[6]**.25/x_casadi[4]**(2.0/3.0)-.4*x_casadi[2]*x_casadi[4]*x_casadi[6]**.75/ \\\n",
" (x_casadi[0]**3*x_casadi[1]**2)-.5*x_casadi[3]*x_casadi[6]**.5/x_casadi[2]**2#>=0\n",
"\n",
"f5_casadi = 10.0*x_casadi[0]*x_casadi[3]**2*x_casadi[6]**.5/(x_casadi[1]*x_casadi[5]**3)+15.0*x_casadi[2] \\\n",
" *x_casadi[3]/(x_casadi[0]*x_casadi[1]**2*x_casadi[4]*x_casadi[6]**0.5)+20.0*x_casadi[1]*x_casadi[5]/(x_casadi[0]**2*x_casadi[3]*x_casadi[4]**2) \\\n",
" +25.0*x_casadi[0]**2*x_casadi[1]**2*x_casadi[4]**0.5*x_casadi[6]/(x_casadi[2]*x_casadi[5]**2)#>=100\n",
"\n",
"f6_casadi = 10.0*x_casadi[0]*x_casadi[3]**2*x_casadi[6]**.5/(x_casadi[1]*x_casadi[5]**3)+15.0*x_casadi[2] \\\n",
" *x_casadi[3]/(x_casadi[0]*x_casadi[1]**2*x_casadi[4]*x_casadi[6]**0.5)+20.0*x_casadi[1]*x_casadi[5]/(x_casadi[0]**2*x_casadi[3]*x_casadi[4]**2) \\\n",
" +25.0*x_casadi[0]**2*x_casadi[1]**2*x_casadi[4]**0.5*x_casadi[6]/(x_casadi[2]*x_casadi[5]**2)#<=3000\n",
"\n",
"g_casadi=vertcat(f1_casadi,f2_casadi,f3_casadi,f4_casadi,f5_casadi,f6_casadi)\n",
"solver = nlpsol(\"solver\", \"ipopt\", {'x':x_casadi,'f':f0_casadi,'g':g_casadi}, {\"print_time\":False,\"ipopt\":{\"print_level\":0}})\n",
"res=solver(x0=x0,lbx=lb,ubx=ub,lbg=g_lb,ubg=g_ub)\n",
"\n",
"print \"Optimal x:\" , res[\"x\"]\n",
"print \"Objective value:\" , res[\"f\"]\n",
"print \"Constraint multipliers:\",res['lam_g']"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"10 loops, best of 3: 62 ms per loop\n"
]
}
],
"source": [
"%%timeit\n",
"solver(x0=x0,lbx=lb,ubx=ub,lbg=g_lb,ubg=g_ub)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter objective inf_pr inf_du ||d|| lg(rg) ls \n",
" 0 2.208886e+03 3.70e+02 7.37e+02 0.00e+00 - 0 \n",
" 1 1.718131e+02 3.81e+01 7.26e+02 5.99e+00 - 2 \n",
" 2 1.777915e+02 1.80e+00 1.84e+03 5.16e+00 - 1 \n",
" 3 7.070532e+02 2.30e-01 3.95e+03 1.93e+00 - 2 \n",
" 4 4.724446e+02 1.44e+00 6.00e+02 5.01e+00 - 2 \n",
" 5 3.154566e+02 9.39e-01 2.33e+02 7.25e-01 - 1 \n",
" 6 5.295164e+02 3.74e-01 3.30e+04 1.13e+00 - 1 \n",
" 7 4.593004e+02 1.05e-01 1.20e+04 1.14e+00 - 1 \n",
" 8 5.041100e+02 2.48e-02 2.43e+03 5.46e-01 - 1 \n",
" 9 5.377245e+02 5.83e-03 4.72e+02 2.88e-01 - 1 \n",
"iter objective inf_pr inf_du ||d|| lg(rg) ls \n",
" 10 5.435357e+02 9.07e-05 8.75e+00 1.95e-02 - 1 \n",
" 11 5.436678e+02 1.33e-07 9.37e-03 1.47e-03 - 1 \n",
" 12 5.436680e+02 6.46e-14 1.57e-09 1.31e-06 - 1 \n",
"\n",
"casadi::SQPMethod: Convergence achieved after 12 iterations.\n",
"CPU times: user 20 ms, sys: 4 ms, total: 24 ms\n",
"Wall time: 25.4 ms\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"CasADi warning: \"Indefinite Hessian detected...\" issued on line 474 of file \"/feedstock_root/build_artefacts/casadi_1493990261301/work/casadi-3.1.1/casadi/solvers/nlpsol/sqpmethod.cpp\". \n",
"CasADi warning: \"Indefinite Hessian detected...\" issued on line 474 of file \"/feedstock_root/build_artefacts/casadi_1493990261301/work/casadi-3.1.1/casadi/solvers/nlpsol/sqpmethod.cpp\". \n"
]
}
],
"source": [
"%%time\n",
"solver2 = nlpsol(\"solver\", \"sqpmethod\", {'x':x_casadi,'f':f0_casadi,'g':g_casadi},\n",
" {\"print_time\":False,'print_header':False,\"qpsol\":\"qpoases\",\"qpsol_options\":{\"printLevel\":None}})\n",
"solver2(x0=x0,lbx=lb,ubx=ub,lbg=g_lb,ubg=g_ub)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 2.43530278 0.55516093 1.45423705 1.51679216 4.30055254 2.87853302\n",
" 0.02576996]\n",
"Objective value: 64.9269246535\n",
"Result code: 6\n"
]
}
],
"source": [
"# # NLopt again, using the barrier method with a subsolver for the steps\n",
"opt = nlopt.opt(nlopt.LD_AUGLAG,7)\n",
"opt.set_min_objective(f0_nlopt)\n",
"opt.add_inequality_constraint(f1_nlopt)\n",
"opt.add_inequality_constraint(f2_nlopt)\n",
"opt.add_inequality_constraint(f3_nlopt)\n",
"opt.add_inequality_constraint(f4_nlopt)\n",
"opt.add_inequality_constraint(f5_nlopt)\n",
"opt.add_inequality_constraint(f6_nlopt)\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"opt.set_xtol_rel(1e-10)\n",
"\n",
"local_opt=nlopt.opt(nlopt.LN_PRAXIS,7)\n",
"opt.set_local_optimizer(local_opt)\n",
"\n",
"# without timeouts, the problem keeps running forever, even though it finds a very close solution\n",
"opt.set_maxtime(100)\n",
"\n",
"# running the solver\n",
"x_opt = opt.optimize(x0)\n",
"print \"Optimal x:\", x_opt\n",
"print \"Objective value:\", opt.last_optimum_value()\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal x: [ 2.95109739 1.38182442 4.86691815 2.84396433 2.00098308 1.75\n",
" 0.17469822]\n",
"Objective value: 64.9064586174\n",
"Result code: 6\n"
]
}
],
"source": [
"# again AUGLAG, this time with the same COBYLA used at the beginning\n",
"opt = nlopt.opt(nlopt.LD_AUGLAG,7)\n",
"opt.set_min_objective(f0_nlopt)\n",
"opt.add_inequality_constraint(f1_nlopt)\n",
"opt.add_inequality_constraint(f2_nlopt)\n",
"opt.add_inequality_constraint(f3_nlopt)\n",
"opt.add_inequality_constraint(f4_nlopt)\n",
"opt.add_inequality_constraint(f5_nlopt)\n",
"opt.add_inequality_constraint(f6_nlopt)\n",
"opt.set_lower_bounds(lb)\n",
"opt.set_upper_bounds(ub)\n",
"\n",
"opt.set_xtol_rel(1e-10)\n",
"\n",
"local_opt=nlopt.opt(nlopt.GN_DIRECT_L,7)\n",
"opt.set_local_optimizer(local_opt)\n",
"\n",
"# without timeouts, the problem keeps running forever, even though it finds a very close solution\n",
"opt.set_maxtime(100)\n",
"# local_opt.set_maxtime(15)\n",
"\n",
"# running the solver\n",
"x_opt = opt.optimize(x0)\n",
"print \"Optimal x:\", x_opt\n",
"print \"Objective value:\", opt.last_optimum_value()\n",
"print \"Result code:\", opt.last_optimize_result()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Graphical comparison of solvability and time it took"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAGMCAYAAAChyhdyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xe8HFX9//HXmwBJIBgpCgktlIAQqgQrJSA2QJooTSEK\nYtcvfKWjRr8WEBX0y1f5UUPvvYigEJogJBBIAoQaSmgJkEAgFJPP749zFiaT3Xv33tx79ybzfj4e\n+7g755w558zc2f3sOTOzq4jAzMxsUbdYqztgZmbWExzwzMysEhzwzMysEhzwzMysEhzwzMysEhzw\nzMysEhzwrFMkDZEUkka1ui+LAkmj8v4c0uq+9LS83aObLDsylx/Rvb3qnI5si/U8BzwD3nuhNvsY\n0ur+tie/Mf5Xq/thtqAkfVvSuZIeljRHUps3T0saLOksSdMkzZY0VtJX6pT7sKQzJD0g6RVJb0l6\nTNJpktbuvi1qncVb3QHrNb5eWt4SOBA4GbitlDcNeBPoD/yn+7vWKSOBIcAJre2G2QI7AlgeuA9Y\nGlilUUFJywG3Ax8G/gg8C+wNXCTpmxFxRqH4ssA6wA3AU8BsYCjwTeArkj4REQ92/ea0jgOeARAR\n5xSXJS1OCnh3lvMK3ur2jtlCQ9ISQJ+I8HHRtUYAT0fEXEnX0EbAAw4H1gB2ioirASSdBtwJ/F7S\nxRExCyAiJgOfLlcg6RLgbuAHwPe6ckNazVOa1in1zuEV0yR9VdL4PKXymKRv5DKrSbokT6G8Lukc\nScvUqX+QpL9KelrSO5Kek3SypA830bcpwNbA6qWp2BGFMltJulHSzNzHeyXt34Ht30HSLZKm5/Wf\nlnSZpHVK5TaSdLmkl/OU0YOSDpXUp536v5v7vFOdvMUkPStpfCl9eG5ruqS3JU2WdFT+8NLMNoWk\n0ZK2k3SXpDclvSDpT5IGlMrWzjkOk/RHSc+SPgB9olDmgLxfZ+f9fIOkLdpov91221i3r6QjJU3K\n+3mGpKslbVoqNyL3e6Sk7+V99JakCZJ2zGU2lHS9pNfy/+3POZg3TdIn8/HxRq7j1Dr7cFVJp0t6\nKv+/XpL0L0n7FctFxJSImNtk03sDj9eCXV5/DvC/wHLA9k3U8VT+u2yTbS40PMKz7rAj8B3gL8Ar\nwP7A6ZLeAX4D3AQcCWxOmj55CzigtrKk1UifSJcETgMeB9YGvgtsI2l4RMxso/3/An4LrAAcVEh/\nKNf/JeBy4AXgD8DrwJ7AqZLWjIij2to4SVsDVwETczszgMHAdrmfj+Ryw4FbgHeB/8vtfQk4FtgY\n2KeNZi4Ajgf2zW0VfQZYOfe91qcdgMuAx3L6K8AngV8CmwDzncNp4KPA7sApwFnANsCPgA0kfbbO\nG++5pKmwPwABPJ/7cyxwKGmkcCSwDGnG4GZJO0fEdQvY7ntyMLoe+BRwNnAiMBD4FnCHpK0iYmxp\nte+T3tBPJR1/PwIuVzrXdQpwPnAF8Dngh8BLwK8a77Z5bAJcA5wBnEcaoe0PzM37oDaDciPp//gX\n0jEzENiIdDrhzCbbeo+kQbm+c+tk35X/bg5cVFpvidz2EqTjd1TOKv+PFn4R4Ycf8z1I58ACGNkg\nf0jOH1Un7Q1g9UL6h0hvKnOBg0v1XAa8AwwopF1JeoNZpVR2OOmc4agm+j8GmFInvQ/pE+wMYHAh\nfUngDmAOMLSduv+Yt/PD7ZS7I/d3o0KaSG84AXymkD4qpw0ppF2c99uypXrPJgXRD+flfqRgeiuw\neKnsQbneEU3ss8iPXUrpf8rpe9bp75g6ba6b/9e3A0sW0gfn/T6FNPXZmXZrx+WIOtv4+dL6HwCe\nBsYU0kbkslOBgYX0jXL6XGC3Uj3jgOebfN3U6vh4Kf3a/D8bUGrv0A6+Lq8BokHeZrnOY+vkLZXz\nzquTt2PhfxD5WDq4I/1aWB6e0rTucEVE1KZFiIhpwGTSG8H/lcreRvpkOQRA0kDSC/Aq4C1JK9Qe\npDfKx0ifujtrM2A14PSIeK7Qx3eA35Gm+Xdup47a6PLLjaYL89Trp4CrIuKBQjsB/Dov7tpOO2cC\nfYE9CvUOyOtdHxEv5eTPAiuSRhQfLO2z2qf0ZvfZ5Ii4opR2TBv9PSEiyhcu7UwK7L/L+xWAvL/P\nAFYHNi2t09F2i74GPAyMK237kqRR1BaS+pfWGR2FWYL8P3oNeC4iLiuVvR1YqdnpVdJ573+X0m4i\nzagNycu1trdRE9P0TVoq/327Tt5bpTJFd5GOoZ1I5wCfB5Ztdip8YbLIbZD1Ck/USXuV9Cm5/GJ8\nNf9dPv9dlxR09s+PZutv1hr576Q6ebW0Ndup40TSm/pfgGMl3U6aUjs/B/f22nmIFPzba+d60kh3\nX+CknPZl0pV6ZxXKrZf/nt5GXSu201axb/OIiOclzaB+fx+pk9bsPi5OM3a03aL1SFcMT2ujzArA\nM4XlRsfoMw3SIR2js9rpS6O6Xy7UQUQ8JenXpCswn8/nY/8JXBwR9zTRRj1v5r996+T1K5V5T0RM\nB/6RF6+WdDbwAOlKz293si+9kgOedYc5HUyHNCIo/j2HxucxZnemU10lIl6WtDnpXMtnga1I59t+\nIWn7iLizi9r5j6TzgP+StHZEPEYKfq8y73m92j47BBhPfc81SF9Q872BtoCACcDBbZQpB8MFOUbb\n01QdEXG0pNOBHUjH0gHAIZJ+FxGHNdlWUe1/vHKdvFra1PYqiYjnJP0D2F/Sj+p8SF1oOeBZb/MY\n6TzCkhHxj/YKt6HRzbm1T9/D6uStXyrTuPJ05duY/EDSRqRzPUeT3sCebKOdj5BGsc2MVM8kXYSz\nr6RTSOegTi69CT2a/76xgPsM3h8tvidfDPFBmh9ZF/fx46W8Rvt4Qdp9lHSe+KZo/mrGXiEiniBd\nQfm/kvoBfwcOlfSHwpR1s3U9L2kqhStlC2pp5Yt3GulPOt/9AdoeOS9UfA7PepWIeJl03mk3SfO9\ncJV8qImqZpHOQ5Q/ld9LupDhG5JWKtS7BGmEFKSLZhrK54fKHiaNPJfL2/ES8C/gS5I2KPafNI0F\n6UrRNkXEeNL00tdIXw6wGPOPfP9Omvo8XOnG43J/+6vOrR8NrCtpl1JabbRRPsfWyFWk/XhI8XL+\nHMC+Qbpo6L4ubPcsYCUajPAkNTud22MkDSzf6hDp/sXa1G5nbwk4H1grX4lca6sP6UrTGRSuvGy0\nXyStT7oS+PHCFP0iwSM8642+S7pQ4FZJZ5HeHBcjncvZmfQGN6qdOu4iXfxyoqR/kaaZboqIlyT9\ngBRs7pF0Mum2hD1In4J/ExGPNqo0O0XSKrz/DRX98/rLMO+5tR+Tbku4TVLttoQdgc+Trpb7Z3s7\nIjuTdNn/YcAjEXFXMTMi3pC0LykwTM7TZI+RRkcfAXYjXfgxpom2JgDn5NHko6TbA3bP23FhM52N\niMmSjiPdlnCrpAt5/7aEAcA+eYTcVe3+iTS1fJykbUkXiLxGujjpM6QLNrZppu89aBvgZEmXki7o\nmkW6oOoA4N+RbgoH3ruNZuO8uHZOOzovz4iIEwv1HkO6BeU8SX8kTWHuRbod4YCIeL1Q9ghJnyVd\nQTqFNN26AemD1RKkWzcWLa2+TNSP3vlgwW5LGFWn/Bjq3yZQa2dEKX0F4DjSRRFvkT6dTiC9ua3f\nRP+XIt3D9yIp2JUvZd+adAXfa7n++4D9m9w3u5FGMc+SroibRnpj/nKdshuTAtEruexDpEDQp1Ru\nFKXbEgp5K5IuaQ/gqDb6tQHp3OdU0q0eL5JGmT8FlmtiuwIYTbqf8N+kEeuLpCm3ZZrtb6HMt/J+\nfSvv5xuBLRew3UbHy+Kke+nuId0W8wYpcJ4LfK5QbkSj45r0pj+mTnq721relvaOc9KFPSfl4+G1\n3N+HSPdNDiytO5p5bxsoPuq9plYm3boyPe/7e4E96pTbDrgkb/eb+fh8gnQl7bDOvG/09ofyhptZ\nxSl9KfGZETGy1X0x6w4+h2dmZpXggGdmZpXggGdmZpXgc3hmZlYJHuGZmVklOOCZmVklOOCZdQFJ\nG0j6T76Rt5b23o+NtrBrCz1JUySN6cR6x0t6pPyNJlZdDnhmXeOPwB0RcWOrO2LvORZYhfTNPWb+\najGzBSXpk6Svtip/F+StpK8de7fHO2VExAuSLiB9x+hfYv7f7bOK8QjPbMF9j/Q1TtcVEyNibkS8\nFfN/b6T1nLOBQbT/o75WAQ54Zgsg/yr0LsA/IuLdUt585/CKaZK+IWmSpLclPSXp0A60u6+kuyXN\nkPSGpCcknVv+JQlJQyWdLel5Se/k82HHSVq6Tp0rSfpzruttSS9JurF4XjKX2yqnz5Q0W9K9kub7\nsV5JY3J7gyWdL+lVSW9K+rukdeqUX1XSRbne1yRdLWmtBtu/g6RbJE3PfXha0mV16r2V9D2VX2li\nt9oizlOaZgtmM9IvANzdwfW+Q/pS6NNIX4z9NdKvpz8bEee1taKkr5N+QeE24GekL1teFdie9CvV\n03K5zUi/HDAD+H+kL5XemPQly5+WtHUtSEsaAtyR+3QW6XfTlib9gsR2pC9+rn1z/+WkX374A+mX\nJvYETpW0ZkQcVeru0qSgcxdwJOlLk38MXClpg9roV9IHc7lVSV+q/CDpC75vJk0LF7d/a9KXd08E\nfpu3b3Du59oUfoU9IuZIuifXZVXX6m+v9sOPhflB+n23AHaqkzeC0jfzF9Keo/Ct+KRfd5gG3NlE\nm5eRvmF/8XbK3U/6nb7yLw7sWqdf1+W0z9epZ7H8tw/p55BmAIML+UuSguUcYGghfUyu89BSfYeU\n2wJ+k9O+USp7Qk4fU0j7Y077cJP/o1Nz+eVbfbz40dqHpzTNFkxtCvGVDq53RkTMrC1ExJukUdDQ\nJtadSQqQO+QflJ2PpA2BjYDzgL6SVqg9SL81+AbwuVx2OeALwPUR8fdyXfH+r4hvRvqNudMj4rlC\n/jvA70inSMrnyuYCfy6l3ZT/Frd1F9JPAp1VKntsnc2r7bcv5ynl9ryc/364ibK2CHPAM1swte/m\nqxt42vBEnbSXgeWbWPc3pJHWFcA0SZdKOkDz/qr5evnvL0gjx+LjJdJUY+0Xr9fO/S//CnnZGvnv\npDp5tbQ1S+nPRfol76JaACpu65rAo1G6wCcinieNKItOzH39C/CKpOsk/ah8/rKg9r/x9yhWnM/h\nmS2Yafnvch1cr9NXbkbEo5LWJ/2a92dI56dOAX4haauIeJz33+T/AFzfoKpXO9uHDmhrOzv6IQGA\niHhZ0ubAlqTbQbYCjidt//YRcWdpldr/ZhpWaQ54ZgtmYv7bzFRkl4mIt0nn3a4DkLQ9cC1wMPB9\n0q99A8yJiH+0U91jpNHPJu2Uq41Kh9XJW79UpqOeAIZK6lMc5UkaBHywXDiXGZMfSNoIGAccDexQ\nKr428EJEvIxVmqc0zRbMfaQLSD7RUw3m83Bl9+a/tdHMfaRg/B1J5WlGJC2ez90REa8AfwO+KGm7\nOmVrI7F7gaeBb0haqZC/BO9fiHJlpzYqrbcisG8p/bA6/am3/Q+TrlZdrlS2DzAcuKWT/bJFiEd4\nZgsg0mXvlwG7SOqbR17d7QZJM0i3JTxDGgGNJAWcs3O/It++cBPwgKTTSefZliKNeHYDjgBG5zp/\nAPwL+JukM0mjpf7Ax4EpwGF5W39Aui3hHkknk25L2IMU8H8TEbWRZUf9DtgbOCXfTjGJdEXrJ0k3\n9RedImkV4AbSucz+uQ/LMP9FL1uTzlde3Ml+2SLEAc9swf2VFHB2BC7tofa+CnybNKJ5mTSi+2FE\n3FwrFBHjJW1KCmw7ke79e50UwEYD/yyUfVLScOCnpPv59iWd47sfOLlQ7mpJnyFNHR5CuiXhIeCA\niDitsxsUEa9K2pJ0y0FtlHcLsE2xn9nZpP29H+kq2ddI9+3tHhHl/f910j2DnR152iLEPwBr1gUk\nXQ8sHRFbtrovluRp1yeAwyOifGuEVZADnlkXkDSMNBraPiJuaHV/DCSdQBqtDovS175ZNTngmZlZ\nJfgqTTMzqwQHPDMzqwRfpdkiK6ywQgwZMqTV3TAzW6iMGzduekQ0+hq5NjngtciQIUMYO3Zsq7th\nZrZQkfRUZ9f1lKaZmVWCA56ZmVWCA56ZmVWCA56ZmVWCA56ZmVWCA56ZmVWCA56ZmVWCA16LTJg6\nkyGHX9vqbpiZVYYDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVULTAU/SUZImSXpA\n0nhJH8/pYyQNL5VdStK5kiZImijpdkkDct4qkq6U9KikJySdKKlvzhshaWau/yFJP2/Ql0GSriks\nf0zSrZImS7pP0qmSlurMDpG0i6SQ9JFC2hBJs3PdD0m6W9LIQv5Okg7Pz38g6ZudadvMzLpPUwFP\n0ieBHYGPRsRGwHbAM22s8mPgxYjYMCI2APYH3pUk4DLgiogYCgwF+gO/K6x7W0RsAgwHvibpo3Xq\nPxg4JfdtReBi4LCIWDciNgWuB5ZpZtvq2Au4Pf8tejwiNo2I9YA9gf+S9A2AiLgqIo7J5U4HftjJ\nts3MrJs0O8IbBEyPiLcBImJ6RDzXTvmptYWImJzX3RZ4KyLOyOlzgIOAfWsjwMI6bwDjgLXr1P9l\nUlAD+D5wZkTcWVj3koh4MY/87swjs39JWhdA0rA8ShufR6xDc/oAYAtSgN6z0cZFxBOkoPujvN5I\nSSfmvDeBKZI+Vl5P0oGSxkoaO+fNmQ13npmZdb1mA94NwKqSHpH0F0lbt1P+dOCwHGx+VQsowDBS\nEHtPRLwGTKEU2CQtD3wCmFRKXwN4tRZ8gQ3KdRY8DGyZR30/A36T078D/Kkwknw2p+8MXB8RjwAv\nS9qsjW28F/hIg7yxwJblxIg4OSKGR8TwPksNbKNqMzPrak0FvIiYBWwGHAhMAy4snsOqU348sCZw\nHLAccI+k9Zrs05aS7iMF2WMiYlIpf1DuQzMGAhdLmggcTwq4AHcCR0o6DFg9Imbn9L2AC/LzC5h/\nWrNIbeS9BAxuso9mZtYDFm+2YJ5+HAOMkTQB2A8Y3Ub5WaTzdZdJmgtsD9wP7F4sJ+kDwErAZODj\npHN4O7bRldlAv8LyJFIwvrJO2f8Bbo6IXSUNyf0nIs6T9G9gB+A6Sd8GxpOmXDeUFEAfICQd0qAf\nmwIPNcjrl/tpZma9RLMXraxbmJYE2AR4qo3yn5a0bH6+JLB+Lv9PYClJ++a8PsAfgBMLo6z2PAIM\nKSyfCOxXu2o017tbvphlIO+fSxxZyF8TeCIi/kwKlBuRAvHZEbF6RAyJiFWBJ6kzNZmD5++B/23Q\nx3WAiU1uj5mZ9YBmz+ENAM6U9KCkB0gBbFQh/1pJz+bHxcBawC15JHgf6ZzWpRERwK7A7pIeBV4G\n5kbEr5vtcL6Y5XFJa+flF0kXmPw+35bwEPB54HXS1Z+/zVOkxdHsV4GJksaTzgGeRZq+vLzU3KW8\nP625Vu22BOAi4M+1i2/q+DRwY7PbZGZm3U8pBrWocelTwPnArhFxbwfW2xXYLCKO7rbOdZKkTYGD\nI+LrbZXrO2hoDNrvBKYcs0MP9czMbOEnaVxEDG+/5PyaPofXHSLiX8DqnVjv8nwVZ2+0AvDTVnfC\nzMzm1dKAtyAi4tRW96GeiPBUpplZL+Tv0jQzs0pwwGuRDVce6PN3ZmY9yAHPzMwqwQHPzMwqwQHP\nzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwq\nwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqwQHPzMwqYfFWd6CqJkydyZDDr50nzb+AbmbWfTzCMzOz\nSnDAMzOzSnDAMzOzSnDAMzOzSljggCdpVld0pFTnSEmD28g/QdJW+fkSko6R9KikeyXdKemLC9D2\neEkXlNJGS3pS0v2SHpF0lqRVCvnXSfqgpCUl3SrJFwOZmfUyvXWENxKoG/AkLQ98IiJuzUn/AwwC\nNoiIjwK7AMt0plFJ6wF9gC0lLV3KPiQiNgbWBe4DbpK0JEBEbB8RMyLiHeCfwB6dad/MzLpPlwU8\nSSPy6OZaSZMlnSRpsZy3l6QJkiZKOrawzixJx0uaJOmfkj4kaXdgOHBuHm31LzX1ZeD6vP5SwLeA\nH0bE2wAR8WJEXJTz/yppbK7/F4V2j5H0oKQHJP2+UPdewNnADcDO9bYzkuOBF4Av5vqmSFohF7kC\n2Kcz+9DMzLpPV4/wPgb8EFgfWAvYLU9NHgtsC2wCbC5pl1x+aWBsRAwDbgF+HhGXAGOBfSJik4iY\nXWrj08C4/Hxt4OmIeK1Bf46KiOHARsDWkjbKI8RdgWERsRHwq0L5PYALgPNJwa8t9wIfqZM+Edi8\nnXXNzKyHdXXAuzsinoiIOaSgsQXpzX9MREyLiP8A5wJb5fJzgQvz83Ny+fYMAqY12Z+vSrqXNAU5\njBSIZwJvAadJ2g14E0DScGB6RDxNmpbcVNJybdSteol529+RNN+0qqQD84hz7Jw3Zza5CWZm1hW6\nOuBFO8sdXb+e2UC//PwxYDVJHygXkrQG8BPgM3kkdy3QLwfdjwGXADuSp0dJI7qPSJoCPA58gDR9\n2simwEMN8vqSguo8IuLkiBgeEcP7LDWwzY00M7Ou1eVTmpLWyOfu9gBuB+4mTSeuIKkPKbDcUmh/\n9/x871we4HUaX3jyEGkqk4h4EzgN+FPtApJ8HvArpID1BjBT0oq8f75tADAwIq4DDgI2zv39KrBh\nRAyJiCGkc3jzTWsq+RFppHl9nfzlSSPFd9vdW2Zm1mO6OuDdA5xICkpPApdHxPPA4cDNwP3AuIi4\nMpd/gxQkJ5LO8f0yp48GTmpw0cq1wIjC8tGkKc4Hcz3XAK9FxP2kqcyHgfOAO3L5ZYBrJD1ACrAH\nA1sCUyPiuUK9twLrSxqUl4+TdD/wCGmadpt8VWbZNrmPZmbWiyiio7OODSqSRgA/iYgdO7DOrIgY\n0Im2bgd2jIgZHV23u0m6DDg8Ih5pq1zfQUNj0H4nzJPmL482M2ubpHH5YsQO66334bXnv4HVWt2J\nsjytekV7wc7MzHpel30jSESMAcZ0cJ0Oj+7yev/uzHrdLU9xntXqfpiZ2fwW1hGemZlZhzjgmZlZ\nJfhLjltkw5UHMtYXqZiZ9RiP8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzM\nrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc8MzMrBIc\n8MzMrBL8i+ctMmHqTIYcfm2ru9GuKf5VdjNbRHiEZ2ZmleCAZ2ZmleCAZ2ZmleCAZ2ZmleCAZ2Zm\nldClAU/SUZImSXpA0nhJH8/pYyQNL5VdStK5kiZImijpdkkDct4qkq6U9KikJySdKKlvzhshaWau\n/yFJP2/Ql0GSrmmvrQbrfkfSvvn5aElP5vbulfTJdvbBBZKGdmS/mZlZ9+uy2xJyINgR+GhEvC1p\nBWDJNlb5MfBiRGyY118XeFeSgMuAv0bEzpL6ACcDv8vrANwWETtKWhoYL+nqiLi3VP/BwClttdWo\nYxFxUinpkIi4RNLngP8HbNTGdv0VOBT4VhtlzMysh3XlCG8QMD0i3gaIiOkR8Vw75afWFiJicl53\nW+CtiDgjp88BDgL2LY/KIuINYBywdp36vwxc305bSNo3j0jvl3R2Thsl6Sd16rwVWFvSWpLeC7CS\nhhaWbwO2k+R7HM3MepGuDHg3AKtKekTSXyRt3U7504HDJN0p6VeFacBhpCD2noh4DZhCKbBJWh74\nBDCplL4G8GotqDVqS9Iw4Ghg24jYmPdHkI18CZgQEY8DMyVtktO/AdQC9FzgMWDj8sqSDpQ0VtLY\nOW/ObKcpMzPrSl0W8CJiFrAZcCAwDbhQ0sg2yo8H1gSOA5YD7pG0XpPNbSnpPlKQPSYiJpXyB+U+\ntNfWtsDFETE9l3ulQXvHSRqft23/nHYq8I085boHcF6h/EvA4DrbfHJEDI+I4X2WGtjkppqZWVfo\n0mm3PP04BhgjaQKwHzC6jfKzSOfrLpM0F9geuB/YvVhO0geAlYDJwMfJ5/Da6MpsoF8Tbb3T5KYd\nEhGXlNIuBX4O3ASMi4iXC3n9ch/MzKyX6LIRnqR1S1cnbgI81Ub5T0taNj9fElg/l/8nsFThKsk+\nwB+AEyOi2SDyCDCkibZuAr6Sp0aRtFyT9RMRbwF/J12kckYpex1gYrN1mZlZ9+vKc3gDgDMlPSjp\nAVJQGVXIv1bSs/lxMbAWcEseCd4HjAUujYgAdgV2l/Qo8DIwNyJ+3WxH8sUsj0uqnfNr1NYk4Nc5\n737gjx3c5nOBuaSpVQAkrQjMjogXOliXmZl1I6X40ntJ+hRwPrBrnVsP2lpvV2CziDi6G/v2E2Bg\nRPy0kHYQ8FpEnNbWun0HDY1B+53QXV3rMv61BDPrTSSNi4jh7ZecX6+/dD4i/gWs3on1Lq9NVXYH\nSZeTRo7blrJmAGd3V7tmZtY5vT7gLYiIOLUb6961QXr5fJ6ZmfUC/i5NMzOrhEV6hNebbbjyQMb6\n/JiZWY/xCM/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/M\nzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBAc/MzCrBv3jeIhOm\nzmTI4de2uhtmPWbKMTu0ugtWcR7hmZlZJTjgmZlZJTjgmZlZJTjgmZlZJTjgmZlZJXR7wJMUks4p\nLC8uaZqka/LyyFxmu0KZXXLa7oW0FSS9K+k77bR3iaQ18/Mpki4t5O0uaXSh3RO7YPtGShpcWL5A\n0tAFrdfMzLpWT4zw3gA2kNQ/L38WmFoqMwHYs7C8F3B/qcxXgLtyXl2ShgF9IuKJQvJmktbvTMeb\nNBIYXFj+K3BoN7ZnZmad0FNTmtcBtZtw9gLOL+XfBnxM0hKSBgBrA+NLZfYC/htYWdIqDdrZB7iy\nlPYH4KiDE4RyAAAgAElEQVRmOyppL0kTJE2UdGwhfZak4yVNkvRPSR/KI9DhwLmSxuegfhuwnSTf\n42hm1ov0VMC7ANhTUj9gI+DfpfwA/gF8HtgZuKqYKWlVYFBE3A1cBOzRoJ1PA+NKaRcBH5W0dnud\nzFOTxwLbApsAm0vaJWcvDYyNiGHALcDPI+ISYCywT0RsEhGzI2Iu8BiwcZ36D5Q0VtLYOW/ObK87\nZmbWhXok4EXEA8AQ0ijtugbFLiBNa+7J/CPAPUiBq1au0bTmIGBaKW0OcBxwRBNd3RwYExHTIuI/\nwLnAVjlvLnBhfn4OsEUb9bzEvNOcAETEyRExPCKG91lqYBPdMTOzrtKT025XAb8HRgDLlzMj4m5J\nGwJvRsQjkorZewErSdonLw+WNDQiHi1VMxvoV6fts0kBb+KCbcK8XW4jr1/ui5mZ9RI9eVvC6cAv\nImJCG2UOB44sJkhaBxgQEStHxJCIGAL8lvqjvIdI5//mERHvAscDB7XTx7uBrfMVoX1yG7fkvMWA\n2lWjewO35+evA8uU6lmHrg2uZma2gHos4EXEsxHx53bK/C0ibi4l7wVcXkq7lPoB71rSCLKe05h/\nRDtS0rO1B9CHFHRvJl0lOi4iahfBvEG6sGYi6RzfL3P6aOCk2kUrklYEZkfEC21sqpmZ9TBFtDUz\nt3DJV0neDHw6IuZ0cd2zImJAE+UOAl6LiNPaKtd30NAYtN8JXdY/s97Ov5ZgXUHSuIgY3pl1F6lv\nWomI2cDPgZVb2I0ZwJktbN/MzOpY5O4Vi4i/d1O97Y7ucrkzuqN9MzNbMIvUCM/MzKyRRW6Et7DY\ncOWBjPU5DTOzHuMRnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZ\nVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYIDnpmZVYJ/\n8bxFJkydyZDDr211N8xsITPlmB1a3YWFlkd4ZmZWCQ54ZmZWCQ54ZmZWCQ54ZmZWCQ54ZmZWCe0G\nPEkrSbpA0uOSxkm6TtI6OW+YpJskTZb0qKSfSlLOGylpmqTxkiZJukTSUpI+K+nOQrk+ku6T9ClJ\noyT9pEE/dpEUkj7SRl/7S7ol1zkkl/9hIf9ESSPz89GSdu/Q3qrf5pGF50tKulWSr341M+tl2gx4\nOShdDoyJiLUiYjPgCGBFSf2Bq4BjImJdYGPgU8D3ClVcGBGbRMQw4B1gj4i4EXgK2D+X+SEwNiL+\n1U5f9wJuz38b+SZwWUTMycsvAT+WtGQ7dS+I9wJeRLwD/BPYoxvbMzOzTmhvhLcN8G5EnFRLiIj7\nI+I2YG/gjoi4Iae/CfwAOLxcSR7xLA28mpMOAo6QNCyvc1hbnZA0ANiCFCT3bKPoPsCVheVppAC0\nX1v1F9qRpOMkTZQ0QdIeOX1EHrldm0ezJ0laTNIxQP88ij03V3NF7oeZmfUi7QW8DYBxDfKGlfMi\n4nFggKQP5KQ9JI0HpgLLAVfncs8DJwB3Ar+KiFfa6cfOwPUR8QjwsqTNygXyKG7NiJhSyjoW+Imk\nPu20AbAbsAlptLodcJykQTnvY6TR6PrAWsBuEXE4MDuPYmtBbiKweb3KJR0oaayksXPenNlEd8zM\nrKt090UrF0bEJsBKwATgkELe/wF9ImJ0E/XsBVyQn19A/WnNFYAZ5cSIeAL4N2lE2p4tgPMjYk5E\nvAjcwvvB6+6IeCJPl56fy84n578jaZk6eSdHxPCIGN5nqYFNdMfMzLpKewFvEjDfaCp7sJwnaU1g\nVkS8VkyPiCCN7rYqpM0For0OSloO2BY4VdIUUtD8au2il4LZQL8G1fyGNG1aXqcjyn1tq+99gbcW\noC0zM+ti7QW8m4C+kg6sJUjaSNKWwLnAFpK2y+n9gT8Dv2tQ1xbA453o4+7A2RGxekQMiYhVgSeB\nLYuFIuJVoI+k+YJeRDxMCtBfaqet20jTsH0kfYgUoO/OeR+TtIakxUgXpdye09+VtEStAknLA9Mj\n4t0Ob6mZmXWbNgNeHpntCmyXb0uYBPwWeCEiZpPOrR0taTJpyvIe4MRCFXvkCzoeADYF/qeJPh0t\n6dnagzR9eXmpzKXUn9a8gQZTjcCvgVVKaf+v0NaduZ0HgPtJwf7QiHghl61t20OkgFvr08nAA4WL\nVrYB/K3QZma9jFJMWzRI+ihwUER8vYvrHQH8JCJ2bKLsZcDh+QKbhvoOGhqD9juhi3poZlVR9V9L\nkDQuIoZ3Zt1F6ptWIuJe4OYmr8jscvlK0SvaC3ZmZtbzFrlvBImI07uhzjHAmCbKvQOc1dXtm5nZ\nglukRnhmZmaNOOCZmVklLHJTmguLDVceyNiKn3w2M+tJHuGZmVklOOCZmVklOOCZmVklOOCZmVkl\nOOCZmVklOOCZmVklOOCZmVklOOCZmVklOOCZmVklOOCZmVklOOCZmVklOOCZmVklOOCZmVklOOCZ\nmVklOOCZmVklOOCZmVklOOCZmVkl+BfPW2TC1JkMOfzaVnfDzKxHTTlmh5a17RGemZlVggOemZlV\nggOemZlVggOemZlVQqcCnqSVJF0g6XFJ4yRdJ2mdnLdOXn5U0r2SLpK0oqSlJJ0raYKkiZJulzSg\nUOcukkLSR0ptDZN0k6TJub1fSKrbb0mbSjotPx8paa6kjQr5EyUNyc+nSFqhM9tfqG+IpL0LyxtK\nGr0gdZqZWffocMCTJOByYExErBURmwFHACtK6gdcC/w1IoZGxEeBvwAfAn4MvBgRG0bEBsD+wLuF\nqvcCbs9/a231B64CjomIdYENgY/luuo5EvhzYflZ4KiObmMHDAHeC3gRMQFYRdJq3dimmZl1QmdG\neNsA70bESbWEiLg/Im4jvfnfGRFXF/LGRMREYBAwtZA+OSLeBsgjvS1IQXDPQlt7A3dExA15nTeB\nHwCHlDslaRlgo4i4v5B8DTBM0rrNbJik5SRdIekBSXfVRoeSRkk6W9KdeeT6rbzKMcCWksZLOiin\nXV3aBjMz6wU6E/A2AMZ1Iu904LAcNH4laWghb2fg+oh4BHhZ0mY5fVi5voh4HOgv6YOl+ocDE0tp\nc4HfkUZ+zfgFcF9EbJTXOauQtxGwLfBJ4GeSBgOHA7dFxCYRcXwuNxbYssn2zMysh/TYRSsRMR5Y\nEzgOWA64R9J6OXsv4IL8/AIK05odMAiYVif9POATktZooo4tgLNzf28Clpf0gZx3ZUTMjojpwM2k\nqdV6XgIG18uQdKCksZLGznlzZhPdMTOzrtKZb1qZBOzeRt7WjVaMiFnAZcBlkuYC20t6kTRy2lBS\nAH2AkHQI8CCwVbEOSWsCL0fEjFL1s4F+ddr8j6Q/AIc1s3FtiHaWa/rlvsxfQcTJwMkAfQcNbbS+\nmZl1g86M8G4C+ko6sJYgaSNJW5JGU5+StEMhbytJG0j6tKRlc9qSwPrAU6TgeXZErB4RQyJiVeBJ\n0rTgucAWkrbL6/UnXZTy8zr9eghYu0GfRwPbkS6eacttwD65rRHA9Ih4LeftLKmfpOWBEcA9wOvA\nMqU61mH+qVUzM2uxDge8iAhgV2C7fJvAJOC3wAsRMRvYEfhhvrjjQeB7pKnGtYBbJE0A7iOd67qU\nNH15eamZS4G9cn07AUdJegSYTrqI5dw6/XoYGJgvXinnvUMKlB8uZT0g6dn8+CMwCthM0gOkC1L2\nK5YlTWXeBfxPRDyX0+ZIur9w0co2pCtVzcysF1GKXwsHSbsAfwS2iYin6uQfBLweEad2cbujgFkR\n8ft2yvUFbgG2iIj/tFW276ChMWi/E7quk2ZmC4EF/fJoSeMiYnhn1l2ovmklIq6IiDXrBbvsr8Db\nPdmnktWAw9sLdmZm1vMWqZ8Hioi3yFdZdnG9o5os9yjwaFe3b2ZmC26hGuGZmZl1lgOemZlVwiI1\npbkw2XDlgYxt4S//mplVjUd4ZmZWCQ54ZmZWCQ54ZmZWCQ54ZmZWCQ54ZmZWCQ54ZmZWCQ54ZmZW\nCQ54ZmZWCQvVryUsSiS9DkxudT/qWIH0M0y9UW/tW2/tF/TevvXWfkHv7Zv7laweEe39tmld/qaV\n1pnc2Z+46E6SxvbGfkHv7Vtv7Rf03r711n5B7+2b+7XgPKVpZmaV4IBnZmaV4IDXOie3ugMN9NZ+\nQe/tW2/tF/TevvXWfkHv7Zv7tYB80YqZmVWCR3hmZlYJDngtIOkLkiZLekzS4a3uT42kKZImSBov\naWyL+3K6pJckTSykLSfpRkmP5r/L9pJ+jZI0Ne+38ZK2b0G/VpV0s6QHJU2S9OOc3hv2WaO+tXS/\nSeon6W5J9+d+/SKnt3SftdGvlh9nuR99JN0n6Zq83PJjrFme0uxhkvoAjwCfBZ4F7gH2iogHW9ox\nUsADhkdEy+/1kbQVMAs4KyI2yGm/A16JiGPyB4VlI+KwXtCvUcCsiPh9T/al1K9BwKCIuFfSMsA4\nYBdgJK3fZ4369lVauN8kCVg6ImZJWgK4HfgxsBst3Gdt9OsLtPg4y/07GBgOfCAiduwNr8tmeYTX\n8z4GPBYRT0TEO8AFwM4t7lOvExG3Aq+UkncGzszPzyS9afaoBv1quYh4PiLuzc9fBx4CVqZ37LNG\nfWupSGblxSXyI2jxPmujXy0naRVgB+DUQnLLj7FmOeD1vJWBZwrLz9ILXvxZAP+QNE7Sga3uTB0r\nRsTz+fkLwIqt7EzJDyU9kKc8WzqlI2kIsCnwb3rZPiv1DVq83/L03HjgJeDGiOgV+6xBv6D1x9kJ\nwKHA3EJay/dXsxzwrGiLiNgE+CLw/Tx91ytFmovvFZ96gb8CawKbAM8Df2hVRyQNAC4F/isiXivm\ntXqf1elby/dbRMzJx/wqwMckbVDKb8k+a9Cvlu4vSTsCL0XEuEZlWn2MtccBr+dNBVYtLK+S01ou\nIqbmvy8Bl5OmX3uTF/P5oNp5oZda3B8AIuLF/AY1FziFFu23fL7nUuDciLgsJ/eKfVavb71lv+W+\nzABuJp0n6xX7rNyvXrC/Pg3slM/1XwBsK+kcetH+ao8DXs+7BxgqaQ1JSwJ7Ale1uE9IWjpfUICk\npYHPARPbXqvHXQXsl5/vB1zZwr68p/Ziz3alBfstX+hwGvBQRPyxkNXyfdaob63eb5I+JOmD+Xl/\n0oVkD9PifdaoX63eXxFxRESsEhFDSO9bN0XE1+gFx1iz/OXRPSwi/iPpB8DfgT7A6RExqcXdgjTv\nfnl6b2Jx4LyIuL5VnZF0PjACWEHSs8DPgWOAiyTtDzxFusqvN/RrhKRNSFM5U4Bv93S/SJ++vw5M\nyOd+AI6kF+yzNvq2V4v32yDgzHzl9GLARRFxjaQ7ae0+a9Svs3vBcVZPbzjGmuLbEszMrBI8pWlm\nZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXg\ngGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdmZpXggGdm\nZpXggNcikkZJOqeL6uov6WpJMyVd3E7ZIZJC0uJ5eYykA7qiH7m+LSVN7qr6FlR5e3uL4n6XtI+k\nG5opa71TTx5nkkZL+lV3t9MZkqZI2q5B3ghJz/Z0n4oc8LqJpFmFx1xJswvL+3Rxc7sDKwLLR8RX\nurjuDomI2yJi3Vb2YWETEedGxOe6oq4cHEfk5132oWph0psDQkdJGinp9lb3Y1HhgNdNImJA7QE8\nDXypkHZuFze3OvBIRPyni+u1BnrbiNGsuy0Kx7wDXmstKeksSa9LmiRpeC1D0mBJl0qaJulJST+q\nV4GkXwA/A/bIo8f9JS0m6WhJT0l6KbcxsL3OtLWepDMl/Xd+vnKevvl+Xl5L0it5/XmmLfIUx08k\nPZCnXC+U1K+Qf6ik5yU9J+mAXO/aDfo3UtITeX89WRspN7u9kvaQNLaUdpCkq/LzvpJ+L+lpSS9K\nOklS/5w3QtKzkg6T9AJwhqQVJF0jaUbe/tsk1X1NSfqspIfzPjgRUGm7bu9I2dzPV/N++GKd9r4A\nHMn7x8X9Dfo1RdIh+f/zhqTTJK0o6W95P/9D0rKF8hdLeiH37VZJwwp5oyX9n6Rr87r/lrRWIf9P\nkp6R9JqkcZK2LOT1z8fYq5IeysdF8Thq9vVwILAPcGje7qtz+npKo98ZSq+1neqtX9jHLTnOSuus\nB5wEfDJvy4xC9rJt7OePSLoxH5OTJX21jW0dI+m3ku7O/5crJS2X82rTtPtLehq4KafvlPfhjLz+\neqVqN5f0YP5fnqHC673UdsP/qdLsxMWSzsnbOEHSOpKOyPv+GUkdnxWJCD+6+QFMAbYrpY0C3gK2\nB/oAvwXuynmLAeNIgWxJYE3gCeDzDeofBZxTWP4m8FhebwBwGXB2zhsCBLB4Xh4DHNDEet8Ers7P\n9wYeBy4s5F2Zn48Ani1t+93AYGA54CHgOznvC8ALwDBgKeCc3Le162zj0sBrwLp5eRAwrCPbm9t4\nHRhaqPceYM/8/HjgqtzPZYCrgd8Wtus/wLFAX6B//p+dBCyRH1sCqtP3FXK7u+dyB+W6avt9JHB7\nB8q+C3yLdNx8F3iuQbujKBwXbRybd5GmxFcGXgLuBTYF+pHe5H5eOraWyfvgBGB8IW808DLwsby/\nzwUuKOR/DVg+5/13/t/3y3nHALcAywKrAA+QjyM6/noYDfyqsLxEPj6OzOtvm/fxur3tOKvTn/eO\njWb2c+7/M8A3ct6mwHRg/Qb1jwGmAhvkdS+tHTOFbTor5/UH1gHeAD6b9+uheZ8sWTieJgKr5u27\no/a/oPDe0N7/lPffHz+ft+Ms4EngqNzut4AnC9txOHBNu+/F3fUm78d8byr1At4/CsvrA7Pz848D\nT5fKHwGc0aD+Ucwb8P4JfK+wvC7pTXJx2g54ba23FvBqPlBPAr5dOHjPBA4uH9SFbf9aYfl3wEn5\n+ekUXujA2rQd8GYAXwb6l/I6sr3nAD/Lz4eS3piWIo2i3gDWKtTzydqLKm/XO+Q36Jz2S+DKev0t\n9W9f8oeZvCzgWeoHvGbKPlbIXypv30rtHRdtHJv7FJYvBf5aWP4hcEWDdT+Y2x6Yl0cDpxbytwce\nbqPtV4GN8/N5AhhwQOH46ujrYTTzBrwtScF1sULa+cCo3nac1enPe8dGafvq7mdgD+C2Uvn/R+FD\nSylvDHBMYXl90nHep7BNaxbyfwpcVFhejBQwRxSOp++U+vZ44TXU1P+UdOzeWMj7EjAL6JOXl8l9\n+2Bbx3f54SnN1nqh8PxNoJ/SPPnqwOA8ZTAjT2UcSfoU3ozBwFOF5adIL8r21m+4XkQ8TnqhbkJ6\nA7kGeE7SusDWpE/njZS3c0ChvWcKecXn84iIN0gv5u8Az+fpnI+01+86VZ0H7JWf7016M38T+BDp\nDWlcYZ9fn9NrpkXEW4Xl40ifbm/IU2CHN+j+PNsZ6RXbaFubKftCIf/N/HQAnfdi4fnsOssDACT1\nkXSMpMclvUZ6c4M0Kp2vb8z7v0ZpavuhPB06AxhYWLetY6ErXg/PRMTcQtpTpBHtPHrJcdaMRvt5\ndeDjpX21D7BSG3UV9/VTpBHUCg3y59kHeZ8+w7z7slzf4DptNvM/LR+H0yNiTmEZOnjcO+D1Ts+Q\nPvF9sPBYJiK2b3L950gHVM1qpGmxF+sXb3q9W0hTbUtGxNS8vB9pGmp8k30rep40fVWzaluFI+Lv\nEfFZ0jTTw8ApTfa76EbgQ5I2Ib0hnZfTp5NeRMMK+3xgpIuO3utCqT+vR8R/R8SawE7AwZI+02A7\n39s2SWpjWztStj3RfpEO2RvYGdiOFKyG5HQ1WqEmn687FPgqsGxEfBCYWVi3rWOho6+H8nY/B6yq\nec+vrkYamcy/cuuPs7a2pT3PALeU9tWAiPhuG+sU9/VqpFHr9AZ9mGcfFI7P4r4s1/dcg34uyHtc\npzjg9U53A68rXSDRP3+y3kDS5k2ufz5wkKQ1JA0AfkM639beVZztrXcL8APg1rw8Ji/fXvjk1REX\nAd/IFxQsRZouqUvpQoqdJS0NvE2a3qh9Ym96eyPiXeBi0uhsOdIbU+2T6inA8ZI+nNtcWdLn2+jT\njpLWzi/6mcCcQp+KrgWGSdotj+B/RONP3B0p254XgSFqcCFNJyxD2vcvk0Ypv+nguv8BpgGLS/oZ\n8IFC/kXAEZKWlbQy6biq6ejr4UXSOaGaf5NGQYdKWkLpto0vAReUV+yFx9mLwCqSlmyQX3YNsI6k\nr+dtXULS5nUuLCn6mqT182vwl8AlbbyeLwJ2kPQZSUuQzsW+DfyrUOb7klZRuvjlKODCOvUs6Htc\npzjg9UL5YNuRNH34JOnT1qmkT9XNOB04mxSYniSd/P1hF6x3C+mNqxbwbie98d1KJ0TE34A/AzeT\npgbvyllv1ym+GHAw6dPiK6Rp1Nqn1o5u73mkUcrFpTerw2r9yFN2/yCdp2lkaC4zC7gT+EtE3Fxn\nO6cDXyFdmPFyXu+OehV2pGwTal9C8LKkeztZR9FZpCmqqcCDvP//asbfSVN3j+Q63mLeqa9fks5V\nPknap5eQj4NOvB5OA9bPU2VXRMQ7pAD3xbzuX4B9I+LhOuv2tuPsJmAS8IKk6Q3KvCciXgc+B+yZ\nt+EF3r/QqpGzSecFXyBdqFT3Cthc/2TSxUf/S9qXXyLdcvVOodh5wA2k87KPA/PdE9kF73HzkHSk\npL+1Wy6fADRrufwpdCLQt4nRqC3CJH2XdFXj1q3uy6JM0hjShU2ntrovPcEjPGspSbsq3Ze0LOmT\n6NUOdtUjaZCkTyvd67Yuaars8lb3yxYtDnjWat8m3fv1OOkcWFsn123RtSTp8vnXSdN4V5KmHs26\njKc0zcysEjzCMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDA\nMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAMzOzSnDAsx4laZSkc1rd\nDwBJIyXd3gPtrCZplqQ+3d1WMyR9XtIVre5HT5L0XUkv5v/D8gtY1xRJ23Wg/CRJIxakzYWBpI0k\n/avV/WiLA94iqKMvyO5qR9IISc92dz96u4h4OiIGRMSc7mynA/v718Ax3dmXntTEcbgE8Efgc/n/\n8HIXtt3uB7iIGBYRYzpZf0h6SdLihbQlcloU0ubbB8UPdJKG5LruK5VZQdI7kqaU6pqdPxy8KGm0\npAE5b5ikGyS9ImmGpHGSts/b+QAwQ9KXOrOtPcEBz3qt4ovcuoakzYGBEXFXq/vSg1YE+gGTWt2R\nTnoV+GJh+Ys5rTOWkrRBYXlv4Mk65b4UEQOAjwLDgaNz+tXAjcBKwIeBHwGvFdY7F/h2J/vW7Rzw\nFnH5U94dkk6UNFPSw5I+U8gfLOmq/IntMUnfKuSNknSJpAslvS7pXkkb57yzgdWAq/MnwUNL7S4N\n/A0YnPNnSRqcs5eUdFauc5Kk4YX1pkg6TNIDwBuSFpe0nqQx+RPlJEk7FcqPkXRAaXtvLyx/TtLk\nvO1/kXRLsXwu83tJr0p6UlLxjaW8L6dIOkLSg7n8GZL65byJxU+2+VP4dEmbFj5dL17o828l3S3p\nNUlXSlqusO5OeTtn5LLrtdeHdvZ30ReBWwr1SdLxecTwmqQJtTdESX3zvnk6f9I/SVL/wrqHSnpe\n0nOSDsjbuHbOG533999yX+6QtJKkE3K/H5a0aaGuwZIulTQt/x9+VMgbJemiesdME8fhOsDkvDhD\n0k05/U+SnsnbPE7SloV1Rkv6VWG57shZ0heAI4E9ctv319nf84y+2tqWNpwN7FtY3hc4q5112qpr\nv2brioippONqA+n/t3fm4X5N5x7/fA0NqUQQRWIsQqNKH1Vibs2UIho0gtIHNcZwaUvdcIuWqqG4\nxhKJIRIxtampGmOoeWq1tyKcSKJNDI2ISnjvH+/a5+zfzm88OUlOzu/9PM95zm9Pa6299trrXdN+\nv+oNrANca2afpr8nzCw/LTAe2FFSt3amb4ESBq852AJ4A+gN/DcwNlfB3gZMBvoA+wPnSfp27trv\nAqOBFYFbgLskLW1mQ4C3SS1BM7sgH6GZzcIr1ynp+HJmNiUd3jvF2wu4B7i8kN6DgD3TceGtygfw\nFuXxwM2SNqh10+kFHQP8BFgJr/i2KpM3f0t5cwFwvSRVCXYwsCuwLtCPtpbvTcDBufP2AKaa2QuU\n5xDgcGA1YC5wWUpzP+BWYCiwMjAOr8y/UC0NNfI7z8a0GQCAXYDtUjjLA4OAbMjvF2n/psB6QF/g\nrJTO3YCTgZ3SsR3KxDUIz5/ewH+ACcDzaXsMPsyIpCXwZ/xSimNHYKikXXNhlS0zdZTDvwMbpc1e\nZpaV7WfSfWXlenTWeKkXM7sPOA8YleLepM5La5X/IncB20nqJWkFYFvg7kbSmmMkcKCkJSX1B5YD\nnq50sqQ18LL8Al4u/gGMlLSPpFWK5ycDOQeo+X4uCsLgNQf/BC4xszlmNgqv8PZMhXlr4HQz+8TM\nXgSuo7Q1+ZyZjTGzOXgFtQyw5Xym53EzG5fmtEYAxYriMjNrMbPZKa7lgF+kFuXDwO9wo1iLPYDX\nzGysmWVGZVrhnLfM7NqUluG4AZrnRc5xeUrbe/hcWJaOkcAeknqm7SHp3ioxwsxeTYbqZ8Ag+aKW\nA4Dfm9mDKc9/BSxLqaGulIZ66AXMzG3PAXoAGwIys7+a2dRk9I8ETjKz98xsJl65H5iuGwTcYGav\nmdnHwLAycd1pZs+Z2SfAncAnZnZTyutRQNbD2xxY2czOSc94InBtLi6oXWYawsxGmtkMM5trZhcB\n3Vh4lXSj9/IJ3iA4IP3dk/a1h8n4+78T/p5XKqN3SfoAeBwfETjPzAz4FjAJuAiYKulRSesXrp2J\nlx79bpcAABSZSURBVLNOR8yRNAfvpMKa8Rbeo+sDZJVZ/lh+iKUl+2Fmn6ehnXJDZY2QNzofA8tI\nWioZpZI4U1wtZvZ5IY1964inD6XptzJDU9Nyxz9OnbvlqoSZT1uWj5jZFElPAAMl3Yn3tk5sIJyl\n8Z5Pn7SdpelzSS2U3m/ZNNTJ+7iBy8J/WNLlwBXAWpLGAqfiDZvuwHO5Dq+AbKVpH+DZCmnKeDf3\ne3aZ7Syf18KHYj/IHV8SeCy3XavMNISkU4Ej8PswoCee/wuD9tzLTcD5+DM4vczxuXgZyrM03qAp\nF9ZheCNqW7wXX2QfM3uouNPMJgPHQWvv75oU3oDcaT2AD4rXdgaih9cc9C0M060JTEl/K0rqUTj2\nTm57jexHGnpaPV0HXlFUo9bxeq6bAqyR4i6Xxll4xZyxau73VDy9gM9X5bfbyRq531k+ZgzHhzW/\nB0xIwzv1hjMHmJ7CW6uQ5jWo8EwKaagnv1+mUMGZ2WVmthnQPx37r5SW2cBGZtYr/S2fFjJAIW8L\naWqUFuDNXDy9zKyHme1R5/UNlbM0X3ca3ktdwcx6AR/ixgSql6n5ins+eIy20Ydyn9K8Daxd2LcO\nucZTjjvwKYOJZvZ2exNkZi14Q6l1EYykvsAXKB027zSEwWsOvgScIF9I8T3gK8C4VGCfBM6XL3z4\nGt7qzS+z3kzSfvIFF0PxuZhshd+7wJerxPsusJKk5ecj7U/jreDTUvp3APbC50AAXgT2k9Q9LZg4\nInft74GN03zDUsCxVK+86uFYSaunOdAz8KG5jLvwVW0nUntRwcGS+kvqDpwDjElDXLfjw807ypfT\nn4Lnef77pkppqCe/xwHbZxuSNpe0RYprFj5U9nnqUV8LXCzpS+ncvrl5tduBH8gXFHXHh2Xby5+B\nmfLFSsum+aWvyleU1kOtclikB94j+hewlKSz8B5exov48PSKklbFy321uNcuNMg6nDRCsxewd2G0\nJmMUPu+5oZxv4HPEtxVPTMPo3wZ+WDxWDUkrSDpb0nqSlkhz5IfTVh+Al62Hzew/jYS9sAiD1xw8\nDayPt9rPBfbPfYt0EN4ynILPs/x3YSjjbnze4H18Xmq/NLcEPsRypnw14anFSM3sdXwBxsR0TsND\noWb2Kf6i757SfyVwSAob4GLgU7ziGY4vi86unY73ti7AJ9z748Nw8/My3oIvoJmILwRqXc2X5hzv\nwFvWY2uEMwK4ER/eWgZf3o2Z/Q3vJf4Gv9+98AUZn9ZKQz35bWbPAx9K2iLt6okbtvfx3sAM4MJ0\n7HR8kcJTkv4NPESa5zKzP+Bzon/KzknXNJy3ydB/B19E8ma67+vwRTT1ULUcluF+4D7g7/g9f0Lp\nkOwIfAHNJDyfR1GZ0en/DEnP15nedpHmSyt9WnEtcAM+1/ch3uA6Iy2sKRfWs2b2RoNJ+BSvKx7C\nP0V4FX/eh+XOGQxc1WC4Cw2VbywEXQVJhwE/NLNt2nHtMGA9Mzu41rmLA6kVPhkYbGZ/asf1k/C8\nnGduI3fOWUC/ankmaTww0syuWxBpqCOMXYBjzGyf9oZRJsyv4BVgt/bOqwWLN2mE6GozG1Dz5EVE\n9PCCLo3cjVYv+XdBP8XnaRbIR9dpiPEIfCK/02JmD3SEsZO0r/xbvRWAXwL3hrFrXszs5c5s7CAM\nXtD1GYAP+2XDg/ukoccORf7BfgvwBzN7tKPD76QchX/y8gbwGfCjRZucIKhODGkGQRAETUH08IIg\nCIKmIAxeEAQlSLpVUoctaFnUSLpIUgy3BmHwgmB+kfR9Sc/KHQhPlTtM3iYdGyZpTjr2gaQnJTU0\nsa8qun1y59KfyB0RZ46Qf6yc8141oEGYVtptQs5Xo6TVJF2f7m2m3PHz2ZK+mH4fXiacE1OeLCl3\nkn1G7tiSkp7Jf0KQvkm8R+7ke6akP0naKnc8c8CdOcaelO5zzdy+j9I5s3Lb2+Lu2X6qUn+kQRMS\nBi8I5gNJJwOX4H4mV8E9n1yBOwjOGJU8lKyMe8kYK83roFo5fbMGOc7MeuCeOE7BfVCOKxdHHRwF\n3Jx93JxWnk7A/XkOSPHsjH8jty7+7eMhZcIZAgxP39gdDpwuacN07FTcQ8nFKY51gSeAV/BvGPvg\n34Q+UKZx0Cvl5UG4I+v+1uYsO/MCs0lu32NmNhV4ndJnEjQhYfCCoJ3IPZqcAxxr7qB6lrmD7t+Z\n2WnF89MH+8Nxby/zpbpdjhT/eLxiH4C7j2qUEvkgXBFhJnCwmU1K8bSY2VBzwc8RwDaS8u7Q+gNf\nwz+Cx8xexR2PX5e+1/spcIS1CeIOw12xnWHJUbWZXZbC/mWFe52A69t9tdzxMoynffkRdCHC4AVB\n+xmAe0m5s56T0zDjYbgz7OkLKlHJP+KzuGPgupFr6q1DqR/EnYCxBefd+bgm495WhuR2D8Fd1+Xv\n8Ty8V/gYrtzxSu7YzrR5LMlzO7C1chp8KZ2StDUu+1NJfqnIX5lPhYVg8ScMXhC0n5WA6XV8bD1I\nrgTQAmwG7LvAU5Ycgzd4TSbpklfPWAl3FF2N4SSDl7zZDE77Wkmu0Z5O4d1cuL53hTim4nVU/j6m\nA+/hrsd+bGZ/rJG2jE4rWRMsPEIeKAjazwygt2pLu9xeztVYWtjyu8K+vKzKd6xUTboR+lLqcLoe\nsrh70Ka3NgOfG6zGWOBKSVviKgPdccfdraTFI/vghvBSfOg0Y3qFOFYDPsf9fH4p7evdTm8unVay\nJlh4RA8vCNrPBNx5bruW8JvZ43lJnLQvL5HTLmMn1ynbjFI9uXrSMwv3mpKXD3oI2FdV1ADMBWDH\n4ItXhgC35Z1dpyHJ6/HFKscCG0jKNwAewp18FxmEz+193Mh9VOAruEPooIkJgxcE7cTMPsRXCl4h\nlyDqLpcw2l3SBR0cneQSTq1/ZU7oLml7/JOCP+NSQBlLFK7vVrw+USIfhC826QkMzxamyGWCfp0+\nYcgYjqtqDKQwnAmcDUwysxuTUT0Klx3qnTu+laRz5ZI8PSQdjxvQcmKn7WF74A8dFFawmBIGLwjm\nAzO7CF/JeCaur9aCK0Lf1cFRbYULsrb+yTX+AC6XNBOXSLoElyjarbDQ5KDC9ZWkYa4BBmefNJjZ\neynuOcDTKZ4/4hI0/8hd92jaN9nMnsl2ynXZjkp/pDAfxIdyL03b/wdsgy8qmYTP3Q0EdjWzJ+rL\nnspIWg2XhuroZxIsZoQvzSAISpB0Cz7v2CUMhKSLgDfM7MpFnZZg0RIGLwiCIGgKYkgzCIIgaArC\n4AVBEARNQRi8IAiCoCkIgxcsFkg6StIl6XfmOb9hxwnzc20DcQyW9ECV4ztImryg4u+qSLpR0s8X\ndToWFJK6SfpLWlVa6ZzXJO2wEJPVLiTdIWn32mcuXMLgdUHkkjHvF7+1Svt/WNhXUvkmP4XHSXpZ\n0seSpqXrDqwWTiHM5ZI0S9nvniQdKOnpJOPyz/T7mEre/ZOsy5nAhfXlwKLFzG42s12y7WRg12tv\neGqTAPpI0nRJY7NKMRmBT9Ox9yQ9mFMlyIexQ0rH6YX9X5fLCq2X27eZXMpo7bQ9SdJO6fcX5Ppy\nk9Um03NJe++to5BLKH2mNlmgNyXdIKlf7pySxk6tvCsTZvbXJ3dOWWkoSVflzv9UbRJRFd8L4Ejg\n0aTuUNbAm9lGyUF4p0Hl5ad+CXS6xkkYvC5GqqS2xeVX2iOHchkwFJeZWQl3UXUmsFsDYQzEPZDs\nLGnVQvpOwb+/uhBXDVgFOBrYGqikV/Zd4HUze6eBNHQ1jkvyN/1wn5AX545dkI71Bd7BvZoUORT3\nQVki5WNmLwCXA9emxs7SwG+BszJ1hAI/Ab4BfBN317UD8Hz7b6tDmZDyYXnc6fVs4DlJ1RQVauXd\nhLz8UPqbAtWloczs6Jxk0Xkkiaj0V6nnczSuELHYY2Z/BnrKv8PsNITB63ocAjwF3IhXcnWTWsPH\nAAea2YNmNtvMPksusA5rIKhDgauAl4FWF1Jqk9M5xszGJBkYM7MXzGywmf2nQnhFyZqMwZLeTr2e\nvMDoEnJx0DckzZB0u1zXrdw9j5d0vlyk9N+S7q5y7iOSBqbfW6fewp5pe0dJL6bfrYKtkh5Nl7+U\nWvcH5MI7JfVwp0r6QYV7LyF9CH4HZWRxzGw2rjCwaSHdXwT2x916rV+mEjob91t5JC7d8xFuBMux\nOXCnmU1Jz26Smd1UKb2SLpXUojZx2m1zx4alZ3OTXPT1tXzaUu/z+XRsFK5MUZNUZt8ws2PwcjOs\njmvK5l2V+2pIGqqO8NYEvow72EbSkbgT7tNSubk37c/3todJGi1pZMqjVyT1k/STVK5aJOVHGpZX\nm5DvO5J+LmnJdGwJSWdKeitde1O6x7JD8Fk6JO2Gl5kDUjrz7tvG08kkmcLgdT0Owb3R3wzsKmmV\nBq79Ni5d82x7I5e7n9ohl4Z8j2IA0I2cmnadbEypZE3GNsAGwI7AWXKtNYDjcf+W2+Niou/jLe9K\nHIKLlK4GzMV7ueV4BL83UtgTge1y2/MYZTPLjmeipKPS9qp4T6QvcATunmyFKmkEQO6OayBlZHGS\nYTuIUg8oAPvhRmw0cD+FhlBqaByBD0OdgmvVlZUDwhtTJ8uHoDeWaorMPoMbkRWBW4DRKnWLtjdw\nG95rvYdkaOXD2HfhPZ4VU9oH1oirHGOpQyapSt5VoiFpqDrYGJiYOcY2s2vw9+eCVG72qnDdXnge\nrYCXifvxer0vbpCvzp17I16+1wO+DuwCZFMTh6W/b+GGdzkqN3paMbP7KO3B5iWYOp0kUxi8LoTc\n+/5auJeM53D3Ud9vIIjewLRCmJPl8zmfKCfyWYUhwMtm9he8IttI0tdz4ZfI6Uh6MoU/W9J2ZcID\nrwxnltl/duqFvoQ7Bs5erqOBM8xscqrMhwH7q/JClRFm9mry8/gzXM5nyTLnPUKbn8ntgPNz22UN\nXhXmAOekXsE43CBtUOX8y+RKCi/hrrdOzh07NR2biTcChhSuPRSvkD7Djc6Baegyz6t4ZfiKmb1e\nJR3n44ZxMK65946kiiMJZjbSzGaY2dzkhq1b4T4fN7NxKW0jaHuGWwJL49p5c8xsDG48G6WWTFKt\nvNsylc/sL3PJVq80VL1UKuO1eMzM7k/pGA2sDPzCXGz4NmBtSb1Sw3cPYGjqjf4THxbP5uYHA782\ns4lm9hE+dH1glXemHjqdJFMYvK7FocAD1ia8eQulrfm5eCWSZ2m88oUyUjBmtjpuqLoBtVrz0NbD\nJM25PZJLQ6ucTi78rZJSwAwql8f38fmiInnj/DHeKgU3+ndmlRTe0vwMn2cpR0vu91t4nvQuc94E\noF+qPDYFbgLWSL2ub+L+JOtlRqGyzKe/HCckBYW+afj3X7ljv0p5uDY+b9VqUOTKCd+iTYPubrxn\nUhxqugh/Vqsrt0CpSBouvMLMtsYrs3OB3+Z61yVIOlXSXyV9mJ7F8pTmbfEZLpPKRx/gHSt1BfVW\npXRVoS8+d1mJinmXeMpKFSzWTfvnKcvzSaUyXot3c79n40b4s9w2eLlaCy/XU3PvxdW0yS71oTR/\n38Ll4xoZISrS6SSZwuB1EeQSLIOA7eUrK6cBJwGbSMpazW/jL3aedWgr6A/jFV67JpolbQWsD/wk\nl4YtgO+niiGT0/lug0G/TKlkTS1agN0LFdUyVRa9rJH7vSbeAJhHkdxcpuY54ETgVXMJnCfx3tYb\ntgBVzOvBXOn8ROBStamED8Hf83vT85iIG7zWhlCaE9obd/D8o3R9TfHY1Lu+Aq+s+xePp/m60/By\nuUIyLB9SX8NpKtC3MGS6Zh3XFdmXOmSSKuRdNeZLGqoMLwPrFAxoR/p9bMHT2zv3TvQ0s43S8Sm4\nUcxYE28gvwvMwjUOAUijHyvXkc5OJ8kUBq/rsA/ei+mP9z42xQvcY7TNo40CfiDpm3L64UbxNgAz\n+xve6rtN0s6Slk2Fe6sy8S2lUrmZpfFK9MFCGr4KLIsboA/wBRJXStpfLgOzhKRNgS9WubeiZE0t\nrgLOVZuczcqSqhnZgyX1l9Qdn/cYk2slF3kEV0PIhi/HF7bL8S4+L7LAMVcimIIvQAF/JmfT9jw2\nxefC9pC0Upq7ugY4ycymp+HVByldBdqKpKFpEcOykpZKw5k9KDOnmPbPxVUklpJ0Fi41VA8T0rUn\nyCWX9sN70TWRtKSkdST9Bp9zPbue68rkXbVzO1Qayswm4/OH+XvssHJj/qnDA8BFknqm925duZwU\nwK3ASSnf8itL5wJ/x3vee6b3/Ex8xCefzrU1r2Zip5NkCoPXdTgUuMHM3jazadkfPvE8WK7KfT/w\nY+AGvKU9DtcuuyYXzrH4oo1f40NBk4H/wbXO3s6d97+Uys3cirfkf5OP38zexOdmDgUwswvwHtFp\n+IvyLm5kT6eyQve9wIbKff9Ug0vxBRAPyOVsnsJ7mpUYgU/oT8N7PydUOfcRvCJ/tMJ2OYbhenIf\nSBpUR/rnlwvx1X3b4632KwrP5B68cj0Ir9heN7Obc9cPBXaXtHOZsD/Ghz+n4b3gY4GBZjaxzLn3\nA/fhFeZbuIp6S5nz5iH1nvfDF1K8h5e/sTUuGyDpI+DfeEOkJ7C5mb1ST5yJLO+yCn2A5v0Ob/OU\nxo6Whrqa0jnE64H+qdx0hHLFIfinP3/Be+VjaJvC+C3+HjwKvIk/q+Oh1bgfA1yHf7oxC68XMkan\n/zMkPQ+Q8ugj888TOg2hlhAsFsiXafc3s6EdHO54YKSZXdeR4QZBoyQj+wKwY+qRLbZIugO4Po0Y\ndBoWmHulIOhI0jLtIOiypBXF88yFLo6YWXs+IVngxJBmEARB0BTEkGYQBEHQFEQPLwiCIGgKwuAF\nQRAETUEYvCAIgqApCIMXBEEQNAVh8IIgCIKmIAxeEARB0BSEwQuCIAiagjB4QRAEQVMQBi8IgiBo\nCsLgBUEQBE1BGLwgCIKgKQiDFwRBEDQFYfCCIAiCpiAMXhAEQdAUhMELgiAImoIweEEQBEFTEAYv\nCIIgaArC4AVBEARNQRi8IAiCoCkIgxcEQRA0BWHwgiAIgqbg/wEhzCjD0ElOswAAAABJRU5ErkJg\ngg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt, pandas as pd\n",
"%matplotlib inline\n",
"\n",
"times=pd.DataFrame.from_dict({\"SLSQP (SciPy)\":5.4,\"COBYLA (NLopt)\":35.3, \"CCSAQ (NLopt)\":41.5, \"Ipopt (CasADi)\":0.6,\n",
" 'SLSQP (CasADi)':0.3,'MMA (NLopt)':24.6},orient='index')\n",
"times=times.rename(columns={0:'time'}).sort_values('time',ascending=False)\n",
"times.plot(kind='barh',legend=False)\n",
"_=plt.title(\"Time to solve problem hs103\\n(in seconds)\",size=18)\n",
"_=plt.xlabel(\"\\nThe following solvers didn't manage to solve the problem:\\n\\n\"+\n",
" \"Ipopt through pyipopt (segment fault in MUMPS)\\n\"+\n",
" \"CP+LDL (CVXOPT)\\n\"+\n",
" \"AUGLAG (help with PRAXIS and DIRECT (timeout)\\n\",size=12)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_Note: the comparison between CasADi and the others is not entirely fair, because the other solvers had all the functions and derivatives calculated in sympy objects, which at each evaluation had to be transformed into text and then evaluated, adding some extra overhead._"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 0
}