From 00297d6881aa423c5dd30143e4b539691984d7c0 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Mon, 28 Dec 2020 08:07:10 -0800 Subject: [PATCH 01/18] draft unit test --- control/tests/obc_test.py | 164 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100644 control/tests/obc_test.py diff --git a/control/tests/obc_test.py b/control/tests/obc_test.py new file mode 100644 index 000000000..9b3260d44 --- /dev/null +++ b/control/tests/obc_test.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python +# +# obc_test.py - tests for optimization based control +# RMM, 17 Apr 2019 +# +# This test suite checks the functionality for optimization based control. + +import unittest +import warnings +import numpy as np +import scipy as sp +import control as ct +import control.pwa as pwa +import polytope as pc + +class TestOBC(unittest.TestCase): + def setUp(self): + # Turn off numpy matrix warnings + import warnings + warnings.simplefilter('ignore', category=PendingDeprecationWarning) + + def test_finite_horizon_mpc_simple(self): + # Define a linear system with constraints + # Source: https://www.mpt3.org/UI/RegulationProblem + + # LTI prediction model + model = pwa.ConstrainedAffineSystem( + A = [[1, 1], [0, 1]], B = [[1], [0.5]], C = np.eye(2)) + + # state and input constraints + model.add_state_constraints(pc.box2poly([[-5, 5], [-5, 5]])) + model.add_input_constraints(pc.box2poly([-1, 1])) + + # quadratic state and input penalty + Q = [[1, 0], [0, 1]] + R = [[1]] + + # Create a model predictive controller system + mpc = obc.ModelPredictiveController( + model, + obc.QuadraticCost(Q, R), + horizon=5) + + # Optimal control input for a given value of the initial state: + x0 = [4, 0] + u = mpc.compute_input(x0) + self.assertEqual(u, -1) + + # retrieve the full open-loop predictions + (u_openloop, feasible, openloop) = mpc.compute_trajectory(x0) + np.testing.assert_array_almost_equal( + u_openloop, [-1, -1, 0.1393, 0.3361, -5.2042e-16]) + + # convert it to an explicit form + mpc_explicit = mpc.explicit(); + + # Test explicit controller + (u_explicit, feasible, openloop) = mpc_explicit(x0) + np.testing.assert_array_almost_equal(u_openloop, u_explicit) + + @unittest.skipIf(True, "Not yet implemented.") + def test_finite_horizon_mpc_oscillator(self): + # oscillator model defined in 2D + # Source: https://www.mpt3.org/UI/RegulationProblem + A = [[0.5403, -0.8415], [0.8415, 0.5403]] + B = [[-0.4597], [0.8415]] + C = [[1, 0]] + D = [[0]] + + # Linear discrete-time model with sample time 1 + sys = ss(A, B, C, D, 1); + model = LTISystem(sys); + + # state and input constraints + model.add_state_constraints(pc.box2poly([[-10, 10]])) + model.add_input_constraints(pc.box2poly([-1, 1])) + + # Include weights on states/inputs + model.x.penalty = QuadFunction(np.eye(2)); + model.u.penalty = QuadFunction(1); + + # Compute terminal set + Tset = model.LQRSet; + + # Compute terminal weight + PN = model.LQRPenalty; + + # Add terminal set and terminal penalty + # model.x.with('terminalSet'); + model.x.terminalSet = Tset; + # model.x.with('terminalPenalty'); + model.x.terminalPenalty = PN; + + # Formulate finite horizon MPC problem + ctrl = MPCController(model,5); + + # Add tests to make sure everything works + + # TODO: move this to examples? + @unittest.skipIf(True, "Not yet implemented.") + def test_finite_horizon_mpc_oscillator(self): + # model of an aircraft discretized with 0.2s sampling time + # Source: https://www.mpt3.org/UI/RegulationProblem + A = [[0.99, 0.01, 0.18, -0.09, 0], + [ 0, 0.94, 0, 0.29, 0], + [ 0, 0.14, 0.81, -0.9, 0] + [ 0, -0.2, 0, 0.95, 0], + [ 0, 0.09, 0, 0, 0.9]] + B = [[ 0.01, -0.02], + [-0.14, 0], + [ 0.05, -0.2], + [ 0.02, 0], + [-0.01, 0]] + C = [[0, 1, 0, 0, -1], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [1, 0, 0, 0, 0]] + model = LTISystem('A', A, 'B', B, 'C', C, 'Ts', 0.2); + + # compute the new steady state values for a particular value + # of the input + us = [0.8, -0.3]; + # ys = C*( (eye(5)-A)\B*us ); + + # computed values will be used as references for the desired + # steady state which can be added using "reference" filter + # model.u.with('reference'); + model.u.reference = us; + # model.y.with('reference'); + model.y.reference = ys; + + # provide constraints and penalties on the system signals + model.u.min = [-5, -6]; + model.u.max = [5, 6]; + + # penalties on outputs and inputs are provided as quadratic functions + model.u.penalty = QuadFunction( diag([3, 2]) ); + model.y.penalty = QuadFunction( diag([10, 10, 10, 10]) ); + + # online MPC controller object is constructed with a horizon 6 + ctrl = MPCController(model, 6) + + # loop = ClosedLoop(ctrl, model); + x0 = [0, 0, 0, 0, 0]; + Nsim = 30; + data = loop.simulate(x0, Nsim); + + # Plot the results + subplot(2,1,1) + plot(np.range(Nsim), data.Y); + plot(np.range(Nsim), ys*ones(1, Nsim), 'k--') + title('outputs') + subplot(2,1,2) + plot(np.range(Nsim), data.U); + plot(np.range(Nsim), us*ones(1, Nsim), 'k--') + title('inputs') + + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp) + + +if __name__ == '__main__': + unittest.main() From 23cb793909b7bcb1cd040c8d4fb0971109d7a0c9 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Wed, 30 Dec 2020 13:50:46 -0800 Subject: [PATCH 02/18] convert to pytest --- control/tests/obc_test.py | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/control/tests/obc_test.py b/control/tests/obc_test.py index 9b3260d44..858989d2a 100644 --- a/control/tests/obc_test.py +++ b/control/tests/obc_test.py @@ -1,11 +1,10 @@ -#!/usr/bin/env python -# -# obc_test.py - tests for optimization based control -# RMM, 17 Apr 2019 -# -# This test suite checks the functionality for optimization based control. - -import unittest +"""obc_test.py - tests for optimization based control + +RMM, 17 Apr 2019 check the functionality for optimization based control. +RMM, 30 Dec 2020 convert to pytest +""" + +import pytest import warnings import numpy as np import scipy as sp @@ -13,7 +12,7 @@ import control.pwa as pwa import polytope as pc -class TestOBC(unittest.TestCase): +class TestOBC: def setUp(self): # Turn off numpy matrix warnings import warnings @@ -154,11 +153,3 @@ def test_finite_horizon_mpc_oscillator(self): plot(np.range(Nsim), data.U); plot(np.range(Nsim), us*ones(1, Nsim), 'k--') title('inputs') - - -def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp) - - -if __name__ == '__main__': - unittest.main() From 8711900e309e1c9a9f589383aa2982df52c82669 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Fri, 12 Feb 2021 19:51:09 -0800 Subject: [PATCH 03/18] initial minimal implementation (working) --- control/obc.py | 217 ++++++++++++++++++++++++++++++++++++++ control/tests/obc_test.py | 214 +++++++++++++------------------------ 2 files changed, 289 insertions(+), 142 deletions(-) create mode 100644 control/obc.py diff --git a/control/obc.py b/control/obc.py new file mode 100644 index 000000000..3c502d2d2 --- /dev/null +++ b/control/obc.py @@ -0,0 +1,217 @@ +# obc.py - optimization based control module +# +# RMM, 11 Feb 2021 +# + +"""The "mod:`~control.obc` module provides support for optimization-based +controllers for nonlinear systems with state and input constraints. + +""" + +import numpy as np +import scipy as sp +import scipy.optimize as opt +import control as ct +import warnings + +from .timeresp import _process_time_response + +class ModelPredictiveController(): + """The :class:`ModelPredictiveController` class is a front end for computing + an optimal control input for a nonilinear system with a user-defined cost + function and state and input constraints. + + """ + def __init__( + self, sys, time, integral_cost, trajectory_constraints=[], + terminal_cost=None, terminal_constraints=[]): + + self.system = sys + self.time_vector = time + self.integral_cost = integral_cost + self.trajectory_constraints = trajectory_constraints + self.terminal_cost = terminal_cost + self.terminal_constraints = terminal_constraints + + # + # The approach that we use here is to set up an optimization over the + # inputs at each point in time, using the integral and terminal costs + # as well as the trajectory and terminal constraints. The main work + # of this method is to create the optimization problem that can be + # solved with scipy.optimize.minimize(). + # + + # Gather together all of the constraints + constraint_lb, constraint_ub = [], [] + for time in self.time_vector: + for constraint in self.trajectory_constraints: + type, fun, lb, ub = constraint + constraint_lb.append(lb) + constraint_ub.append(ub) + for constraint in self.terminal_constraints: + type, fun, lb, ub = constraint + constraint_lb.append(lb) + constraint_ub.append(ub) + + # Turn constraint vectors into 1D arrays + self.constraint_lb = np.hstack(constraint_lb) + self.constraint_ub = np.hstack(constraint_ub) + + # Create the new constraint + self.constraints = sp.optimize.NonlinearConstraint( + self.constraint_function, self.constraint_lb, self.constraint_ub) + + # Initial guess + self.initial_guess = np.zeros( + self.system.ninputs * self.time_vector.size) + + # + # Cost function + # + # Given the input U = [u[0], ... u[N]], we need to compute the cost of + # the trajectory generated by that input. This means we have to + # simulate the system to get the state trajectory X = [x[0], ..., + # x[N]] and then compute the cost at each point: + # + # Cost = sum_k integral_cost(x[k], u[k]) + terminal_cost(x[N], u[N]) + # + def cost_function(self, inputs): + # Reshape the input vector + inputs = inputs.reshape( + (self.system.ninputs, self.time_vector.size)) + + # Simulate the system to get the state + _, _, states = ct.input_output_response( + self.system, self.time_vector, inputs, self.x, return_x=True) + + # Trajectory cost + # TODO: vectorize + cost = 0 + for i, time in enumerate(self.time_vector): + cost += self.integral_cost(states[:,i], inputs[:,i]) + + # Terminal cost + if self.terminal_cost is not None: + cost += self.terminal_cost(states[:,-1], inputs[:,-1]) + + # Return the total cost for this input sequence + return cost + + # + # Constraints + # + # We are given the constraints along the trajectory and the terminal + # constraints, which each take inputs [x, u] and evaluate the + # constraint. How we handle these depends on the type of constraint: + # + # We have stored the form of the constraint at a single point, but we + # now need to extend this to apply to each point in the trajectory. + # This means that we need to create N constraints, each of which holds + # at a specific point in time, and implements the original constraint. + # + # To do this, we basically create a function that simulates the system + # dynamics and returns a vector of values corresponding to the value + # of the function at each time. We also replicate the upper and lower + # bounds for each point in time. + # + + # Define a function to evaluate all of the constraints + def constraint_function(self, inputs): + # Reshape the input vector + inputs = inputs.reshape( + (self.system.ninputs, self.time_vector.size)) + + # Simulate the system to get the state + _, _, states = ct.input_output_response( + self.system, self.time_vector, inputs, self.x, return_x=True) + + value = [] + for i, time in enumerate(self.time_vector): + for constraint in self.trajectory_constraints: + type, fun, lb, ub = constraint + if type == opt.LinearConstraint: + value.append( + np.dot(fun, np.hstack([states[:,i], inputs[:,i]]))) + else: + raise TypeError("unknown constraint type %s" % + constraint[0]) + + for constraint in self.terminal_constraints: + type, fun, lb, ub = constraint + if type == opt.LinearConstraint: + value.append( + np.dot(fun, np.hstack([states[:,i], inputs[:,i]]))) + else: + raise TypeError("unknown constraint type %s" % + constraint[0]) + + # Return the value of the constraint function + return np.hstack(value) + + def __call__(self, x): + """Compute the optimal input at state x""" + # Store the starting point + # TODO: call compute_trajectory? + self.x = x + + # Call ScipPy optimizer + res = sp.optimize.minimize( + self.cost_function, self.initial_guess, + constraints=self.constraints) + + # Return the result + if res.success: + return res.x[0] + else: + warnings.warn(res.message) + return None + + def compute_trajectory( + self, x, squeeze=None, transpose=None, return_x=None): + """Compute the optimal input at state x""" + # Store the starting point + self.x = x + + # Call ScipPy optimizer + res = sp.optimize.minimize( + self.cost_function, self.initial_guess, + constraints=self.constraints) + + # See if we got an answer + if not res.success: + warnings.warn(res.message) + return None + + # Reshape the input vector + inputs = res.x.reshape( + (self.system.ninputs, self.time_vector.size)) + + return _process_time_response( + self.system, self.time_vector, inputs, None, + transpose=transpose, return_x=return_x, squeeze=squeeze) + +def state_poly_constraint(sys, polytope): + """Create state constraint from polytope""" + # TODO: make sure the system and constraints are compatible + + # Return a linear constraint object based on the polynomial + return (opt.LinearConstraint, + np.hstack( + [polytope.A, np.zeros((polytope.A.shape[0], sys.ninputs))]), + np.full(polytope.A.shape[0], -np.inf), polytope.b) + + +def input_poly_constraint(sys, polytope): + """Create input constraint from polytope""" + # TODO: make sure the system and constraints are compatible + + # Return a linear constraint object based on the polynomial + return (opt.LinearConstraint, + np.hstack( + [np.zeros((polytope.A.shape[0], sys.nstates)), polytope.A]), + np.full(polytope.A.shape[0], -np.inf), polytope.b) + + +def quadratic_cost(sys, Q, R): + """Create quadratic cost function""" + return lambda x, u: x @ Q @ x + u @ R @ u diff --git a/control/tests/obc_test.py b/control/tests/obc_test.py index 858989d2a..b340d2c12 100644 --- a/control/tests/obc_test.py +++ b/control/tests/obc_test.py @@ -9,147 +9,77 @@ import numpy as np import scipy as sp import control as ct -import control.pwa as pwa +import control.obc as obc import polytope as pc -class TestOBC: - def setUp(self): - # Turn off numpy matrix warnings - import warnings - warnings.simplefilter('ignore', category=PendingDeprecationWarning) - - def test_finite_horizon_mpc_simple(self): - # Define a linear system with constraints - # Source: https://www.mpt3.org/UI/RegulationProblem - - # LTI prediction model - model = pwa.ConstrainedAffineSystem( - A = [[1, 1], [0, 1]], B = [[1], [0.5]], C = np.eye(2)) - - # state and input constraints - model.add_state_constraints(pc.box2poly([[-5, 5], [-5, 5]])) - model.add_input_constraints(pc.box2poly([-1, 1])) - - # quadratic state and input penalty - Q = [[1, 0], [0, 1]] - R = [[1]] - - # Create a model predictive controller system - mpc = obc.ModelPredictiveController( - model, - obc.QuadraticCost(Q, R), - horizon=5) - - # Optimal control input for a given value of the initial state: - x0 = [4, 0] - u = mpc.compute_input(x0) - self.assertEqual(u, -1) - - # retrieve the full open-loop predictions - (u_openloop, feasible, openloop) = mpc.compute_trajectory(x0) - np.testing.assert_array_almost_equal( - u_openloop, [-1, -1, 0.1393, 0.3361, -5.2042e-16]) - - # convert it to an explicit form - mpc_explicit = mpc.explicit(); - - # Test explicit controller - (u_explicit, feasible, openloop) = mpc_explicit(x0) - np.testing.assert_array_almost_equal(u_openloop, u_explicit) - - @unittest.skipIf(True, "Not yet implemented.") - def test_finite_horizon_mpc_oscillator(self): - # oscillator model defined in 2D - # Source: https://www.mpt3.org/UI/RegulationProblem - A = [[0.5403, -0.8415], [0.8415, 0.5403]] - B = [[-0.4597], [0.8415]] - C = [[1, 0]] - D = [[0]] - - # Linear discrete-time model with sample time 1 - sys = ss(A, B, C, D, 1); - model = LTISystem(sys); - - # state and input constraints - model.add_state_constraints(pc.box2poly([[-10, 10]])) - model.add_input_constraints(pc.box2poly([-1, 1])) - - # Include weights on states/inputs - model.x.penalty = QuadFunction(np.eye(2)); - model.u.penalty = QuadFunction(1); - - # Compute terminal set - Tset = model.LQRSet; - - # Compute terminal weight - PN = model.LQRPenalty; - - # Add terminal set and terminal penalty - # model.x.with('terminalSet'); - model.x.terminalSet = Tset; - # model.x.with('terminalPenalty'); - model.x.terminalPenalty = PN; - - # Formulate finite horizon MPC problem - ctrl = MPCController(model,5); - - # Add tests to make sure everything works - - # TODO: move this to examples? - @unittest.skipIf(True, "Not yet implemented.") - def test_finite_horizon_mpc_oscillator(self): - # model of an aircraft discretized with 0.2s sampling time - # Source: https://www.mpt3.org/UI/RegulationProblem - A = [[0.99, 0.01, 0.18, -0.09, 0], - [ 0, 0.94, 0, 0.29, 0], - [ 0, 0.14, 0.81, -0.9, 0] - [ 0, -0.2, 0, 0.95, 0], - [ 0, 0.09, 0, 0, 0.9]] - B = [[ 0.01, -0.02], - [-0.14, 0], - [ 0.05, -0.2], - [ 0.02, 0], - [-0.01, 0]] - C = [[0, 1, 0, 0, -1], - [0, 0, 1, 0, 0], - [0, 0, 0, 1, 0], - [1, 0, 0, 0, 0]] - model = LTISystem('A', A, 'B', B, 'C', C, 'Ts', 0.2); - - # compute the new steady state values for a particular value - # of the input - us = [0.8, -0.3]; - # ys = C*( (eye(5)-A)\B*us ); - - # computed values will be used as references for the desired - # steady state which can be added using "reference" filter - # model.u.with('reference'); - model.u.reference = us; - # model.y.with('reference'); - model.y.reference = ys; - - # provide constraints and penalties on the system signals - model.u.min = [-5, -6]; - model.u.max = [5, 6]; - - # penalties on outputs and inputs are provided as quadratic functions - model.u.penalty = QuadFunction( diag([3, 2]) ); - model.y.penalty = QuadFunction( diag([10, 10, 10, 10]) ); - - # online MPC controller object is constructed with a horizon 6 - ctrl = MPCController(model, 6) - - # loop = ClosedLoop(ctrl, model); - x0 = [0, 0, 0, 0, 0]; - Nsim = 30; - data = loop.simulate(x0, Nsim); - - # Plot the results - subplot(2,1,1) - plot(np.range(Nsim), data.Y); - plot(np.range(Nsim), ys*ones(1, Nsim), 'k--') - title('outputs') - subplot(2,1,2) - plot(np.range(Nsim), data.U); - plot(np.range(Nsim), us*ones(1, Nsim), 'k--') - title('inputs') +def test_finite_horizon_mpc_simple(): + # Define a linear system with constraints + # Source: https://www.mpt3.org/UI/RegulationProblem + + # LTI prediction model + sys = ct.ss2io(ct.ss([[1, 1], [0, 1]], [[1], [0.5]], np.eye(2), 0, 1)) + + # State and input constraints + constraints = [ + obc.state_poly_constraint(sys, pc.box2poly([[-5, 5], [-5, 5]])), + obc.input_poly_constraint(sys, pc.box2poly([[-1, 1]])), + ] + + # Quadratic state and input penalty + Q = [[1, 0], [0, 1]] + R = [[1]] + cost = obc.quadratic_cost(sys, Q, R) + + # Create a model predictive controller system + time = np.arange(0, 5, 1) + mpc = obc.ModelPredictiveController(sys, time, cost, constraints) + + # Optimal control input for a given value of the initial state + x0 = [4, 0] + u = mpc(x0) + np.testing.assert_almost_equal(u, -1) + + # Retrieve the full open-loop predictions + t, u_openloop = mpc.compute_trajectory(x0, squeeze=True) + np.testing.assert_almost_equal( + u_openloop, [-1, -1, 0.1393, 0.3361, -5.204e-16], decimal=4) + + # Convert controller to an explicit form (not implemented yet) + # mpc_explicit = mpc.explicit(); + + # Test explicit controller + # u_explicit = mpc_explicit(x0) + # np.testing.assert_array_almost_equal(u_openloop, u_explicit) + +def test_finite_horizon_mpc_oscillator(): + # oscillator model defined in 2D + # Source: https://www.mpt3.org/UI/RegulationProblem + A = [[0.5403, -0.8415], [0.8415, 0.5403]] + B = [[-0.4597], [0.8415]] + C = [[1, 0]] + D = [[0]] + + # Linear discrete-time model with sample time 1 + sys = ct.ss2io(ct.ss(A, B, C, D, 1)) + + # state and input constraints + trajectory_constraints = [ + obc.state_poly_constraint(sys, pc.box2poly([[-10, 10]])), + obc.input_poly_constraint(sys, pc.box2poly([[-1, 1]])) + ] + + # Include weights on states/inputs + Q = np.eye(2) + R = 1 + K, S, E = ct.lqr(A, B, Q, R) + + # Compute the integral and terminal cost + integral_cost = obc.quadratic_cost(sys, Q, R) + terminal_cost = obc.quadratic_cost(sys, S, 0) + + # Formulate finite horizon MPC problem + time = np.arange(0, 5, 5) + mpc = obc.ModelPredictiveController( + sys, time, integral_cost, trajectory_constraints, terminal_cost) + + # Add tests to make sure everything works From 6f9c092b26f0167c5f5f47ccdeb182475a835f5d Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 13 Feb 2021 13:57:37 -0800 Subject: [PATCH 04/18] minor refactoring plus additional comments on method --- control/obc.py | 221 ++++++++++++++++++++++++++++++-------- control/tests/obc_test.py | 14 +-- 2 files changed, 186 insertions(+), 49 deletions(-) diff --git a/control/obc.py b/control/obc.py index 3c502d2d2..e96339c39 100644 --- a/control/obc.py +++ b/control/obc.py @@ -16,16 +16,46 @@ from .timeresp import _process_time_response -class ModelPredictiveController(): - """The :class:`ModelPredictiveController` class is a front end for computing - an optimal control input for a nonilinear system with a user-defined cost +# +# OptimalControlProblem class +# +# The OptimalControlProblem class holds all of the information required to +# specify and optimal control problem: the system dynamics, cost function, +# and constraints. As much as possible, the information used to specify an +# optimal control problem matches the notation and terminology of the SciPy +# `optimize.minimize` module, with the hope that this makes it easier to +# remember how to describe a problem. +# +# The approach that we use here is to set up an optimization over the +# inputs at each point in time, using the integral and terminal costs as +# well as the trajectory and terminal constraints. The main function of +# this class is to create an optimization problem that can be solved using +# scipy.optimize.minimize(). +# +# The `cost_function` method takes the information stored here and computes +# the cost of the trajectory generated by the proposed input. It does this +# by calling a user-defined function for the integral_cost given the +# current states and inputs at each point along the trajetory and then +# adding the value of a user-defined terminal cost at the final pint in the +# trajectory. +# +# The `constraint_function` method evaluates the constraint functions along +# the trajectory generated by the proposed input. As in the case of the +# cost function, the constraints are evaluated at the state and input along +# each point on the trjectory. This information is compared against the +# constraint upper and lower bounds. The constraint function is processed +# in the class initializer, so that it only needs to be computed once. +# +class OptimalControlProblem(): + """The :class:`OptimalControlProblem` class is a front end for computing an + optimal control input for a nonilinear system with a user-defined cost function and state and input constraints. """ def __init__( self, sys, time, integral_cost, trajectory_constraints=[], terminal_cost=None, terminal_constraints=[]): - + # Save the basic information for use later self.system = sys self.time_vector = time self.integral_cost = integral_cost @@ -34,20 +64,28 @@ def __init__( self.terminal_constraints = terminal_constraints # - # The approach that we use here is to set up an optimization over the - # inputs at each point in time, using the integral and terminal costs - # as well as the trajectory and terminal constraints. The main work - # of this method is to create the optimization problem that can be - # solved with scipy.optimize.minimize(). + # Compute and store constraints + # + # While the constraints are evaluated during the execution of the + # SciPy optimization method itself, we go ahead and pre-compute the + # `scipy.optimize.NonlinearConstraint` function that will be passed to + # the optimizer on initialization, since it doesn't change. This is + # mainly a matter of computing the lower and upper bound vectors, + # which we need to "stack" to account for the evaluation at each + # trajectory time point plus any terminal constraints (in a way that + # is consistent with the `constraint_function` that is used at + # evaluation time. # - - # Gather together all of the constraints constraint_lb, constraint_ub = [], [] + + # Go through each time point and stack the bounds for time in self.time_vector: for constraint in self.trajectory_constraints: type, fun, lb, ub = constraint constraint_lb.append(lb) constraint_ub.append(ub) + + # Add on the terminal constraints for constraint in self.terminal_constraints: type, fun, lb, ub = constraint constraint_lb.append(lb) @@ -60,8 +98,15 @@ def __init__( # Create the new constraint self.constraints = sp.optimize.NonlinearConstraint( self.constraint_function, self.constraint_lb, self.constraint_ub) - + + # # Initial guess + # + # We store an initial guess (zero input) in case it is not specified + # later. + # + # TODO: add the ability to overwride this when calling the optimizer. + # self.initial_guess = np.zeros( self.system.ninputs * self.time_vector.size) @@ -75,14 +120,18 @@ def __init__( # # Cost = sum_k integral_cost(x[k], u[k]) + terminal_cost(x[N], u[N]) # + # The initial state is for generating the simulation is store in the class + # parameter `x` prior to calling the optimization algorithm. + # def cost_function(self, inputs): - # Reshape the input vector + # Retrieve the initial state and reshape the input vector + x = self.x inputs = inputs.reshape( (self.system.ninputs, self.time_vector.size)) # Simulate the system to get the state _, _, states = ct.input_output_response( - self.system, self.time_vector, inputs, self.x, return_x=True) + self.system, self.time_vector, inputs, x, return_x=True) # Trajectory cost # TODO: vectorize @@ -104,43 +153,71 @@ def cost_function(self, inputs): # constraints, which each take inputs [x, u] and evaluate the # constraint. How we handle these depends on the type of constraint: # - # We have stored the form of the constraint at a single point, but we - # now need to extend this to apply to each point in the trajectory. - # This means that we need to create N constraints, each of which holds - # at a specific point in time, and implements the original constraint. + # * For linear constraints (LinearConstraint), a combined vector of the + # state and input is multiplied by the polytope A matrix for + # comparison against the upper and lower bounds. + # + # * For nonlinear constraints (NonlinearConstraint), a user-specific + # constraint function having the form + # + # constraint_fun(x, u) + # + # is called at each point along the trajectory and compared against the + # upper and lower bounds. + # + # In both cases, the constraint is specified at a single point, but we + # extend this to apply to each point in the trajectory. This means + # that for N time points with m trajectory constraints and p terminal + # constraints we need to compute N*m + p constraints, each of which + # holds at a specific point in time, and implements the original + # constraint. # # To do this, we basically create a function that simulates the system - # dynamics and returns a vector of values corresponding to the value - # of the function at each time. We also replicate the upper and lower - # bounds for each point in time. + # dynamics and returns a vector of values corresponding to the value of + # the function at each time. The class initialization methods takes + # care of replicating the upper and lower bounds for each point in time + # so that the SciPy optimization algorithm can do the proper + # evaluation. + # + # In addition, since SciPy's optimization function does not allow us to + # pass arguments to the constraint function, we have to store the initial + # state prior to optimization and retrieve it here. # - - # Define a function to evaluate all of the constraints def constraint_function(self, inputs): - # Reshape the input vector + # Retrieve the initial state and reshape the input vector + x = self.x inputs = inputs.reshape( (self.system.ninputs, self.time_vector.size)) # Simulate the system to get the state _, _, states = ct.input_output_response( - self.system, self.time_vector, inputs, self.x, return_x=True) + self.system, self.time_vector, inputs, x, return_x=True) + # Evaluate the constraint function along the trajectory value = [] for i, time in enumerate(self.time_vector): for constraint in self.trajectory_constraints: type, fun, lb, ub = constraint if type == opt.LinearConstraint: + # `fun` is the A matrix associated with the polytope... value.append( np.dot(fun, np.hstack([states[:,i], inputs[:,i]]))) + elif type == opt.NonlinearConstraint: + value.append( + fun(np.hstack([states[:,i], inputs[:,i]]))) else: raise TypeError("unknown constraint type %s" % constraint[0]) + # Evaluate the terminal constraint functions for constraint in self.terminal_constraints: type, fun, lb, ub = constraint if type == opt.LinearConstraint: value.append( np.dot(fun, np.hstack([states[:,i], inputs[:,i]]))) + elif type == opt.NonlinearConstraint: + value.append( + fun(np.hstack([states[:,i], inputs[:,i]]))) else: raise TypeError("unknown constraint type %s" % constraint[0]) @@ -148,28 +225,22 @@ def constraint_function(self, inputs): # Return the value of the constraint function return np.hstack(value) - def __call__(self, x): + # Allow optctrl(x) as a replacement for optctrl.mpc(x) + def __call__(self, x, squeeze=None): """Compute the optimal input at state x""" - # Store the starting point - # TODO: call compute_trajectory? - self.x = x - - # Call ScipPy optimizer - res = sp.optimize.minimize( - self.cost_function, self.initial_guess, - constraints=self.constraints) + return self.mpc(x, squeeze=squeeze) - # Return the result - if res.success: - return res.x[0] - else: - warnings.warn(res.message) - return None + # Compute the current input to apply from the current state (MPC style) + def mpc(self, x, squeeze=None): + """Compute the optimal input at state x""" + _, inputs = self.compute_trajectory(x, squeeze=squeeze) + return None if inputs is None else inputs.transpose()[0] + # Compute the optimal trajectory from the current state def compute_trajectory( self, x, squeeze=None, transpose=None, return_x=None): """Compute the optimal input at state x""" - # Store the starting point + # Store the initial state (for use in constraint_function) self.x = x # Call ScipPy optimizer @@ -185,11 +256,33 @@ def compute_trajectory( # Reshape the input vector inputs = res.x.reshape( (self.system.ninputs, self.time_vector.size)) + + if return_x: + # Simulate the system if we need the state back + _, _, states = ct.input_output_response( + self.system, self.time_vector, inputs, x, return_x=True) + else: + states=None return _process_time_response( - self.system, self.time_vector, inputs, None, + self.system, self.time_vector, inputs, states, transpose=transpose, return_x=return_x, squeeze=squeeze) + +# +# Create a polytope constraint on the system state +# +# As in the cost function evaluation, the main "trick" in creating a constrain +# on the state or input is to properly evaluate the constraint on the stacked +# state and input vector at the current time point. The constraint itself +# will be called at each poing along the trajectory (or the endpoint) via the +# constrain_function() method. +# +# Note that these functions to not actually evaluate the constraint, they +# simply return the information required to do so. We use the SciPy +# optimization methods LinearConstraint and NonlinearConstraint as "types" to +# keep things consistent with the terminology in scipy.optimize. +# def state_poly_constraint(sys, polytope): """Create state constraint from polytope""" # TODO: make sure the system and constraints are compatible @@ -200,7 +293,7 @@ def state_poly_constraint(sys, polytope): [polytope.A, np.zeros((polytope.A.shape[0], sys.ninputs))]), np.full(polytope.A.shape[0], -np.inf), polytope.b) - +# Create a constraint polytope on the system input def input_poly_constraint(sys, polytope): """Create input constraint from polytope""" # TODO: make sure the system and constraints are compatible @@ -212,6 +305,48 @@ def input_poly_constraint(sys, polytope): np.full(polytope.A.shape[0], -np.inf), polytope.b) +# +# Create a constraint polytope on the system output +# +# Unlike the state and input constraints, for the output constraint we need to +# do a function evaluation before applying the constraints. +# +# TODO: for the special case of an LTI system, we can avoid the extra function +# call by multiplying the state by the C matrix for the system and then +# imposing a linear constraint: +# +# np.hstack( +# [polytope.A @ sys.C, np.zeros((polytope.A.shape[0], sys.ninputs))]) +# +def output_poly_constraint(sys, polytope): + """Create output constraint from polytope""" + # TODO: make sure the system and constraints are compatible + + # + # Function to create the output + def _evaluate_output_constraint(x): + # Separate the constraint into states and inputs + states = x[:sys.nstates] + inputs = x[sys.nstates:] + outputs = sys._out(0, states, inputs) + return polytope.A @ outputs + + # Return a nonlinear constraint object based on the polynomial + return (opt.NonlinearConstraint, + _evaluate_output_constraint, + np.full(polytope.A.shape[0], -np.inf), polytope.b) + + +# +# Quadratic cost function +# +# Since a quadratic function is common as a cost function, we provide a +# function that will take a Q and R matrix and return a callable that +# evaluates to associted quadratic cost. This is compatible with the way that +# the `cost_function` evaluates the cost at each point in the trajectory. +# def quadratic_cost(sys, Q, R): """Create quadratic cost function""" + Q = np.atleast_2d(Q) + R = np.atleast_2d(R) return lambda x, u: x @ Q @ x + u @ R @ u diff --git a/control/tests/obc_test.py b/control/tests/obc_test.py index b340d2c12..f04e9cc91 100644 --- a/control/tests/obc_test.py +++ b/control/tests/obc_test.py @@ -32,7 +32,8 @@ def test_finite_horizon_mpc_simple(): # Create a model predictive controller system time = np.arange(0, 5, 1) - mpc = obc.ModelPredictiveController(sys, time, cost, constraints) + optctrl = obc.OptimalControlProblem(sys, time, cost, constraints) + mpc = optctrl.mpc # Optimal control input for a given value of the initial state x0 = [4, 0] @@ -40,12 +41,12 @@ def test_finite_horizon_mpc_simple(): np.testing.assert_almost_equal(u, -1) # Retrieve the full open-loop predictions - t, u_openloop = mpc.compute_trajectory(x0, squeeze=True) + t, u_openloop = optctrl.compute_trajectory(x0, squeeze=True) np.testing.assert_almost_equal( u_openloop, [-1, -1, 0.1393, 0.3361, -5.204e-16], decimal=4) # Convert controller to an explicit form (not implemented yet) - # mpc_explicit = mpc.explicit(); + # mpc_explicit = obc.explicit_mpc(); # Test explicit controller # u_explicit = mpc_explicit(x0) @@ -64,7 +65,7 @@ def test_finite_horizon_mpc_oscillator(): # state and input constraints trajectory_constraints = [ - obc.state_poly_constraint(sys, pc.box2poly([[-10, 10]])), + obc.output_poly_constraint(sys, pc.box2poly([[-10, 10]])), obc.input_poly_constraint(sys, pc.box2poly([[-1, 1]])) ] @@ -78,8 +79,9 @@ def test_finite_horizon_mpc_oscillator(): terminal_cost = obc.quadratic_cost(sys, S, 0) # Formulate finite horizon MPC problem - time = np.arange(0, 5, 5) - mpc = obc.ModelPredictiveController( + time = np.arange(0, 5, 1) + optctrl = obc.OptimalControlProblem( sys, time, integral_cost, trajectory_constraints, terminal_cost) # Add tests to make sure everything works + t, u_openloop = optctrl.compute_trajectory([1, 1]) From bd322e7befd85e67492d87f87e5291bfcf800357 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sun, 14 Feb 2021 21:59:49 -0800 Subject: [PATCH 05/18] remove polytope dependence; implement MPC iosys w/ notebook example --- control/obc.py | 91 +++++++++++----- control/tests/obc_test.py | 66 ++++++++++-- examples/mpc_aircraft.ipynb | 202 ++++++++++++++++++++++++++++++++++++ 3 files changed, 328 insertions(+), 31 deletions(-) create mode 100644 examples/mpc_aircraft.ipynb diff --git a/control/obc.py b/control/obc.py index e96339c39..ff56fbb43 100644 --- a/control/obc.py +++ b/control/obc.py @@ -103,9 +103,8 @@ def __init__( # Initial guess # # We store an initial guess (zero input) in case it is not specified - # later. - # - # TODO: add the ability to overwride this when calling the optimizer. + # later. Note that create_mpc_iosystem() will reset the initial guess + # based on the current state of the MPC controller. # self.initial_guess = np.zeros( self.system.ninputs * self.time_vector.size) @@ -128,24 +127,25 @@ def cost_function(self, inputs): x = self.x inputs = inputs.reshape( (self.system.ninputs, self.time_vector.size)) - + # Simulate the system to get the state _, _, states = ct.input_output_response( self.system, self.time_vector, inputs, x, return_x=True) - + # Trajectory cost # TODO: vectorize cost = 0 for i, time in enumerate(self.time_vector): - cost += self.integral_cost(states[:,i], inputs[:,i]) - + cost += self.integral_cost(states[:,i], inputs[:,i]) + # Terminal cost if self.terminal_cost is not None: cost += self.terminal_cost(states[:,-1], inputs[:,-1]) - + # Return the total cost for this input sequence return cost + # # Constraints # @@ -188,7 +188,7 @@ def constraint_function(self, inputs): x = self.x inputs = inputs.reshape( (self.system.ninputs, self.time_vector.size)) - + # Simulate the system to get the state _, _, states = ct.input_output_response( self.system, self.time_vector, inputs, x, return_x=True) @@ -234,7 +234,7 @@ def __call__(self, x, squeeze=None): def mpc(self, x, squeeze=None): """Compute the optimal input at state x""" _, inputs = self.compute_trajectory(x, squeeze=squeeze) - return None if inputs is None else inputs.transpose()[0] + return None if inputs is None else inputs[:,0] # Compute the optimal trajectory from the current state def compute_trajectory( @@ -242,7 +242,7 @@ def compute_trajectory( """Compute the optimal input at state x""" # Store the initial state (for use in constraint_function) self.x = x - + # Call ScipPy optimizer res = sp.optimize.minimize( self.cost_function, self.initial_guess, @@ -263,14 +263,33 @@ def compute_trajectory( self.system, self.time_vector, inputs, x, return_x=True) else: states=None - + return _process_time_response( self.system, self.time_vector, inputs, states, transpose=transpose, return_x=return_x, squeeze=squeeze) + # Create an input/output system implementing an MPC controller + def create_mpc_iosystem(self, dt=True): + def _update(t, x, u, params={}): + inputs = x.reshape((self.system.ninputs, self.time_vector.size)) + self.initial_guess = np.hstack( + [inputs[:,1:], inputs[:,-1:]]).reshape(-1) + _, inputs = self.compute_trajectory(u) + return inputs.reshape(-1) + + def _output(t, x, u, params={}): + inputs = x.reshape((self.system.ninputs, self.time_vector.size)) + return inputs[:,0] + + return ct.NonlinearIOSystem( + _update, _output, dt=dt, + inputs=self.system.nstates, + outputs=self.system.ninputs, + states=self.system.ninputs * self.time_vector.size) + # -# Create a polytope constraint on the system state +# Create a polytope constraint on the system state: A x <= b # # As in the cost function evaluation, the main "trick" in creating a constrain # on the state or input is to properly evaluate the constraint on the stacked @@ -283,26 +302,48 @@ def compute_trajectory( # optimization methods LinearConstraint and NonlinearConstraint as "types" to # keep things consistent with the terminology in scipy.optimize. # -def state_poly_constraint(sys, polytope): +def state_poly_constraint(sys, A, b): + """Create state constraint from polytope""" + # TODO: make sure the system and constraints are compatible + + # Return a linear constraint object based on the polynomial + return (opt.LinearConstraint, + np.hstack([A, np.zeros((A.shape[0], sys.ninputs))]), + np.full(A.shape[0], -np.inf), polytope.b) + + +def state_range_constraint(sys, lb, ub): """Create state constraint from polytope""" # TODO: make sure the system and constraints are compatible # Return a linear constraint object based on the polynomial return (opt.LinearConstraint, np.hstack( - [polytope.A, np.zeros((polytope.A.shape[0], sys.ninputs))]), - np.full(polytope.A.shape[0], -np.inf), polytope.b) + [np.eye(sys.nstates), np.zeros((sys.nstates, sys.ninputs))]), + np.array(lb), np.array(ub)) + # Create a constraint polytope on the system input -def input_poly_constraint(sys, polytope): +def input_poly_constraint(sys, A, b): + """Create input constraint from polytope""" + # TODO: make sure the system and constraints are compatible + + # Return a linear constraint object based on the polynomial + return (opt.LinearConstraint, + np.hstack( + [np.zeros((A.shape[0], sys.nstates)), A]), + np.full(A.shape[0], -np.inf), b) + + +def input_range_constraint(sys, lb, ub): """Create input constraint from polytope""" # TODO: make sure the system and constraints are compatible # Return a linear constraint object based on the polynomial return (opt.LinearConstraint, np.hstack( - [np.zeros((polytope.A.shape[0], sys.nstates)), polytope.A]), - np.full(polytope.A.shape[0], -np.inf), polytope.b) + [np.zeros((sys.ninputs, sys.nstates)), np.eye(sys.ninputs)]), + np.array(lb), np.array(ub)) # @@ -316,9 +357,9 @@ def input_poly_constraint(sys, polytope): # imposing a linear constraint: # # np.hstack( -# [polytope.A @ sys.C, np.zeros((polytope.A.shape[0], sys.ninputs))]) +# [A @ sys.C, np.zeros((A.shape[0], sys.ninputs))]) # -def output_poly_constraint(sys, polytope): +def output_poly_constraint(sys, A, b): """Create output constraint from polytope""" # TODO: make sure the system and constraints are compatible @@ -329,12 +370,12 @@ def _evaluate_output_constraint(x): states = x[:sys.nstates] inputs = x[sys.nstates:] outputs = sys._out(0, states, inputs) - return polytope.A @ outputs + return A @ outputs # Return a nonlinear constraint object based on the polynomial return (opt.NonlinearConstraint, _evaluate_output_constraint, - np.full(polytope.A.shape[0], -np.inf), polytope.b) + np.full(A.shape[0], -np.inf), b) # @@ -345,8 +386,8 @@ def _evaluate_output_constraint(x): # evaluates to associted quadratic cost. This is compatible with the way that # the `cost_function` evaluates the cost at each point in the trajectory. # -def quadratic_cost(sys, Q, R): +def quadratic_cost(sys, Q, R, x0=0, u0=0): """Create quadratic cost function""" Q = np.atleast_2d(Q) R = np.atleast_2d(R) - return lambda x, u: x @ Q @ x + u @ R @ u + return lambda x, u: ((x-x0) @ Q @ (x-x0) + (u-u0) @ R @ (u-u0)).item() diff --git a/control/tests/obc_test.py b/control/tests/obc_test.py index f04e9cc91..a98283034 100644 --- a/control/tests/obc_test.py +++ b/control/tests/obc_test.py @@ -10,7 +10,7 @@ import scipy as sp import control as ct import control.obc as obc -import polytope as pc +from control.tests.conftest import slycotonly def test_finite_horizon_mpc_simple(): # Define a linear system with constraints @@ -21,8 +21,7 @@ def test_finite_horizon_mpc_simple(): # State and input constraints constraints = [ - obc.state_poly_constraint(sys, pc.box2poly([[-5, 5], [-5, 5]])), - obc.input_poly_constraint(sys, pc.box2poly([[-1, 1]])), + (sp.optimize.LinearConstraint, np.eye(3), [-5, -5, -1], [5, 5, 1]), ] # Quadratic state and input penalty @@ -48,10 +47,12 @@ def test_finite_horizon_mpc_simple(): # Convert controller to an explicit form (not implemented yet) # mpc_explicit = obc.explicit_mpc(); - # Test explicit controller + # Test explicit controller # u_explicit = mpc_explicit(x0) # np.testing.assert_array_almost_equal(u_openloop, u_explicit) + +@slycotonly def test_finite_horizon_mpc_oscillator(): # oscillator model defined in 2D # Source: https://www.mpt3.org/UI/RegulationProblem @@ -65,8 +66,7 @@ def test_finite_horizon_mpc_oscillator(): # state and input constraints trajectory_constraints = [ - obc.output_poly_constraint(sys, pc.box2poly([[-10, 10]])), - obc.input_poly_constraint(sys, pc.box2poly([[-1, 1]])) + (sp.optimize.LinearConstraint, np.eye(3), [-10, -10, -1], [10, 10, 1]), ] # Include weights on states/inputs @@ -85,3 +85,57 @@ def test_finite_horizon_mpc_oscillator(): # Add tests to make sure everything works t, u_openloop = optctrl.compute_trajectory([1, 1]) + + +def test_mpc_iosystem(): + # model of an aircraft discretized with 0.2s sampling time + # Source: https://www.mpt3.org/UI/RegulationProblem + A = [[0.99, 0.01, 0.18, -0.09, 0], + [ 0, 0.94, 0, 0.29, 0], + [ 0, 0.14, 0.81, -0.9, 0], + [ 0, -0.2, 0, 0.95, 0], + [ 0, 0.09, 0, 0, 0.9]] + B = [[ 0.01, -0.02], + [-0.14, 0], + [ 0.05, -0.2], + [ 0.02, 0], + [-0.01, 0]] + C = [[0, 1, 0, 0, -1], + [0, 0, 1, 0, 0], + [0, 0, 0, 1, 0], + [1, 0, 0, 0, 0]] + model = ct.ss2io(ct.ss(A, B, C, 0, 0.2)) + + # For the simulation we need the full state output + sys = ct.ss2io(ct.ss(A, B, np.eye(5), 0, 0.2)) + + # compute the steady state values for a particular value of the input + ud = np.array([0.8, -0.3]) + xd = np.linalg.inv(np.eye(5) - A) @ B @ ud + yd = C @ xd + + # provide constraints on the system signals + constraints = [obc.input_range_constraint(sys, [-5, -6], [5, 6])] + + # provide penalties on the system signals + Q = model.C.transpose() @ np.diag([10, 10, 10, 10]) @ model.C + R = np.diag([3, 2]) + cost = obc.quadratic_cost(model, Q, R, x0=xd, u0=ud) + + # online MPC controller object is constructed with a horizon 6 + optctrl = obc.OptimalControlProblem( + model, np.arange(0, 6) * 0.2, cost, constraints) + + # Define an I/O system implementing model predictive control + ctrl = optctrl.create_mpc_iosystem() + loop = ct.feedback(sys, ctrl, 1) + + # Choose a nearby initial condition to speed up computation + X0 = np.hstack([xd, np.kron(ud, np.ones(6))]) * 0.99 + + Nsim = 10 + tout, xout = ct.input_output_response( + loop, np.arange(0, Nsim) * 0.2, 0, X0) + + # Make sure the system converged to the desired state + np.testing.assert_almost_equal(xout[0:sys.nstates, -1], xd, decimal=1) diff --git a/examples/mpc_aircraft.ipynb b/examples/mpc_aircraft.ipynb new file mode 100644 index 000000000..53b8bb13b --- /dev/null +++ b/examples/mpc_aircraft.ipynb @@ -0,0 +1,202 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Model Predictive Control: Aircraft Model\n", + "\n", + "RMM, 13 Feb 2021\n", + "\n", + "This example replicates the [MPT3 regulation problem example](https://www.mpt3.org/UI/RegulationProblem)." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import control as ct\n", + "import numpy as np\n", + "import control.obc as obc\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# model of an aircraft discretized with 0.2s sampling time\n", + "# Source: https://www.mpt3.org/UI/RegulationProblem\n", + "A = [[0.99, 0.01, 0.18, -0.09, 0],\n", + " [ 0, 0.94, 0, 0.29, 0],\n", + " [ 0, 0.14, 0.81, -0.9, 0],\n", + " [ 0, -0.2, 0, 0.95, 0],\n", + " [ 0, 0.09, 0, 0, 0.9]]\n", + "B = [[ 0.01, -0.02],\n", + " [-0.14, 0],\n", + " [ 0.05, -0.2],\n", + " [ 0.02, 0],\n", + " [-0.01, 0]]\n", + "C = [[0, 1, 0, 0, -1],\n", + " [0, 0, 1, 0, 0],\n", + " [0, 0, 0, 1, 0],\n", + " [1, 0, 0, 0, 0]]\n", + "model = ct.ss2io(ct.ss(A, B, C, 0, 0.2))\n", + "\n", + "# For the simulation we need the full state output\n", + "sys = ct.ss2io(ct.ss(A, B, np.eye(5), 0, 0.2))\n", + "\n", + "# compute the steady state values for a particular value of the input\n", + "ud = np.array([0.8, -0.3])\n", + "xd = np.linalg.inv(np.eye(5) - A) @ B @ ud\n", + "yd = C @ xd" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# computed values will be used as references for the desired\n", + "# steady state which can be added using \"reference\" filter\n", + "# model.u.with('reference');\n", + "# model.u.reference = us;\n", + "# model.y.with('reference');\n", + "# model.y.reference = ys;\n", + "\n", + "# provide constraints on the system signals\n", + "constraints = [obc.input_range_constraint(sys, [-5, -6], [5, 6])]\n", + "\n", + "# provide penalties on the system signals\n", + "Q = model.C.transpose() @ np.diag([10, 10, 10, 10]) @ model.C\n", + "R = np.diag([3, 2])\n", + "cost = obc.quadratic_cost(model, Q, R, x0=xd, u0=ud)\n", + "\n", + "# online MPC controller object is constructed with a horizon 6\n", + "optctrl = obc.OptimalControlProblem(model, np.arange(0, 6) * 0.2, cost, constraints)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "System: sys[7]\n", + "Inputs (2): u[0], u[1], \n", + "Outputs (5): y[0], y[1], y[2], y[3], y[4], \n", + "States (17): sys[1]_x[0], sys[1]_x[1], sys[1]_x[2], sys[1]_x[3], sys[1]_x[4], sys[6]_x[0], sys[6]_x[1], sys[6]_x[2], sys[6]_x[3], sys[6]_x[4], sys[6]_x[5], sys[6]_x[6], sys[6]_x[7], sys[6]_x[8], sys[6]_x[9], sys[6]_x[10], sys[6]_x[11], \n" + ] + } + ], + "source": [ + "# Define an I/O system implementing model predictive control\n", + "ctrl = optctrl.create_mpc_iosystem()\n", + "loop = ct.feedback(sys, ctrl, 1)\n", + "print(loop)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Computation time = 8.28132 seconds\n" + ] + } + ], + "source": [ + "import time\n", + "\n", + "# loop = ClosedLoop(ctrl, model);\n", + "# x0 = [0, 0, 0, 0, 0]\n", + "Nsim = 60\n", + "\n", + "start = time.time()\n", + "tout, xout = ct.input_output_response(loop, np.arange(0, Nsim) * 0.2, 0, 0)\n", + "end = time.time()\n", + "print(\"Computation time = %g seconds\" % (end-start))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([-0.15441833, 0.00362039, 0.07760278, 0.00675162, 0.00698118])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9UUlEQVR4nO3deXhU1fnA8e87a/YEkrAl7LK4oyKKirUVKypqXUCrUqkLFreiaAXF2v7EBRUES8GtSutStahVUFsF676ComyCIEsSIAshezLr+f1xh5CQsIRkuJPk/TzPfe5y7r3zThjmnXPuueeKMQallFIq1jjsDkAppZRqjCYopZRSMUkTlFJKqZikCUoppVRM0gSllFIqJmmCUkopFZM0QSmllIpJmqCUOkhExIjIIbF6PqVijSYopZRSMUkTlFJNJCKHisgHIlIiIitF5LzI9g9E5Jo6+40VkU8iyx9FNn8nIhUicomInCYiuSJyp4gUichGEbm8zvFNPV+GiCyMxFUsIh+LiP4fV62Wy+4AlGpNRMQNLACeAX4JnAK8ISKD93acMeZUETHA0caYdZFznQZ0ATKALOBE4G0RWWKMWXMA53sAyAUyI7udCOhYZqrV0l9XSjXNiUAS8KAxxm+MeR9YCPy6Gee82xjjM8Z8CLwFjD7A8wSArkBPY0zAGPOx0cE2VSumCUqppukG5BhjwnW2bcKqAR2IHcaYyt3O1e0Az/UwsA54V0R+EpFJB3gepWKCJiilmmYL0H23azs9gDygEkios73Lfpyvg4gk7nauLZHlJp3PGFNujJlojOkDnAvcKiKn70cMSsUkTVBKNc2XWInjDyLijlxHOhd4CVgGXCgiCZHu31fvdmw+0KeRc/5ZRDwiMgwYCfwrsr1J5xORkSJyiIgIUAaEIpNSrZImKKWawBjjB84DzgKKgDnAb4wxPwCPAn6sxPF34IXdDv8T8PdIL7ud15m2ATuwak0vAL+LnIsDOF8/YBFQAXwOzDHGfND8d62UPUSvoSplj0jt63ljTLbNoSgVk7QGpZRSKiZpglJKKRWTtIlPKaVUTNIalFJKqZgUU0MdZWRkmF69etkdhlJKqYNo6dKlRcaYzN23x1SC6tWrF0uWLLE7DKWUUgeRiGxqbLs28SmllIpJmqCUUkrFJE1QSimlYpImKKWUUjFJE5RSSqmYpAlKKaVUTNIEpZRSKiZpglJKKRWTYupG3ZZy2mmnNdg2evRorr/+eqqqqjj77LMblI8dO5axY8dSVFTExRdf3KB8/PjxXHLJJeTk5DBmzJgG5RMnTuTcc89lzZo1XHfddQ3Kp0yZwvDhw1m2bBkTJkxoUH7//fdz0kkn8dlnn3HnnXc2KJ85cyaDBg1i0aJFTJ06tUH5E088wYABA1iwYAHTp09vUP7cc8/RvXt3Xn75ZebOndugfP78+WRkZDBv3jzmzZvXoPztt98mISGBOXPm8MorrzQo/+CDDwB45JFHWLhwYb2y+Ph43nnnHQDuvfdeFi9eXK88PT2dV199FYDJkyfz+eef1yvPzs7m+eefB2DChAksW7asXnn//v158sknARg3bhxr166tVz5o0CBmzpwJwBVXXEFubm698qFDh/LAAw8AcNFFF7F9+/Z65aeffjp33303AGeddRbV1dX1ykeOHMltt90GtK/PnsFgjOGxvz5G3/59WbhgIbNnzbbKjMFgjfM5/fHpdMnqwpuvvsmLz74Ihtoyg+Ghpx8itWMqb770JgteXrDz5LX7PPj3B/HEeXj9H6/z4cIPa4/bOb//hfsxGF5/6nW+/t/X9co9Xg+TnpxEmDCvzXmNFV+sqHf+pLQkbnj0BgyG+TPns/679bXxA6R1TmPsfWOt8ofmk7s2d9f5DWT2zGTUXaPAwCv3vULh5sJ6f59u/box8taRGGN4+Z6XKS0orfe37X5Ed4aPH47B8PLkl6kqq6qND6DX4F6ccuUpAPzztn8S9AXrHd93aF9O+PUJGAwv3vxi/X84A/1/3p9jfnUMgZoAr97xaoN/28NGHMbhIw6nurSahfcsrP277XTUeUfR/+f9Kcsv490H321w/OiZo5lx2owG21tK1BOUiIwAZgFO4GljzIPRfs3WwNT9T2qgrDpAUYWP0io/wbDBISAIIjYHqmwXNmGqAlWU1JTgD/kJm7A1Yc1XFK7As9HDT5t/orC6kLAJY4yx5hheWfMKXyV/xebVm9lQuqF2+84kct8X99GhpANbV25l5faVtduNMYQJc9V/riJudRzF3xeTl59XW7bT6IWj8Xb1UvZtGUUFRQ3i/+1/f4sn3UPp0lK2F21vUH7D4htwJbvYsXwHO4p3NCi/9YNbcXgdbF+7ndKS0gbld35iJdWijUWUlZbVK3N4HEz90kqqBbkFVJRV1Ct3GRczllpfsNu2baOqvKpeeYGngLnfWT/othZtpaaypl55yY4SXlj1AiLCppJN+Kp8uwoFqkureXPdmyCQV5GHv8YfKbL+Y4fKQizevBhBKKwubJCAcspy+HLblwhCqa+UsD9crzyvIo9lBcsAqAxUNvjb5Ffm80PxD4R8IaqCVbWvW/v+qgqIL4nHV+ajOljd4Pii6iLiyuOoqqzCF/I1KN9SsaXBtpYU1dHMRcQJrAXOAHKBr4FfG2NWNbb/4MGDTWsY6igcNpTXBCmp9lNSFaCkOkBJlZ+y6kDtemlkubwmQJU/RKU/SJUvMveHCIX3/+/udAgJbicJXieJXheJHheJXidJXhcp8W5S492kxXtIS3CTlmCtpyd6yUj20DHRg9fljOJfQ+20M5FUBauoDFTWLu9tXh2spjpYTVVw13J1sJqaYE3tvCZUs+8X3wOnOPE4PbgdbjxODx6Hx1p3unE7dk0uh6ve8s6p7rpTnDgdTlziqp2DA4c4cOBExAHGgVN2LTvEATiQyH6CA+vreddy7T5GEHEgCCayh0TKapcNOMT6PFvHya7ziUS+gK1jjdQ5h5HIDz7rC1pqXzdyxM59hF0x1lu2jqv7dbmn7849/c/e01ft7rWWvZ2kqeeONhE4+ZCMFjiPLDXGDN59e7RrUEOAdcaYnyJBvAScDzSaoJqrZPtW3ht/Dv7qhr8kcHvY1rUHeRnJfJm/msI4Q231xECHbr3o2K0XQb+fjd99Zm0OW//w4bAhqUtfPOk9qKmsYseqLzAhMEEgBCYE3h7H4skYQLC0jPKv38IRFCTkwBEQHAFh8NlXc/yxwyjN+5H3/jYNt/Xpj/znEC6+7jYGHD2Ytd8v5V+PPxz5FRupaRnDyVffQWpWP1Z/8ylfvfY04bAhFDYEw4ZgOEzaGTfgTs+mat2XlH31er237nQIR10xhezsbEpWfMiPH7yG2+nA7XLgdgpup4O//eNF+vXsxisvPt9umviMGIzbcOwJx3LTxJuoDFRyx5Q7KK0uxbitMuMy9Ozfk+NOPI6qYBXvLH6HoCOIcZnafTxJHsQjjf4C3aMQSFCQgJCWmEbXjK64jIsfv/sRggIBICAQhGOPGMphhx5DZWkVb770MgSsz144AAThlOEXcOgxp5Cfs4X5c6YTDoIJQDhoMGHhlFHXkn3ECeSuW8378yKfLauFirAxHH7edaT0OoLCdd/xw4KnGnz2ss4aj7dLH0p/XEr+R89Zfzuz62s1/cwbcad3afSzB5AxciKulEwqV39E+bdvNyjP/NVknAmpVCxfRMXyRQ3KO436Ew53HOXfvEXlDx83KO9ymdUoU/rla1Sv/6pembi8dB79ZwBKPv0nNZu+q1fujE8h8wKrBrbjw3n48n6oV+5KziDjXKv5tnjRk/gLfqpX7u6YRfqImwDY/p+/ECjOq1fu6dSHjsPHAVC04BGC5fVrmN6sgXT42VgACl+/n1B1/RpgXM+jSTv51wDkv3IPJli/FhPfdwipJ1wIwLYXJ7G7xIHDSD72HMKBGgr+9acG5UlHDifpyOGEqkop/PcDDcqTjzmbxENPJVhWSNHChpcOeox5iLX3ndVge0uJdoLKAnLqrOcCJ9TdQUTGAeMAevTo0awX89dU0X1dNY31/XCGgxy7yvpwXQtUuyGnI2xKFxYf7mB9Wi7VgQIMBveAcG3yQKzTiWsjQedmnIkh0js39urfWFM2JB8O1n//UG3pBh5nm+/veDO8JF3tw+F3IH5BfIL4hWDW15BSRYeehaQf4bC2+wSHz9rvdz/ry6BBR7EorYCpnyU3ePXZt5xK1x59ePXflTyR+z7BUJhA2BAIhQmEwvTvnITPIWwrraGo0k8wVL+pYORfPsGZkIpv1XKqNpfgdgoupwOXQ3A5hdnv/0inDil8l1PCjio/TofgFMERmZfXBIhzO/f4q7K5jDHUhGrwOX0Ek0OEXWFwGcIuQ35aIc+vmE9FoJINXTZTllCJcYcxLjBuw1cdlnLhvy+jOljF1qGbCTmC4DG1n/7FLGbxG5GkeVLD115tVrPuh424HfFUJfsI+wz4BVNhJQJvSgZdex2FuLwsf/dNjF8I+w3GJ4T90Knv8WQNOhNfjWHpU1MJ+Q0mJNaPH2NIO/p4Co4cTnX5Dra99l2D189Z34EP1/awviTeb5gE8z3VJGwpI7C9hO2bd5WLCA6Bz9ZvJ91ZgD+/lEpf0Go+Fqn9ceRyCqnxbkIJXuLcjjplVv3i5EPSyerbjY2Sw8ffxVnnZtc5fn1Kb7r27MuKz/NYtD6x9v/NzsakG84aSGaXbnz+3kbey01q0Gw9+cKjSO3YkcVv/Mjigoaf7Xt/fQzx8QkscC/no5KG5Y+OHYwgvBz+mi8q6ycYb1wc0357PADPVX/MUv+G2vgBUtI6cO9VQwB4cse7rAzXb7LK7JLO3VcPQRAe2/YmP66un2B69O7EHddYX2nTNr5MzobyeuWHHNaFCdda5X9enUHhNn+98iOOzmL8tSciApOXdaS0pP57O35wd64adyIAEz5Pw+errtdMd/LQnlxx7VAAxn+Q2uBvM3xYby4eM5Sa6iomfNqwfOTP+zLy4qGUFG9n0lcNyy86oz9njBxK/pY87vmmYfkT157QYFtLinYT3yjgTGPMNZH1McAQY8xNje0f1Sa+Dx8m9J/78A19CF9VKr516/CtX0fN8hWEKypIGDKE9GuvJfGUk2ur83tijCFoggRCAQJha/KFfLXNMnWbaSoDlVQEKqzJb83L/eVU+Cso9ZdS5iujzF9Gub+88ap+RLwrnhRPCqneVFI8KdbkTSHZk0ySO4kkd5K17LGW413xJLgTiHfFW8uuBOJccZFmDfAHw2yv9FFYbk3bK/1sr/CzvcJHcaWfoko/xZW+Ok2VwT3GFvmrgIRAAnjcITyuEG5XCJcriDgCIH5wBDDijyz7Mfgx4rMm/BiHr942HD5rf4cPxI/I/n1WjXFA2IMJezFhL0TmJuyp3V5/m3fXviEPxngxoV37YKymrJ1EwO1w4HIKTodVA3U6BLcjktSdUlvucjpwOwSPy2HVWp0OPC6ps+zAE5m7nYLH6cTtEjxOB15XpNzlwON01i7Xbm+wT/3lfX2OlYoVdjXx5QLd66xnA9G9qrYnJ92Ec9kLJGx6koTxn4HLA0CoopKSf/2L4nnzyLn2WryHHkr6NVeTcuaZiKvxP4+I4Barvb6lhE2Ycn85Zf4yynxllPpKKfWXWnNfKWX+stp5mb+MzeWbKdteRoW/gqpg1b5fIMLlcO26DhG5LuF2uK1f2zvb8eMFR7yDuAzwECY9bF2MD4RDBCNTIBwgEPITNEFCJkDINExgwci0Nw7cOPHiEi8uvDglDrfE4ZQ03BKH2xFXO/c64vE44/E64vE6rSnOEU+cK5KInYkkuBLwOD21icLlsGqBTqdY8zoJZee6y+HA6SAyl9rE4xQrwdTdd2etUSkVfdGuQbmwOkmcDuRhdZK4zBizsrH9o95JYu278OIoGP5nOGVCvSLj91O6YCHbn34a/4YNuHv2oNv995Nw3HHRi6eFhMKherW0cn95vZpc3ckf8hMIBxrMDaa291ftMganOK2L4OLAgQOHw7oIXnuB3enG4/DgcrjwOr14nV7iXHHWsstLnDOOOFcccc642tpcvCueOJe17nK0yTsdlFJNsKcaVFQTVOSFzwZmYnUzf8YYc9+e9j0ovfhevBQ2fAQ3fg2pWQ2KTThM+eLFFDz8CIHcXNKvG0fm9dcj7parLSmllNplTwkq6iNJGGPeNsb0N8b03VtyOmhGPADhILw7pdFicThIOeMMer/2Gqm/+hXb5z7OxsuvwL9x48GNUyml2rn2N9RRx95wyi2w8jWrJrUHzqREut1/H1kzZ+LftImfLryIkvnzo9ZLTSmlVH3tL0GBdf0prSe8fTuEAnvdNWXEmfR549/EH3kkW6fcTd6EWwhX7X+nBKWUUgemfSYodzyMeBAKf4Avn9j37l260OPZZ+h020TK33uPzb+9iuCOhkOyKKWUajntM0EBDDgL+v0SPngQSvP2ubs4HKRfcw1Zs2ZSs3o1my6/gsAWe3rMK6VUe9Am+/ju94jSwRrIK4R/HMXYPzzE2Kuu3ueI0qUDB/JIUiLjN2wgf/gZzO7QgS1ul45m3gqHOtpJRzPXzx60o8+eMVjDuRvGXn4JYy+/hKLCfC4ec23t9p3z8b8ZxSXnDicnbwtjfv/HXYP+Rfb54Pnp0G94g9dtKW0yQe03Vxyk94WitbDmHeDq/TpsncfD9I4duGlHCROLi5nTIS2qYSqlWhET3jXVlEHhGghUQ2UR1JREBjmMlOevhK//BiE/bF8HJRXWdiIDga7xwxvbrGvl25aDLxApjySSL7fAkwshFIS873cdtzPJ/Pc7KLkPqv2wsYQGw82+swS23gVVYchtZCzJRSth231QGoatjZS/9Gu4u7Dh9hYS9fugmsK20cxfHw/f/ROufBN6n7rfh/lz88i55hoCW7eSNWM6yaefHsUglVItJhQAX7k1+SvAVwH+cvBX1pkqIvMqCOycV1nbAtXWcqDamoLVu5b3MmTZ/hFwecHptUa8cXrA6bbmDre1zRFZd7p2bXe6weHabb5zu9NadrisY8S5az+Hq065c9e67L7s3FVedz2r+YMZ2HajblPYlqB8FfDkadaH9XefQFLmfh8a3LGDnOt+R82qVWT/5TGSf/7z6MWplLKEw1ZtpHoHVBVHlkvqz2tKoKbUqsX4ynbNfeVW8/7+EAe4E8GTAO7I5EmwOlq568xdcZHl+F3LrrjIcmTu8u7a5vTU2eaNrEeSksNJgxF12zhNUPuybQU89QurBnXZK+DY//4joYoKNo/9Lb61a+n+5BMknnhiFANVqg0Kh6BqO1TkQ0WB1RxWWQhVRdZy1fZd8+piKwntrabiToC4VGvypkBcSp15MnhTwZsEnqTIehJ4ds4Tre2eRCuJtLNkYQdNUPvj66fhrYlwxr1w8s1NOjS4Ywebf3Ml/rw8evztaRKOOSZKQSrVioTDUFkAZXlQvg3Kt+4232YlpKqiyLWV3ThckJABiRmQ0NFaTugI8R0gvmNkeed62q6k5PIe9LeqDpwmqP1hDLzyG1jzNlz1X8hu8Pfaq0BBAZvGjCFUvIOe//g7cYceGqVAlYoRgRoozYGSTVCy2ZpKc61bN8pyoWxr5MmKdYgDkjpDchdI6gJJnaz1pE7WlBiZJ6RbyUZrMG2eJqj9VV0CTwyzlq/9ABLTm3R4IC+PjVeMwfh89Hz+Obx9+rR4iEodNMZYzWrFP0HxBmu+Y4O1XLLJapKry+GGlG6QkmUNxpySBanZ1jylKyR3hcRM6zqLUhGaoJoidynMOxs6Hw5XLrDaopvAt2EDm8b8BnE66fnC83iys6MUqFItJFBjdXPe/iMU7Zz/aG3z1X0MuVgJp0Mva0rrCWk9rKlDT6smpMlHNZEmqKZavRBeGQOHDIdLX7S6ZDZBzZo1bPrNlTjTUun10ku4OnSIUqBKNUHQb933V/gDFKzeNd+xof41oJRsyDgE0vtZ9wp27AMdeltJSK/vqBamCepALHkGFt4Cgy6H8//a5Lbwqm++ZfPYscQdeSQ9nn0Gh8cTpUCVakRVsXVz584pf4WVkMKR5xyL00o8nQZC5qGQOQAy+kH6IU1uNVCqOex65HvrNvgqKM+HDx+0Luie/scmHZ5w7DF0feB+tky8ja1TptBt2jREL/iqaKjeAVuWwZZvYcs31nJpzq7y5K7Q+Qhr/MnOh0PmQCsZaW1IxTBNUPty2iSo2AYfT7d6HJ0wrkmHp55zDoGcHApnzsLToyeZN94QpUBVuxEOWUPk5HwJOV9B3hKr88JOHXpD9vEw5FrociR0PrJJN58rFSs0Qe2LCJw9HSoK4Z0/WN1fD/9Vk06Rft11+Dduomj2bDw9upN63nnRiVW1Tf4qKxlt+sya5y21huEBq1NC9vFwzBXQ7RjoOsi6N0ipNkAT1P5wuuDiv8E/zodXr7Ha8I9sOOr0nogIXf/vzwS2bGHrXVNwd+tGwuCm3WOl2pFADeR+DRs/hg0fWzWkkN+6f6jzEXD0r6H7CdB9iNV7TpuNVRulnSSaoroEXroMNn0KZz4AQ69v0uGh0lI2XvprQsXF9Hr5JTy9ekUlTNXKGAPb18O69+DHd2HjpxDyWQmp69HQa5g1BFePE61heZRqY7QXX0sJ1MBr18LqN+Gkm2H4n5s0bp9/82Y2XnIpzrQ0er3yMs5k/cJpl4I+2PCRlZB+fBd2bLS2p/ezbm3o8zPoeZI1koJSbZwmqJYUDlnXo75+Go66xOqC3oT7pKq+/ppNv72KpGHDyP7rbKQJCU61YoFqWLcYVr0Ba/9j3QDrirdqR/3OsKYOveyOUqmDTruZtySHE85+xOp6/v5Ua5Tl0f+wRkLeDwnHH0/nyZPIv3cqRbNnk3lz0wamVa1I0Gclo5X/hrX/tZ4rFN8BDjvfmnoNsx7HoJRqQBPUgRKBU2+3up4v+D08OwIuegYy++/X4R0uu4yaVasomjMX78CBpPzyl1EOWB00xsDWZfDtC7D8X9ZziRIy4KjRkaR0SpNHJlGqPdIE1VzHjrFqUq+NgydOhTOnwuCr99mzSkTo8sc/4lu3ji2TJuPp1Yu4/vuX3FSMqiyC71+2ElPBSuvhc4eOtEYi6XOajlGnVBPpNaiWUr4N3rgB1i2y7tY/bzYkd97nYYH8fDZcfDGO+AR6/+sVnKl6UbzV2bIMvnwcVrxqdQfPOs5KSkdcaDXnKaX2SjtJHAzGWB0n3p1ijWV23l9g4Dn7PKzqm2/ZdOWVJJ5wAt2feBxx6i/tmBcKwpq34Iu5sPlz6wmsgy6zhsfqpM8BU6op9pSgtPtYSxKxhpe57iPr+TcvXQavXWc9xG0vEo49hi53T6Hyk08onDnrIAWrDoivHD77Czw2yHq4ZdkWOPN+uHUVnP2wJielWpBeg4qGzAFwzWL4cJr1ZbbiVThuLJx6m3W9qhEdRo+mZsVKtj/1FPGDjib59NMPbsxq76qK4csnrKa8mhKr991Z06D/CL22pFSUaBNftJXmwUcPwbfPW08bHXItnHJLo+OlhX0+Nl1+Bf6NG+k9/1860kQsKN8Gn8+GJc9a498NHAmn3ArZx9kdmVJthl6DslvxT/DBg/D9K9b1isFjrZt8Ox9Rr8dfIC+PDRdehKtzZ3q99E8cCQn2xdye7dgEn86yfliEA3DERVZi6nyY3ZEp1eZogooV+aus50v98JY16Gynw+DIUdaU1h2Aik8+Jefaa0kZOZJuD+kzpA6qwrXwyQzrh4Q4rI4Pp0ywHuynlIqKg56gRORPwLVAYWTTncaYt/d2TLtIUDtVboeVr1k3cuZ8aW3rebLVRT37eAoXLKHor4/T+e4pdLz8cntjbQ+2LLOe+bV6AbjiYPBvYeiNkJpld2RKtXl2JagKY8wj+3tMu0pQdRVvgOXzYcV865HcgMFJ7pfZVGwO0utPY4k/8TTrqajJXVp+FAJjIFBl9VDzlVtjxPnKwV9pjR8XqIZgza55KAAmDBhrbsLWORxO6+ZUV2Ryeqwve2+SNehpXCrEpUWmVOsxJnYK1MAPC+Gbv1sDt3pTrGuEJ14PiRn2xqZUO6IJqrWoLLIeSJf7NaEfv2DDUz9iQobeZxbiigsDYj2kLqUrJHcDT0IkKXjqzD3WDaNBvzUP+azlYI11od9XDr6KXQnJXx5JOE0gjl0TYl1HC4es6zX7Ky7NegBkYifria+Jnaz1lG5WN/3UbGvZHd+02PYlfxV88w/4/iXrUempPaxrgsdfo6OHK2UDuxLUWKAMWAJMNMbsaGS/ccA4gB49ehy3adOmqMTTWtWsWMHGyy4n/tBe9LjtfKRiG5TlWffflG+zaj6hgDUoachvzcMBq8dgbS2mztybbE2eJKvG4E2yluNS6mzbuU+ilRxccdbcHW+Nvr23mk84HEmIO+OpsZJhTanVPbum1HquVk2JlYwrC6AiMlUWWrW33SWkWwkrrYc12nfdKbX73gdbNQZKc2DbCshfYQ3YmrfE+vscOhKO/Q30Pq1Jj0xRSrWsqCQoEVkENHZjz13AF0ARYIB7ga7GmKv2dj6tQTWu5NVX2XrXFDKuv57Mm2+yO5zoClRbybcsz+qiX5ZrzUtzrRueSzZZSa8ub6rVJJeYGamJZVo1u/xVkL8SfKW79u10OBxzudWDUpvxlIoJUXnchjFm+H6++FPAwua8VnuWdtFFVC1ZStHcucQfM4ikYcPsDil63PGQ3teaGmMMVORb3cB3bITSzZGaWKFVCytaB5s+s2qVnQ6FIy+yuvJ3OdJa1yfSKtVqRO0qtYh0NcZsjaxeAKyI1mu1B13+eDc1K1ey5fY/0Pu1V3F362Z3SPYQsTqKJHeBHifYHY1SKoqi2fD+kIgsF5HvgZ8Dt0Txtdo8R3w8WbNmYgIBcm+5BeP32x2SUkpFVdQSlDFmjDHmSGPMUcaY8+rUptQB8vbuTdf77qPmu+/Jf3i/O0cqpVSrpF2XWpmUEWfS8crfsOO55yh75x27w1FKqajRBNUKdZo4kfhBg9h61xR8P22wOxyllIoKTVCtkHg8ZM18FPF6yfv97wlXV9sdklJKtThNUK2Uu0sXuj38ML5169j2pz8TS4P+KqVUS9AE1YolnXIyGTfcQOkbb1Ayf77d4SilVIvSBNXKZYz/HYknn0z+vVOpWbXK7nCUUqrFaIJq5cTppNvDD+Hs2JHc308gVNbIWHZKKdUKaYJqA1wdO5L16AwCW7eyZfKdej1KKdUmaIJqIxKOOYbOf7idisWLKX7mWbvDUUqpZtME1YZ0GDOG5DPPpGDGDKq+/trucJRSqlk0QbUhIkLX+6bi6d6d3FtvJVBQYHdISil1wDRBtTHOpCSyHptFuKKSvFtvxQSa8IRbpZSKIZqg2qC4/v3p+n9/pnrJUgoenWl3OEopdUA0QbVRqeeeS4fLLqP4mWcoe/ddu8NRSqkm0wTVhnWadAdxRx/F1sl34tugg8oqpVoXTVBtmMPjIXvmTMTtJu/m3xOuqrI7JKWU2m+aoNo4d9eudJv+CL5169h6z5/0Jl6lVKuhCaodSDr5ZDJvvomyBQvY8cKLdoejlFL7RRNUO5F+3XUk/fzn5D/4IFXffGN3OEoptU+aoNoJcTjoNu1B3FndyPv9BL2JVykV8zRBtSPOlBSyH/sLoYoK8m7Rm3iVUrFNE1Q7EzegP12n3kv10qXkP/Sw3eEopdQeuewOQB18qeecQ833yyn++9+JP+ooUs8daXdISinVgNag2qlOt00kYfBgtt59NzVr1tgdjlJKNaAJqp0St5usmY/iTEkh98abCJWW2h2SUkrVowmqHXNlZJD92CwC27aRN/E2TChkd0hKKVVLE1Q7Fz9oEF3unkLlJ59QOHOW3eEopVQt7SSh6DB6NDWrVrH9qaeIO+xQUs46y+6QlFJKa1DK0uXOO4k/9li23HmXdppQSsUETVAKAPF4yJ41E2dyMrk33EiopMTukJRS7ZwmKFXLlZlJ9l8eI5ifr50mlFK20wSl6ok/+mi63PNHKj/9lIIZM+wORynVjjUrQYnIKBFZKSJhERm8W9lkEVknImtE5MzmhakOprSLL6bDZb+m+G/PULpgod3hKKXaqebWoFYAFwIf1d0oIocBlwKHAyOAOSLibOZrqYOo86RJ1kgTU6ZQvXyF3eEopdqhZiUoY8xqY0xjXb7OB14yxviMMRuAdcCQ5ryWOrjE4yHrsVm40tPJvfFGgoWFdoeklGpnonUNKgvIqbOeG9nWgIiME5ElIrKkUL8EY4qrY0ey5/yVUFkZuTfdTNjvtzskpVQ7ss8EJSKLRGRFI9P5ezuskW2msR2NMU8aYwYbYwZnZmbub9zqIIkbOJBuDzxA9bJlbPvznzGm0X9GpZRqcfscScIYM/wAzpsLdK+zng1sOYDzqBiQMuJMfNePp2jOXOIGDKTjb8bYHZJSqh2IVhPfm8ClIuIVkd5AP+CrKL2WOggybryRpNNPJ3/aNCo//9zucJRS7UBzu5lfICK5wFDgLRH5L4AxZiXwCrAK+A9wgzFG7/psxcThoNu0aXj79CZ3wi34N22yOySlVBsnsXRNYfDgwWbJkiV2h6H2wp+Tw8ZRo3Gmp9PrpX/iTE62OySlVCsnIkuNMYN3364jSagm8XTvTtasWfg3bSJv4kQdDkkpFTWaoFSTJZ4whC5TplD50ccUPDLd7nCUUm2UPg9KHZAOl16C78cfKX72Wbz9+pF24QV2h6SUamO0BqUOWOfJk0g8aSjb7rmHqm++sTscpVQbowlKHTBxuciaMQNXt67k3nQzgbw8u0NSSrUhmqBUszjT0ug+dy7G7yfn+hsIV1baHZJSqo3QBKWazdunD1kzZuD78Ufy7rgDEw7bHZJSqg3QBKVaRNKwU+g86Q4qFi2mcOYsu8NRSrUB2otPtZgOY8bgW7ee7U8+ifeQvqSed57dISmlWjGtQakWIyJ0uXsKCUOGsPWuKVR9+63dISmlWjFNUKpFidtN1qyZuLp2JffGmwhs0UHslVIHRhOUanGuDh3oPncOxufTnn1KqQOmCUpFhbdvX7IenYFv7Vrt2aeUOiCaoFTUJA0bRudJk7Rnn1LqgGgvPhVVHcZcgW99pGdf3z6knn++3SEppVoJrUGpqBIRuky5i4QTTmDrlLup+kZ79iml9o8mKBV14naTNfNRa8y+G2/UMfuUUvtFE5Q6KKyefXMxgQA5468nVKE9+5RSe6cJSh003j59yJr5KL7169ly2236NF6l1F5pglIHVdLJJ9P5zslUfPABBTNm2B2OUiqGaS8+ddB1vPxyfOvWUfy3Z/D2PUSfxquUapTWoJQtutx5JwlDT2TrPfdQtXSp3eEopWKQJihlC3G7yZ45E0+3buTeeBP+XO3Zp5SqTxOUso0zNZXsuXMxoRC548drzz6lVD2aoJStvH16kz3zUXw//cSW22/Xnn1KqVqaoJTtEk86yerZ97//Ufjoo3aHo5SKEdqLT8WEjpdfjn/9erY//Tc8fQ8h7YJf2R2SUspmWoNSMaPz5MkkDD2RbX/8o47Zp5TSBKVih7jdZD8aGbPvJn0ar1LtnSYoFVOcaWnWmH1+vz6NV6l2ThOUijnePn3ImjEd39q1bJk0WZ/Gq1Q7pQlKxaSkYcPo9IfbKX/vPYpmz7Y7HKWUDZqVoERklIisFJGwiAyus72XiFSLyLLI9HjzQ1XtTccrryT1ogspmjOXsrfftjscpdRB1txu5iuAC4EnGilbb4wZ1Mzzq3ZMROhyzz34N25iy+Q7cXfvQfyRR9gdllLqIGlWDcoYs9oYs6alglFqdw6Ph+zHZuFKT7eexltQYHdISqmDJJrXoHqLyLci8qGIDNvTTiIyTkSWiMiSwsLCKIajWitXejrZc/5KqKyM3JtuIuzz2R2SUuog2GeCEpFFIrKiken8vRy2FehhjDkGuBV4UURSGtvRGPOkMWawMWZwZmbmgb0L1ebFDRxItwcfpOa779l2z58wxtgdklIqyvZ5DcoYM7ypJzXG+ABfZHmpiKwH+gNLmhyhUhEpZ/4S3403UjR7Nt6BA0gfO9bukJRSURSVJj4RyRQRZ2S5D9AP+Ckar6Xal4zrx5P8y19S8NDDVHz8id3hKKWiqLndzC8QkVxgKPCWiPw3UnQq8L2IfAfMB35njCluXqhKgTgcdHvgfrz9+pF36634NmywOySlVJRILLXlDx482CxZoq2Aat/8uXlsHDUKZ1oavV5+CWdKo5c4lVKtgIgsNcYM3n27jiShWiVPdhbZj83Cn5NDnj7oUKk2SROUarUSjj+eLlOmUPnhRxTOnGl3OEqpFqYPLFStWodLL6Hmh9Vsf+ppvAMGkjryHLtDUkq1EK1BqVavy513Ej/4OLbedRfVK1baHY5SqoVoglKtnng8ZM+ahTO9I7k33kiwqMjukJRSLUATlGoTXOnpdJ89m1BJCbk3/x7j99sdklKqmTRBqTYj7rDD6Hb/fVR/8w3b7p2qwyEp1cppJwnVpqScfTY1a9ay/Ykn8A4YQMcrLrc7JKXUAdIalGpzMn9/M0mnn07+Aw9Q8emndoejlDpAmqBUmyMOB92mTcPbty95t+hwSEq1VpqgVJvkTEoke84cxOkkd/z1hEpL7Q5JKdVEmqBUm+XJziL7L4/hz8sj75ZbMcGg3SEppZpAE5Rq0xIGD6brn+6h8rPPyJ/2kN3hKKWaQHvxqTYv7aKL8P24juJ58/AecggdLhltd0hKqf2gNSjVLnS6/TYSTx3GtnvvpeIT7dmnVGugCUq1C+J0kjVjBt5DDiHv5pupWb3a7pCUUvugCUq1G86kJLo/8QSO1FRyxl1HYMsWu0NSSu2FJijVrrg7d6LHk08Qrqlh87hx2v1cqRimCUq1O95+/ciePZvAps3k3ngTYR1YVqmYpAlKtUuJJwyh6wMPUPX112ydNBkTDtsdklJqN9rNXLVbqSPPIbB1C4XTZ+Dq0oVOt9+GiNgdllIqQhOUatfSr7mG4NZtFD/zDI7EBDJvuMHukJRSEZqgVLsmInSechfhqiqK/jIbR1wc6VdfbXdYSik0QSmFOBx0vW8qxu+j4OFHEG+cPkdKqRigCUoprBt5u02bRtjnJ3/qVMTrocOoUXaHpVS7pr34lIoQt5usR2eQOGwY2/54D6Vvvml3SEq1a5qglKrD4fGQ/ZfHSBgyhC2TJlP2n//YHZJS7ZYmKKV244iLo/ucvxI/aBB5t06k5NXX7A5JqXZJE5RSjXAkJtLj6adIPPFEtt51F9vnzbM7JKXaHU1QSu2BIyGB7MfnkvzLX1Lw4DQKZs3CGGN3WEq1G5qglNoLh8dD1ozppF50IdvnPk7+vVN1WCSlDhLtZq7UPojLRdepU3GmpFL87LOEysvpdv99iNttd2hKtWnNqkGJyMMi8oOIfC8ir4tIWp2yySKyTkTWiMiZzY5UKRuJCJ3+cDuZEyZQtmABOeOvJ1RebndYSrVpzW3iew84whhzFLAWmAwgIocBlwKHAyOAOSLibOZrKWUrESHjd9fR5d7/o/KLL9h4yaX4N22yOyyl2qxmJShjzLvGmGBk9QsgO7J8PvCSMcZnjNkArAOGNOe1lIoVHUaNosff/kZo+3Y2jL6Eyi++tDskpdqkluwkcRXwTmQ5C8ipU5Yb2daAiIwTkSUisqSwsLAFw1EqehJPGEKv+f/ClZnB5muuYcdLL9sdklJtzj4TlIgsEpEVjUzn19nnLiAIvLBzUyOnarR/rjHmSWPMYGPM4MzMzAN5D0rZwtO9O71eeonEk09i25/+xLap92GCwX0fqJTaL/vsxWeMGb63chG5EhgJnG523SSSC3Svs1s2sOVAg1QqVjmTkug+Zw4Fj0yn+Nln8f3wA92mP4K7c2e7Q1Oq1WtuL74RwB3AecaYqjpFbwKXiohXRHoD/YCvmvNaSsUqcTrpfMcf6PbQNKpXrWLD+b+i4sMP7Q5LqVavudegZgPJwHsiskxEHgcwxqwEXgFWAf8BbjDGhJr5WkrFtNTzzqP3/Pm4Oncm57rfkf/Qw5hAwO6wlGq1JJaGbhk8eLBZsmSJ3WEo1Szhmhryp02j5J8vEXf0UWRNn4Enu9E+QkopQESWGmMG775dhzpSqoU54uLoes89ZM18FP/6n9hwwQWULlig4/gp1USaoJSKkpQRI+j9+mt4+/Rhy+1/IPfGmwgUFNgdllKthiYopaLI0707PV98gU63307lJ5/w07nnUfrGG1qbUmo/aIJSKsrE6ST96qvo/frrVm3qjknkjr+eQL7WppTaG01QSh0k3j696fn8c3SadAeVn3/OTyNHUvzii5iQdnBVqjGaoJQ6iMTpJH3sWPq88W/iDjuM/P+7lw0Xj6Lqm2/tDk2pmKMJSikbeHr1ose8Z8maMZ1QcTGbLruMLZMmEywqsjs0pWKGJiilbCIipJx9Nn3ffov0a6+l9K23WD/iLIr/8Q+M3293eErZThOUUjZzJCbSaeKt9HnjDeKPPpr8+x9g/chzKX3rLX28vGrXNEEpFSO8fXrT/emn6P7E4zji49ky8TY2XHwxFZ98qt3SVbukCUqpGCIiJP3sZ/R+/TW6PTSNcGkZOddcw+bfXkX18uV2h6fUQaUJSqkYJA4HqeedR5933qbznXfiW7OGjaNGs3ncOKq++cbu8JQ6KDRBKRXDHB4PHX8zhr7vvUvmLbdQs3wFmy67nE2/uZLKzz/Xpj/VpmmCUqoVcCYlkXHdOA5ZvIjOkyfh37iRzb+9io2XXkr5++9rZwrVJunjNpRqhcJ+P6Wvvc72p54ikJeHu2cPOl5+BakXXoAzKcnu8JRqkj09bkMTlFKtmAkEKF+0iOK//4PqZctwJCaSetGFdLziCjw9etgdnlL7RROUUm1c9fLlFP/jOcreeQdCIZJOPZW00aNIOvVUxO22Ozyl9kgTlFLtRCC/gJKXX6LkX/MJFhbizMwg7VcXkHbxRXh69rQ7PKUa0ASlVDtjgkEqPvqIkn/Np+LDDyEcJmHIEFIvvIDk4WfgTEq0O0SlAE1QSrVrgfx8Sl//NyWvvkogJwfxekn6xc9JHTmSxGHDcHg8doeo2jFNUEopjDFUf7uMsoULKXvnHUI7duBISSHlzF+SPGIEiccfj2iyUgeZJiilVD0mEKDyiy8oW7iQ8vcWEa6qwpGcTNJpp5E8fDhJw07BkZBgd5iqHdAEpZTao3BNDZWffUb5e4uoeP99QqWliNdL4sknk3Taz0gaNgx31652h6naqD0lKJcdwSilYosjLo7kX/yC5F/8AhMMUrVkKeWLF1O+2EpYAN5+h5A47FSSTh1G/LHH6nUrFXVag1JK7ZExBv+6dVR8/AkVH39E1ZKlEAggCQkkDD6OxBNOIGHICcQddijidNodrmqltIlPKdVs4cpKKr/8ispPPqbyiy/x//QTAI7kZBIGDybhhCEkHHcccQMH6s3Bar9pE59SqtkciYkk/+LnJP/i5wAECgqo+uprqr78ksqvvqTif/8DQOLiiD/ySOKPOYb4YwYRP2gQrg4d7AxdtUJag1JKtZhAfj7V335L1TffUP3tMmpWr4ZgEAB3jx7EH3E4cYcfQdwRRxB3+GE6sK0CtIlPKWWDcHU1NStWULVsGTUrVlKzfDmBLVusQhE8vXrhHTiAuAED8Q7oT9zAgbi6dEFE7A1cHVTaxKeUOugc8fEkHH88CccfX7stWFxMzcqVVC9fTs3KVdQsX0H5O//ZdUxqKnH9+uHpdwjePn3xHtIXT9++uDIzNXG1M5qglFIHlatjR5KGDSNp2LDabaHycnxr11KzZg2+H9bgW7uWsoVvES4vr93HkZKCt3dvPL164enV05r37ImnZ08ciTquYFukTXxKqZhkjCFYWIh//Xp863/Ct34d/g0b8W/cSHDbtnr7OjMy8GRn487Oxt09O7LcHXdWFu7OnbRHYYyLShOfiDwMnAv4gfXAb40xJSLSC1gNrIns+oUx5nfNeS2lVPsiIrg7dcLdqROJQ4fWKwtXV+PfvBn/xk34N27En7OZQG4e1d9+W/s8rDonwpWZibtrV1zduuLu2g13l864OnXG1akT7s6drOZDvfE45jS3ie89YLIxJigi04DJwB2RsvXGmEHNPL9SSjXgiI8nbsAA4gYMaFBmAgEC27YRyMkhsHUrgS1brfnWLfhWraZi8fsYv7/Bcc70dFyZmbgyMqwp05o70zNwZaTj7NARV3pHnGlpiEuvjhwMzforG2PerbP6BXBx88JRSqnmEbcbT/fueLp3b7TcGEOopIRgfj7BggIC+fkE8wus9aIigkVF+NavJ1hUBIFAIy8gOFNTcXbsiLNDB5xpaTjTUnHVLqfhSE3FmZKKMzXF2jclBUlI0E4eTdSSPwOuAl6us95bRL4FyoApxpiPW/C1lFLqgIgIrg4drBuHBw7c437GGMKlpQQLCwkW7yC0o5jg9u2EincQLN5OaHsxoZISAjk51CxfTmjHDkxjCW0nlwtncjKO5GRrnpKMM8ladyQl4kxKwpGYhCMpCUdSIo7ERJyJiUhCQv15fDzicEThLxN79pmgRGQR0KWRoruMMW9E9rkLCAIvRMq2Aj2MMdtF5Djg3yJyuDGmrJHzjwPGAfTo0ePA3oVSSrUwEamtEXn3Y39jDKaqilBJCaGyMkKlZYRKSwmVlRIuKyNUUkqoopxweQWh8jLC5RX4izYQKisnXFlJuLIS9rPTmsTH49g5JcQjCQk44uJxxMUhCfHWcnwc4o1D4rw4ds7jrG2OOK9V5vXg8HqRyOTweBCPx1rfuex225YQm92LT0SuBH4HnG6MqdrDPh8Atxlj9tpFT3vxKaXaKxMOY6qrCVVUEq6IJK2qqvrzykrCVdWEq6sJV1Viqqt3rddUY6prGiybmprmB+d243C7rWS1M2l5PEh8PH1ef63Zp49WL74RWJ0iflY3OYlIJlBsjAmJSB+gH/BTc16rKU477bQG20aPHs31119PVVUVZ599doPysWPHMnbsWIqKirj44oaX0saPH88ll1xCTk4OY8aMaVA+ceJEzj33XNasWcN1113XoHzKlCkMHz6cZcuWMWHChAbl999/PyeddBKfffYZd955Z4PymTNnMmjQIBYtWsTUqVMblD/xxBMMGDCABQsWMH369Ablzz33HN27d+fll19m7ty5Dcrnz59PRkYG8+bNY968eQ3K3377bRISEpgzZw6vvPJKg/IPPvgAgEceeYSFCxfWK4uPj+edd94B4N5772Xx4sX1ytPT03n11VcBmDx5Mp9//nm98uzsbJ5//nkAJkyYwLJly+qV9+/fnyeffBKAcePGsXbt2nrlgwYNYubMmQBcccUV5Obm1isfOnQoDzzwAAAXXXQR27dvr1d++umnc/fddwNw1llnUV1dXa985MiR3HbbbYB+9vSzd+CfvYtHjWrZz57Xw+gxVzB+/HgqS0o459xzIWzAhDHhMITDXH7OOVxx5pkUFhQw5u4/Yoy1HWMw4TBjhw3jV0cdTW5BPje88IJVw9t5DmOYP2o00dTca1CzAS/wXuTi387u5KcC/yciQSAE/M4YU9zM11JKKdVEImI149Xpebizq4a3Vy8Shw6luqgIZ8eGg/mmDB9O5iWXUJOTg+fjht0IsmY0/EHSkvRGXaWUUrbaUxNf++gKopRSqtXRBKWUUiomaYJSSikVkzRBKaWUikmaoJRSSsUkTVBKKaVikiYopZRSMUkTlFJKqZikCUoppVRMiqmRJESkENjUAqfKAIpa4Dytgb7Xtqe9vE/Q99pWNfW99jTGZO6+MaYSVEsRkSWNDZvRFul7bXvay/sEfa9tVUu9V23iU0opFZM0QSmllIpJbTVBPWl3AAeRvte2p728T9D32la1yHttk9eglFJKtX5ttQallFKqldMEpZRSKia1qQQlIiNEZI2IrBORSXbHEy0i0l1E/iciq0VkpYj83u6Yok1EnCLyrYgstDuWaBKRNBGZLyI/RP59h9odU7SIyC2Rz+8KEfmniMTZHVNLEZFnRKRARFbU2dZRRN4TkR8j84bPWG+F9vBeH458hr8XkddFJO1Azt1mEpSIOIG/AmcBhwG/FpHD7I0qaoLARGPMocCJwA1t+L3u9Htgtd1BHASzgP8YYwYCR9NG37OIZAE3A4ONMUcATuBSe6NqUfOAEbttmwQsNsb0AxZH1tuCeTR8r+8BRxhjjgLWApMP5MRtJkEBQ4B1xpifjDF+4CXgfJtjigpjzFZjzDeR5XKsL7Ese6OKHhHJBs4BnrY7lmgSkRTgVOBvAMYYvzGmxNagossFxIuIC0gAttgcT4sxxnwEFO+2+Xzg75HlvwO/OpgxRUtj79UY864xJhhZ/QLIPpBzt6UElQXk1FnPpQ1/ae8kIr2AY4AvbQ4lmmYCfwDCNscRbX2AQuDZSHPm0yKSaHdQ0WCMyQMeATYDW4FSY8y79kYVdZ2NMVvB+pEJdLI5noPlKuCdAzmwLSUoaWRbm+5DLyJJwKvABGNMmd3xRIOIjAQKjDFL7Y7lIHABxwJzjTHHAJW0nWageiLXX84HegPdgEQRucLeqFRLE5G7sC5JvHAgx7elBJULdK+znk0bajLYnYi4sZLTC8aY1+yOJ4pOBs4TkY1Yzba/EJHn7Q0panKBXGPMztrwfKyE1RYNBzYYYwqNMQHgNeAkm2OKtnwR6QoQmRfYHE9UiciVwEjgcnOAN9y2pQT1NdBPRHqLiAfrguubNscUFSIiWNcpVhtjZtgdTzQZYyYbY7KNMb2w/k3fN8a0yV/axphtQI6IDIhsOh1YZWNI0bQZOFFEEiKf59Npox1C6ngTuDKyfCXwho2xRJWIjADuAM4zxlQd6HnaTIKKXJC7Efgv1gf9FWPMSnujipqTgTFYtYllkelsu4NSLeIm4AUR+R4YBNxvbzjREaklzge+AZZjfRe1maGAROSfwOfAABHJFZGrgQeBM0TkR+CMyHqrt4f3OhtIBt6LfD89fkDn1qGOlFJKxaI2U4NSSinVtmiCUkopFZM0QSmllIpJmqCUUkrFJE1QSimlYpImKKWUUjFJE5RSSqmY9P/Yemvdpb3iwgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "