From 79d11935799b5bfacb5abab3048dc6dd3b93faaa Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 19 Mar 2022 11:20:01 -0700 Subject: [PATCH 1/4] add not implemented test/fix for continuous MPC --- control/tests/optimal_test.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/control/tests/optimal_test.py b/control/tests/optimal_test.py index f059c4fc6..270a22131 100644 --- a/control/tests/optimal_test.py +++ b/control/tests/optimal_test.py @@ -117,7 +117,7 @@ def test_discrete_lqr(): assert np.any(np.abs(res1.inputs - res2.inputs) > 0.1) -def test_mpc_iosystem(): +def test_mpc_iosystem_aircraft(): # 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], @@ -171,6 +171,21 @@ def test_mpc_iosystem(): xout[0:sys.nstates, -1], xd, atol=0.1, rtol=0.01) +def test_mpc_iosystem_continuous(): + # Create a random state space system + sys = ct.rss(2, 1, 1) + T, _ = ct.step_response(sys) + + # provide penalties on the system signals + Q = np.eye(sys.nstates) + R = np.eye(sys.ninputs) + cost = opt.quadratic_cost(sys, Q, R) + + # Continuous time MPC controller not implemented + with pytest.raises(NotImplementedError): + ctrl = opt.create_mpc_iosystem(sys, T, cost) + + # Test various constraint combinations; need to use a somewhat convoluted # parametrization due to the need to define sys instead the test function @pytest.mark.parametrize("constraint_list", [ From 19c1d582a95b019cdd137b898cf25a61f2efbe2b Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 19 Mar 2022 17:29:16 -0700 Subject: [PATCH 2/4] optimal improvements: defaults, docs, t_eval, cost * Added config.defaults['optimal'] * Improved keyword handling (including unknown check) * Allow t_eval in input_output_response to chnage timepts * Save cost in optimal control result * Fixed doc/optimal.rst example (generates working solution) * Updated optimal docstrings and user documentation * Added optimization tips to doc/optimal.rst * Fixed user documentation warnings/errors (not related to optimal) --- control/config.py | 3 + control/iosys.py | 126 +++++++++++++++++++++------------- control/optimal.py | 69 ++++++++++++++++--- control/statefbk.py | 2 +- control/tests/optimal_test.py | 59 ++++++++++++++++ doc/control.rst | 2 +- doc/examples.rst | 1 + doc/optimal.rst | 76 +++++++++++++++----- doc/steering-optimal.png | Bin 39597 -> 29001 bytes doc/steering-optimal.py | 1 + doc/steering-optimal.rst | 14 ++++ 11 files changed, 276 insertions(+), 77 deletions(-) create mode 120000 doc/steering-optimal.py create mode 100644 doc/steering-optimal.rst diff --git a/control/config.py b/control/config.py index afd7615ca..8acdf28e2 100644 --- a/control/config.py +++ b/control/config.py @@ -103,6 +103,9 @@ def reset_defaults(): from .iosys import _iosys_defaults defaults.update(_iosys_defaults) + from .optimal import _optimal_defaults + defaults.update(_optimal_defaults) + def _get_param(module, param, argval=None, defval=None, pop=False, last=False): """Return the default value for a configuration option. diff --git a/control/iosys.py b/control/iosys.py index 2f63a7e8b..d5a7b755b 100644 --- a/control/iosys.py +++ b/control/iosys.py @@ -1571,7 +1571,7 @@ def __init__(self, io_sys, ss_sys=None): def input_output_response( sys, T, U=0., X0=0, params={}, transpose=False, return_x=False, squeeze=None, - solve_ivp_kwargs={}, **kwargs): + solve_ivp_kwargs={}, t_eval='T', **kwargs): """Compute the output response of a system to a given input. Simulate a dynamical system with a given input and return its output @@ -1654,50 +1654,57 @@ def input_output_response( if kwargs.get('solve_ivp_method', None): if kwargs.get('method', None): raise ValueError("ivp_method specified more than once") - solve_ivp_kwargs['method'] = kwargs['solve_ivp_method'] + solve_ivp_kwargs['method'] = kwargs.pop('solve_ivp_method') # Set the default method to 'RK45' if solve_ivp_kwargs.get('method', None) is None: solve_ivp_kwargs['method'] = 'RK45' + # Make sure all input arguments got parsed + if kwargs: + raise TypeError("unknown parameters %s" % kwargs) + # Sanity checking on the input if not isinstance(sys, InputOutputSystem): raise TypeError("System of type ", type(sys), " not valid") # Compute the time interval and number of steps T0, Tf = T[0], T[-1] - n_steps = len(T) + ntimepts = len(T) + + # Figure out simulation times (t_eval) + if solve_ivp_kwargs.get('t_eval'): + if t_eval == 'T': + # Override the default with the solve_ivp keyword + t_eval = solve_ivp_kwargs.pop('t_eval') + else: + raise ValueError("t_eval specified more than once") + if isinstance(t_eval, str) and t_eval == 'T': + # Use the input time points as the output time points + t_eval = T # Check and convert the input, if needed # TODO: improve MIMO ninputs check (choose from U) if sys.ninputs is None or sys.ninputs == 1: - legal_shapes = [(n_steps,), (1, n_steps)] + legal_shapes = [(ntimepts,), (1, ntimepts)] else: - legal_shapes = [(sys.ninputs, n_steps)] + legal_shapes = [(sys.ninputs, ntimepts)] U = _check_convert_array(U, legal_shapes, 'Parameter ``U``: ', squeeze=False) - U = U.reshape(-1, n_steps) - - # Check to make sure this is not a static function - nstates = _find_size(sys.nstates, X0) - if nstates == 0: - # No states => map input to output - u = U[0] if len(U.shape) == 1 else U[:, 0] - y = np.zeros((np.shape(sys._out(T[0], X0, u))[0], len(T))) - for i in range(len(T)): - u = U[i] if len(U.shape) == 1 else U[:, i] - y[:, i] = sys._out(T[i], [], u) - return TimeResponseData( - T, y, None, U, issiso=sys.issiso(), - output_labels=sys.output_index, input_labels=sys.input_index, - transpose=transpose, return_x=return_x, squeeze=squeeze) + U = U.reshape(-1, ntimepts) + ninputs = U.shape[0] # create X0 if not given, test if X0 has correct shape + nstates = _find_size(sys.nstates, X0) X0 = _check_convert_array(X0, [(nstates,), (nstates, 1)], 'Parameter ``X0``: ', squeeze=True) - # Update the parameter values - sys._update_params(params) + # Figure out the number of outputs + if sys.noutputs is None: + # Evaluate the output function to find number of outputs + noutputs = np.shape(sys._out(T[0], X0, U[:, 0]))[0] + else: + noutputs = sys.noutputs # # Define a function to evaluate the input at an arbitrary time @@ -1714,6 +1721,31 @@ def ufun(t): dt = (t - T[idx-1]) / (T[idx] - T[idx-1]) return U[..., idx-1] * (1. - dt) + U[..., idx] * dt + # Check to make sure this is not a static function + if nstates == 0: # No states => map input to output + # Make sure the user gave a time vector for evaluation (or 'T') + if t_eval is None: + # User overrode t_eval with None, but didn't give us the times... + warn("t_eval set to None, but no dynamics; using T instead") + t_eval = T + + # Allocate space for the inputs and outputs + u = np.zeros((ninputs, len(t_eval))) + y = np.zeros((noutputs, len(t_eval))) + + # Compute the input and output at each point in time + for i, t in enumerate(t_eval): + u[:, i] = ufun(t) + y[:, i] = sys._out(t, [], u[:, i]) + + return TimeResponseData( + t_eval, y, None, u, issiso=sys.issiso(), + output_labels=sys.output_index, input_labels=sys.input_index, + transpose=transpose, return_x=return_x, squeeze=squeeze) + + # Update the parameter values + sys._update_params(params) + # Create a lambda function for the right hand side def ivp_rhs(t, x): return sys._rhs(t, x, ufun(t)) @@ -1724,27 +1756,27 @@ def ivp_rhs(t, x): raise NameError("scipy.integrate.solve_ivp not found; " "use SciPy 1.0 or greater") soln = sp.integrate.solve_ivp( - ivp_rhs, (T0, Tf), X0, t_eval=T, + ivp_rhs, (T0, Tf), X0, t_eval=t_eval, vectorized=False, **solve_ivp_kwargs) + if not soln.success: + raise RuntimeError("solve_ivp failed: " + soln.message) - if not soln.success or soln.status != 0: - # Something went wrong - warn("sp.integrate.solve_ivp failed") - print("Return bunch:", soln) - - # Compute the output associated with the state (and use sys.out to - # figure out the number of outputs just in case it wasn't specified) - u = U[0] if len(U.shape) == 1 else U[:, 0] - y = np.zeros((np.shape(sys._out(T[0], X0, u))[0], len(T))) - for i in range(len(T)): - u = U[i] if len(U.shape) == 1 else U[:, i] - y[:, i] = sys._out(T[i], soln.y[:, i], u) + # Compute inputs and outputs for each time point + u = np.zeros((ninputs, len(soln.t))) + y = np.zeros((noutputs, len(soln.t))) + for i, t in enumerate(soln.t): + u[:, i] = ufun(t) + y[:, i] = sys._out(t, soln.y[:, i], u[:, i]) elif isdtime(sys): + # If t_eval was not specified, use the sampling time + if t_eval is None: + t_eval = np.arange(T[0], T[1] + sys.dt, sys.dt) + # Make sure the time vector is uniformly spaced - dt = T[1] - T[0] - if not np.allclose(T[1:] - T[:-1], dt): - raise ValueError("Parameter ``T``: time values must be " + dt = t_eval[1] - t_eval[0] + if not np.allclose(t_eval[1:] - t_eval[:-1], dt): + raise ValueError("Parameter ``t_eval``: time values must be " "equally spaced.") # Make sure the sample time matches the given time @@ -1764,21 +1796,23 @@ def ivp_rhs(t, x): # Compute the solution soln = sp.optimize.OptimizeResult() - soln.t = T # Store the time vector directly - x = X0 # Initilize state + soln.t = t_eval # Store the time vector directly + x = np.array(X0) # State vector (store as floats) soln.y = [] # Solution, following scipy convention - y = [] # System output - for i in range(len(T)): - # Store the current state and output + u, y = [], [] # System input, output + for t in t_eval: + # Store the current input, state, and output soln.y.append(x) - y.append(sys._out(T[i], x, ufun(T[i]))) + u.append(ufun(t)) + y.append(sys._out(t, x, u[-1])) # Update the state for the next iteration - x = sys._rhs(T[i], x, ufun(T[i])) + x = sys._rhs(t, x, u[-1]) # Convert output to numpy arrays soln.y = np.transpose(np.array(soln.y)) y = np.transpose(np.array(y)) + u = np.transpose(np.array(u)) # Mark solution as successful soln.success = True # No way to fail @@ -1787,7 +1821,7 @@ def ivp_rhs(t, x): raise TypeError("Can't determine system type") return TimeResponseData( - soln.t, y, soln.y, U, issiso=sys.issiso(), + soln.t, y, soln.y, u, issiso=sys.issiso(), output_labels=sys.output_index, input_labels=sys.input_index, state_labels=sys.state_index, transpose=transpose, return_x=return_x, squeeze=squeeze) diff --git a/control/optimal.py b/control/optimal.py index 493b6bc3d..9dc8b225b 100644 --- a/control/optimal.py +++ b/control/optimal.py @@ -17,9 +17,19 @@ import time from .timeresp import TimeResponseData +from . import config __all__ = ['find_optimal_input'] +# Define module default parameter values +_optimal_defaults = { + 'optimal.minimize_method': None, + 'optimal.minimize_options': {}, + 'optimal.minimize_kwargs': {}, + 'optimal.solve_ivp_method': None, + 'optimal.solve_ivp_options': {}, +} + class OptimalControlProblem(): """Description of a finite horizon, optimal control problem. @@ -110,6 +120,10 @@ class OptimalControlProblem(): values of the input at the specified times (using linear interpolation for continuous systems). + The default values for ``minimize_method``, ``minimize_options``, + ``minimize_kwargs``, ``solve_ivp_method``, and ``solve_ivp_options`` can + be set using config.defaults['optimal.']. + """ def __init__( self, sys, timepts, integral_cost, trajectory_constraints=[], @@ -126,13 +140,22 @@ def __init__( # Process keyword arguments self.solve_ivp_kwargs = {} - self.solve_ivp_kwargs['method'] = kwargs.pop('solve_ivp_method', None) - self.solve_ivp_kwargs.update(kwargs.pop('solve_ivp_kwargs', {})) + self.solve_ivp_kwargs['method'] = kwargs.pop( + 'solve_ivp_method', config.defaults['optimal.solve_ivp_method']) + self.solve_ivp_kwargs.update(kwargs.pop( + 'solve_ivp_kwargs', config.defaults['optimal.solve_ivp_options'])) self.minimize_kwargs = {} - self.minimize_kwargs['method'] = kwargs.pop('minimize_method', None) - self.minimize_kwargs['options'] = kwargs.pop('minimize_options', {}) - self.minimize_kwargs.update(kwargs.pop('minimize_kwargs', {})) + self.minimize_kwargs['method'] = kwargs.pop( + 'minimize_method', config.defaults['optimal.minimize_method']) + self.minimize_kwargs['options'] = kwargs.pop( + 'minimize_options', config.defaults['optimal.minimize_options']) + self.minimize_kwargs.update(kwargs.pop( + 'minimize_kwargs', config.defaults['optimal.minimize_kwargs'])) + + # Make sure all input arguments got parsed + if kwargs: + raise TypeError("unknown parameters %s" % kwargs) if len(kwargs) > 0: raise ValueError( @@ -271,9 +294,10 @@ def _cost_function(self, coeffs): logging.debug("initial input[0:3] =\n" + str(inputs[:, 0:3])) # Simulate the system to get the state + # TODO: try calling solve_ivp directly for better speed? _, _, states = ct.input_output_response( self.system, self.timepts, inputs, x, return_x=True, - solve_ivp_kwargs=self.solve_ivp_kwargs) + solve_ivp_kwargs=self.solve_ivp_kwargs, t_eval=self.timepts) self.system_simulations += 1 self.last_x = x self.last_coeffs = coeffs @@ -393,7 +417,7 @@ def _constraint_function(self, coeffs): # Simulate the system to get the state _, _, states = ct.input_output_response( self.system, self.timepts, inputs, x, return_x=True, - solve_ivp_kwargs=self.solve_ivp_kwargs) + solve_ivp_kwargs=self.solve_ivp_kwargs, t_eval=self.timepts) self.system_simulations += 1 self.last_x = x self.last_coeffs = coeffs @@ -475,7 +499,7 @@ def _eqconst_function(self, coeffs): # Simulate the system to get the state _, _, states = ct.input_output_response( self.system, self.timepts, inputs, x, return_x=True, - solve_ivp_kwargs=self.solve_ivp_kwargs) + solve_ivp_kwargs=self.solve_ivp_kwargs, t_eval=self.timepts) self.system_simulations += 1 self.last_x = x self.last_coeffs = coeffs @@ -548,7 +572,7 @@ def _process_initial_guess(self, initial_guess): initial_guess = np.atleast_1d(initial_guess) # See whether we got entire guess or just first time point - if len(initial_guess.shape) == 1: + if initial_guess.ndim == 1: # Broadcast inputs to entire time vector try: initial_guess = np.broadcast_to( @@ -804,6 +828,15 @@ class OptimalControlResult(sp.optimize.OptimizeResult): Whether or not the optimizer exited successful. problem : OptimalControlProblem Optimal control problem that generated this solution. + cost : float + Final cost of the return solution. + system_simulations, {cost, constraint, eqconst}_evaluations : int + Number of system simulations and evaluations of the cost function, + (inequality) constraint function, and equality constraint function + performed during the optimzation. + {cost, constraint, eqconst}_process_time : float + If logging was enabled, the amount of time spent evaluating the cost + and constraint functions. """ def __init__( @@ -833,15 +866,19 @@ def __init__( "unable to solve optimal control problem\n" "scipy.optimize.minimize returned " + res.message, UserWarning) + # Save the final cost + self.cost = res.fun + # Optionally print summary information if print_summary: ocp._print_statistics() + print("* Final cost:", self.cost) if return_states and inputs.shape[1] == ocp.timepts.shape[0]: # Simulate the system if we need the state back _, _, states = ct.input_output_response( ocp.system, ocp.timepts, inputs, ocp.x, return_x=True, - solve_ivp_kwargs=ocp.solve_ivp_kwargs) + solve_ivp_kwargs=ocp.solve_ivp_kwargs, t_eval=ocp.timepts) ocp.system_simulations += 1 else: states = None @@ -858,7 +895,7 @@ def __init__( # Compute the input for a nonlinear, (constrained) optimal control problem def solve_ocp( - sys, horizon, X0, cost, trajectory_constraints=[], terminal_cost=None, + sys, horizon, X0, cost, trajectory_constraints=None, terminal_cost=None, terminal_constraints=[], initial_guess=None, basis=None, squeeze=None, transpose=None, return_states=False, log=False, **kwargs): @@ -966,6 +1003,16 @@ def solve_ocp( return_states = ct.config._get_param( 'optimal', 'return_x', kwargs, return_states, pop=True) + # Process terminal constraints keyword + if constraints is None: + constraints = kwargs.pop('trajectory_constraints', []) + + # Process (legacy) method keyword + if kwargs.get('method'): + if kwargs.get('minimize_method'): + raise ValueError("'minimize_method' specified more than once") + kwargs['minimize_method'] = kwargs.pop('method') + # Set up the optimal control problem ocp = OptimalControlProblem( sys, horizon, cost, trajectory_constraints=trajectory_constraints, diff --git a/control/statefbk.py b/control/statefbk.py index a866af725..6598eeeb8 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -734,7 +734,7 @@ def dlqr(*args, **kwargs): The dlqr() function computes the optimal state feedback controller u[n] = - K x[n] that minimizes the quadratic cost - .. math:: J = \\Sum_0^\\infty (x[n]' Q x[n] + u[n]' R u[n] + 2 x[n]' N u[n]) + .. math:: J = \\sum_0^\\infty (x[n]' Q x[n] + u[n]' R u[n] + 2 x[n]' N u[n]) The function can be called with either 3, 4, or 5 arguments: diff --git a/control/tests/optimal_test.py b/control/tests/optimal_test.py index 270a22131..53d8fe8e4 100644 --- a/control/tests/optimal_test.py +++ b/control/tests/optimal_test.py @@ -45,6 +45,9 @@ def test_finite_horizon_simple(): np.testing.assert_almost_equal( u_openloop, [-1, -1, 0.1393, 0.3361, -5.204e-16], decimal=4) + # Make sure the final cost is correct + assert math.isclose(res.cost, 32.4898, rel_tol=1e-5) + # Convert controller to an explicit form (not implemented yet) # mpc_explicit = opt.explicit_mpc(); @@ -564,3 +567,59 @@ def final_point_eval(x, u): optctrl = opt.OptimalControlProblem( sys, time, cost, terminal_constraints=final_point) res = optctrl.compute_trajectory(x0, squeeze=True, return_x=True) + + +def test_optimal_doc(): + """Test optimal control problem from documentation""" + def vehicle_update(t, x, u, params): + # Get the parameters for the model + l = params.get('wheelbase', 3.) # vehicle wheelbase + phimax = params.get('maxsteer', 0.5) # max steering angle (rad) + + # Saturate the steering input + phi = np.clip(u[1], -phimax, phimax) + + # Return the derivative of the state + return np.array([ + np.cos(x[2]) * u[0], # xdot = cos(theta) v + np.sin(x[2]) * u[0], # ydot = sin(theta) v + (u[0] / l) * np.tan(phi) # thdot = v/l tan(phi) + ]) + + def vehicle_output(t, x, u, params): + return x # return x, y, theta (full state) + + # Define the vehicle steering dynamics as an input/output system + vehicle = ct.NonlinearIOSystem( + vehicle_update, vehicle_output, states=3, name='vehicle', + inputs=('v', 'phi'), outputs=('x', 'y', 'theta')) + + # Define the initial and final points and time interval + x0 = [0., -2., 0.]; u0 = [10., 0.] + xf = [100., 2., 0.]; uf = [10., 0.] + Tf = 10 + + # Define the cost functions + Q = np.diag([0, 0, 0.1]) # don't turn too sharply + R = np.diag([1, 1]) # keep inputs small + P = np.diag([1000, 1000, 1000]) # get close to final point + traj_cost = opt.quadratic_cost(vehicle, Q, R, x0=xf, u0=uf) + term_cost = opt.quadratic_cost(vehicle, P, 0, x0=xf) + + # Define the constraints + constraints = [ opt.input_range_constraint(vehicle, [8, -0.1], [12, 0.1]) ] + + # Solve the optimal control problem + horizon = np.linspace(0, Tf, 3, endpoint=True) + result = opt.solve_ocp( + vehicle, horizon, x0, traj_cost, constraints, + terminal_cost= term_cost, initial_guess=u0) + + # Make sure the resulting trajectory generate a good solution + resp = ct.input_output_response( + vehicle, horizon, result.inputs, x0, + t_eval=np.linspace(0, Tf, 10)) + t, y = resp + assert (y[0, -1] - xf[0]) / xf[0] < 0.01 + assert (y[1, -1] - xf[1]) / xf[1] < 0.01 + assert y[2, -1] < 0.1 diff --git a/doc/control.rst b/doc/control.rst index 87c1151eb..8bd6f7a32 100644 --- a/doc/control.rst +++ b/doc/control.rst @@ -114,7 +114,7 @@ Control system synthesis lqe mixsyn place - rlocus_pid_designer + rootlocus_pid_designer Model simplification tools ========================== diff --git a/doc/examples.rst b/doc/examples.rst index 91476bc9d..89a2b16a1 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -29,6 +29,7 @@ other sources. robust_mimo cruise-control steering-gainsched + steering-optimal kincar-flatsys Jupyter notebooks diff --git a/doc/optimal.rst b/doc/optimal.rst index e173e430b..8da08e7af 100644 --- a/doc/optimal.rst +++ b/doc/optimal.rst @@ -79,7 +79,7 @@ Every :math:`\Delta T` seconds, an optimal control problem is solved over a :math:`T` second horizon, starting from the current state. The first :math:`\Delta T` seconds of the optimal control :math:`u_T^{\*}(\cdot; x(t))` is then applied to the system. If we let :math:`x_T^{\*}(\cdot; -x(t))` represent the optimal trajectory starting from :math:`x(t)`$ then the +x(t))` represent the optimal trajectory starting from :math:`x(t)` then the system state evolves from :math:`x(t)` at current time :math:`t` to :math:`x_T^{*}(\delta T, x(t))` at the next sample time :math:`t + \Delta T`, assuming no model uncertainty. @@ -219,9 +219,11 @@ with a starting and ending velocity of 10 m/s:: To set up the optimal control problem we design a cost function that penalizes the state and input using quadratic cost functions:: - Q = np.diag([0.1, 10, .1]) # keep lateral error low - R = np.eye(2) * 0.1 - cost = opt.quadratic_cost(vehicle, Q, R, x0=xf, u0=uf) + Q = np.diag([0, 0, 0.1]) # don't turn too sharply + R = np.diag([1, 1]) # keep inputs small + P = np.diag([1000, 1000, 1000]) # get close to final point + traj_cost = opt.quadratic_cost(vehicle, Q, R, x0=xf, u0=uf) + term_cost = opt.quadratic_cost(vehicle, P, 0, x0=xf) We also constraint the maximum turning rate to 0.1 radians (about 6 degees) and constrain the velocity to be in the range of 9 m/s to 11 m/s:: @@ -230,20 +232,19 @@ and constrain the velocity to be in the range of 9 m/s to 11 m/s:: Finally, we solve for the optimal inputs:: - horizon = np.linspace(0, Tf, 20, endpoint=True) - bend_left = [10, 0.01] # slight left veer - + horizon = np.linspace(0, Tf, 3, endpoint=True) result = opt.solve_ocp( - vehicle, horizon, x0, cost, constraints, initial_guess=bend_left, - options={'eps': 0.01}) # set step size for gradient calculation - - # Extract the results - u = result.inputs - t, y = ct.input_output_response(vehicle, horizon, u, x0) + vehicle, horizon, x0, traj_cost, constraints, + terminal_cost=term_cost, initial_guess=u0) Plotting the results:: - # Plot the results + # Simulate the system dynamics (open loop) + resp = ct.input_output_response( + vehicle, horizon, result.inputs, x0, + t_eval=np.linspace(0, Tf, 100)) + t, y, u = resp.time, resp.outputs, resp.inputs + plt.subplot(3, 1, 1) plt.plot(y[0], y[1]) plt.plot(x0[0], x0[1], 'ro', xf[0], xf[1], 'ro') @@ -252,15 +253,13 @@ Plotting the results:: plt.subplot(3, 1, 2) plt.plot(t, u[0]) - plt.axis([0, 10, 8.5, 11.5]) - plt.plot([0, 10], [9, 9], 'k--', [0, 10], [11, 11], 'k--') + plt.axis([0, 10, 9.9, 10.1]) plt.xlabel("t [sec]") plt.ylabel("u1 [m/s]") plt.subplot(3, 1, 3) plt.plot(t, u[1]) - plt.axis([0, 10, -0.15, 0.15]) - plt.plot([0, 10], [-0.1, -0.1], 'k--', [0, 10], [0.1, 0.1], 'k--') + plt.axis([0, 10, -0.01, 0.01]) plt.xlabel("t [sec]") plt.ylabel("u2 [rad/s]") @@ -272,6 +271,47 @@ yields .. image:: steering-optimal.png +Optimization Tips +================= + +The python-control optimization module makes use of the SciPy optimization +toolbox and it can sometimes be tricky to get the optimization to converge. +If you are getting errors when solving optimal control problems or your +solutions do not seem close to optimal, here are a few things to try: + +* Less is more: try using a smaller number of time points in your + optimiation. The default optimal control problem formulation uses the + value of the inputs at each time point as a free variable and this can + generate a large number of parameters quickly. Often you can find very + good solutions with a small number of free variables (the example above + uses 3 time points for 2 inputs, so a total of 6 optimization variables). + Note that you can "resample" the optimal trajectory by running a + simulation of the sytem and using the `t_eval` keyword in + `input_output_response` (as done above). + +* Use a smooth basis: as an alternative to parameterizing the optimal + control inputs using the value of the control at the listed time points, + you can specify a set of basis functions using the `basis` keyword in + :func:`~control.solve_ocp` and then parameterize the controller by linear + combination of the basis functions. The :mod:`!control.flatsys` module + defines several sets of basis functions that can be used. + +* Tweak the optimizer: by using the `minimize_method`, `minimize_options`, + and `minimize_kwargs` keywords in :func:`~control.solve_ocp`, you can + choose the SciPy optimization function that you use and set many + parameters. See :func:`scipy.optimize.minimize` for more information on + the optimzers that are available and the options and keywords that they + accept. + +* Walk before you run: try setting up a simpler version of the optimization, + remove constraints or simplifying the cost to get a simple version of the + problem working and then add complexity. Sometimes this can help you find + the right set of options or identify situations in which you are being too + aggressive in what your are trying to get the system to do. + +See :ref:`steering-optimal` for some examples of different problem +formulations. + Module classes and functions ============================ diff --git a/doc/steering-optimal.png b/doc/steering-optimal.png index 6ff50c0f423ca3c58abffeb34f6be372333e1659..f847923b5f8f86f6901d2d60c95d714f36289f18 100644 GIT binary patch literal 29001 zcmbrm1z1(>w=KSCBqXGzRYZ`G25B}SDhL81CEbm5H;5pJh>C!qNDC<45&|M3(jl!N z-68d#%kOv2J?GwY&;Q=%`h0%yxcAz7tvBWzbBrSp8QD99_wD|DIF z`u=@qHwivIhkt(quam14AGy|-YWR>-&MJ4^5QNMG{fm_&lWl__&gIuIa#~(5mw$MB z-PS&qS@V1>qM!2C_xpFlotvDT?+L}YSTm;uRtOi?s1>sYG%%kuZsklyeG;bS;s?yM~0F-4_2p zet1z{UfzbsDU5U^R(uXles1o^58cu+a-sNg`TgudLR8I>WNpd!0*2odV%R?>s!K2-$NvFVrvrjMI4QXk?;Ns((?yXJj{Y<|1f^oN* zOa{}PCY5oE^IL_}y#oD;)4gwRD#*)++`jr@_~D?DFKM0cZiIY>ECz$i9I$6XMNMtE z<@QX(fdLx_XMNPE!Q^mnosf`FR!fV@DDcP~Y2<4-4AFS4#Ky{sNcrxx$Q-YeJxV_q zDlxm_v80<^R{nGBqI%lHy*;lX`*My0MiG0u$BBtd-QncA^?p?|-HIwIge%pHye}o( z^RC4n9qhIZ6zX3wZ6VC*x0T#zKi3ub4Kqt05>;!uz0i*xG*o7H=YzxN&!1aLZF*Zq z%IsQyb)~d=PWWFuKAe*w&eh^NUVjuQR))t)NlDq#f9D?0@zH8vj8px-(nF@JI(z&3 z(TNw;B}SN}{dBF%^sX2;BJ1mIEHO7zM53+I{dH8bG}k|K7a5eUT)O?1FdD;typ=5z z4bQ0Qcd%Sx`Y{A=dbwi4qcr-_Bibuhu4Fn*)KAYmlA)OhJgM88VL9PMq8>gZQLdQ? zIJhZ!yf-u5t1d&WUv5uvFYu_$YpX{&9Hv<2c<=RVLmi!`Glv{Ep0U>NO%t!T(%y^0 z!oiP5y24iGn zV-x2v)Ao0in~eFz>P|W?rc+tWLvPM-EfWO6w8B;+(PttI|*UsJJz_68e{FN75noUJqB ziy%Gwt8`*%$$n~eb(JdQc)KX@6fzAv+|9PeZjcvC_VT^@+EzDMp|**;g98H(qM`_~ zWS4%_43-=n?r+L-Bu4l3X+27W=SnqsbW$&IL=hq> zsFZ9-ajK^mc>ZPGZ?Z8MnH-ucF>AYZ>sC#|$4oeysfY zRR5)PBpM8h2F*I3^3X~tM`J>w&s>y9JJ{dP@ZDXBdhnoW@?$8|sDgq*YCVDa%Y!WI z4a7OsT^0#~LgU4Gq1orQ2azMD^7*G$_mk=(l{16zv*xd5~iB4oYW4W;Gq=w!O#l^75>S z9sCxbpIG{?Q=U1K;vWb>yFiB0yhEyGbn zO+_UpMMg%ZtgA~`f3R{j`M`E4apBa8v}faZr{2Wv!CNzk$ECS(ND&79+x@zAKElkt zOIHev>KJ-@dgw$QU-qKY6Y%HvMP+`6ob~dEz^BpeSJ4^akH%nE8LVuPS<1-;ciwF* z%+Seq`QbijGEwhOC+5sL>$-oT?&NKXqr9FTgHC~N@c!mpMQ%jIdCF^v&d!sKSW81C zI3*<|$+Znebw049W+qOKMKy2VW*T*KbuIQk@W3c5UoF4+Lg-X_z@JO-xUhZ(b`q`6 zd8S`Az87>LbFxo%gCo&wVy-*AZD)Cuf{JR^+B0|V7{+-9LVa7g!{{0$1~y^gej72L z?{BVo9c>RLBqtA=@1;gbaG)V-wANctEnWJD`9VwcnJYHm&)W|bW2sB7Ho_$lVO!pL#+1#Y-}FobHkLY801CuAInrvNFqjZCyBoLdh(*$wpCn~Fs9#Y1K(21r8f-QDQO z@vcT zQsM?-xj|H=Evo;+b-a<>@T$0&90Cjky~y6Mh+imJ>*(l!j72?Aa2F>@HN4*cU>fpf z?X2Ffuc%KFVq!-wE%=A>O`RSqG^ zH!c=6#Em-EkY2lX&1?T>@_Lu(gxK`merxOEBd zO@@*c`|Z09*^9suHk)lvsQj~#KYCvpGb>KX7PjU=fx)C{!F;Ak`I>IbcIZ~QW>s{u zv2RH1Ae0)(nh=5Uh-oIP7vCwswtfHM!}@@}Q$}qq;|I&G4G(V&h6x)mtNSI0cG>Fb zpU$XD20wmG2`~!%axYk^Jw^9qoRXC>eQeHz&;CMWXEy8q4L1IN0@nZU7q3IyW72>K zA!0ui1o?}UD~afw(!<6~#h7pTcZzgYn=L;~!RSIBj8{3&T<{bgr<#%t7q%S*6BB7} zpP;`xI;QPuS4&Qv;g}rvU2WCAV~7{i8dBvpHx0>1%ypI!8yj1%?2(|=S5^^`f#xJA z7c}l9rKH5b#6Ma}!78SH@ZbSCJ$*DFI#}h^+N)RSNC;1#R&sRY4G9U6G+=`V^?uk- zw6#k1QfkyM#qmc_C;`^QYM&x^Huf01FjEv{g&K$Uv&EV^i19l-MVe3c(ag{BZj(6 zCs;3^v`pvh7n-xSh1fbNq%J8ampo`DG>?+O~x78)LjSIxf0E z;Iro1E|X`;UZPzxn=flbYjYQR#`1rj3$K z5`C6-8NdG6naT5%vOoDk^>{;@)bX<#`qpo+Se^X1#}PiJj$L(wNai`!vB0S}l&22g zP`1C(JZ~ElRHWxpoWohsH^xkd=&$wNRyToqX$bb6B43jhC)F)Zs%0&Dn9nwJ{DL#^ z_ahiy3XZUZ|KvDOe|;+ar!Xg%XvTN@x>^5D3r*|(V?v#Mh5QZW;t>`()j-_b!wii- zY$`S<^f_A!ubPeta}IK!cD{OgAorWE?l7Y-x-RG^-*BY2;my=|m7r&IlOuL~t76Mh zSl{s2aMOGiF8Lkr&ken|AcnMDXNaz=nh&}Tocsfyjf)!(dg z-Z^X^KRMRKSW^*_mC7_>&k_*SfVw8`G|6a6-Uvu4}3|r>NJbjIV+Ja@r-L3ncuQ}lk zQr&pt1_qNFJ#4C=c8&V|18x`ihU|raH0gWMehN2ZR=1O{edEy$AG#%MD&H4%N~2$D zsX-Z+i#6$wN~Im?^Dl+h>%553T3_QLYnZG{IQ_*KF6Z|-*x3)bT=LK#+zI-c;ZK!p zg;izFM-Z^`YQ^sH3r{QHke|&{GI+-DhZXjj7fK97Mdjd>4C8Km&U48u#%`c}I)gF6 zQ=`(EG}o|*Im3%RD{x*3m-UJyd9;b``7i?b4LN0kMW}Od@RgjWj0nAl1lm4=x(VV= zmx5Q=4e8u%tL`R6)dSm)FHxQ%Ka<%+O%bd@MAkx`fVGoOffe-XxKxFZEJ!Al*&FNc;1b)51TQ`2 zVwARPO!%M#KUxy}&OP{SxVvDZeLR0AOU}uk{13+F+jW+(v})O$pNAKdB~BmP0I`#2 zRx;w8nS&j<@lfoE6`7y%Mr;oM-4~@vDnVhz$X4uof;wppx%o5HVVs4J6R+eno1t`2 zx_SAK^Pz$G>7Vq@^95K0n-BMRCI}N-2i_DmsXrL^A^oFd*CQ%H=>Iz8n|_g8P{_$c z&AAqtDh#IYXTP4q-jM;`>0e1xjr2I^CkQn#TH5s%XKUe5fw zNNz*PsMM$AJFQQ4=NOzPZ_o?4bvko)PrCX`#eOHlFHdusG?FNLdLzDb~C4hp`aNAiI5Q&P#t`Lon5qpv9pGf&eZP}cF~ z)&cc3Ugt}}!m^qH6>|FzuA4M48e~2PuAk$d>o2%+iei>m34KBrm_OQ9(QXRm)d?3na{?ad0oTiy>Z|>siozvOU?B! zHJ|A8i`Wc8G>|L#NXJO~@S_XcJZ)EoLFUf)DePBW16w|io<2!V_RXPWXRoGX!*#vB zG`O(fMRb?l5i5kUY0z znQZFmacajrIis?v?~P%pPW0>^`T4eRR&i-uSa7VnGX(~jahc5ojkUFshliWpoR=@# zf4{m1Wub2_7d!jJEiFgat-s&=Q80o-(Cx^Aaw}%QRK<@lqAy&!?rEJbhUsqP=^%hY zP$vkoCU;WM(EPg3R*Yds$xzpGDK{rE8h9S*<-04(UA*p9D@M9OV-!0PEs?KaAr}2O zw_tu=p3r68avX&z^n#YgUtfyHt6m-a7W{8VHCitYt@NXB{w9VJW1 zGji%kw`d8vbQWx*N@uqElVjhlg??h>#SQHB(x08lGqIxMvH=0of3_DRv#uvkkJuG$ z@3}-r@qW6c9U~=vA+E9F1#dt zUR-*<^><#$Z_-zH^KgAvjT)?Z2FK5_CrfOR7aJ9|smmONF-!TB9+X6qo{gBDH80RD zX(}lM0{WSxw*XW#yng~tW1Z2^zw3alycs*>UyWPyVR>+kn$GS;>2)G9=dBgPB=XP0 zRT#stZZfRy7sA9)`SFf>&!u`uy&Jaw1Qp4P8`AWmj*R6_6XWm7$TbG zxL&e&|J2UUCmBuVF19lllh6hBE|1txVOoKlIcH6?71KBDK95wo&8gXs*U|%-DGZ}H zXinJ?e}g3W9pr$hSJGh6E_8lwmudY}8EL)mo3kPYRdhqbwuuz_ar^jo0pQ z#zaX_u;BteC^;yt^K|m`xCFnal=c!_a?v);sf*;3rmEXwXn(`_^#XhM{5LFTV4ID^ z9lgB7%mwZyj@LIeHF4C1Yy*iF3S3{C+v}bvU^gEPT*ELu8a8d7`hedy;#;GhJnffN zb*<*lPu+!!_Riuo?6&Gm5{OTX0dHG8n?ihAS{zX5KT4#eq?jsxc%*dha#9g-Qq?j? zwhz+f*_DX4Nk6GMxPJ+ffw3(->6rQhxho&(9^VpubC%_EjU9Qy9rsG)*~^y?SpxR3 z*x1+#ZDnlxI1BDIFaz)FX2+D9oBK?_f+*rN(}StzNVA`xliE8vK0=5X^I>FPkX%SM z8@?@i*G%2xQd}yx@=O_(eBUvwf4#wG+K}TPrMLFh|46Yfh(Tcf^5shllnfLM3?mCW zq!;g1c1rDzg`)HTFcMJgE1w)6p>;2-wDkCL8+&G4hiMrK31AFN)yw@=c3+nK4s@by z>VuD8c$VR-h(XYeps?t?hyaVLa-OoZ{4(%_z=APB8O{7|>AN$NZ06TC-PNTJE@{sNvl^ zH3patf}^p@5A#^&gy48rn7Rd17zPNYkxd&6OI^IpkObgIYp^8FmsrcWKb=dSgLHbd zdyovsb2WT&P1!SdP>m@yt<|>K37}YU1%yiR3)(kF!tvS@3~15lYb-0u^W+#^WD|ei z`*ok_jDXvs@9yW~lN6mMotJ29Ad_vaBYhY=xEKcR!TP;s3t)5GY-O&22aBHZ z5g<{MGVzm(DqL7T3qo`TgEh&JsNXp6tWYGZT0bJ`nZhwvIM`WkR%Y=JYdAStB_ttf zf#?~9A_pm-tu+suoBYUQzBBT#8g)aDsw{>AKaDqg@^=RjC1zb4C5**})U)WbE6>lb z6oV*6GQVVH-{4HIJ)`0VTx^L=FNgg|=_!HdB*34yz!vHOPU^BzEnTNkM7PV{MX9ntcT{Bd|Ed#7*dF-bkx*%DP0jLJ$Kb}rR9}e=!;T%nE*+U z5J&4o*zeskc|Yl#2H=+jkaHs}Ts3O#YWo!Ox1Sp6e81HQD8&Ty=#2{|N9BYre8XBR ztQc=%IoVU@OzN6=aq0YdTp;b1b*VK5&2!$h{)5`d{NTgdx!#54=ZDqc<$r6?kC$3& z#WjOvZS|9q!krx+Snat$>*xh?_iKtjIZ(Vm1$bZFt_58;HM$1$NM@<(*6Q2VMn^pu zG*h^;37dM^;0gsr)oqC3nET#PDvW*4t4@=8eQt}o`q?Y2A>W-piUcO`)<1*cbu>64Ton>B*#_&vb8FkNMDi*@W5h{kC^USb zELbpyj~2}8-z#j4f`zs$v)P&8KEA!}<0^B>XHgg+ZCHJe+&nE)*!UmTcV)_$+2%by zM1w#_)8>;Co)RNBXWP2l&UoDn-cO)J49m&Lv|gz{ky_f0`~K=)Db2MhECZkLsHm5K zyxkECZOcO=-QbzI&+7-`2i@>f!@$I2*oZd~YB5%T-`|v3!i=C0!QipUH!p5SMa@q` zEjJQPj;#nB31>lLM-vhHW7VDPj%OFLG`P?jp4W3bTaGUkziw+s+0+X;9L`t9U>sTy zm03#&5AVD0OlfK0K}j}^{-7FGNH2N&JM5uPkKP5J1!1&2gB9fSTX>$Cp0Mm?;LYf$ z_PRYjVvVyy_Skcu?c%i!Yyzo<1Jk(PVW65W5TX5W1yQ*mYoEh#2ZNgl*Fv<&%vKBK zX2+puzeW$mlG+~MX$`>wkzk-R+ArnG@9cA#Cu>dLce$WOI~rJ&TPA#9RAvzICJH$p zPB8r-O7^@W?kjkH%4VEbTXV}jP$<;UAstqy*o+IxSu0A3sxD}wpW_JhglZmIdJAT0 z=af~ibl@P*wL*r-sUNd1edLAO8AcmQOj%$nqA}m&tQ6p8d`N=&EQal|Ta!3ulXPb6 z?~+0>?p$7gZ>5^={1%uO)qF$V$7AsZ9Evw`h@hx#iNs@NyAiv$hMbSUd!--s;LH`M z;%}QGOAiZlCbQFn4*giUh$?sVWPj;5KcLAcZKbxPqO-eYobyrKjC6xkxLs81YCxLc zGvjg*nH&0ilDFbeqCufAidg%ch?d2P)QwoqLLz}nz8mnlb|#4ZEJ(3lAPSUtZ(4#{ z#3~`dfKu6@o12vJ?ppL!02 zA69V!%B6~=kh1u#VO4vs#wkReHUT3-PC-G-#PN>kDouq6j`N_QwQECBt>nSZ$#I23 zxNKr8k%Bu9iFO%h{yRmREFkSHpJ)vX9F@v^z4Oeu9_92bwFax}o5Q&Ua&dt*x-X|X zANu5*T~gF{E(+>`_a}W^8N`~d&w!=Pp@=4bx5BJSTDSh`tZs=Zo?5Cz|H3X^tP-g4 zbR_t|h8hM?=y1o!PmXO4#$!xZPks4ZT`*;waS^@>YiSl-@|GrIL60t<@6f3K%kpPlrQdvtVM@VkvUu#+)P#eajBwWZwDz4`DiR_ZEYRnaD(H_jkm;}Zu|_fD?HQL=`Y9vgX-?^ zi_)$(lruE!RuKXDs$)WSzf4E}heD&lZq2Z1Xou!WuPvFs<$=2AV-=+EGmBo5%>N+o za`fWj5%Iswo)+unwWd`pS9+^&zzlkm<4h{Pl=PFgno<9^*{Of1WP;uPDd9H<>i@}Ef2R^3?5I8mKE_3uolk+Gm{@_~1 zTB6+4#qcI|q9g-+DDFO26GYwtlg-FIUlOc?8^~o|UL{bCV;sPe@k)yeYp8hw)T>C7 z<_MFn7sBzXd=BivZSSzYblORr!gU@MMr}N^hPC+LVw;+qDVUkZ^UA%wy(>X`<2P?7 z2NCkm<7O7)hECkcUv^l%--9mQesPtx*qQCt-QV>yz&UHwrWP#2%H8qud?}Xo|1jS+R@RYT{QAmPzMMz8>s#o|4Kt{1uH*?$K;c3Y>5Yzbr_feJ`DCYPlPr^XV z1!YEjWNp(k9V{BtPxT$;?%x*ztr#`AWDn|Z-c7l|A-Uamn-J7zSy{xfZd(VY<1JJ7 zG$rAJmh{oqZE#`mNgu8$xXiXwpb{npIXNq+U3G?f$$|lgUL`);wl4GE*b<&UzvNiI zZ%IK*dy$Kai_rnKBW3hhw}6d@<217*=StNaGX*vEg@;P7X}MQ@JXXdop`M_|V4PNv z5*bJL_kDG)U5kzaCl4r_XBv+GNueGYZ<&3$n*aF z`Dy3qVZfAt@$so1 zs_`sBS8q$w9ivN^X|7*Q?`pk*cnE;hTLw0ph~9GG&0=Vcygpv1-j8c)lE;!YH#a9_ z3o9PBmoFNOf$#p&^Y-mqG(@7Zrf>ci-t5|1gRS`$$a$zgmq*{x;h#C0A{!0BjE;FA@0COAVefES9jX;d#L|D1GiBY4^L?kR)%aCP5rSV?4gMI1a?+_#2D6i2z$q- zH_W9f6rjqOk>9%F(pc3Q7zZIR!?!;@~U{5T7x02otIzOzPoEL|HMc6Rp`kH9T2a(nF75w^_2NYWcc$3W_4ymZCdyb?5_d#(^$tySS-F4u3)BPyV zir6=Tc`xDFv-vF_Q&aXWpSs>%6Oj~%Q6E>EN#8Xb`!Gs+kw1U_JgyB!YP}JMj8#zZ z43Gtfhf?++pWl2cEfH`i297wawbh_(aJl9Kd$1ze|x^#)1l~wkowExIxmb*9%Smig{GE=ZJ zvCGIztmowAxt4l^Mo85hn zSsJbAF8QJd3k|4QboSPv)Kv#>#2wWBG%3J8E4^B`n}&ui@Rpswc=77f3v9dClN)Go zR6d*1#bFe%Z;>&$7%rbd{#yAF4czjUQIDa)4AKtUvi{EkcTYoLAhC1~rFjVq8ZF>%3J0mWbZTx@8nrEUp8QfzmCz>asqU_>uGW*> zpG^bGtn<-M|3(d-&k8~ktS3ssK(OCRq4ln^2N|7R1g_t5X& zY_u;|*$L}50Ly61hdk-_4g34bb9Vf@TGT>_@ORh^lJoLM4SiQAD3?fh~ z@G*UA5WmT3Xgr1wdKR}~VjXiJ=~Dxx(sO~6iz@{9)247@rqO_-?K`liq5S{-ct;r6 zlN)K0^u)BhVQSe5&d&U>QVGC4Fqh_g8VrU}(a|4)ftX$CWo2iXn_I??ijqhKTc95} zJkbREDfVn61oauz--a47Ac7ih9E&)OKZMD&|4~gXb8;xi?7yXfx}9KxK7pC;0kEQG zK<8HiG6q(DF<6kD{ay>REa(_QF4LGB>C#MOl52?NmuFbu8-%1#ck_UVl(hBgYHcW3 zrKS>;ndOR%>TYRkhtuBkqX06b5i%=i1M85KYW6QjwI`!SgR+k(($~FAa|ns<6S$?+ zIg>t%EN9 zGU(pwzCan>=^&C)*wx0$Zu=0OZ$ox)6B96dA-iL4WInqqu`pOvun1u=H$7}Hn8K!i zo(mEu#M#=-uUEZ5G4IJ%XOa9f|9YG^Ch8wyllS9O9{u9`zcmXCYiQ);P3Mbr=V_`>yK`(c5lkuqYRI%d&vuXA9cN%>)=txb(OwEip>DHIL*=G{j( zZ&5s;l8P8-TDXO7nG;5^QsW!19>v7PDS?<_*_}p;@ISxt?e~4%QfmsZsg=a62(qzZ zL4Z}cc5U7GC%UzlvwL@u^7H?|2NYkJrTl{+P;VcfjuO0ORw`hew=ofXj0hr;`qice zwxEYcx&9*yR2X>9*zF*4R#jCMHZ^`UW_rL3Mp$F;LYS<#w|M63*KA2XZ;RetCusDE z#d${Fw48DIAJdE5qBImG?Vdk>{bJ|f_;l`I%Tkk3jSMc-s)VHHQVfJh`3(B0XI zLIQB*_kd%6M@kNzAJN6F&PMy;_jrSbo`#9b;QnrXA#BGQ)>!MkxdyT5RcR?za|CI# z?`mBnko#9hpb(%rb?Ow#&mxmRO@P(V{(Uz9N@K`=g?XLULh||jV81RVKzl0m3PI${ zCJSOhqbE3)6n)L3c>3eeW)~Zp&@c{4!ykU zUJYL3g-eylGkV25}rJn{g?;e8s(>6Gfun*8ITX$w0t9KguwvM%M}#S53JAr zfo<&PC#m#^IkvdCSf@-~RW$_6y^;Cocx*^5(;y$-zvs<80vIZ(bqc*Dr$zn`L7ksJ z*BbsU#njhS5d&DR14Z{%z0iJ(bC)lZfSS?_15~GhR$K|)@Rmcs&*W#%5+i_S(D|-E z+R-1o-+R*w%8c57z0G>?ZsRrp4JNELa0q?^I2{XtA7mOV81zl0Ykx8A$3|aoKslT^g zWQa4l(4X(JGA52D4<$&8V9bRI9d5P~$8MlN4$6Grr$EJFW)nx;35Sv(wN2+2)i*(N zO!SW$PqU$7lgjM@eU}OU-3A08gQw~$|i+A$4$m&^9R z%eu7uc;w;Vrk;_eIc(!QO@CLKLk}~36L)7L*t}6kcPol>p>~4Ils&k+@?$af{ol0{ zatC|z+#|`5t*!g0XBM?w0+guLz(*%meFG*o@MQI53W}wV;B2jP@j{0-{a2DNmgnLC zDU|+aBRSBOCq}q##cSNFLq801KaB{wz&5D59 zow;d)>zfDm$9ufe4VGyJOCM^BPu{lK$L0Hr_Xj3)w^a?iod}OfPC0p}FJ{?{-xWA` zb7E$}KyP#3;IKO>Ur$dux%}`h#;vGEEV=yXu1i<{ezT%K?Lm!vTB`+Ljv&9bw$~vY zJP8V4H#Rnkq57D=U)}iv_r3d~PAveBp$}GMs3{%V7R79(ISKDaPwRX)U;G|W4W1`WcXF0G&I zW?mN%$BBBTs~$^@PzHMsALIe&2){)G+d}6#?roJgNXuvR5L9wwAcdee;i8f6vhbMK zOsp965q%hBTLGkmq6s7aZSApugXKMd&o1yBXxRu({h{0H7cSs^u13^5xMA z(F1hs#p(!ic9B^-DygsH*IpdaW#w`QO zRy5$G4LS}Y%wo=ni2h=-f;ixs%QT%9OoCUgPykf<{rh*EgXIT$_(ONa6!T&Ath8Tc zuPlRvR}HC+v-70O)cg!9X8{wqeCMg zCox49yhk zsz*Z{=#ha#4%#}c&w zr%hqGqu^J=&@dK6%%&3Zt+zK*UJf8%Iy!>2d@=;@w?t9=RSAtjP$3+Bt zC!ADY)5vICEr}fdasQnMDQd;w87IXMx^#FFv%yVfUhjQ}b?bGHG4!8wjo;pd8oABv zk<>ui1-OD@|0P@k-hLzqphu;jl$aRZ@_jo~xQ;!&|6QW&O2bJ2R4JDC1{b$c!~JLp z`V#xZIf^CMO~QZ6WkIiAu{tHA*(VhP$mDbx2?_e)-bOV>IJ!lh^^=zN5R1Kx>EXc{3#{7^Xk_kWJw zpx%5p8_ww@!@_CoNW8e-e(sSmj8A-rIkfI70o_*OxoX4+J#*I%78Y#3yc9nJeQx38 z*kE@;wLcn2dv0@I*#rdMcla{wB2SZ&*rcVIVWRotg3M6JPUP zV0@$5%`-541~3=#KOMlT*RQVwNS%TBXIE0!&>#ak3iPC{f$(4ojaAG?e|N!1sQMOy zwn~x9zq(Uo1Thyk&38o{HQ15K9P%UR8fx$CbSv2QJlu0Z*;s(f3N|*E&<H_j;ypkd+Bf|ot7PkXD9SYT;xWNw$335D!^nt{~OSAg?4KHEarD-JWP`y2vBEYX?XoLc;fclH|B{%^&8lK!xgcY+}Z!M zjsBPBa-YVHpPfKn5&{{+rVs&YalnShwLV$kL*G>|a9Im6z@kC>cmTgVtLV&(QMrCS z6vHm3s`{YfdWI69uQiwlb|^M%%NRoONt=N0t~I0)HZIX!LnA%R80a!8s;W%@|4lv& z*BZW1AFX!hcC7kIgO>cz&O`9-T9 zSjo6Byl>ka?Kx6wB?uCGID~2-HK>6%Vh15C8Wt-9^ogNjLq+Fn3BaKdpnU@u z%$C>z5w*nlBlZ|H)sOSWzR*ZNYhdwUX=$l*Tn6VA0VspnbMjFfqd8{(q4B2+zULe; zB`7EYmlJN5=I>?K9C~mX!2WFoeiaJ|B_R6*hN=?NRwDey)8-&)?9C>s7d!rt*jgMU zgUS{tB2+7UtCqU3ype;K_L+(#$dQVf`SkM}18*vdpt`>na&gsnAvYT>4p9;n8t~EQ z1r?i&h?w{ah{J6n-~$HX(RBE|1u#(%f{hB|5mxAq1gg({WB>?8KS0bzc~J)5P`Ll? zKZM>eJcNXf9~mAVhCOPAEWO(DkwSdcjkUC%nAOw>;1yxrKUZZC2fLhAO6rG=FKie{ zS5sMeIai?F$8X~r+8&xB>ZD|9nptt3o%VktT_jEwg+lAJd7{M=nDsbF$9&qDcUu3* z78mbvB)ZT)JAL^txE;TA`&pMI2s4!zZWj}L>jE&1Hjk1pN&G!)0tkDc+CxDc@Tw8Y z6LkA^167Ic1?*oCVu*d+rjy?ft>xZT5_#D@)-Mc_>B-51przWVXgFAM^NS#$Wn3Cw zgDBXeXHnl6d|HX?Pju&#j=6U7=oTAeO;1mMgKk04A;_*jz5Edh=xKB}|6ej!d?%A3 zW6+U0erA1rU7j~N2KFm#gQ?}mIamlU3%J*JJZ$b$plSGT48V#h4p)d*Dn^BWVd0ZSA8)l8=WV$k^1 z8G*Oc2QdWR?RMGw3>m|Df^PF4=5=g-=YZ(lAaL zV~B&w3W{*j1P#!6g|=4%kNsxn<`nY#eadqr_V@R7tMA*?FvdhxpO$f=mBsWWL7JG~ z{^M~q&8Wz!qWBudkiZv!yBHRJSd{-ix8(kq0USWFB%_+DK>~t=y!FZYTK5F?civOScaEs?4cY-CNFyI7ItBJt#V&DT%^$q<9_+aEB zdNmR~5n%lo2vg&pOhvEXS;;6iM3*9jsO7n9nEG&qMj+h*$HaI31Ul1!VWbQE$RIXG z!&w&y0t7Q-Q6Z?FXkTxobB5U{YU>083ZqN_BdCyGaG*g6kkdC(C8!bSR_Z&(eK{J* z-a44sMO62ZoCl!s32B*(r&ag_*?ScJzklS{Zri6%rch!7uZ5H5zjGH=w~-)N=8V#Qbg)cNwyyUadNhj9 zR*-5|6X2c6cemNtu^WNIx~~E6*n_Iu<p^U+Gb#>J{GJ z%1oo6BL<_c3pg?3=E%jrQ-ud>k2yCT3ap&Eeear}PEg}Z!0dze%BzjG0YC&Msm8Ax zfHlA4m{>M?K1*Akc0SlNhyw#pMEXtt3w`GMTYnhpgjEXpSJ0~tZlicrDM8)QTpVK; zmn(35!@mfo$iM8;B~DJ~;c?TFSAcY2|2yPe_*cVnuP<7dGIVxzDc-p8PPc<1x3?EPrzm-qUfRzq zaVJH@p$V|wsNWF_IvZBsb8%RaH*eoguLho6mHsn-j#0{o3T)&gu<`iKK9L|h8N-px zKJ7gL5I^7$7bS3`pa>fL@r2+^Xo3@9td!4!j0zGFP*zR$BIo${P@Dwy6Q(8|{gQ+l z_dgcK$btMf^w0yeJ!yCKQN_3qG0^V#$c+>cQuH7I;C>jF1qB4wz6f?tLd1`Xi~9u8 zf+?5+fio^*TUuH;+5kyB>4{*Z++@8mT3@;p)0^p)6VD{7|LKIcPr~7*P1Do185bM=UKZr(w;)NmY}IqeN_| zo_O^KS1e3sRPD@rnvej1VbEg_ZF^8t{QAiVt%JZByy?P%#tve5#zMzH$SF|pRBSt^ z4y62UBQGM3a=q{&K#HR0jx9wBOQevEdjqP zsHj9;J>7dr%Il^iiH-~vgPrqN>n8C2-k>>4(qG1}H8{Oj^?+uArC`T%i$TeDX{uQI zveH0!D>NDq$>1x(NZB$~T*j3$ZbVRq7wswlJ`!b1;QSG46f>qv`xk?2!?SNBSEwTD zH97VtTVkTl<>2DtT6}2BrW!ahvKE>$_+e)mQ`p}|tZ^_%D#)GD%gM0CWkeWD!GL#Y zVXzNuXU;jhN8zu2e8+hPnFDSg8$r*mxdG>-fQg|6%m`Ot-FPx)b^>nsmUon)eegO` zt=Il48!kLZT)ggtuh%;{T}dmUjGQst5B(9}^Yox!a=Id9gXR$n4Xg4Sv&t$Lr%b8K z5d;F9Fum>;j;DfLs~Ev{g9M?(rx&P{2m)}5eoo=~TgqTuQ4i0cFQe^MCFNYfs zy`n8?)McTv1J86wE5|JZ6D8>HR_3LoYQ~|y$x$Sf_KlJzC+G&3%Si_+f~8bmKzQ}b zNzj8{jzQl-PtJ*}@kbB-g5(%-^|V~_HLf~xs)Nk^#Z%Bp@R&=Gh-8O<$a$}-cx;v~ zgg1}>qwN_^t_m^*7U@wx^`cUEp0wKU{6b3jm&j=vUjgv!W9GG^YzKnZ8P3G{I59Di zD(rr0*##?G6ngV6B>f^Dyf8p2oX%b3^2Lkyaz&(yKKDsj#~Vrv?nu$ihxXrM=J)Q> z(llew#-o#Yhs?0c ziPD)cmy zF^^FPRJ=<7Z}wY(t1S49*P%5VM%7VZ)?nyWVBPRh$7erbN>us?q_`0lz`4{NesgEg0*YplzWD)tx?4(X)2W?|9FP7|4|5)Z({N z|I%5#&6~)c>3H|$pSxcZ)Xv;p8ZG5QH`PY18ZB}bRgge$iitfdIt}_A96RyX&Zo|p zr(lX~@qTEH$RtXf@o8bF6-*=Uyvdnc8}@4cqJzLp)zncc#DaMtRq#vyy~QpGcf`IO z&J+Qz1>(=zR)z!#e&g4V_9eo)d?jx?nnt$=xu5Fvm;0m>$=D=zaFa&9Jh$jWCrpIi zj%}X8#Hm?G^ETR7!5IP0yU*E8Ebh$ER`SoG0~JH>-VQA2(zJQ+RUP{%`dSNV4K|sV zuw}m?f(9|3mD|7h+_xA+UesL7aDZGeE$6Zb#^AuRBwM9KcvVIk|7~OC<&QDnhd$QQWQmp3?au*DV4}f$_-s1 z8ImURJe^~TNJ5h-AqpkKrAdYIK3n(Rb^q`G{jc@D?^>;u@i^yqetYll{yv}Q*&ef1 z*maY}3Dy}MAVU-iihz5N|KC5ikF@fQ(Mt+`b})3_$#CEglPvWEqPmZ#$_rJK-@9^$ z^UH7`G$e8J!I0Ls_%3HuQ(|bMBImY^eewG1oOZ90{NU+MUo?p>!@}bBnYPulR8Gyc2}Vav8Xd&6!9LN9 z@pN5RqeLt&-l_O|sbw~DlM*os2cz3OwyX-Ck&w_%fr=v&zDGme)n}%PaSg_RhdTu> z3tAk_NGqs6yJ*dK5tr5KiSNgT9kRM_e<$S+>AB>%Dly5IK7jltT>e$<_XR9AR z)0Q!t(37*P^RSSX&OMBwwo7FWY-Q+J5dZ@f``a?byn4GDfUfzDEu<}Iee!w|*E#?6 z!lfvRstodLn7c1JXNcbPQJSTmr?`56p_jD%sf}q< zNO_SSqA$KgpP3xP@_VOs;Z}w}Q`7@vHYQ!q?AA@mqq=tY8s$T_^EG#Lu-4~Myn5Ho zDp2q2U&%9C8p0~quV}w}RqCT#twzjD{^)g)j7s=R1G>h1ytAR)mW?XuwKby7Wyu`E z)c^F?T<4M`9%cu#VvT*LHrZ_pirT+rnU!CXUA9@j&#mQe*S37kPN@j7@5Zz-x-_YC z?_Npm-Al-_`9T7U3u%uLJ{&~*Yb589T?h&mE(sRI)%D^%8hK6TeT1bf{+D_G<-P5t z6PX7U?qD4Y- zL!QX~8R4zDYRq%5QWyuEStCQt^Hfq;Zm>ggSm1WDB=KQdg1g8`K0T zV9H>kv{6^TYI}q9zHJjJqer@rK6lE$U>UP5=n2+UQZL^HqgEIrWw=Z|-&e>oQ?JN;VWs(bap{SV zMtUz_#fojMw2o{^*l+{cnZupAE8)>TI+|T-2SYKsW2@xkEUgX3JP2w)x&rxjA@&DQ zBPpRx3;s#Xc?;NVvN9P&ig~&wthZk@HPIIdZJR*%QM>sv_}2$jd$tv6ece36h{e@MCJkg{g-ii<5@mcXVg6>T_H z+RR~!$PnCM^5!dWrITcG+qg9aiS;aRPChe3|5n$9X#E(Cu4XF>xTEzJuI1&++Fr`= z7vi(yKA$f^8-@ebLnc-I*sm07!J+dv;wHdaY(xJJm{TOAwE>ttJsKpaa{mU?f7nDk z_}vHO1<;a^l9S8mQg0dOqx|WlS5IEd>&fq{?A)ss@KB)JEPmZ;r8e=t-9Z%vG6F^1E<&v-toaA5a_v%T-7*7^ z_T$mz_t8OlIdF(bY?mZgt*R4ME4KX8@SjYth2ekG>Cb@Do2CGE!~Notzw`JE3+3nD zFc#o1O-Ck1oRK3qkrDhvCbxUucjzsY#c^(`FH>Y&{zF&&V@&+=scqbrGwt_#*>{iRMfhmoX}(#txMNopIJ6u(oL;$ zA_>00q7qBoU{6L*-Lt3RG`>dAdd8BYARf@v4Q7257&iQD$E~F^1uDYQiyG9vUz?Z? zDIWFad0fVJ;l41(!wedtzQT3Aa&ixJte#B_jagxA`VM#fe*Z+KOKE%_xLeTuB~iKq zAjRsjHuCWB^nYT=I@Q2{V)yX!8W_h`hjOosh*`T&QdYJhKc^3BX;fU?IaGE+{%7w` z6gV`I84<)SX)PeFFQ#?urA48mE$k`;3nkufA87br`{sU}gTU)Gri~3>`m5R(EjumO z^N>U3N>Gw)hMsdwMLD~P)c7R{ zR2WznEPc*4Us;dmVobuJ7+TTw7ceI7*bb>g^ktgZWYLGM%x^B;y7H#0auEH8*NMEP z>{}iklIeTq%q%`nQ&v4fqym+`*vXLkMLc<3W&Uk{LrGKSj%vk^ zysCZqjE<4A?=PzomE&d>b46_svNKoPuiE`GJk9Kz_oV!Jtve_{I~1$2*?XEX1^W6y zbF&^iC6Aj~db?5*v*~ z4N=G~)9D8)5onq1ZsuLs)zhJiJv;aB*dlVk54k%aV1123HzF{SCQ_f}vF2w9OUF1~ z?kv)dwir@$$=f;=l6=n{j@HmL9eTbY1%_$D zS{DOV;?rBvPUI#;xk;a85N-Pe~<)|H8AYdKVjTq$C ze%J=DH4?z511P*^9t%;$3^G3)IjZ=ELn8l2iMJ)clXB}V0>?{!;K3FKZrR!e{AJ4O zWsP_gye}(0i>R|*aA;qhnGtgR#Se|&kDTyIq|lp(Y|MSlmKNnPDM+!(+4luXIO8<(VC#r;MF_$}92S!nz!F zEPchb*RhntgdhfAYO=rLEkVbZ#y2YE@zV8y!Y@}!DZ9C&UESP{fyx!$E%}2~Ha2ea8>-H4| z6?=U1=7ed@QsaCfdDQId&kzB`mjO7g(<&hV-9l{xO-tJ7ty+u~sMKru6*UMA{* z>BXMkX}E_BDOOPhPvsqpTkH65ZE;dQB@@R_Ie6xs^jKf@N?;!P78c71$Vkcqc!&I{ z3amEORRQ=L6*tsW;vX+}NsOf*vz)q}CF~f^g_Od<6|}ta`1s~m0T(!2B-9zb)}=Q5 z%*!f9vZUx)At?g`1IRa@U4XfDA%#GH8~S`yiY(U=rybHJ(I=ZZ#4RouC<9ybT@3de zRn0!f+^6GT(wol0R_w{Qr;W=|`5@Z`_{nc3Hulis+7}QCjtEtA-yZv2tuTkPcG%U)=pCZE)j0@aK*cbq)nb?9VMJ|Q zo4p-V7E&vh#U>zmt9i&5Q)+~Ro3P?>KE4c{sAqjaPVcPU-!Va}%lq#j8$6HR^mD@} zo)g8hO%s!~rP(M3O<(l>>g=9F^kH9B!z~SYA3_$(35Le&w*APy+SucKhVo^+yYo;2 zoV19l(%>UWv0ADa)K4gyXtD_j@r@N+yv`_ToIV%v>P-Rd^4_Z^*7jtLMy$qXBgbOu z8zD>guZ5eHwzHJ<897E;F90DeK^l?hMo+a#XvxnTnZHFK>C~wS84oW$3D;FTqv0hg zFu3VP7lv)5PZ$RHZ8>|kVS>rk7$b;tQIIV6Hgu#1J!$ln3+g8O6(seDn3TQlF5QiTmH#$OQS1FHKOOl(UKi0y8J_# zggLI)zB@FE;0(0!yp`vY4}_hIJyRl9%E`2$f{%>?ODz6mKeyl41)%a4T)gio{Kvze zopHCP&N2aD2I=-2@FRL{OSBp`nczZ&J#Zt!Bql4PJ4^gw02%=C3Xr7zsQCEv@Lm>u zO$KXY-}t*UL2piuzgUP}6&VggLxKWsB!-P_UtfnSo|#@op4f{d1Hr(ln<#b3nFZuD z2(q#Jd3*!b&jB>jbkPdBT4BHz@E=i_MZoNyV~@vPa8}yC2dJ z0!Kw(aZ_^WZ-il3cI~U`9rDHtCkgTdo+(+gj2-6a+g)D>Eb6HyR>>D8i$HR5!oN*D zc<>9bW7Qck3S30S-JV*auzo9Mih;<9epl;eJu=}BkwVd>8h>=~z(`?`EUWWID#d3g zW*-Uq66hvACjakmS}rb@pt>JPdq0}(^|&4zM|##qb@i!&sd+W%9k$K{uIoB=j>Dh0 zbqS;(&<#ZriT4Z8W}-%SBK!C1r^)u#(a)5zO#u+2BF5oOXsNnW{8tFrM@+bN4GmH# zgY6HBS#ACs)B9lAvgEhEJq-;LX~+A&v}wZP21}&`AR&lPN9V<-SFc_Px@`cBjPDWR z4cvnor5~m~B&OaC9Z#VD5Ng}{1YvQU4L1S&;n0+dRWBmrQh7yJEPx-!cHL> zN%;Dk8g&2Wo4leVTJk+|{fnO8e$Oe4Nr)sCp^Y_UW$!HRaY_5T2QB-?9Kr{;m7Ysh=!w`R6g_^_)`g#4nn1Y znton-WMxj}lqE48q3QPir1M@`+zt}wJRnjhHzW?IJzQ>N?*Ewb;+lb zYubLnnup%dVgQk(K|p+ z5e%_;nWd-S+rgvE0)yDS9INt_&e3Ba{!%L{W(y))Ebh9)ilGEYvON(`+|Evas==ir zm#qF_nx?+KJ}xHuYoYcn1|u4XZHPja9m!XO_I-~yGGX<^74)<|{c0fz#uRxs+z&>W z(0DliLbIe7k%YJ{;Q`GSTzlgTBHhuM(UjSh%z5|aj3B)6L{5jwO@U1?rS)AN=tR=h zjq}t-RPd#rNroaEdpaR&kSN1>fNw;Q{mbAc!Sr#*M~Ua5BP`$0b2JE%X@dVDy#Vf1 z;pGqueZRh9C7cR{^b8Sx%>~)zlZQiFVpf(s<}eKhcCe;|XP4e{Vk9JGD2lp7s1hv!ZNwX1%ju+n8?PuMO#;hT8sP`PFzZg2*iOC!iS6y#Jyuz>|(oUXR zH*QRwi~o{QfZ+M-j<{n95rbKQg@QJ|5MdcE@UIo3Y>a~cX8`9M)j;eeoEWiGVE({8 z4TL=dM4`ZGxjjl0hRMI=CoWE|n7}~?wIEINffqq|K8@G=k30gAZtwtPfH*Y8q9wLX zCfX}shuJ?ZZS^OVB+jG_I0-G;D_zoKi2gSD9*1IuKdg*`6mseS{K!)4)?H-E3wrjC zDfU2M^u|X7U0n8C3f?v115i_4pi>trJI_WElh9wcFiZ==)QxER$g75-tgx$cZ`1K z7`%!#g50s;*PT1-ffy;B7~ajfa9t%TJGb#H?sim+wU zt@iz%2X2In)JW|QVWq!*bpAj3U8-k;I6TS`_pl3#90#4MfVR6nhHmMvsf)nBAcqHW z34q;}t87PTKvTEFZmEN3KgXa=Uyben&+g*xmB2Ri!#PhC0x38jXT?Xv@3N4keofd z)Z-rVk|SEP6;>1+7LfY%Fy_DR@zZB2t9ol?cs)hOnBXt3LE#()7zj00_xheQ2-q^O z-m-Dy3c{4Z9D$aaKeiskb9k)67j%PBq6RMr1XMhvT)Yq-E{blVI^hN2nMA(IRw`9Q zKz@J!qdlu6S0e09uPqP{y@H%N2GI^rF%qE%b@%U=On=965ZrjE7B(Pq3zLVy15?KC zl>)M?ArmZgE>cO7BS^>{)G+fFk*0)nMt7zb1`Gn^0+=qnCd>Ox09pla@mA2J3UN4B zsMPuJ2Bd+&tS5^A5DEosx<^F`cwREy!!i}&wi3sB!0&Yo7(cO)62uVgKqcfyQttBa zqXt06w|Zm2dnuP>bbm-86f)X|8f2eHq9NxWd5H)cg}4xXXWttbMB3=$^e{ba0T9Ur ztaf@zAt8*RL7tb_lWRXd$>LBZ_C2?aYH*jV$M7QOs39ruOD{N9bgNFBnFUIR2uv)j zF>#r)sW}(0ksX%AbdY84lF+yE&+uaqQ^ZL39%m>qQ*fJ;FcuVpDX5Q(BmMm;Xk9MA z2^0{dg366U!q`Y4=Tm)Yl#8g?>QPUbP^p{=Egd6YZr{Awk8|D(KyZ`D86tB4kfV(8 z=|P_-lQJ>MkSU$KW&<%V7Vz9CJRpwW2X7Ziny7fe^D82IW5K`N3+qghB`(lw6#E`v6JY@SFQVxNPh}Zko>c zC6n%i_RLrz>3&?mulB=?vCY=@y1}F9FrRaJY!Ve#EV#F$a2g)Y+A{A`tFt1hSorv4 z`F9asPhT=h+<(Wnan~V32v}0C)gMogio_OP&wR_Y5jSWyCPQ>}SD`>DE-t2wl?J;I zbO{b-qYi}zsnnc+O?VLE2an!Lb|>UGzCY`+|L3?S2+({3>OW2jSLV3GoK?EEF{@kP QSVPJVePg{G9h<=a0(3Y8y8r+H literal 39597 zcmce;WmuJ66!&>(5TrY$l@4hTB_u_V?hfe^kVaZcq@`3qx_F>=KQ6nQ%b)iNo+}ee)@P{a#t~o5LFKTSqbHDRv?Jr-y}C% zF!B+d2%KB+YhsEysv`0U*69dh_+6js|M|~0AN|kW6$1;5s80u1B`>S0BCKq`78m#2M%dcgHZ(TA^E=(2sWO+%(2JIP zk!Bi5#i>(vUr`ZTi-T;g-X-Bnrbt?JtNX@{hY9pIA3uJ4J4HDxENmsyrLFi$_8(uV zw(GwRp4(qhyu7?z-P|73J1|<}w@j)&7 z<=3z3qN1W5@9${f};%M}L7t8++&$n`!MMWv`@$t8pS|#uq?ity% zF9)2`8GN!sDR>eMuP!rcBqm`~y}#ZUr_xJCMz)&t_}+t3o$|L@9AiZq%&J+UM8(C$ zr~V9Vs_6(UEG(lIA0h(jOTM3bd#mf~p^v4nND;1Qhvt8N&XSy*oY;-M#|ZqpQ)m#h z7jSVRzuXouRPOM|r@eh3JLV(d~AY(nyTr156bka zWqvL+YF_>xm0pVzSt5K=tabTk{BeN=N)#+)43BYBl#s(z6t`hrR9nEsjYcD$3gc#* zSPAd_ELSfta;Nz^GU#d@mnE&UE{>B02^J)w-v6Hl?WR+p79{iW=Kd^7%R^lI5U( zt%tvVle~NPj!9a&{oASe<=K(UJN2knuLzMtehHV4jwWHhfmmDgUMAx<_(Ujuc^LZN zz+W-Sq9fdX&Qu!&>;?S&Lj>P<>MMkv@6iZ55+Zu)oaU9C4!Xi|QG?HZlvc)CC$T>! z7Zw&?JDxVB0OL&1YjiKLn=DB={Pm0Y^XJb6KejhF*-%kYd!|i;xZ!rSf4IDe`hiJu zpGmM-8P!LmT!-^6b)3)1tgU>1Je`UZ=6#>%&!d0-bhOlWcmMcGxxiwn<$2%7NQOYY zvx|$Eva+%#qe15zY2oGNig(pQ)ED>g)w^z~Vu*3Yp#1q?(r7`3wL`q|O?kfdqeSk%KfY8=E`Pc%v2 zuRni08IwNu#8L^`@73aJ9HUp2R2t`vVe|6x^85bTe&XiRdYr0$o0N&0`w=!{&pX*j zbi`56ztfI8nQBXg@B<8@$$MTXB>D}t4%19bOsKGNUkz$8?0RLLZ?({36p?6 z7JP<#vzE%`)#)-KIL*X&C+B$o>pL`5G_<%=Cnu+|_wu-w#>)N|d-d43xS=q&N5^|J z)p-lkrvuCh_t}!D4(c7I@e#4}EuK5XualCFz9&6S^))dwy9Gm2PH|dlq2l4;p%nMD zD?Wug!#_MZ>DymyJ~*dlCB}gPIJWSn4ULuQ3!a^wy|wI*?{gOssejV{PT0innbvXb z(1ttA8|Lm*xvY*3IpTidxOwkOGD=}b%7w+nxcMYS7D_=I)#AP9etvo2ORSevbaizl z4?EB%D@2Gg@WGmJ>zo(GCC|Snf#G-awqL%0L9Bw~NWc4(_y#WSP*0lwt!1TRt@m%= zA`q`~a`4M_V!&flgY6bLvqq~L7#M^fD?NG?uF4i#Uw5e2=hl|37Gb?U*SqD29Se@u7o;E=GZan(vOqPMiP z3_2cwn}qvhV92$f&WKV~ao)c{sL>mEzCnSGg=KmCWAc?EOP_nKfpOkmV^fn1Tar0! zTC%z|3@j|&vrDkc2ZiH|Vq&u_3~TL|zi|Ja?8WrPP^58;e*RoMcB1?`KU2gd^*NWA zSmQ(h0RcfnbMty`(DfZ}L+`r$DNN0Fua2*S2H&1|7dDer{b)IX@Zh=g4Pz*kyBn9c zCD?i>wR3u!0u2-IcC8ddI9L=al5vSeW3R1c0ip~nHxgXNdzetz3X{IqZy$v6xhzR?nY2(c1pLv6`{BwSE%$^?Of=x{U3@z6u;Ace zdzVUhJETRu5pv+O?8CEsZe&D@PcMmWYHDgZ3m&TqT8h^oL-j zsHk{*{7OwOH5%tuV#eM=BT032we=vpK9*57!fdhWS-SW$j)(EI(Pw8q{+CDnhet=9 zDV)0YpW4Cdh{5B_RhqOOuEp@ahk!ls`E%!Fsc!YtaW;rdHC+?l6#GeQc zv$WAI4IKfH@!?HpJE2H3zIJsc#MIo%KFi;gEA1K-NP5 zRtstATa`2(AzvvuEk70+Z1a6(LZZnJ!^5hd?8Z>R%k+OMQ2Am^$HY`)iF=Dk=r*{M zbU|Bu*gja^-rh#-fpn(+M;Pwz?tD&jcRh3yq6V1wP*~n$$`@-fLC&XHs3!mLAtu;V z)3ex{{=RRpL!}-0cHqm_sPeW*8^B*p93QCCM)DjHO=5nSf$9^{7VbO5DTwFAX z?xR1!KB4B7Q6ldlpP13xMJHY?#2wlN!KH(k&L=k-N=PGe z;Yabypt<>E4sZ2=#L6+o()*!h{&B0eA(*om-j^mm7pKh()7$8ud_!U-Wtj}f)au;P zMcnUKHQ(r@$G7|WPNphSirR0c`}5R?xysCUX}w*P&a5!7n}a z{O;pe$)5a8Mm8ye>Zi}6l!jzg4FcULJc+)!$GHCdTP`qq?DbKIE#gru&iflE&K$Ql zRPl>fUtu$bHF9b!gQaxO&Qj^?>t|u8%01(e(v3^%jb1|DdH^%-xqN})DGTx2#==m> zSx&XHH$I~G4ML}yBUlR6LuWquTUtH7|Quphyu{zNV ztvgn(+Eh(Np$ppB?)3_9GhZ5yr^*$go$R5#evM#Yc&p2MU*}(Z!X#_3-^tHT72ejc zk!%Us(n*WW;S5SXGsH-)G=ZO=AJP$(W$6_5zuWetjaTcB-pvqw#$2PiOOLqw+w;EO zTNSt69L_gF7Q0JV?%eddl;qn6B}UjG6GlQ71x;J+JAO0=Yieq)LZ%FH`mOiDlD02CxWj?uFhWtH zwCFeuoEd*{oE`5kp3k%f#~8mB_*~}Tieol(>Gu*`{gtbvxIOVBHrg_TUuhy$iAB=` zFG;_ma0~4FAZO?RwgqhFxBNfr!2%-lCCbw(kTF@qS^@NAbHDHWJ!8I(sD1o*9f zd+qS28sk+|2`6I>B2$gkLRO23hsBy^Nc#pAz4Kn2g2dT-TS0=EfwNWB)6~+bHD6!) zwQs+PGIJ6dV%J#sVXx1IZ)|b~o6}nKBXm8T)*^U1ow0E}!^`@vG1bT{ca#{%ONF<_xOcGn}`CsFf|hZcJ60 zCM-64=fMq6E@85V5NbTXD=u|ns?;4(gzhL6R5S5&QJsMu?_fKIanFC%@Jp#V;UpU) z`kKHmHx@QDh4+GlRGB6=jWp5db{l(P6`ju8)ITFeI-L())!!8N*}NS`iK9Wo!P!Va zOIF0I39+x+iD)=~Zja){S^Id&SIhDOZ zq_PbPT0v5WuQ=|^58{rOAG&!B>+nEVD*7f2F>S2(yQ5>2a!+589E;VGGxw3Jx|yPc zspMQI1s$4E13H!vcCnbfT%Xk~0X8}X3jtStu|ic<7=&@?%_QX}A-Bl(QsdC`3*L$C zSFg)4r~ZX>%pRneVDa5T=v8&(?;-?HKB+@aQ>IyA5F@0pUcIU3UrcZxtc0Ncq^}CU zZF_8j&V6L10Fy_sy-1Yt?kD0LmEjlH6^=4@(%vWL5$D7YuI1L%2YEkvpBTz5nVsxf zubaLX^bkAbrwq5V4F%JkQbLg$Jn5tw+{7&H;O~w?uWN3^a$j6mc+LEKm}}5ZjBSfT zLZlxYHjK(8O&F@3cr>6vL%yLskiR|uktpKXq3l1x$2{*I==91LbOcW=QGT7H30|ch z>aSexN?gl7{1c~IRO`FYpdCz8&(ZFmBZt~B&-rzG1pR0ve@rR*A$~-DwfOGsAAeAH z8tqxqSiEP$6-PMQte3C#mr}e2^O9+JR@w8PD5~;}buY}I>Y;%pf5g`~$~X%&y?cEx z6VKS~_rs)zC+U<#u?q9&2B&3;@0E*cBY4J>vt;r&ixOYw314#E84|k8NiF3Wch{rx z3&+b*tF!aeeSbSm`BhfILwBzb{?96(HPYS$m{LJX;3F~HH9vcDGAa?|{V_ad@WV*c zr&&lNwE8`LL&~>aA6;Ej(zXW<>Cul3(~X|sJzOTwY1!l%le)bjfZk9{xbmUou390* zzyPO=vf4X-C9cMT@1#tuhEA&FxOc_BlkE^DrADTI`k3@&kU0erN1Ni!>g{#trPDS4 zc|^}9&cRgxveg&wa2>u>)93~A^);YJfTK{DzgV5k>BK&^c92di{+v&2>0q*}yat|B z+d-wCZF9ldd}2qwk=rMZHkmfb5CxQ6t4ZgD!T!fsmMU=v9@?{Js`a;G_sG#77h<-3 z!qAYqm8b}T6-~jA=I+z)X8S+xVuPg(aAe0ehaa5wujESfoeDWwCD{52cSgkr>ddmr z$xVeqPNRRbTrHt|%VjOnz3)_Nx(4U(C%c}b068nNA9qnaJv>TvTA(@sNsf?|R7*vO ziOJ!WW>M#lyd$`osRovC8MfX$I)DBVfgVA+geu6laQN8H_?fk?kl`~Vdx%0dv-FPPCa zd0vp}AS-us47umcgYd5JC0Q9+mRb2PTbc(NeaeU`Gs$KP_HAu5I1&@Ht`veL&3Szk zuzV2H2U&>a^HCdmh&z0L?{1ylV|q`Vg)Sp^yi}ob=gjK&b0ZZYF8tD0N=RqTc8vEh zilG0Lh-mRk8IGBZ3@SdYXmE0JvOsxgW!vA9410OG^R`&Nze*E<%jt+x))!x=IAtHBdk2Q_&c~D4Em8itW&qT>9 z{S{BGuH_yt(O4_A2__rxYIe}k(ai>v@A#Z9dP6BuIrM3wSN=_O5f|jwrue4)FZOg_ zX}?85NRs{t&j^+MI@7&0T+7TP5Z5#yI6Q@#o|~Icm>bx%#2O@LMc>TxVP(C2gt^Hl z%aqea!Jm(Xl<14(u#q6s`@j8;ILMqW21L8HZ0{w0U9zYm`UbzO+~!u)N&>s+LUFA7 zefq6(azYk$JZ9Nyk)SHP1x*Ty%<#BJJ8lvn!ZeDNXtVxAWh)7w-BOU7DtX|IGMX(x z@$dSA*PzZ32Vu_}-z`v1j?J#>D^wJx0Ru94m zX9Dij4vV(`!^4p$i2n^YbO$C=xd7P@*_@Li(tpJ%;&FInN>P=zei;?@BW#wqk1CS#0D(x>v9Udn_aO!BPtJ=?!z{0!?M`JB+sJ9nL)_A# zZHmX3=VfkI#|{4PVy|b3?^94vNS@4CBHt?mrT-kdFqkzaH%79c6b}1$ecG<5tSkc= zZ?#z`>YpEzLuJpH><0#Ewod!HxsR+wkuR7Of)Q_i@;$a!&AmqlO#&G?IrLO6{TB~H zF;}6U$Hv3E3xzFpa%>Be|MO|lb|s=tpJ?O&hCirLyub&+Af8TQ<;M?$g0mbcKMPRK zKq3@;D%Iq*_p+*r50s}FJ#u2=7pvWomY_AC`rm`fxJjf|Co-47q+eH}@>V4BY9(cz z+~{;Gwv_9OP9=T>(S_dtmHy!SDAEGR^VLW?($mxJ zrTdqc+sW;p2se_-Xr;ThUFqI$y-UE2S1$@(;Pj+U~QNkvQX{Yvy-$3@a8(gG5AQQu@ZI*a%SQIzSUA zqNk4sP1$niTR~Jz4ALD8Lsb;xeP}?9oRqC3yrg%ii&{S3>YR9v>&BexUeQqlf64?< z{C^Y|g5(>YaPKd+l9H0$_=lTyj9O6P0$hC3WazRDVrnt9!Yp3cvwZkhg zvia{8o=xnC2j9e)MqzPyonacr>uuEeQg-^?B?ikMnc%zjKaJO2L(B~({4u` zGtf9N@hCy-L?D}N>-`BL=Va1RMSo!fdtI|_UIh^GKYfk}PBn>7KUc)|t+%s64leOd zdM&OfsGCTEH=Uz!P$YW+-G{Rt)~?W*CCh{Z7sAYaClz%Mj#b4lYJERAIc61R#sMX*woNRLbe|~ zwniB5-X;CLHQMh`*=n=TK`f+bOY-s1(28~{Tk(F#Lo9^Pr&iu2>KWqe3J$^P_^2rR z7;qnP4@2W>pLIIoyiG_T0Feo_@$i8GMQE?SQym615s!i*O}vsu zLzX$L5(Rh+jCB19YP|&Fauy zL8#xH6VeD}3%hTnim5{v*Pkk=r0#HxKol6gRG2$A(@MJ^hSm777G%Ab`aLX!Oie+M z>Rr06A)l<^X`fTwl2Xt7(EDZ&6e6y5>#-(GOctn$T9rG;s92Q+>-BkVGrl}j)Y<0R z>?tx4kO!^p@PybM91Kt0V`Ab6$GN2;mHgT`s%Z3Ms3QN%uv;bUuti(bG>>4Vyo^*s)LiWLJhTS+Tk;rk_KAb1+D2?DSiwnSK_0|q zWMP9cMH|!qChtgbnB#7$_YL>n@ocn1MG1x?(o0=7ruQZ5g_wyLcfS1~%1(^|U8~_$ zR-B;7M^{mCy56_|SB=?~>iBD)CKZH%kJy{~k8*NaF8nAg0%(lrPdp|Y(vw4DKF{{z znBrG__)f|0zEEsPA`<(9I4aW6uY2mDO$0Rwk%nv+`m3T;GJDU4H?L!Gm1Y~#{Uuga zMmae|mNf1ta55--`sF19KmOSL2r(nL5U&s-A`-HK;j7AFhxC-m7;&lnQ6_Fn|0KB{ z7AK*^df$@3FU&`R@q2vO=xA^JP#qzWMmch_(Vp9L^DGZHW^a=v%nsd|NToS=|M7vM z+Y8-aHI0SLKS_U%urj@G>PLT77D0$=UrrCt|1}~0y`*-4^BZ*_`Lg?DlSE*9>CgC? zWLhD)*@sdC1;j|11X7AmzuaWtN7$KY=eWja9VyUhDOC~0Cf~OIC=PPTU6%Xz*nVe* z0$m^Djl^rQMdRpf@^g7{yP=(z<;4&fERlnF%8YgU`;(--(5R>)VUTlUv;n*1rsDPf5rx>vHd()J{=m3~haafAKn3S&TMT+PX?NnWsi~aOG{XYyn^LdBaRs z_k*Whx~i}q-syVnpxcWr2?vYUc{ zt>T(9yDMw^5Xhp8NgpLrL3d8mW(|gx3bMmd?Xjf>L9QDbNtS)FozTLC-r2CR6k0a^ zOg(|3GOE9Sut_UI|0?tL+i~3Qgyv60xZ^i)$br;mwYK9?uen_m8qPz{QR2LZmNuGBu6Uda}qlB(BHn%&my251uC=mn4c%UE*(?6>H~wkHyzV0o!ntM=w#Jgt!`-9hf7cW+hKHV> zo+-IW+1WRN53#bg?gT~uE}4EkKBXA3^@r5dVcM&!Qx4V;8EnwMVo4a><7)J8T8nwB z#tz(>Sn@{o=EYYQyENWgIWTBCuE?LCGI^@K`~YrpzgXv2Y_r}~p$l|Q7Zw(V4RNwZ zzW&lR;o(6?XBo);B|y$;$fWCFm(koB7LDH+4!}a1ZhklpdsElQLriuFu`EZY6)Rgb zcEjgKk29MG588Mm^U819_CS)Tf2&jtQEA*9fozpRtNNer)Bkj$a^a=Xkr81h`I9}) z#XUa?9CTXEB-gd+I8oY}3gygXMMUK5jIuiKN3m5xs!^`$I?2OSgLFVvz>XRL4 zx^&}52<{4w*tL2@6Dt(NxfK)p70_~o$Kn!@K)TUL-9y>&>C@f(f&y*o@dh``nYlSs zK9-=XzsIuE6~@&_wx=M)TQhYUm6DLqtf}pMu{k90YM-v&f#PfZ!PP+m%qJIlKSC!l zG7<~&o)Xf`SFeHtul~M+q)!iW2B{=z1oZ3FIUL-^mQ7mmY63`z2N|$#CsNWMOqQPf zxDQlL>u zJW2f2cXbJ&gs5U-Vu#1a<=@n6-&B>5-`7l3F0eKrbO=pKBKgZbtP#-_H)aKP2fsya4)5hW~ec zrcd{gS-q<(9|Y8~Vomm+^Yx<(+V)kDctk?E(G6XloSYoWd;X$8@Mv*>5Ca`R!nU@y zs+t;Epo-k~7sLc7XT*m=JqOrnyzsF|>XGO~Tv!-75-Ql=-&Zfzf?$sp&misWpVj~F zrF2K7D~q4SideD7)9qd_e*1{rXOJb$1sRI@_QL)=aExKOdd6 z7lH;YnU9Z8$I=o#yx_L;LW4$R<$G#${bcP@9UK}Onj_$8q9MzPfg%bA1aY^kzt22$ z%OCijtw#2Kv>xh##;v@TR^C!xK>=4LwxVv zJthf>`B)*B1}0`^c|fw{o4xjWpqm@JA&Fnh%gihb_)j#HdQr>@@lomNl+^(jUT()f zv@P)xWcMzB&j@R7Zk}m*UK14^O$5DgGRWyHI$vP)QPN|L&wsRzNf-0@0A5%c9Z?M? z{v<(nLn>Fm`X!PaLQ+KRYFUJ^b4d2)>nRM%>*7tnwDO~T;xlgYB+~q#dPjSzM4R&n z+PCS&52qmGlSdFjZ-7V!``z?zczu2S8IV_T{Ndr@U9+?L$X@JbiA`uqQWEnU5f>i7 zh}i!A*|AhbJg5E9H}BHU?R)Z5e|23F|Jiz?0zDCQhB>sMZcK>gx9r1a^r8LHwTNS; zz4RM8q~?~E%)q``$0{i*{)c!$+51fAw8Vm<3ggcIsK)rsI)a~U0G!a#8H)8qi*xsN zN(u(T)x!e~(bLm&3&}k2YjC7&j~5aoC~^Dj6r>#j=CT5*Gne+ry1Iz4^Fq|&p{wLk z56R)5KXL#jAs3d4Vau(E208~lP*A#_|FTAO6l)e%0Y-pFFG&geMNDf$N;>dT^e*+I zVoK=p!eaF?CE!+KI6nCH?c3Y&ot2FZ?UIF|AzX0&HYF*?gfeD*u~gf$HSBYro}wZ8 zfa9@rhdCFow3Mv_rXfqz4HXbK?%$LNvWtMk^{w^Bz=YtnoG$rAK`~xh*>(faJpc{B z-V^JW18Jat1@sjeq}`t4bkLj$g~mO~%D?}H_E=3VLgrm?X{kd7rOG9qEvj7z3$ar1 zg-n!ywUsZ>3j%7n_u{TDPSc*p_UC)`_V-)Du#k!Of0^#OTamGEEB)%~YPy6s4+2T? zh_YWn+E42;wk(;v*D9 z%&e>!+2YS$81Y0J@$^F>e>kaY`hX{mX4!#0a&80g!>?bzKCxzGWQ5uft>a&0)341; z7o-H!iKHRG;4C2R_SvhQM&c$2N&HRql+%Uzc^Y|ndFX-Yg!89{;k?DrxJ|k}QA`3) z*D3LbT>!98WQP4@6%OO!uOVJbz&~ahJ=kA_hN?)>)6(L(xEws`&o3&%2Gd{V=YMy7 zb+HaRIO^p~fu*kQ?xSJB(OVoGHy*vYyXsbi1}@$@^+A*Wxd0%8h{J#HD!y-R$#QUC z9dRXnjveYa6rE@-Y1u;FW?bB(;`ID{mMkg1DS?ET*$=tujmoJZ=c|znJ?Fo-dI1ax zZ*6Ujii;ESo}|u)B1i+mKn!#fV&GNl+|ts^l)nK}xB=NlHIPwmo5NK6_7k^x>eaFK zvoN^O6J(9Sf=JyPjbOYAMvd( z@6dlP5<blZvX)Sxa@_bDsHGs`+vHS(^HE1dKU_$3qi7|X(8`QE?_z2fX_n` za9mu&4J=Xg6+lTdQp$f*mDji)_t}AJj91bnR1!J zPGL`qNr9Lm26Te^SQOmc(8vh+_=iU%At51MQ&T#3pMH}Y&JxosnzRQm7`g*IQ0c;G zDCje|j3p;I6r`k4fa2;y^9@HrLV_QVA0i40+X!6`7h=XRzQ^;-utG!`X#M}Md0}L@ z|G%P*GxQuCxmbvC;0uso*X4Q6S*Ul00$gIEuz`t@({iYDSq3lSwQ%b7`t|E7v(6V6 zl~7plg9I&C+q1Yx2T|*^mCnV*1(r4Q?L#>^Oej&JI1=Byxd9IjXG=`$d3_%sx=nh+ zaj~h^VEzrCIcfv}LwuvhRu?!JfbjsS5&mSrz2yJ<5!`U+@-hRGV=OksB+|I!?x7nT zHSk)^~0}Petv+)0@^BlvcDJuY1zuGZEi2neMoJCR?Ho} zRJXDZRG7oBvILB`ZoMTpf!RYw`TeJdE1l_*&joUxZ&Imd2w>FJ)gi-y!zXS4GX0My zb?s(9vI5Q$|p@*;#@4XyAaU&cR3!j$nsPq<7QydZpAiF#W`WLEG+6(m8 zAqbPO&~*>zpMrbI|2AJ}0zw%OI#WKTY)PFkZ{SN#jcu9UzeMsYB`>zO?O52@pa570 zbVg_?9=WeBh@Or$ zI}X$_%Qg*~bi~wsrTR+9QvV~CO2}%(hGYL%lKJ05|F0pYyA9~Bz?TEIdK-xtB`O{Q z*Ouz*4{qKEGByZxD3AfdeAY}PA+ug1Y9rM2w6sWQm|L!Pg_GT*2NVvl_fPQ-utVjn zsHhiG?qlAh4g+mk`|J!fOB4jux~Hv`EiIA*Rv!g61TyeoGZd_GFr_F1;-Q;2b9w2H z%{UmB0ys;tme5P?+qZ*FDl4n#0icK9xVXa93I&zQlnVqZURCInzMHgvMC<@gz2(pg zOcu~-B@5tp+S*m{9%+gYGPY=2czyaT| zbaHfrMoI}?>GZCKy1LnGy+$Maf+si~uARp)nxr0`t$#UNQNJvo$FGRJ*uoASG^zg_ zg^3r8r0E(HdV6}X5lGK}x}aum{tj>kfh7Q#F@kTn0B^{1?4P2nZSI(RFF?81*T+YT zl7^n%ghTlD?V_$3WQeW?W{<2769XgVY?bi?s?Yz{rhuF1!j8 znH6-XGcz+CV(Fv<{R^sLWB1KKzNJG=LqcL18?YWH!Vw3_SQ)?9RGEHC-@}IwsV%og zb4{QA1gbx_8UX3CyI#LHzx2%jJ*rLJZ1aWM@+bHZjyU5w#|(btpHljGj2gaDHbZu- zjC`$2Y+P(Cyb3PZHv?j-4)sMJjSHov%5^Xx;xS2zOcLtm(zunhc zE$5q=HWhSmEsbsqqO=z(%E|(4!GP%Y!6BGvWN?HecN_8v5bck%2~oqTr2sSm)+-c{ zC94ew=nQ|Yw`QK1OKydB8F(xkBAmCSrJS%chzJRTgRal9a{YhdAiA-TP zdp}RsrVX?H*G(lOn~=?uE7Y3;aY`mT%Q2rec;$9WK6Y+bgiy|x^S5L>4aI@)9wN$3 zOgK1|XD(F7Z6*AQd!H9w?r8M203%*54Ows0TKc!qJ2Zv@nroj&P<-latDc*-;ONVt zoB@}6HPMfV^BMar$nKY zx328WxV95M_4J(?;Icp`)zYRKdhcaWYVz#t{eo91O@Cw6c&qjI7)_c!s`D=e^GaFM>0Ie?v_WY% zE@%Aua2%Jv^gUWSG~+MEtRbi-u4r)IahVYgqu-j?x2U0WhvJZsQRYYQn1_d6PrBp_f<4@1whXdb#3!IoW(P{1QFWGxX1)q%Am}k0FP5|Oj0nOpwKa_09ix2gjUODNZ2s+P>`KGoS^DK~ zx1X(9H$cH`^KWAI6+Y$juKmHyh`wYPgxN0NM>FIr{JZ00bT*BuN(LdS;8a%E-y|G3 z;fIAW^*HW+aOMyCQ5KbI*DbP3mtAWRnJ@1ZNX(r#Pt3MG^9{inVyZq+rhelyW7$1K z4@GX_qlkq^YDokm**xs_Zty_tR1AO~%jQPr2?TV4MO_*y*t4^xd z;X|}-3xOD>j_r6ES}zc~5`b{#g{Pqf232Ea*&l^wsalS1XsPOv<1hw~(!0b!1Apo{ zoxrb`eBVn0pv}?KF}0*8$yPRXEkf;&icQ&;vbWsM0L4ld5ao#=A4byqVk6a}2s|oO z0)n9nol*2E8|bQRU|Ajo8S*l=D`G*Tsrc&6T5sE2I#dR=17qV3-)vg~t^Vzk8)peU`!JDCJF3nm}bx;VX37@uE^xA&6-0r#p{t(%j8ulwIpZ(!3uuq<2|f zxc4t@&WsFaWf5X&@S{v&XguB&_4-MP?61SQ1=75SaW8Unt?^o{Q1t?~NNSGuo5q0( z#YKE=tmfR5nj-%q$?ZoQPKX+rchmBHooU0Jg z{CxFaj-AV;MGbVRhZo(6p5PpIpN(u!4~`GlmSP4l8T#jKEB(E^o_GAI)Y3F@i04Qln^Wb%;*4%G2epr!R^^0qp^v zbGSw}8e_Sw@iPFOt1ivBrs$hPjB3ddfboWE7^Cu6_z^-^Vc4iULIpcyHgm-Cyi-q+ zwPQN33)*VAPxvX4YWR<*Nut!ggsO6MN}w68Hkx!^s#o4udxyzPPS-P;Y~+dU+F7#I zoMJ$ti5=p7&Bw9;t%PfzYAgd^!xM_SX5+`B-BI#biHW*lDDO=4BtK{tGot>>QRz!I zjTcEz9f<(NRU}9{ZaHB&7$}xZlbxqTTl$RXCZ--YcGcKAx*tTho9yB6vX76vssD@q z^$9;p*&8su!==B14LLak2&B)giL4bqmsU20Mw!{(e~=cr8&$+KTaq$gQJe_Ux+!hR z-y{>F{{Gh**4tc(g1A12KJjTu3TmuwwC7$j_Qr%Qfr_BNaoh##Nnk4xf7Fzu#>;oU zWjXVlAHzgV+0nO=U2^EWo$K2V{{4PtQQB-2X$lEPA2C$=<9t6UPSa4L7=l~kLIh4% zS=E2Toz?FQ^7y;is%&|$ca2H+4B^!l{JL~M6JPrGn=7AqCC#?^y@KEn6|oMmDOJA%44s!9vs2RGJjPBr#ELMwnv>& z$3d@3Gw^X1Csk8E4YlF+=~Glij=KtaOX?=`rO$OBQ_j!tVZA@}M?Sja7=rjCiuR&`ge0 zX6!!D^u0Baw)MZ8|AY>_5DlLbG+{ zY$&u+zIQ<+i-0U0aZZG(vu1ACr_{vL?ga09&El^Bf2 z&#GolHLyeK|B1#4^*eW%Fm}&zuy|5uGnz~Ke0jy^>~rV+0=Y&B4S=S2SH9A(YjF7+ z%0=3!jQwhy4yLC(KKlgrx*=1k}@?_J$h9MVQZI6CE1j*T%SF;V9|N7Ix1 zRB`G1UEHQZ8W)j%q>kG%*Z#b=J~5@TbJyXq`2|MzYyO0~)83Z;xh|8ZL=G#{jhaUp zKJK3xfeBUJn8QQ2Mc?oy;z;_TOeuxW)jM*Z0?htA9Yr0f z)}0|L7W0PE?b86r7hGoOOiagvaM_W|E8o%C_0LqrU;>bC;Wk?H*Pr8_CEDuA{Y%cx zh4ZA~VHkG)pd!~oHus=+U!_=mARjr@%VT@qQk84OTcBje=^f+I4*qoilNzH3{93(e zY`uO!RPs-+$*@t-hpt`Xai?7%doF^gU$@^0c(ay*W$&e=Xq$4gjb<7Lq7eoJTL@rB z5LX(Kvn&UHOyP>u;38uafDk@>JtQu8cKYR83U_~ZD0I&u5Z5mA){ghG^=KwrcnwEb z%&(=>5Qv(Hx{9e`Yq)A))>kU$lN(t|16>ZKd2gpPTLXH2`2RhXr;Yv`9Qq{QrQh@8 z*>?hjNl+0+lZcx_VB)H>AZrpjceXrBFtV{3iB!*#=Z3wti>OD~Z7h4I1>dK*yJS-9r~~bIR+N6>qsx(L^Jf0WkLy*B95$n#8ycy6>8^4)Sx>Mx3ao15gz$3t za=js={2?wcmj97fB=wgrhD9O5(UZHw$pHah?#T7pzG?U%q6568t_%56+ED(PucyNL`NrPz<8D-O z-Gk3e_X@hZyEXLmqCks79-eY`b}n$X(!U`jj+{C?;|7L87pu-Oo%rVowY?BRCNqO( znwX%{sg0IsqnpQ(jwOhM$z|rJqCEM6V$}9Dpj`d*-s^f004e4Zuy?2B>f$M>y(imjTLDMLcUO*kBs?56{vte=ZGGR{APmDf zh%f_WMxMaPU(*o=q53S=U-f$SYt#9C&d2_I`J16f-89YS*8^UNU{=5P)qm}>ic6GI z8vF&ZT>iOH2fv#PGNZ}$lkEQJBYOMji=a8sp@-8MESKt80|o*ylxSxwp#)4kJor#E zIyNTV2f5pH!)b4^nd0BoL6FitaV^&Y(8Hi8X?8x{5;Qnugy6%;_4h6E9=c&iveX8SIrj z^KwXzER&or`pW(B{%s=QnAf3M6h&(6p4Qg1`PBBG^iP+!PmGm93-gE}=$U^Ee-TcY zUN`7E`+(h0Z0b>*8G!=#*If*gvklR@>~>_OW<2V<*vLwr-tT9%KmraUX!oelb0A5OI-!JiS!izHU!3#vbL8hB zD4i^S#byL54TT`E4J#F$^Bl@SX+$@GpKbpc>PFb0bm@~o6EKb9?b znlVl!6N8{fl13*V@-T9Vo7#MJJO0_nK#S2$`m^0hlkJ5WGDwL&T>GLjI7?2OvF#W= z`rKLCo*hQmQN$BM>6dU!QaZG(9bT~Bg)9)Rx*iSxD|GH3!eQidSQq=*{2Eh&hANP9 zc{X1+8cva?OMj#j$7p0~R*z8{N?;1L^nT=gl;eOT4@m{?kSFYDRVGP;yZODd?ljgl zgmOpBzNO_C%U$sy$wEEQ7og3A{P>4{U^5eevIITriHK%&bo2vvcLAUp&A2rGN2?W1 zZ0oK5p0IybKKc1IvZm-ww4?Y)>K)k7K*oy6xP>p2WBqZ)O<;iVGNb^V;-?E zbX>}5FKL0_F^3C}`(qDwe=0a^iPysk!ZVl8K8>x&yg;_w&8qWi|ihP8j0BxsMJ zH8nLyQ2Qr-hRr zK2WR9JU_C@p6(@uLxd(K&xxyd&JR7RyVq>H?4Rje14o#y)55kgd~EBsZQxtJSK0$< zPz$%n zqQs1Kd4SWMmUnQ{Yc`a@A~;>d1uBpC-ds2w`ZwPi$%fBcu>Rn6Dq_wkSFjupA3c+^{tkavIMF%Dfk9Q|bCbB@BZ_rws`*RW$v=C3ox;c&i2 zvhD0vT~6s4u97Svh-0=d4JBm*vJSzy*_a;ecB=CtQ3a8WzKYY|&@cRepnmcVp$L+SPyf_nFh89^qmMLD8Oqdx=M}oH zFe`E%yb4@_eSCdxSD(@STYuY^AJi}+WA|ERe_}@E1BoD8|1-NucJnEAwE44@^QDMV ziq_+crx5IgJmQpVqLgmc$}=SZznI$|6oc70@#`{CyP z)#DM49ugINrh9ov8Kvs)xih)YzC!t=`MxQ|;W%S=Q$BJ=kBj1uEk0cQ$aC|#NDxE> zLxR-RA^gny>VGTkP2h5F_kPj4ga(PENlFoqN&|@|Q5hN(X`Yl)RGLSnLW5>WB?(P5 zkJ3EQglM3V<|IlJN|ZX^%Ub(AXP^DPd!MuS`mFUl&wB2<{qO6(uHW?=zQZJ%59%Va zp=k@M;^ndeNs}RBlT_E2xu-hZLA9O}jMLM&;lL1kRQfe` z|3-Nd$joqP?F)5G-U-SMSw7E?kNMOr7%whTRiLD0@%q4LrEQPZ9|RbmowMu2j<~Mv z>62*@Kbe-7j{y zrz9|Wdps4KuV`aZclqRUX5v za%4(i>9?Uyi1IR}maiX$>9ek6splvXkRY#HwvtB`qSk%xHV0Et+PjZT z|8R>O$l2tTa0$rOKsd-F(4wdA`n7xGmoHj8f^RE5 z-RkQ5JbfK(ZG{1#AkFj|(b?^Q5ulS!!u6pmQlFX=LgpsAM!xnG-ofA(Jb&m&8NG|u zu4K|m9ah%_2ULhqKW&yZr;mP5fi>*bI|e@ zlXQadV8?*ku)KzZ3>V9SgS_RzxQl%6rLW7+oJwO0esWo6<+~W5t?{Ycb8Rm{NwOc2 zq}0SS70dfYXJvo>av{X-KHYUgfoG4-1TB0q$tumHsR=cz!IZT{Zbe2@Kr4Ef4a^cC zu#%Occ1AZ$_E&B}1}AvtkoS_i(Alcqk_|r~b)me6rVnRNpsVtb3)RFNI#2KN$1lh% zOU8#L4cmQuDxo3swv#QnjOO8)n}PZ4Ez0g^Vyk~UdYIbG5yXu{eD^p8(a%)JZy;SpB zhlsPnxFi8^OUuC+W#^YorlK;U)iEyhpp&1(L; z;P+LeUJ^Z2(V$bYyR?9s;-Ln~jIOq6 zV?E2~;PMQv+m577@K}0n=PP|Z&GG_zare>UDDqR3nfKTC%vh~h%ng5VN<0?$QjVv# zO7NonljxD+mwn@x1b3R;@tc_#I=IBj5%^>H=+#T+%;1__cj&vm!MHK|%URp%Ya&n2 z286w~?@6C&h;f>9{*R`5Eg}rh8PA8&Zgb7CghGaOl^0nTm1I0_c|T_|_}Ql z%OmG&Gg>8P=^lT1)<>`FlWnPQs%Nrk?oCs4(=U9|TdAVlCBz%k0j?U--LoDJID3Cm z-&UP^(_x157wHBG|Ig{`nojPJ;~Jfl`R5X@8P0;ku~UjAHCN+2!|qeN%}HwNOmD?LaU!EgYU#>O4mqCS`oEk% z0MkRK+t@)0?)<&dD`y!cIvTr);`=M2xHc+bkKGjwaO%*T7#=KgUMTGzYKi3&u$e9S z;GbecN*9l#QlVczMy!O}SVq3=_}h=BS8)4DqEo4HG}D;^3Gj}es|YgsJVr5KZ69lp z6_mMKlKtrnCjZO+M!7sk4?MLf4`Wy&S)#nSJO_{V-b%Sooase_z7hlbl$o} zji0?d5NKuhq%)KTSEh{U{BY8@M)w!p08z0dF$>0`UYB+fBhR%RsY@A6G<(QE640j- z(LN4JI8ZZo>$izzZa_bHXiUjy&*iOerOV4>!x_iU?+r03jk}KT?l7bcM1I&_b6XpG z(mAUmA;`f)MFIfJQ=e2R$!DDMf9uLDAM{>zv9vf>r6)OEnjA{-96G+FK1rL}8l={U z?(B$m)fyBksG56|gcDP`7%Ju(FD`j6t+6p@7TU@d?!8C$jZwCKESgL2M`m>$7A~hd z6j$E$#O>JN{xEu%`OZ5=4`MPK^8DA{ys77F{6;p&d+D zXFGeSFY$S2#S-71`}&UTm=?Cy)pAcfn#WA(bmRUUZWPDM9;|43J;P`Q1!&hq;`&E~ z!ism^fD3c#WN-d(k1*?sn|{}!elKphuMt<1HNwMke+bQ;kEDC(FCZxz6*Co(P`Ul1 z1|8B@{1rv7jk(5)lv~4&_~-OlH-%6ZiM=l;+0&mK%vo?tr`a1vBk^{4^2Lr42SJ`p zz7cTZYjq>QevoN<$)oPF{_suDKX*!dIDOWWvwJ> z;G=mz5IS|^>ZPpArsKJY7W|t9*Bdo9tL$#*^gst$=G~qR>i%on-6w1Lvbo8n2PFb= zR~^$@71@(j-Hw+e-S;q#c;E2T79rt7bn@A!^XGol>oEWN*?SVI5bJbar0*)M-!CPL zY$5oH&f7tZD&)=v=e_9sRkqEq=JwfjjcvGA&{(HqVlbF~$9%#0a_I7uXm~@vyo$*2 zET8m1htIUw@LG;0)Q@tnYU&l8w^2PRK=u3Sk2(-+T63d!StjV^-;-y2lpFj!JY+Y; zSBLvTKo_@->C^oftt|YQ4$&<*z==GuK5DCIMe+l^@pHR}_kEjf8k*3(BE%i{i<9Du zn(OA%L+oA2Vp5vPqw~#5o5PhQJFDB1`xzgd-OLf_*2$5RGK<8K{c=}qjE8}+$;t7?V0G!gcRr}7n@E2S1|6ZgH`)}7iB!?8+Tfk zia^10cz5i%j_-Y~`R*sH2lz&2Y#kg+2Gf5pj`E^Zs;!pv{ZYXu!#?cKZzdg=GErUC z?Zc0qBST$D`4oO3q|L~KyCimd+ZXmfe>P%h$HQbDmfyE!AC-*DI17+F58vr^rsfHm zGSjvtc*}XOn$DKK&X>SVG6`?5gI?Bo-1pQ(Zr)8kCT~H-tjzDe`os$%MR2lk!A0BF zvM{Vjm=hI%U#Em59;nGjc8XrCJ9_^5V60+amsK^MYM1izmEG5^SQjT?G7q1gsCj=MCrSwFF0MP?<6D` zeHlSna13A|;M%mF%}#K*o7Ov|jIut%j0tNNt>?tC*$8eTSdSp*~V}8Gj%3YPqIh>gh^&?%i4p6i0!%L@C|v+GAfn9yukkIa6XWcdck|2A6Ck~xxc zsAog8@5Rg{itce0U1SNns$~84(dk{77z;)08z|I@x7{RfGM>9@=4X%j!od&2#j}$e zzK-?SZ#tZ=ssBZ;Yy5GHA>B zO!I4whoC-QXW&k`SZYN%?0!?w+at&Ny2lmDOh>rjg4c)hAh0^qLzP7y-Jsnzm(24D zBRRYU%=4m~Odj8iIkU5TyEeUa{p8s0(P$sVX*tJ5`cM@w z=t(Bs2N31&1`2yW)mz+10fHxZ*aj;2n9YA#-c&KC~Lx){dKxJI^)T{b$K@ z%Azzvy608c%X<$p`?kIWlkB_mlcEeUm6#o9YE*ND*W$*)dvi8b-n-rX!%^K`YJZkq zDm0-?dzm6(*F#jx-d_=zfyyBNr=iy*S@rICa6 z{fNSBE#vEjRWFeyLq3NvS=0aD8)?h=TD5I9b3J_o?TKT)3 zn|csYS*+$7cD%hK4!JHX7i#JN8co{6l-7Y+PbCkJ^VHOrsYXWr@y<%g`si@`_K}DF zvB^iXj=tZ31Uj&mZPo0>J=2XzmQ)4;PYGuF^d+V0Xga0?UHWhKaazPn6Eu~q78{#K z-i>rGBJY?Qu`(vNlAe4Qpw5l^y;vhDIcKY@+OL4;3kJ)JQAhDh?{EJ_bVU;AlKZ!h zFZZ4YJdG5&P}VaI;Kc83+;7M#y`~lDTD;}PLkokeqU@x>jRNuK_^DGH(0^Yi+o{mDg!m%!)_VwVyDD{iluBRh87@y=HROpH&( zY-3Xm0r{p*G4S5JecKN|N#vpypj@c>H22yy>Vr~J=YVq{Jf6@|;PYyrLOJtt*MnZC zE&NS2fE>q2%NhWY+&MoV{Kvz?W3#P_K=kwG=A2{sXbz`+o>AJtoTizd++5@z`quDv3{Ep zM8I!EKI5pX>rYZ5 z-+XBfJxm3#=j;JHxQ7Q;t68B2Z~9Yh4rtyqJ)xG3Yh;dV z%5^d~Z**x}Z4*sr5(H>65x;9}T&pfS$KvSZGz+QbzV99tEL=Ei|3X$?X0^XmSX>+k z0zD>P{yr;z*s%Ne{Ts)ca5R)E%bLiO_-!D_{v6}V-~$6*o1lwXhxR7Xt^?Rp-!giS zFLno^TeX5a7P{rYdP1hcI4#VKg#cnX_hxdh(mlhqipaJX`J4SE48HhCFz)}~+BEF>0?tP<1+qw@pjj4O*pBx&x9m9B6M@gK?0IO|?9&c>bld%K_@1NbHlUdbR;XHliMi#osU z>4{rWoVx)`dW4SbGc&E~8TSqm8RNO*N>!%?Kjakc0X-G;PSD2jh(wUS)AwJbNl2@| zG1{C#jrYCYnBhG%`}dEDC_AXld^+AKlqhq)7$)3<|DJRCZ^h0UD~Pvm1RKfJ+B!rn zLCO?tNMMbkiH^zV-qL+sTelLiJ1CosKu`2j4y#sMg~xY5RS7H0Nyq;c!KZ;xRs-UM zKr}sxV}RhkG1?;#2)+eTUT=~34EE$Xaxkn%>2qBdhO-<;*e&!D6mJPMWBNZ zW?&t7Xn%!1bOHr##>a;cIT`&Yf5(HM{rEzl=}x)6bD&$l$69b&Os0>^63!6>&Ozhy z7LAczUKj!BUpN08#37EkLD3Qjds+N3>uXwnk4;R3VAV!&<*cEA&_sRED1*9vpmj-6CVNmKKrFZXb8q{AIcF$6FP0N6lc=eYadd^D3e_ zh}K+0zCuP;whcn#UI^)AXFRlQxevm^TH&)*>nI0NaLi8+M1VdSUy6?l@!k_^YQoWO z@F@WGtZ|3O+T7KKM=xHy*!sm8Qttv|#VV$mnQ)@H~&~ zjDnK2Y4bfv_Wf43OEZuCJax<9b}Z4r{Zmy8Ealz*o5ZMR?EUNfy+c;IIz_$h~5>&T_it*VsF z7qzvk_QqwYTAx)`R*w66PpFY-2NInU#AI+-4!F6bQR7ZHuB=Qi*eZ+bCJmG3pd~cR^2P2q@NJ7=4$u9=J z;!a=IW{gIB6Zu;k8yj_PGb*l|Q8G^?o_dI6wEO8LgY<9s=|pfv$p(Os1#z4pN}}LT zZ{*^#nG!h7|NG~3#6HVr@DX3RGcYkV?rV{`4r{6&KtY%%B0r;w?wMcgI4QnmXsL1_4MfqSsmt%Eionh95q zRRzzcP2>|HjPrtakBp?Tx6}y{&q>^_p9^2KzjEqxOq8a5#;#?I!fIn9a%N(SfMc6) z*~0Cl6ZDGeVK1oBO%N8!NIG=o_YLs4Gs>W_00r+EtjjXVys|yMX z>%sCRRBqxs0G{VYEHMj|p?bBp^VZ?6e)T7dQ_j0$gFJ1eu?WPX`70_b{V}ho^!&+c z*?dpzWmW0=f4P+mu|qTX|44v2iPNgLL~O@j+98DKC(Ix4FFH9mq(}V$se#BIe-?&L zJK>IIvMjb8=i!B9cXM-dkwN-7$Rg?HnH)*@8gheExU^L8u%hMzNhzu4P%LPkADf&s zvV6G{c{h{<`O0ej=r-kHy}Y2>*)Y`!WpM$y70~Mljmpr-$kfg*0%6(=%+7|!mIrbq zh}1rNMum)$Y!)z~Ln9{UE_(q-LKHcWg~Ek^xM`fJQsmr#By}3M; z?F5MR2Pm4GnGP*%?13%5KnfWc|=6(iJuAKuv1ZRquEdON)pD!o4fl3U69rt2ONRu4TFz*XFPSo z#%MySX*1$HU>gc;fo;$eCxL!qPRn}=duNuo73|!(2ISrv5?+{@scS@!>;TgTNN;&_ z^Bqaa$@ySZtKJo-L-_%^n>(WUY^Q-5iJyf0}VW)(2Y;b{j|J<>eC$6Ui3V&t}EP0TDN?p=8}U{M00p7tXLXzC*X{ zU3f8}r9dfd+}ua8U6e@ZlAO(+nUik zATqut+VOl=hoN=(^5xasT#=Nw?=NPWK@1eY9FqUK?Z@M#QH4ZA1Xfr3{44)OQA;h# z#lCCTu6G-OG z@;@kkuUE$fTX>pqY~+F*G=6-=fP)bLh{HGbEWr08f88 zCD2}EtpmgxNemha$oAto#Dz1G`JkKe`ZK%vW_b3!Mxqt&WtXHXG#vax$PL(B4ciO&cX*r-|#Z`x+ zYE@S5nBHMUjBVf~pd4L|JS_>Gke^=O=+Y87k|_=yUEP212BF1p#Grck`~%WZ71^_? z8d~fh(9PAz=jP_#+z=2L7#qw;PoE(wAuj$z)8@j33;y}Ie^bu^f+`p~H#76kcW!N# zOpSEqdH}^r-u|3aL-5HGmpE{si(T4#X+hcSjD|*NJ`$d6afzKf8`udqQ zfj$=q19?-t0aE1EZ@W=1aG0?P8|q&CimM0D>sL~$889@~aL>wYy2A9wBk_WU~ zI-K1`mH!9D$9NmL{|+1m7C4IVB2G6dgKKqliVMGf4PkxZ;8mXQg`60*!)h2z>3;d1 z`F=IJKPvm}kY&O1z}KjzBlP^w^@JXJ<`4QhR;*1mUnDI5 z8sotgN|JS8fU172A&y}33JSAt4h+LLE8VJg%Y9}_Q zE32M7H!Z^(zjwZwX~UU@{$GK1ueyW7TmCc%u+Wdt)f(o*W2a%$<(L^EDdFE@9{C_DfNuf+vEebO;GZS6aEe|#UheLs0&r!Io0hR*- zb!hHY5<_|SDN(24|NZHQ!l_whMfN8-*=WGqcCXlx@LM^tQZu8%8{VGL=T1i#44*a)#hmm819s?_s znmp6m8T1f{&X;Dnm^7@dNW(aPL3ptI7X)G;&f8sZo&o9M40L!%(m^nlNZ?PF3X(_8 zNR%pN&BGmc`LYz}ek;un2y%LWHe#|l)B z{46WU3s7Pv$`)qm(nBqQhn&KHc_*$DshFWAE%h^BpY#1Gk>iZwzQDN4V4{r0_{vVd zn!9GR3m4vbeBAuywNZk==@%&hb*FY;zy5a3dS$&YG#aKSxTaNwRHoLJrVoBGx1NeU zY_0`H#>gMJmO(|sofJGOpKLCzNUluKo9rlpUAyoxP27%uf|r+uj70%%smWd{2vh4w zys6Z5a{cDiTPNaX#Zp!26E-xYoHV+Zwz)6z2}e|EId`s-T&8h?&)M8u0oAvw@7=u{ zjQ$t#!4flO)6MJ{#;0mPu8$^~$H3;I9sDN?sor%FrUcNgb#*ryX56&m%B_Cx!5F*a z_{eNakNFSVp{tm^uQ%8dSDJaA8Ka$w?~NQb)g~9BZAZrZ8m8uND)HY`+JVN|O(m*y zMQX5x(`jyU_{3u5cUVJ1+t7Alt6Lq5T=l`_=EfZqddQ!r7ri)b9r?JEulQ7mVbS_h zE5?86NCp@rvaaTQw<9gPVryGw>DtTt?fl@i!1&3XS??~u0EY%#IY55(H>S+fb3Bsk zrfzX&s4FhdHEo-+ZXI83s41~vf594~oTknifwGvpA#oEegZ9pAC=jPz1^VumTYKJ2 zDyn98aBM2L)u^oBDRfow+S?$D*=Y$k{*+c0lQ)8Ytcw`BA-Zo|<;B4!Vhm4CC@{ez z&p%z_(rIPp8gHi+dFPm+$tWd$GM2|=I60v@#6wW`q&zJ&{s&{XJ1mUZFpHOc@L>+? z)4H>$YXD(R(+1A|YAEs!C&;O7$5##DfSpV)b<4aBbzHlndKC~%9!CRGMAfs8s+l^x zbw7SfcYHa);i$ej+l7rjg4ONC>QGX}Qhe^K*m${U<}gqc-jt1fyLf3N04~}U-?E+m zI3vY9C7Y?AhnbqkQrfJyv`(m^bY8r9*ZEZxSWDF1+jwcHUi&wFdXbb+{}sj%S0jZ; zH70{CyNx~bR)en?5SBGwYv#8hRe$2u&hXTz#&Uy&cKxlGDBd33VV&x&tt}iSZ(J|c zRnJ}{rAJ$qxLZ;*E-vn1Tvo%+#dCQpR(EG-;n;AeHRs$T{+wBIk-ldm z3>!z|OB=p{XUe`|-2VYn{jkN;A@g}?$pZUC)12|Uqb~u@8zU#j2uz5gVBO88?sE7> zgo|qb%NK&ruV!8(V7A4a6kf^ly!H`5gnn4{(^eX|Z%c|>E5p^j=cwjZoimQJ4;U>S zX7r^fwXvTf0*a1IvySFo$!z;FJ2_>&U^TY+f-`maQ%53iuI~J8EH=NKt1Kow&gcFS zXDe>XL%XJ93+}pGXi7R(`k#8H>i~GRhQ)=2<%?T_t!Ks<}(U!P0*r!o3#x3k)`@*x&7;y|i?X zik-Qx1p>`yBkCZ}&UJeMDExI+BSs=>nuI!H_GUNaYeB!dKixsIZ_N{54O_{^Bw;#U z`~CcEp1JJI&$o+x`KSH%6zfGCKjSxB58!OprR6ZA1tVF7@Y02`!H zZ2=l^lt99?2>;mNn*sHqAB?poD~G_rm@>NrcO$n5O9`)HkiXrURugP3rS-+>AUy?z zY%A)Ml}Z%1uiaP4o&D_Dvl1935skuqGF*in^`k3-)*Zw-a2B z^hFd3X!wAcJ@2+XVt?Qag9;|VZ$SA6H0MTei+IcdaQ5!xTy@Wi^6`(-<1^bk(_Hqh z!OoP)|MJj?>7~I*o47vKqx!N$2H)vO5%rg4{agOrR#@(ju)w3 zbh+$REV2d!d4lZ#b+~_PFMV=3jsED$-3`upaW73Pf18S7I#2ex6l{L$HF1!GJ(470`%M+NM*5RO-xV(Y7Iybh;Jba_-7XSI`D5T=>X>WdoPLg4dn+RB2a%M) zuI*j<)2}A;gRrI8K}J*$01S@q8dR9w$0GJMXxBloU2Hnv6|mF;beE3)xP0|IR1JFM z*r0k|t9AZ*>Gt>MP8skNc4k(+omefrbIfEqS`!NyoGey+Dy6*7NGO^zxn$H@j?1g) zg5P>7;SSE5Cp90W6-jwjW6w=jT$2yuc1)*<@;vZNkV*Wyvgxgkc8CEBNpHBma;ee} zuL}6DvtzRH&+oe?G_LYAoy2Phl=zqDB%QFn^HD6@he4(*^Kz67=iSRul*<4}fN%h> z-lNE(No!r^_vd7?xX9=*!ttCAYd?irkaYau!nLId%^y)HRX|%bS9J-ar7#4!Va$&j z3eZ^_#C8T(?iAKDYSo#&@8-*|v^Cg8P+I#3ILH%28ba&-R?SGDn6sof1Ht5Ic(%+; z*IImZJE+}$ReNMl?w!zj%tjztM2;WQzYlLJ&3;?;Ml@CZOkpU8%4?w-FbO!SCFd{o z(0=E!Yn+=_O?!AE5J9cmv!_+coj>)*$`e2-r%v{@WggxB!)}!BY`9s{8nR;naMV+L z%n?vSFv1amMMs=O8ufe5s_JtSC=Z|6x~2qYI@j8z#t(`4fmgD>w%jAp1@a6f+AYEo$B{P$uT54MWhc5Z9!t$yf{DkD9xG!wF((%!9|y# zHH02u_JUdo&EoreTfX#_A0QJK2u;4#ed_(^GU@XpBR`O1=DrJaks|#|6%YOnu6;Rq z?qu#R{p&|^(4n?OB9bVq%&luUldx*=Vb<>h%9X~0Ko=9T{!^n8I-^+{7R1o|U%7nS zjHs=%w+w+2h|{I_wdyuK1p+I+>O6LVc{ODh3)tjM0LVGiZO;ls7c^5uzejG6dDkJ; z*wUVwl12Ts{9Zt%jz&IiX9xEqK>#s+<{9=6l&@%QXbEb(KAu2eJYWTPI}2>|HD&Hy zs~|sfVe$Rha`jv`z-y%rEspVh3vnin3Sm@x#KiczlVWo+GF~t7z5d|XRw8vh@H+Df z2WN2NWR~Tv#>)2izvUE^q%L&jZ2k8AJG3I~Ojefvkj8iwbaBuoh)c>j54gbdmKG+= z5D#NU5;{ssSn!1W5L}}@)X+;4HGNJ)L_r%0Ai;&5OpPwz%1`&WDgSb!13RxYKY3lsg!f zP)J{B(20WA`I@~oskxE4(emTSq@3Co&jWC{6ZJ9oEMzK8ProHH33fxokKvXIXk^Sb z*J#h#2Iark7GXXNjWZH4(f_e`$^P){E$DJ3RNlaQemc7roo#Aj08s6kFxHt{WefZ{pGKQ5!F{ z=B9qhr7ryGDA2$3hsD++lNah2fKYUr&+wi><-hPlli@JWPvN@2r(r+7~(lYSz&)BWc4 z{zBC`J}MdtCZ?qRMK1Fw?M!rs7TZjp)kMwbIcjN;7f7U=<)O3|#;v_vJW~^prg%@K zqEhpCTEr&YvCome6}^A2z136X@o?pwrTs{Z>8J1bmkTX7Pelp&oB>1TvCr8W^}w7} zqFIFtKqa>xsu`^+K3SA}D!F#(v$IXL{U)r#w%rr+^6kt)L#yVW{oOVGtj6nEBUMw` zc6q`;!hR<&gQcxGuUKYYukLXJNi0SZSrZTGB>z}mG3TPObzgq(*q%%KVEGD_bvQ?7 zf)n@wlUK86G~IpcRcny#)Ra#&n@pdty%|5h+zcvgWlz-)9j!?1h*Zzqq0EE;C2r&v zx$}16FDvibC$uT!%;w?z5k1Q`uszSTr|T1?8-l)oi*`c+CK#|c$$GA|jc5R2;uLls zh7o)SgH|rgCtMPOTk6~vlk!JZnSZ!(@%kGE^Ip0#1Dh$o)n&SMYd`gY=d*C`!+FBs zxl9H0_{_7LrvMDq8Toy`Dl+IxW^$V9L~?iAG1?Z6B3o!+`Q zzbuW2#q#$YQO|^*xfIgQWmoVi+|s%CIYH^gUf~bU#XG~kwsp2}K!Rj6UJ|k~`-Fdc zA5cHeuOG9+uRhx!-hrwTX17%p5YSx2xuM2#fvfSiraPsdyPe*a^-diFZ{GX@qgTZ* z8RV{nQa`@xjaaJPnUH&8p5hV&lWcZRyc5;vA@0scN7OEe72);(XlOHBXeACUFPJ%x zx6`3@$Ygs-i@T(9-D?Ohg#6pqWnrp(nhd89_J(jSkHRL(d-<0or8xGg$B;i8NoTmE zoM}l4SajZGO2_Nip_r{V>VKsOK@Dh=x-)&+foujy77S0yHM6s^Pj&RJaytGbyEBC) zTP?!Q)u^^hI`?oU+YGpO)#^`*;0)pxbW&rJgKe!1Ajj|6=P zxa54`fKgalJ4K2eAJQ z%i^vtgnldD9Z9)rD_|yQs$y5W2d+|laur+33hDIzo;&_*?9_A3HzYo@nEx`m^-|9t z=XV_*ChMQJ_g!>^m~+tZMQ_2k{R=-HFf8nOEBdmw&-f}-ShAzJsG;z)h_VOOp+ zQx1U?(?})cR6LG*jVYO^&$1Wf3lhq|;`hNO5qqoVt*o@L7(dULv^v>VrpV(l>BB8d+oHMHOmMNVZgRs*4 zQgD>4+~Y;guMbefzB;KO@1DOyOttx_eok4h`<6JGlUpQjTkWy=GD!>eMaUQLECF7) z*HR@I;|UcW=f3~*>X5%c6Ptjin}D*o@=m@buYXuCo3bnXd9ScXw94~TCFdEPZ&HVP zo!5ic08kvw6TY^5LT$wUiK&elA@%;&EqAQWr{|rSh&=cj>IUe`fQ>ynnL_X z2AA!dKD2DnUn|o0diR|ByJPbKIm4$!aIi_^rZA}3V7_^WIj`{c>fTESTy^?73dRgA zWA7-pXpOwZhEbDv!osP~=oFdR8MDPM_hw>g>~`20-8@0zrBKzb1kOJGXE&Rh%=(Az zS4@RON|`9X>F+a9U6}i++)v|txpKBjd)LE`iPJ|SO6oa5*_oDOrp~=wVR1oOjd%I` znB->bKONY3`@Q!^~RZ+qwCX~)^CpcCGaPTEsOW8fW+(&~rXluy1gTH)(-5-bdB~ML$?wnc}h??}c>MY3mAkyxh zqT@lS%JvfJXnAd`nAwDD%C}_YET7Cd9kZ@3Q<=|Dro9SOx|!fm>PH)#kPz$LrKvd; z8Urs=nE^-devyLaSn8BYs_&EMon^0ztxsGMbkDjqc3Dbb_P1}$jSHVnM~@rG@X{cZ$Xv-+Nj61L5Z7zqg=sSNZ zPWliviaKP@0p_?pzMTCUm1S$ed6gHVb3PlyH`Odg!4`+hbal`}9C~6ssJM zTXd|D z;3frje9#PPFOgJ(9mAmoihx#!g){5s&3)20vR>6@SGHBAOCN%C2~=^`LXq$7}5ukfFo>YEs^?gYKd6*=7C^_&)*Rsrnyd z5M>I{VSro0bC@~+qH<@yHLG*khpti8lPBQ=J99UoT=XK}tR{bOXL- z7I^PL0Pu8+ZMU+}*QUlu|Jgyk0IHe{ar?1eol)|DFMdo3QgH#v4oPk8N4)?m5Iin~`=>3xlLgXb!0S+x zmpBCEg(DAn52zT!fv0(~=V&Sc^#VFGyz-%Z@!>z?a_a}DNBlj2 zg{!@eV{=%F!?Aruw1q;OqEz_V#@ zSn2`~GSA&K^YOX4_}x!hvc&9}XSKzA#D^=>Q)7k;-yD9rs%L8}WNUAS5Oh&dkwh2F zF_@r4L<11%3-i6~>()H(fOD5(`(~Fks;{UXv9wvjKj-I<+1S_|cQt(nB$}xxI{=&4 z;evu6N162A@NrcYHzt4$Yk9l|pOlfMDJxhnyXI34rrGZetAH6e{t!j&tsA00hHemR!*4u5mIa;38467Xe5 z(lwt7dVc?OGBctl@})Uy-pN<}WiCzQwUtzuh-Dh}U>0Z~v#VgcUbHfQzp3rsxYq9O z&o>E9Ki=T`N(d(>r~O#-)`EMrh1PkC|vKzl20Tx-q&UfO@rOWiyJU~ zg%-T4HU9nN9fGMsM2S$1#nhh|4{C)%nJ1-hA>` z01|#!)dfLyHVQ2qJ(yniS6wj>qPt-`IWC z;`ciI^WB3ENtDA_ife@jKdv(YFet>s#l;0WNa63MXdv|x@ZiDOkYC>3-f1~GuAi%* zvpd`GIW6fNPtKl3 z04vcc_#`19AVB2hZ{1>6Ja+8;>$|EuL*9MBp9%0z4VDPsVGFB=asNfxrY7l%irOaW zc}vM-PekzK9F4S#g+`{KX}P(9$C9tMU9c%?v8M^VNA={<)2FP*lG}1z07bnKx6ktZ zYnrudM-akFV2u=VUMI)!>apZ}>jz9pr~av1<%OLf5<>Oz?o;#kQI&3oUCrnF<+g>t ze*ax}n`c31xcc!zKs$zgd9^zG`i!^^ZjGCSTgJgw>AixnO&^x_*{FAd6tP<`PW97m zp}~p)&+v`p8M_~hvbU*bh&}K0vE{1gWgOh~!0X(s{stONEiF@6t-v^-!05xbxjFyb z-04+L7{_AW{}B=*rchLrgO`8&6MAerJp9gg9U?Gz!ypCD8Zwm~E{c!ef_Ziza70_Nqp}=L>yA|h{i#_|Wp*55lfkeg!M zyZbM>p_E+$xJN}r;Q+?JV9*fLzZ{5<=l)+WH>3}&Udk(5{9K3k3=#ypDJk$x-}3Xg z2rG+K-Rse~Ttxw3>x3a4ZVlNi;7`{S=}+FiR7hnLi0qa^5fBtia5Z`1Fg`vGCn<8D z%)}d%9puiO@^8%a9GV?JgvA9Ov);Rwqnq;dI)c^lU2s0~2`Y;O*!`Qk%O)A^K86UL zpU+1+uJleGxL0uL@uj%GVRpxVdU!sbxj2@ln%G8U_P^lxXMGl)C{836$rh+;(t zhy3>Qrv*7-)PzctZ^sTASm6?zAUI35wYTqZxky92coa-5; zne1xK$=jed_Glh(>PDbuek}7{!EYkz;HBfobZewV)PlgoQ0EEF6q=b=k?%VKM0hlh>AlqzSW^n=m>?BA{%vUJrWb zj7T9(0Y`FQ`SWul2M6(sKVih$M>N_GdM$|NdOL4IEsmJ+4EadC2Z3P%VhvxD9a*7M zyP{ssKttXPJSocEmgQfIiweN;gKI^hM7^W|<1aLHUB&!+>#gs3S}0FXD=4fA4hbQC zEMjsH<*QynbI$H%QgsuU7W0b|n1du44kLrECm z%xwiiJeVudstScpOlfdhvrBPFoB`*rW`Tkh)Ml@vuU|iH=(*O{iD4cHarLHacVDO~V zjj*Xj`uy4NI@OElNQ|*C=KsB)w^TGQOiie_BE7L6mrbp2V)r(6aO*~`u$#|0AL*8)T#vzFr1-U^RR zwze@oU!&;X$s zUVr_WFJkmIP!|B3##c**R7eSkb(luFC~%J9(}vs1iCdgJXc~h zI{pL(aM)N9!C_&WH8eC*gsrjj|P&PZ&`)o_6~*GO{XL`@_Y>cUDCS zoi@yEwIP@p z1#wPP6~u{g_}imLI}wOY25YXu!wO4B5V2gpSc7hmOS;b*?+d5BXUH{9X=(XqWC-9Q ze&X>5q^LWpJ6Tv*9zJ=(0xa3g{QR0$uU|;-StBDOiT@rcV47k@V1ddp{2OjThmRdw z1L6}g8Ygmjtrx`C`#ezm_}-oF4lIbk6Tqn<6-U_<1EeirsF=dskm8Fk@kSGHxDX8U zhY6-bKeBmp$0H0r#MvUoYPc4FL!?lM96kQDchB68Dw+o&<%nDbgdQ!#`K6GMieZ?f z8ZW^JhUY2@3mY5EQq^U)yXPX~OhN*=Gl~x2*XgeXK`-iOm_bxnMSR&=@Q$wfG!Q|9#*>u#vz~2N5IrY_)e zyftk4D>x{F)6)^Sk!o!u%wV&^zSo=<2NHXB=u@W*SFU~$J^7s2RAkHypVRgG_wQHp zz;|3`er+4A2HO$E_N}YeuEl*>^vghjYcxcp%2)02#FJd`69SoVzLH$&*RMaOrL`Fe zg9A|NYQF;3W4DS+NIZY{E(%|0(H0)KxwNYxI+Qq#t@r(u>o>7(&2=iuw%8`=^FE3{ y?1%qb+7 Date: Sat, 19 Mar 2022 22:09:03 -0700 Subject: [PATCH 3/4] unit test for nonuniform sampling in simulations --- control/tests/iosys_test.py | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py index bdc5ba31a..e76a6be55 100644 --- a/control/tests/iosys_test.py +++ b/control/tests/iosys_test.py @@ -1648,7 +1648,9 @@ def test_interconnect_unused_output(): outputs=['u'], name='k') - with pytest.warns(UserWarning, match=r"Unused output\(s\) in InterconnectedSystem:") as record: + with pytest.warns( + UserWarning, + match=r"Unused output\(s\) in InterconnectedSystem:") as record: h = ct.interconnect([g,s,k], inputs=['r'], outputs=['y']) @@ -1679,13 +1681,17 @@ def test_interconnect_unused_output(): pytest.fail(f'Unexpected warning: {r.message}') # warn if explicity ignored output in fact used - with pytest.warns(UserWarning, match=r"Output\(s\) specified as ignored is \(are\) used:"): + with pytest.warns( + UserWarning, + match=r"Output\(s\) specified as ignored is \(are\) used:"): h = ct.interconnect([g,s,k], inputs=['r'], outputs=['y'], ignore_outputs=['dy','u']) - with pytest.warns(UserWarning, match=r"Output\(s\) specified as ignored is \(are\) used:"): + with pytest.warns( + UserWarning, + match=r"Output\(s\) specified as ignored is \(are\) used:"): h = ct.interconnect([g,s,k], inputs=['r'], outputs=['y'], @@ -1697,3 +1703,25 @@ def test_interconnect_unused_output(): inputs=['r'], outputs=['y'], ignore_outputs=['v']) + +def test_nonuniform_timepts(): + """Test non-uniform time points for simulations""" + sys = ct.LinearIOSystem(ct.rss(2, 1, 1)) + + # Start with a uniform set of times + unifpts = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] + uniform = [1, 2, 3, 2, 1, -1, -3, -5, -7, -3, 1] + t_unif, y_unif = ct.input_output_response(sys, unifpts, uniform) + + # Create a non-uniform set of inputs + noufpts = [0, 2, 4, 8, 10] + nonunif = [1, 3, 1, -7, 1] + t_nouf, y_nouf = ct.input_output_response(sys, noufpts, nonunif) + + # Make sure the outputs agree at common times + np.testing.assert_almost_equal(y_unif[noufpts], y_nouf, decimal=6) + + # Resimulate using a new set of evaluation points + t_even, y_even = ct.input_output_response( + sys, noufpts, nonunif, t_eval=unifpts) + np.testing.assert_almost_equal(y_unif, y_even, decimal=6) From f3d46bcc9c7be4d34357d227d592c2afdbe2caa0 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Tue, 22 Mar 2022 22:18:52 -0700 Subject: [PATCH 4/4] rebase cleanup --- control/optimal.py | 15 ++++----------- control/tests/optimal_test.py | 2 +- 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/control/optimal.py b/control/optimal.py index 9dc8b225b..cb9f57ded 100644 --- a/control/optimal.py +++ b/control/optimal.py @@ -16,8 +16,9 @@ import logging import time -from .timeresp import TimeResponseData from . import config +from .exception import ControlNotImplemented +from .timeresp import TimeResponseData __all__ = ['find_optimal_input'] @@ -155,11 +156,7 @@ def __init__( # Make sure all input arguments got parsed if kwargs: - raise TypeError("unknown parameters %s" % kwargs) - - if len(kwargs) > 0: - raise ValueError( - f'unrecognized keyword(s): {list(kwargs.keys())}') + raise TypeError("unrecognized keyword(s): ", str(kwargs)) # Process trajectory constraints if isinstance(trajectory_constraints, tuple): @@ -997,16 +994,12 @@ def solve_ocp( # Process keyword arguments if trajectory_constraints is None: # Backwards compatibility - trajectory_constraints = kwargs.pop('constraints', None) + trajectory_constraints = kwargs.pop('constraints', []) # Allow 'return_x` as a synonym for 'return_states' return_states = ct.config._get_param( 'optimal', 'return_x', kwargs, return_states, pop=True) - # Process terminal constraints keyword - if constraints is None: - constraints = kwargs.pop('trajectory_constraints', []) - # Process (legacy) method keyword if kwargs.get('method'): if kwargs.get('minimize_method'): diff --git a/control/tests/optimal_test.py b/control/tests/optimal_test.py index 53d8fe8e4..1aa307b60 100644 --- a/control/tests/optimal_test.py +++ b/control/tests/optimal_test.py @@ -455,7 +455,7 @@ def test_ocp_argument_errors(): sys, time, x0, cost, constraints, initial_guess=np.zeros((4,1,1))) # Unrecognized arguments - with pytest.raises(ValueError, match="unrecognized keyword"): + with pytest.raises(TypeError, match="unrecognized keyword"): res = opt.solve_ocp( sys, time, x0, cost, constraints, terminal_constraint=None)