From 5acf66ab558a9b97a0ff427e141eb681ccced129 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Wed, 25 Sep 2024 07:53:18 +0200 Subject: [PATCH 01/12] variable measurements and repeats --- ...ontiguous_parameter_identification_wang.py | 323 ++++++++++++++++++ 1 file changed, 323 insertions(+) create mode 100644 examples-gallery/plot_non_contiguous_parameter_identification_wang.py diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py new file mode 100644 index 00000000..f3341232 --- /dev/null +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -0,0 +1,323 @@ +# %% +""" +Parameter Identification from Noncontiguous Measurements with Bias +================================================================== + +In parameter estimation it is common to collect measurements of a system's +trajectories from distinct experiments. For example, if you are identifying the +parameters of a mass-spring-damper system, you may excite the system with +different initial conditions multiple times and measure the position of the +mass. The data cannot simply be stacked end-to-end in time and the +identification run because the measurement data would be discontinuous between +each measurement trial. + +A workaround in opty is to create a set of differential equations with unique +state variables for each measurement trial that all share the same constant +parameters. You can then identify the parameters from all measurement trials +simultaneously by passing the uncoupled differential equations to opty. + +Mass-spring-damper-friction Example +=================================== + +The position of a simple system consisting of a mass connected to a fixed point +by a spring and a damper and subject to friction is simulated and recorded as +noisy measurements with bias. The spring constant, the damping coefficient +and the coefficient of friction will be identified. + +**State Variables** + +- :math:`x_i`: position of the mass of the i-th measurement trial [m] +- :math:`u_i`: speed of the mass of the i-th measurement trial [m/s] + +**Parameters** + +- :math:`m`: mass [kg] +- :math:`c`: linear damping coefficient [Ns/m] +- :math:`k`: linear spring constant [N/m] +- :math:`l_0`: natural length of the spring [m] +- :math:`friction`: friction coefficient [N] + +""" +import sympy as sm +import sympy.physics.mechanics as me +import numpy as np +from scipy.integrate import solve_ivp +from opty import Problem +from opty.utils import parse_free +import matplotlib.pyplot as plt + +# %% +# Equations of Motion. +# -------------------- +# +number_of_measurements = 4 +number_of_repeats = 3 +t0, tf = 0.0, 10.0 +num_nodes = 500 + +m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction') +bias_ident = sm.symbols(f'bias:{number_of_measurements}') + +t = me.dynamicsymbols._t +x = sm.Matrix([me.dynamicsymbols([f'x{i}_{j}' for j in range(number_of_repeats)]) + for i in range(number_of_measurements)]) +u = sm.Matrix([me.dynamicsymbols([f'u{i}_{j}' for j in range(number_of_repeats)]) + for i in range(number_of_measurements)]) + +n = sm.Matrix([me.dynamicsymbols([f'n{i}_{j}' for j in range(number_of_repeats)]) + for i in range(number_of_measurements)]) + +xh, uh = me.dynamicsymbols('xh, uh') + +# %% +# Form the equations of motion, such that the first number_of_repeats equations +# belong to the first measurement, the second number_of_repeats equations belong +# to the second measurement, and so on. +def kd_eom_rh(xh, uh, m, c, k, l0, friction): + """" sets up the eoms for the system""" + N = me.ReferenceFrame('N') + O, P = sm.symbols('O, P', cls=me.Point) + O.set_vel(N, 0) + + P.set_pos(O, xh*N.x) + P.set_vel(N, uh*N.x) + bodies = [me.Particle('part', P, m)] + # :math:`\tanh(\alpha \cdot x) \approx sign(x)` for large :math:`\alpha`, and + # is differentiale everywhere. + forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(30*uh)*N.x)] + kd = sm.Matrix([uh - xh.diff(t)]) + KM = me.KanesMethod(N, + q_ind=[xh], + u_ind=[uh], + kd_eqs=kd + ) + fr, frstar = KM.kanes_equations(bodies, forces) + rhs = KM.rhs() + return (kd, fr + frstar, rhs) + +kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction) + +kd_total = sm.Matrix([]) +eom_total = sm.Matrix([]) +rhs_total = sm.Matrix([]) +for i in range(number_of_measurements): + for j in range(number_of_repeats): + eom_h = me.msubs(fr_frstar, {xh: x[i, j], uh: u[i, j], uh.diff(t): u[i, j].diff(t)}) + eom_h = eom_h + sm.Matrix([n[i, j]]) + kd_h = kd.subs({xh: x[i, j], uh: u[i, j]}) + kd_total = kd_total.col_join(kd_h) + eom_total = eom_total.col_join(eom_h) + rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})]) + rhs_total = rhs_total.col_join(rhs_h) + +eom = kd_total.col_join(eom_total) +uh =sm.Matrix([u[i, j] for i in range(number_of_measurements) + for j in range(number_of_repeats)]) +rhs = uh.col_join(rhs_total) + +for i in range(number_of_repeats): + sm.pprint(eom[i]) +for i in range(3): + print('.........') +for i in range(number_of_repeats): + sm.pprint(eom[-(number_of_repeats-i)]) + +# %% +# Generate Noisy Measurement Data with Bias +# ----------------------------------------- +# Create 'number_of_measurements' sets of measurements with different initial +# conditions. To get the measurements, the equations of motion are integrated, +# and then noise is added to each point in time of the solution. Finally a bias +# is added. +# + +states = [x[i, j] for i in range(number_of_measurements) + for j in range(number_of_repeats)] + \ + [u[i, j] for i in range(number_of_measurements) + for j in range(number_of_repeats)] +parameters = [m, c, k, l0, friction] +par_vals = [1.0, 0.25, 2.0, 1.0, 0.5] + +eval_rhs = sm.lambdify(states + parameters, rhs) + +times = np.linspace(t0, tf, num=num_nodes) + +measurements = [] +# %% + +for i in range(number_of_measurements): + for j in range(number_of_repeats): + seed = 1234*(i+1)*(j+1) + np.random.seed(seed) + x0 = 4*np.random.randn(2*number_of_measurements * number_of_repeats) + sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(), + (t0, tf), x0, t_eval=times, args=(par_vals,)) + seed = 10000 + 12345*i + np.random.seed(seed) + measurements.append(sol.y[1, :] + 1.0*np.random.randn(num_nodes)) +measurements = np.array(measurements) +print(measurements.shape) +print(sol.y.shape) + +fig, ax = plt.subplots(2*number_of_measurements*number_of_repeats, 1, + figsize=(6.4, 1.5*number_of_measurements*number_of_repeats), + sharex=True, constrained_layout=True) + +for i in range(2*number_of_measurements*number_of_repeats): + ax[i].plot(times, sol.y[i, :], color='blue') + ax[i].set_ylabel(f'{states[i]}') + +# Add bias to the measurements +bias = 0.0 *np.random.uniform(0., 1., size=number_of_measurements) +measurements = np.array(measurements) + bias[:, np.newaxis] + +print(measurements.shape) + +# %% +# Setup the Identification Problem +# -------------------------------- +# +# The goal is to identify the damping coefficient :math:`c`, the spring +# constant :math:`k` and the coefficient of friction :math:`friction`. +# The objective :math:`J` is to minimize the least square +# difference in the optimal simulation as compared to the measurements. If +# some measurement is considered more reliable, its weight :math:`w` may be +# increased relative to the other measurements. +# +# +# objective = :math:`\int_{t_0}^{t_f} \left[ \sum_{i=1}^{i = \text{number_of_measurements}} (w_i (x_i - x_i^m)^2 \right] \hspace{2pt} dt` +# +interval_value = (tf - t0) / (num_nodes - 1) + +w =[1 for _ in range(number_of_measurements)] + +nm = number_of_measurements +nr = number_of_repeats +def obj(free): + sum = 0.0 + for i in range(nm): + for j in range(nr): + sum += (w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - + measurements[i])**2) + + +def obj_grad(free): + grad = np.zeros_like(free) + for i in range(nm): + for j in range(nr): + grad[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] = 2*w[i]*interval_value*( + free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - measurements[i] + ) + return grad + + +# %% +# By not including :math:`c`, :math:`k`, :math:`friction`, :math:`bias_i` +# in the parameter map, they will be treated as unknown parameters. +par_map = {m: par_vals[0], l0: par_vals[3]} + +# %% +bounds = { + c: (0.01, 2.0), + k: (0.1, 10.0), + friction: (0.1, 10.0), +} +noise_scale = 2.0 +known_trajectory_map = {} +for i in range(number_of_measurements): + for j in range(number_of_repeats): + seed = 10000 + 12345*(i+1)*(j+1) + np.random.seed(seed) + known_trajectory_map = (known_trajectory_map | + {n[i, j]: noise_scale*np.random.randn(num_nodes)}) +print(len(known_trajectory_map)) +print('states', len(states)) +problem = Problem( + obj, + obj_grad, + eom, + states, + num_nodes, + interval_value, + known_parameter_map=par_map, + time_symbol=me.dynamicsymbols._t, + integration_method='midpoint', + bounds=bounds, + known_trajectory_map=known_trajectory_map +) +print(problem.num_free) +# %% +# This give the sequence of the unknown parameters at the tail of solution. +print(problem.collocator.unknown_parameters) +print(problem.collocator.unknown_input_trajectories) +# %% +# Create an Initial Guess +# ----------------------- +# +# It is reasonable to use the measurements as initial guess for the states +# because they would be available. Here, only the measurements of the position +# are used and the speeds are set to zero. The last values are the guesses +# for :math:`bias_i`, :math:`c` and :math:`k`, respectively. +# +initial_guess = np.array([]) +for i in range(number_of_measurements): + for j in range(number_of_repeats): + initial_guess = np.concatenate((initial_guess, measurements[i, :])) +initial_guess = np.hstack((initial_guess, + np.zeros(number_of_measurements*number_of_repeats*num_nodes), + [0.1, 3.0, 2.0])) + +print(len(initial_guess)) + +# %% +# Solve the Optimization Problem +# ------------------------------ +# +# Here, the initial guess is the result of some previous run, stored in +# 'non_contiguous_parameter_identification_bias_solution.npy'. +#initial_guess = np.load('non_contiguous_parameter_identification_bias_solution.npy') +solution, info = problem.solve(initial_guess) +# %% +# This is how the solution may be saved for a future initial guess +# ```np.save('non_contiguous_parameter_identification_bias_solution', solution)``` + +print(info['status_msg']) +print(f'final value of the objective function is {info['obj_val']:.2f}' ) + +# %% +problem.plot_objective_value() + +# %% +problem.plot_constraint_violations(solution) + + +# %% +# The identified parameters are: +# ------------------------------ +# +print(f'Estimate of damping coefficient is {solution[-3]: 1.2f}' + + f' Percentage error is {(solution[-3]-par_vals[1])/solution[-3]*100:1.2f} %') +print(f'Estimate of the spring constant is {solution[-1]: 1.2f}' + + f' Percentage error is {(solution[-1]-par_vals[2])/solution[-1]*100:1.2f} %') +print(f'Estimate of the friction coefficient is {solution[-2]: 1.2f}' + f' Percentage error is {(solution[-2]-par_vals[-1])/solution[-2]*100:1.2f} %') + +# %% +# Plot the measurements and the trajectories calculated. +# +fig, ax = plt. subplots(number_of_measurements, 1, + figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True) +for i in range(number_of_measurements): + ax[i].plot(times, solution[i*num_nodes:(i+1)*num_nodes], color='red') + ax[i].plot(times, measurements[i], color='blue', lw=0.5) + ax[i].set_ylabel(f'{states[i]}') +ax[-1].set_xlabel('Time [s]') +ax[0].set_title('Trajectories') +prevent_output = True + +# %% +# +# sphinx_gallery_thumbnail_number = 3 +problem.plot_trajectories(solution) + +plt.show() From 231908bcb6cd7555168e2030f73c84c31494b572 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Wed, 25 Sep 2024 13:15:11 +0200 Subject: [PATCH 02/12] including noise for every measurement --- ...ontiguous_parameter_identification_wang.py | 136 ++++++++++-------- 1 file changed, 75 insertions(+), 61 deletions(-) diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py index f3341232..8aeb22b7 100644 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -1,7 +1,6 @@ -# %% """ -Parameter Identification from Noncontiguous Measurements with Bias -================================================================== +Parameter Identification from Noncontiguous Measurements with Bias and added Noise +================================================================================== In parameter estimation it is common to collect measurements of a system's trajectories from distinct experiments. For example, if you are identifying the @@ -15,6 +14,18 @@ state variables for each measurement trial that all share the same constant parameters. You can then identify the parameters from all measurement trials simultaneously by passing the uncoupled differential equations to opty. +This idea is due to Jason Moore. + +In addition, for each measurement, one may create a set of differential +equations with unique state variables that again share the same constant +parameters. In addition, one adds noise to these equations. (Intuitively this +noise corresponds to a force acting on the system). +This idea is due to Huawei Wang and A. J. van den Bogert. + +So, if one has 'number_of_measurements' measurements, and one creates +'number_of_repeats' differential equations for each measurement, one will have +'number_of_measurements' * 'number_of_repeats' uncoupled differential equations, +all sharing the same constant parameters. Mass-spring-damper-friction Example =================================== @@ -26,8 +37,12 @@ **State Variables** -- :math:`x_i`: position of the mass of the i-th measurement trial [m] -- :math:`u_i`: speed of the mass of the i-th measurement trial [m/s] +- :math:`x_{i, j}`: position of the mass of the i-th measurement trial [m] +- :math:`u_{i, j}`: speed of the mass of the i-th measurement trial [m/s] + +**Noise Variables** + +- :math:`n_{i, j}`: noise added [N] **Parameters** @@ -49,14 +64,15 @@ # %% # Equations of Motion. # -------------------- -# -number_of_measurements = 4 -number_of_repeats = 3 +#========================================= +# Basic data may be set here +number_of_measurements = 3 +number_of_repeats = 5 t0, tf = 0.0, 10.0 num_nodes = 500 +#========================================= m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction') -bias_ident = sm.symbols(f'bias:{number_of_measurements}') t = me.dynamicsymbols._t x = sm.Matrix([me.dynamicsymbols([f'x{i}_{j}' for j in range(number_of_repeats)]) @@ -97,6 +113,8 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction) +# %% +# Stack the equations appropriately. kd_total = sm.Matrix([]) eom_total = sm.Matrix([]) rhs_total = sm.Matrix([]) @@ -111,6 +129,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): rhs_total = rhs_total.col_join(rhs_h) eom = kd_total.col_join(eom_total) + uh =sm.Matrix([u[i, j] for i in range(number_of_measurements) for j in range(number_of_repeats)]) rhs = uh.col_join(rhs_total) @@ -123,12 +142,11 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): sm.pprint(eom[-(number_of_repeats-i)]) # %% -# Generate Noisy Measurement Data with Bias -# ----------------------------------------- +# Generate Noisy Measurement +# -------------------------- # Create 'number_of_measurements' sets of measurements with different initial # conditions. To get the measurements, the equations of motion are integrated, -# and then noise is added to each point in time of the solution. Finally a bias -# is added. +# and then noise is added to each point in time of the solution. # states = [x[i, j] for i in range(number_of_measurements) @@ -144,7 +162,8 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): measurements = [] # %% - +# Integrate the differential equations. If np.random.seed(seed) is used, it +# the seed must be changed for every measurement to ensure they are independent. for i in range(number_of_measurements): for j in range(number_of_repeats): seed = 1234*(i+1)*(j+1) @@ -154,24 +173,10 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): (t0, tf), x0, t_eval=times, args=(par_vals,)) seed = 10000 + 12345*i np.random.seed(seed) - measurements.append(sol.y[1, :] + 1.0*np.random.randn(num_nodes)) + measurements.append(sol.y[0, :] + 1.0*np.random.randn(num_nodes)) measurements = np.array(measurements) -print(measurements.shape) -print(sol.y.shape) - -fig, ax = plt.subplots(2*number_of_measurements*number_of_repeats, 1, - figsize=(6.4, 1.5*number_of_measurements*number_of_repeats), - sharex=True, constrained_layout=True) - -for i in range(2*number_of_measurements*number_of_repeats): - ax[i].plot(times, sol.y[i, :], color='blue') - ax[i].set_ylabel(f'{states[i]}') - -# Add bias to the measurements -bias = 0.0 *np.random.uniform(0., 1., size=number_of_measurements) -measurements = np.array(measurements) + bias[:, np.newaxis] - -print(measurements.shape) +print('shsape of measurement array', measurements.shape) +print('shape of solution array', sol.y.shape) # %% # Setup the Identification Problem @@ -185,7 +190,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): # increased relative to the other measurements. # # -# objective = :math:`\int_{t_0}^{t_f} \left[ \sum_{i=1}^{i = \text{number_of_measurements}} (w_i (x_i - x_i^m)^2 \right] \hspace{2pt} dt` +# objective = :math:`\int_{t_0}^{t_f} \sum_{i=1}^{\text{number_of_measurements}} \left[ \sum_{s=1}^{\text{number_of_repeats}} (w_i (x_{i, s} - x_{i,s}^m)^2 \right] \hspace{2pt} dt` # interval_value = (tf - t0) / (num_nodes - 1) @@ -197,8 +202,9 @@ def obj(free): sum = 0.0 for i in range(nm): for j in range(nr): - sum += (w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - - measurements[i])**2) + sum += np.sum((w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - + measurements[i])**2)) + return sum def obj_grad(free): @@ -212,17 +218,20 @@ def obj_grad(free): # %% -# By not including :math:`c`, :math:`k`, :math:`friction`, :math:`bias_i` +# By not including :math:`c`, :math:`k`, :math:`friction` # in the parameter map, they will be treated as unknown parameters. par_map = {m: par_vals[0], l0: par_vals[3]} -# %% bounds = { c: (0.01, 2.0), k: (0.1, 10.0), friction: (0.1, 10.0), } -noise_scale = 2.0 +# %% +# Set up the known trajectory map. If np.random.seed(seed) is used, the +# seed must be changed for every map to ensure they are idependent. +# noise_scale gives the 'strength' of the noise. +noise_scale = 1.0 known_trajectory_map = {} for i in range(number_of_measurements): for j in range(number_of_repeats): @@ -230,8 +239,9 @@ def obj_grad(free): np.random.seed(seed) known_trajectory_map = (known_trajectory_map | {n[i, j]: noise_scale*np.random.randn(num_nodes)}) -print(len(known_trajectory_map)) -print('states', len(states)) + +# %% +# Set up the problem. problem = Problem( obj, obj_grad, @@ -245,11 +255,9 @@ def obj_grad(free): bounds=bounds, known_trajectory_map=known_trajectory_map ) -print(problem.num_free) # %% # This give the sequence of the unknown parameters at the tail of solution. print(problem.collocator.unknown_parameters) -print(problem.collocator.unknown_input_trajectories) # %% # Create an Initial Guess # ----------------------- @@ -257,7 +265,7 @@ def obj_grad(free): # It is reasonable to use the measurements as initial guess for the states # because they would be available. Here, only the measurements of the position # are used and the speeds are set to zero. The last values are the guesses -# for :math:`bias_i`, :math:`c` and :math:`k`, respectively. +# for :math:`c`, :math:`friction` and :math:`k`, respectively. # initial_guess = np.array([]) for i in range(number_of_measurements): @@ -267,20 +275,11 @@ def obj_grad(free): np.zeros(number_of_measurements*number_of_repeats*num_nodes), [0.1, 3.0, 2.0])) -print(len(initial_guess)) - # %% # Solve the Optimization Problem # ------------------------------ # -# Here, the initial guess is the result of some previous run, stored in -# 'non_contiguous_parameter_identification_bias_solution.npy'. -#initial_guess = np.load('non_contiguous_parameter_identification_bias_solution.npy') solution, info = problem.solve(initial_guess) -# %% -# This is how the solution may be saved for a future initial guess -# ```np.save('non_contiguous_parameter_identification_bias_solution', solution)``` - print(info['status_msg']) print(f'final value of the objective function is {info['obj_val']:.2f}' ) @@ -304,20 +303,35 @@ def obj_grad(free): # %% # Plot the measurements and the trajectories calculated. +#------------------------------------------------------- # -fig, ax = plt. subplots(number_of_measurements, 1, - figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True) -for i in range(number_of_measurements): - ax[i].plot(times, solution[i*num_nodes:(i+1)*num_nodes], color='red') - ax[i].plot(times, measurements[i], color='blue', lw=0.5) - ax[i].set_ylabel(f'{states[i]}') -ax[-1].set_xlabel('Time [s]') -ax[0].set_title('Trajectories') +sol_parsed, _, _ = parse_free(solution, len(states), 0, num_nodes ) +if number_of_measurements > 1: + fig, ax = plt. subplots(number_of_measurements, 1, + figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True) + + for i in range(number_of_measurements): + ax[i].plot(times, sol_parsed[i*number_of_repeats, :], color='red') + ax[i].plot(times, measurements[i], color='blue', lw=0.5) + ax[i].set_ylabel(f'{states[i]}') + ax[-1].set_xlabel('Time [s]') + ax[0].set_title('Trajectories') + +else: + fig, ax = plt. subplots(1, 1, figsize=(6.4, 1.25)) + ax.plot(times, sol_parsed[0, :], color='red') + ax.plot(times, measurements[0], color='blue', lw=0.5) + ax.set_ylabel(f'{states[0]}') + ax.set_xlabel('Time [s]') + ax.set_title('Trajectories') + + prevent_output = True +# %% +#Plot the Trajectories +problem.plot_trajectories(solution) # %% # # sphinx_gallery_thumbnail_number = 3 -problem.plot_trajectories(solution) - plt.show() From 97b0da6693edecccf2754455e2c5722988efbbc6 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Wed, 25 Sep 2024 13:35:00 +0200 Subject: [PATCH 03/12] shortened title --- .../plot_non_contiguous_parameter_identification_wang.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py index 8aeb22b7..e1cac73e 100644 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -1,6 +1,6 @@ """ -Parameter Identification from Noncontiguous Measurements with Bias and added Noise -================================================================================== +Parameter Identification from Noncontiguous Measurements with Bias and Noise +============================================================================ In parameter estimation it is common to collect measurements of a system's trajectories from distinct experiments. For example, if you are identifying the @@ -302,8 +302,8 @@ def obj_grad(free): f' Percentage error is {(solution[-2]-par_vals[-1])/solution[-2]*100:1.2f} %') # %% -# Plot the measurements and the trajectories calculated. -#------------------------------------------------------- +# Plot the measurements and the trajectories calculated +#------------------------------------------------------ # sol_parsed, _, _ = parse_free(solution, len(states), 0, num_nodes ) if number_of_measurements > 1: From e39b75a07bcfbf269e8c28469d28ff4ed356833e Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Wed, 25 Sep 2024 13:51:34 +0200 Subject: [PATCH 04/12] no idea which title is wrong --- ...t_non_contiguous_parameter_identification_wang.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py index e1cac73e..4a35ad4c 100644 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -64,13 +64,13 @@ # %% # Equations of Motion. # -------------------- -#========================================= +# # Basic data may be set here number_of_measurements = 3 number_of_repeats = 5 t0, tf = 0.0, 10.0 num_nodes = 500 -#========================================= +# m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction') @@ -89,6 +89,7 @@ # Form the equations of motion, such that the first number_of_repeats equations # belong to the first measurement, the second number_of_repeats equations belong # to the second measurement, and so on. +# def kd_eom_rh(xh, uh, m, c, k, l0, friction): """" sets up the eoms for the system""" N = me.ReferenceFrame('N') @@ -115,6 +116,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): # %% # Stack the equations appropriately. +# kd_total = sm.Matrix([]) eom_total = sm.Matrix([]) rhs_total = sm.Matrix([]) @@ -164,6 +166,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): # %% # Integrate the differential equations. If np.random.seed(seed) is used, it # the seed must be changed for every measurement to ensure they are independent. +# for i in range(number_of_measurements): for j in range(number_of_repeats): seed = 1234*(i+1)*(j+1) @@ -220,6 +223,7 @@ def obj_grad(free): # %% # By not including :math:`c`, :math:`k`, :math:`friction` # in the parameter map, they will be treated as unknown parameters. +# par_map = {m: par_vals[0], l0: par_vals[3]} bounds = { @@ -231,6 +235,7 @@ def obj_grad(free): # Set up the known trajectory map. If np.random.seed(seed) is used, the # seed must be changed for every map to ensure they are idependent. # noise_scale gives the 'strength' of the noise. +# noise_scale = 1.0 known_trajectory_map = {} for i in range(number_of_measurements): @@ -242,6 +247,7 @@ def obj_grad(free): # %% # Set up the problem. +# problem = Problem( obj, obj_grad, @@ -330,8 +336,10 @@ def obj_grad(free): # %% #Plot the Trajectories +# problem.plot_trajectories(solution) # %% # # sphinx_gallery_thumbnail_number = 3 +# plt.show() From f2cb12d6e038227ddc97149d9c2031fdf815e366 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Thu, 26 Sep 2024 07:59:09 +0200 Subject: [PATCH 05/12] cosmetic changes, bit faster execution --- ...ontiguous_parameter_identification_wang.py | 67 ++++++++++--------- 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py index 4a35ad4c..4878433f 100644 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -1,3 +1,4 @@ +# %% """ Parameter Identification from Noncontiguous Measurements with Bias and Noise ============================================================================ @@ -22,9 +23,9 @@ noise corresponds to a force acting on the system). This idea is due to Huawei Wang and A. J. van den Bogert. -So, if one has 'number_of_measurements' measurements, and one creates -'number_of_repeats' differential equations for each measurement, one will have -'number_of_measurements' * 'number_of_repeats' uncoupled differential equations, +So, if one has *number_of_measurements* measurements, and one creates +*number_of_repeats* differential equations for each measurement, one will have +*number_of_measurements* * *number_of_repeats* uncoupled differential equations, all sharing the same constant parameters. Mass-spring-damper-friction Example @@ -32,7 +33,7 @@ The position of a simple system consisting of a mass connected to a fixed point by a spring and a damper and subject to friction is simulated and recorded as -noisy measurements with bias. The spring constant, the damping coefficient +noisy measurements. The spring constant, the damping coefficient and the coefficient of friction will be identified. **State Variables** @@ -40,9 +41,6 @@ - :math:`x_{i, j}`: position of the mass of the i-th measurement trial [m] - :math:`u_{i, j}`: speed of the mass of the i-th measurement trial [m/s] -**Noise Variables** - -- :math:`n_{i, j}`: noise added [N] **Parameters** @@ -51,6 +49,7 @@ - :math:`k`: linear spring constant [N/m] - :math:`l_0`: natural length of the spring [m] - :math:`friction`: friction coefficient [N] +- :math:`n_{i, j}`: noise added [N] """ import sympy as sm @@ -66,12 +65,12 @@ # -------------------- # # Basic data may be set here -number_of_measurements = 3 -number_of_repeats = 5 +number_of_measurements = 4 +number_of_repeats = 3 t0, tf = 0.0, 10.0 num_nodes = 500 +noise_scale = 1.0 # - m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction') t = me.dynamicsymbols._t @@ -86,9 +85,7 @@ xh, uh = me.dynamicsymbols('xh, uh') # %% -# Form the equations of motion, such that the first number_of_repeats equations -# belong to the first measurement, the second number_of_repeats equations belong -# to the second measurement, and so on. +# Form the equations of motion, I use Kane's method here # def kd_eom_rh(xh, uh, m, c, k, l0, friction): """" sets up the eoms for the system""" @@ -101,7 +98,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): bodies = [me.Particle('part', P, m)] # :math:`\tanh(\alpha \cdot x) \approx sign(x)` for large :math:`\alpha`, and # is differentiale everywhere. - forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(30*uh)*N.x)] + forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(60*uh)*N.x)] kd = sm.Matrix([uh - xh.diff(t)]) KM = me.KanesMethod(N, q_ind=[xh], @@ -115,19 +112,24 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction) # %% -# Stack the equations appropriately. +# Stack the equations such that the first number_of_repeats equations +# belong to the first measurement, the second number_of_repeats equations belong +# to the second measurement, and so on. # kd_total = sm.Matrix([]) eom_total = sm.Matrix([]) rhs_total = sm.Matrix([]) + for i in range(number_of_measurements): for j in range(number_of_repeats): - eom_h = me.msubs(fr_frstar, {xh: x[i, j], uh: u[i, j], uh.diff(t): u[i, j].diff(t)}) - eom_h = eom_h + sm.Matrix([n[i, j]]) + eom_h = fr_frstar.subs({xh: x[i, j], uh: u[i, j], + uh.diff(t): u[i, j].diff(t), xh.diff(t): u[i, j]}) + \ + sm.Matrix([n[i, j]]) kd_h = kd.subs({xh: x[i, j], uh: u[i, j]}) + rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})]) + kd_total = kd_total.col_join(kd_h) eom_total = eom_total.col_join(eom_h) - rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})]) rhs_total = rhs_total.col_join(rhs_h) eom = kd_total.col_join(eom_total) @@ -155,6 +157,7 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): for j in range(number_of_repeats)] + \ [u[i, j] for i in range(number_of_measurements) for j in range(number_of_repeats)] + parameters = [m, c, k, l0, friction] par_vals = [1.0, 0.25, 2.0, 1.0, 0.5] @@ -164,22 +167,24 @@ def kd_eom_rh(xh, uh, m, c, k, l0, friction): measurements = [] # %% -# Integrate the differential equations. If np.random.seed(seed) is used, it -# the seed must be changed for every measurement to ensure they are independent. +# Integrate the differential equations. If np.random.seed(seed) is used, the +# seed must be changed for every measurement to ensure they are independent. # for i in range(number_of_measurements): - for j in range(number_of_repeats): - seed = 1234*(i+1)*(j+1) - np.random.seed(seed) - x0 = 4*np.random.randn(2*number_of_measurements * number_of_repeats) - sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(), - (t0, tf), x0, t_eval=times, args=(par_vals,)) - seed = 10000 + 12345*i + seed = 1234*(i+1) + np.random.seed(seed) + + x0 = 3.0*np.random.randn(2*number_of_measurements * number_of_repeats) + sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(), + (t0, tf), x0, t_eval=times, args=(par_vals,)) + + seed = 10000 + 12345*(i+1) np.random.seed(seed) measurements.append(sol.y[0, :] + 1.0*np.random.randn(num_nodes)) + measurements = np.array(measurements) -print('shsape of measurement array', measurements.shape) -print('shape of solution array', sol.y.shape) +print('shape of measurement array', measurements.shape) +print('shape of solution array ', sol.y.shape) # %% # Setup the Identification Problem @@ -209,7 +214,6 @@ def obj(free): measurements[i])**2)) return sum - def obj_grad(free): grad = np.zeros_like(free) for i in range(nm): @@ -236,7 +240,6 @@ def obj_grad(free): # seed must be changed for every map to ensure they are idependent. # noise_scale gives the 'strength' of the noise. # -noise_scale = 1.0 known_trajectory_map = {} for i in range(number_of_measurements): for j in range(number_of_repeats): @@ -319,7 +322,7 @@ def obj_grad(free): for i in range(number_of_measurements): ax[i].plot(times, sol_parsed[i*number_of_repeats, :], color='red') ax[i].plot(times, measurements[i], color='blue', lw=0.5) - ax[i].set_ylabel(f'{states[i]}') + ax[i].set_ylabel(f'{states[number_of_repeats* i]}') ax[-1].set_xlabel('Time [s]') ax[0].set_title('Trajectories') From 3523f76afe1cdad6476a8a2cdb9c5e0be43c730d Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Thu, 26 Sep 2024 08:08:52 +0200 Subject: [PATCH 06/12] try again --- .../plot_non_contiguous_parameter_identification_wang.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py index 4878433f..2597502f 100644 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -85,7 +85,7 @@ xh, uh = me.dynamicsymbols('xh, uh') # %% -# Form the equations of motion, I use Kane's method here +# Form the equations of motion, I use Kane's method here. # def kd_eom_rh(xh, uh, m, c, k, l0, friction): """" sets up the eoms for the system""" From c0cd80599d81bb3c7ee0faa31c28bf6b9b0855c3 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Sat, 16 Nov 2024 14:44:38 +0100 Subject: [PATCH 07/12] just to see if it wasses all test now --- .../plot_non_contiguous_parameter_identification_wang.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py index 2597502f..fd96d8c9 100644 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py @@ -1,4 +1,3 @@ -# %% """ Parameter Identification from Noncontiguous Measurements with Bias and Noise ============================================================================ From 3f6139d5bf02f62dd0637e47b856d7b1cdb62320 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Tue, 19 Nov 2024 10:49:12 +0100 Subject: [PATCH 08/12] betts examples 10.144 and 10.145 --- examples-gallery/plot_betts_10_144_145.py | 220 ++++++++++++++++++++++ 1 file changed, 220 insertions(+) create mode 100644 examples-gallery/plot_betts_10_144_145.py diff --git a/examples-gallery/plot_betts_10_144_145.py b/examples-gallery/plot_betts_10_144_145.py new file mode 100644 index 00000000..176bdbe3 --- /dev/null +++ b/examples-gallery/plot_betts_10_144_145.py @@ -0,0 +1,220 @@ +""" +Van der Pol Oscillator +====================== + +https://en.wikipedia.org/wiki/Van_der_Pol_oscillator + +These are example 10.144 / 145 from Betts' book "Practical Methods for Optimal Control +Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems. +It is described in more detail in section 4.14. example 4.11 of the book. + +""" +import numpy as np +import sympy as sm +import sympy.physics.mechanics as me +from opty.direct_collocation import Problem +from opty.utils import create_objective_function +import matplotlib.pyplot as plt + +# %% +# First version of the problem (10.144). +# -------------------------------------- +# **States** +# :math:`y_1, y_2` : state variables +# +# **Controls** +# :math:`u` : control variables +# +# Equations of Motion. +# -------------------- +t = me.dynamicsymbols._t + +y1, y2 = me.dynamicsymbols('y1, y2') +u = me.dynamicsymbols('u') + +eom = sm.Matrix([ + -y1.diff() + y2, + -y2.diff(t) + (1 - y1**2)*y2 - y1 + u, +]) +sm.pprint(eom) + +# %% +# Define and Solve the Optimization Problem. +# ----------------------------------------- +num_nodes = 2001 +iterations = 1000 + +t0, tf = 0.0, 5.0 +interval_value = (tf - t0)/(num_nodes - 1) + +state_symbols = (y1, y2) +unkonwn_input_trajectories = (u, ) + +objective = sm.integrate(u**2 + y1**2 + y2**2, t) +obj, obj_grad = create_objective_function(objective, + state_symbols, + unkonwn_input_trajectories, + tuple(), + num_nodes, + interval_value +) + +instance_constraints = ( + y1.func(t0) - 1.0, + y2.func(t0), +) + +bounds = { + y2: (-0.4, np.inf), +} + +prob = Problem( + obj, + obj_grad, + eom, + state_symbols, + num_nodes, + interval_value, + instance_constraints= instance_constraints, + bounds=bounds, +) + +# %% +# Solve the optimization problem. +# Give some rough estimates for the trajectories. + +initial_guess = np.ones(prob.num_free) + +# Find the optimal solution. +for i in range(1): + solution, info = prob.solve(initial_guess) + initial_guess = solution + print(info['status_msg']) + print(f'Objectve is: {info['obj_val']:.8f}, ' + + f'as per the book it is {2.95369916}, so the deviation is: ' + + f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %') +solution1 = solution +# %% +# Plot the optimal state and input trajectories. +prob.plot_trajectories(solution) + +# %% +# Plot the constraint violations. +prob.plot_constraint_violations(solution) + +# %% +# Plot the objective function as a function of optimizer iteration. +prob.plot_objective_value() + +# %% +# Second version of the problem (10.145). +# --------------------------------------- +# This is example 10.145 from Betts' book "Practical Methods for Optimal Control +# Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems. +# It is same as problem 10.144 but a bit reformulated. +# It has two control variables and one additional algebraic equation of motion. +# As opty needs as many state variables as equations of motion, I simply call +# the additional control variable v a state variable. +# It is described in more detail in section 4.14, example 4.11 of the book. +# +# This formulation seems to be more accurate compared to the above, when +# considering the violations of the constraints. +# +# **States** +# :math:`y_1, y_2, v` : state variables +# +# **Controls** +# :math:`u` : control variables +# +# Equations of Motion. +# -------------------- +y1, y2, v = me.dynamicsymbols('y1, y2, v') +u = me.dynamicsymbols('u') + +eom = sm.Matrix([ + -y1.diff() + y2, + -y2.diff(t) + v - y1 + u, + v - (1-y1**2)*y2, +]) +sm.pprint(eom) + +# %% +# Define and Solve the Optimization Problem. +# ------------------------------------------ +state_symbols = (y1, y2, v) +unkonwn_input_trajectories = (u,) + +objective = sm.integrate(u**2 + y1**2 + y2**2, t) +obj, obj_grad = create_objective_function(objective, + state_symbols, + unkonwn_input_trajectories, + tuple(), + num_nodes, + interval_value +) + +instance_constraints = ( + y1.func(t0) - 1.0, + y2.func(t0), +) + +bounds = { + y2: (-0.4, np.inf), +} + +prob = Problem( + obj, + obj_grad, + eom, + state_symbols, + num_nodes, + interval_value, + instance_constraints= instance_constraints, + bounds=bounds, +) + +# %% +# Solve the optimization problem. +# Give some rough estimates for the trajectories. + +initial_guess = np.ones(prob.num_free) + +# Find the optimal solution. +for i in range(1): + solution, info = prob.solve(initial_guess) + initial_guess = solution + print(info['status_msg']) + print(f'Objectve is: {info['obj_val']:.8f}, ' + + f'as per the book it is {2.95369916}, so the deviation is: ' + + f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %') +# %% +# Plot the optimal state and input trajectories. +prob.plot_trajectories(solution) + +# %% +# Plot the constraint violations. +prob.plot_constraint_violations(solution) + +# %% +# Plot the objective function as a function of optimizer iteration. +prob.plot_objective_value() + +# %% +# Plot the Difference between the two Solutions. +# ---------------------------------------------- +diffy1 = solution1[: num_nodes] - solution[: num_nodes] +diffy2 = solution1[num_nodes: 2*num_nodes] - solution[num_nodes: 2*num_nodes] +diffu = solution1[2*num_nodes:] - solution[3*num_nodes:] +times = np.linspace(t0, tf, num_nodes) + +fig, ax = plt.subplots(3, 1, figsize=(7, 4), sharex=True, constrained_layout=True) +ax[0].plot(times, diffy1, label='Delta y1') +ax[0].legend() +ax[1].plot(times, diffy2, label='Delta y2') +ax[1].legend() +ax[2].plot(times, diffu, label='Delta u') +ax[2].legend() +ax[2].set_xlabel('Time') +ax[0].set_title('Differences between the two solutions') +avoid_printing = True +# sphinx_gallery_thumbnail_number = 2 From 4c60b6794ba9c72e122a7dfba0103fd6faf59701 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker <83544146+Peter230655@users.noreply.github.com> Date: Tue, 19 Nov 2024 11:02:36 +0100 Subject: [PATCH 09/12] Delete examples-gallery/plot_non_contiguous_parameter_identification_wang.py --- ...ontiguous_parameter_identification_wang.py | 347 ------------------ 1 file changed, 347 deletions(-) delete mode 100644 examples-gallery/plot_non_contiguous_parameter_identification_wang.py diff --git a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py b/examples-gallery/plot_non_contiguous_parameter_identification_wang.py deleted file mode 100644 index fd96d8c9..00000000 --- a/examples-gallery/plot_non_contiguous_parameter_identification_wang.py +++ /dev/null @@ -1,347 +0,0 @@ -""" -Parameter Identification from Noncontiguous Measurements with Bias and Noise -============================================================================ - -In parameter estimation it is common to collect measurements of a system's -trajectories from distinct experiments. For example, if you are identifying the -parameters of a mass-spring-damper system, you may excite the system with -different initial conditions multiple times and measure the position of the -mass. The data cannot simply be stacked end-to-end in time and the -identification run because the measurement data would be discontinuous between -each measurement trial. - -A workaround in opty is to create a set of differential equations with unique -state variables for each measurement trial that all share the same constant -parameters. You can then identify the parameters from all measurement trials -simultaneously by passing the uncoupled differential equations to opty. -This idea is due to Jason Moore. - -In addition, for each measurement, one may create a set of differential -equations with unique state variables that again share the same constant -parameters. In addition, one adds noise to these equations. (Intuitively this -noise corresponds to a force acting on the system). -This idea is due to Huawei Wang and A. J. van den Bogert. - -So, if one has *number_of_measurements* measurements, and one creates -*number_of_repeats* differential equations for each measurement, one will have -*number_of_measurements* * *number_of_repeats* uncoupled differential equations, -all sharing the same constant parameters. - -Mass-spring-damper-friction Example -=================================== - -The position of a simple system consisting of a mass connected to a fixed point -by a spring and a damper and subject to friction is simulated and recorded as -noisy measurements. The spring constant, the damping coefficient -and the coefficient of friction will be identified. - -**State Variables** - -- :math:`x_{i, j}`: position of the mass of the i-th measurement trial [m] -- :math:`u_{i, j}`: speed of the mass of the i-th measurement trial [m/s] - - -**Parameters** - -- :math:`m`: mass [kg] -- :math:`c`: linear damping coefficient [Ns/m] -- :math:`k`: linear spring constant [N/m] -- :math:`l_0`: natural length of the spring [m] -- :math:`friction`: friction coefficient [N] -- :math:`n_{i, j}`: noise added [N] - -""" -import sympy as sm -import sympy.physics.mechanics as me -import numpy as np -from scipy.integrate import solve_ivp -from opty import Problem -from opty.utils import parse_free -import matplotlib.pyplot as plt - -# %% -# Equations of Motion. -# -------------------- -# -# Basic data may be set here -number_of_measurements = 4 -number_of_repeats = 3 -t0, tf = 0.0, 10.0 -num_nodes = 500 -noise_scale = 1.0 -# -m, c, k, l0, friction = sm.symbols('m, c, k, l0, friction') - -t = me.dynamicsymbols._t -x = sm.Matrix([me.dynamicsymbols([f'x{i}_{j}' for j in range(number_of_repeats)]) - for i in range(number_of_measurements)]) -u = sm.Matrix([me.dynamicsymbols([f'u{i}_{j}' for j in range(number_of_repeats)]) - for i in range(number_of_measurements)]) - -n = sm.Matrix([me.dynamicsymbols([f'n{i}_{j}' for j in range(number_of_repeats)]) - for i in range(number_of_measurements)]) - -xh, uh = me.dynamicsymbols('xh, uh') - -# %% -# Form the equations of motion, I use Kane's method here. -# -def kd_eom_rh(xh, uh, m, c, k, l0, friction): - """" sets up the eoms for the system""" - N = me.ReferenceFrame('N') - O, P = sm.symbols('O, P', cls=me.Point) - O.set_vel(N, 0) - - P.set_pos(O, xh*N.x) - P.set_vel(N, uh*N.x) - bodies = [me.Particle('part', P, m)] - # :math:`\tanh(\alpha \cdot x) \approx sign(x)` for large :math:`\alpha`, and - # is differentiale everywhere. - forces = [(P, -c*uh*N.x - k*(xh-l0)*N.x - friction * sm.tanh(60*uh)*N.x)] - kd = sm.Matrix([uh - xh.diff(t)]) - KM = me.KanesMethod(N, - q_ind=[xh], - u_ind=[uh], - kd_eqs=kd - ) - fr, frstar = KM.kanes_equations(bodies, forces) - rhs = KM.rhs() - return (kd, fr + frstar, rhs) - -kd, fr_frstar, rhs = kd_eom_rh(xh, uh, m, c, k, l0, friction) - -# %% -# Stack the equations such that the first number_of_repeats equations -# belong to the first measurement, the second number_of_repeats equations belong -# to the second measurement, and so on. -# -kd_total = sm.Matrix([]) -eom_total = sm.Matrix([]) -rhs_total = sm.Matrix([]) - -for i in range(number_of_measurements): - for j in range(number_of_repeats): - eom_h = fr_frstar.subs({xh: x[i, j], uh: u[i, j], - uh.diff(t): u[i, j].diff(t), xh.diff(t): u[i, j]}) + \ - sm.Matrix([n[i, j]]) - kd_h = kd.subs({xh: x[i, j], uh: u[i, j]}) - rhs_h = sm.Matrix([rhs[1].subs({xh: x[i, j], uh: u[i, j]})]) - - kd_total = kd_total.col_join(kd_h) - eom_total = eom_total.col_join(eom_h) - rhs_total = rhs_total.col_join(rhs_h) - -eom = kd_total.col_join(eom_total) - -uh =sm.Matrix([u[i, j] for i in range(number_of_measurements) - for j in range(number_of_repeats)]) -rhs = uh.col_join(rhs_total) - -for i in range(number_of_repeats): - sm.pprint(eom[i]) -for i in range(3): - print('.........') -for i in range(number_of_repeats): - sm.pprint(eom[-(number_of_repeats-i)]) - -# %% -# Generate Noisy Measurement -# -------------------------- -# Create 'number_of_measurements' sets of measurements with different initial -# conditions. To get the measurements, the equations of motion are integrated, -# and then noise is added to each point in time of the solution. -# - -states = [x[i, j] for i in range(number_of_measurements) - for j in range(number_of_repeats)] + \ - [u[i, j] for i in range(number_of_measurements) - for j in range(number_of_repeats)] - -parameters = [m, c, k, l0, friction] -par_vals = [1.0, 0.25, 2.0, 1.0, 0.5] - -eval_rhs = sm.lambdify(states + parameters, rhs) - -times = np.linspace(t0, tf, num=num_nodes) - -measurements = [] -# %% -# Integrate the differential equations. If np.random.seed(seed) is used, the -# seed must be changed for every measurement to ensure they are independent. -# -for i in range(number_of_measurements): - seed = 1234*(i+1) - np.random.seed(seed) - - x0 = 3.0*np.random.randn(2*number_of_measurements * number_of_repeats) - sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(), - (t0, tf), x0, t_eval=times, args=(par_vals,)) - - seed = 10000 + 12345*(i+1) - np.random.seed(seed) - measurements.append(sol.y[0, :] + 1.0*np.random.randn(num_nodes)) - -measurements = np.array(measurements) -print('shape of measurement array', measurements.shape) -print('shape of solution array ', sol.y.shape) - -# %% -# Setup the Identification Problem -# -------------------------------- -# -# The goal is to identify the damping coefficient :math:`c`, the spring -# constant :math:`k` and the coefficient of friction :math:`friction`. -# The objective :math:`J` is to minimize the least square -# difference in the optimal simulation as compared to the measurements. If -# some measurement is considered more reliable, its weight :math:`w` may be -# increased relative to the other measurements. -# -# -# objective = :math:`\int_{t_0}^{t_f} \sum_{i=1}^{\text{number_of_measurements}} \left[ \sum_{s=1}^{\text{number_of_repeats}} (w_i (x_{i, s} - x_{i,s}^m)^2 \right] \hspace{2pt} dt` -# -interval_value = (tf - t0) / (num_nodes - 1) - -w =[1 for _ in range(number_of_measurements)] - -nm = number_of_measurements -nr = number_of_repeats -def obj(free): - sum = 0.0 - for i in range(nm): - for j in range(nr): - sum += np.sum((w[i]*(free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - - measurements[i])**2)) - return sum - -def obj_grad(free): - grad = np.zeros_like(free) - for i in range(nm): - for j in range(nr): - grad[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] = 2*w[i]*interval_value*( - free[(i*nr+j)*num_nodes:(i*nr+j+1)*num_nodes] - measurements[i] - ) - return grad - - -# %% -# By not including :math:`c`, :math:`k`, :math:`friction` -# in the parameter map, they will be treated as unknown parameters. -# -par_map = {m: par_vals[0], l0: par_vals[3]} - -bounds = { - c: (0.01, 2.0), - k: (0.1, 10.0), - friction: (0.1, 10.0), -} -# %% -# Set up the known trajectory map. If np.random.seed(seed) is used, the -# seed must be changed for every map to ensure they are idependent. -# noise_scale gives the 'strength' of the noise. -# -known_trajectory_map = {} -for i in range(number_of_measurements): - for j in range(number_of_repeats): - seed = 10000 + 12345*(i+1)*(j+1) - np.random.seed(seed) - known_trajectory_map = (known_trajectory_map | - {n[i, j]: noise_scale*np.random.randn(num_nodes)}) - -# %% -# Set up the problem. -# -problem = Problem( - obj, - obj_grad, - eom, - states, - num_nodes, - interval_value, - known_parameter_map=par_map, - time_symbol=me.dynamicsymbols._t, - integration_method='midpoint', - bounds=bounds, - known_trajectory_map=known_trajectory_map -) -# %% -# This give the sequence of the unknown parameters at the tail of solution. -print(problem.collocator.unknown_parameters) -# %% -# Create an Initial Guess -# ----------------------- -# -# It is reasonable to use the measurements as initial guess for the states -# because they would be available. Here, only the measurements of the position -# are used and the speeds are set to zero. The last values are the guesses -# for :math:`c`, :math:`friction` and :math:`k`, respectively. -# -initial_guess = np.array([]) -for i in range(number_of_measurements): - for j in range(number_of_repeats): - initial_guess = np.concatenate((initial_guess, measurements[i, :])) -initial_guess = np.hstack((initial_guess, - np.zeros(number_of_measurements*number_of_repeats*num_nodes), - [0.1, 3.0, 2.0])) - -# %% -# Solve the Optimization Problem -# ------------------------------ -# -solution, info = problem.solve(initial_guess) -print(info['status_msg']) -print(f'final value of the objective function is {info['obj_val']:.2f}' ) - -# %% -problem.plot_objective_value() - -# %% -problem.plot_constraint_violations(solution) - - -# %% -# The identified parameters are: -# ------------------------------ -# -print(f'Estimate of damping coefficient is {solution[-3]: 1.2f}' + - f' Percentage error is {(solution[-3]-par_vals[1])/solution[-3]*100:1.2f} %') -print(f'Estimate of the spring constant is {solution[-1]: 1.2f}' + - f' Percentage error is {(solution[-1]-par_vals[2])/solution[-1]*100:1.2f} %') -print(f'Estimate of the friction coefficient is {solution[-2]: 1.2f}' - f' Percentage error is {(solution[-2]-par_vals[-1])/solution[-2]*100:1.2f} %') - -# %% -# Plot the measurements and the trajectories calculated -#------------------------------------------------------ -# -sol_parsed, _, _ = parse_free(solution, len(states), 0, num_nodes ) -if number_of_measurements > 1: - fig, ax = plt. subplots(number_of_measurements, 1, - figsize=(6.4, 1.25*number_of_measurements), sharex=True, constrained_layout=True) - - for i in range(number_of_measurements): - ax[i].plot(times, sol_parsed[i*number_of_repeats, :], color='red') - ax[i].plot(times, measurements[i], color='blue', lw=0.5) - ax[i].set_ylabel(f'{states[number_of_repeats* i]}') - ax[-1].set_xlabel('Time [s]') - ax[0].set_title('Trajectories') - -else: - fig, ax = plt. subplots(1, 1, figsize=(6.4, 1.25)) - ax.plot(times, sol_parsed[0, :], color='red') - ax.plot(times, measurements[0], color='blue', lw=0.5) - ax.set_ylabel(f'{states[0]}') - ax.set_xlabel('Time [s]') - ax.set_title('Trajectories') - - -prevent_output = True - -# %% -#Plot the Trajectories -# -problem.plot_trajectories(solution) -# %% -# -# sphinx_gallery_thumbnail_number = 3 -# -plt.show() From 936aebe6a887bdd5c4d4add11f5f935f88d7d229 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Tue, 19 Nov 2024 11:42:27 +0100 Subject: [PATCH 10/12] changed underline of a title --- examples-gallery/plot_betts_10_144_145.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples-gallery/plot_betts_10_144_145.py b/examples-gallery/plot_betts_10_144_145.py index 176bdbe3..6c626b63 100644 --- a/examples-gallery/plot_betts_10_144_145.py +++ b/examples-gallery/plot_betts_10_144_145.py @@ -40,7 +40,7 @@ # %% # Define and Solve the Optimization Problem. -# ----------------------------------------- +# ------------------------------------------ num_nodes = 2001 iterations = 1000 From 11e0768f85efb229a315bc284323527f6235f987 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker Date: Sun, 24 Nov 2024 09:40:40 +0100 Subject: [PATCH 11/12] example 10.1 of Betts' book --- examples-gallery/betts_10_1_opty_solution.npy | Bin 0 -> 33784 bytes examples-gallery/plot_betts_10_1.py | 165 ++++++++++++++++++ 2 files changed, 165 insertions(+) create mode 100644 examples-gallery/betts_10_1_opty_solution.npy create mode 100644 examples-gallery/plot_betts_10_1.py diff --git a/examples-gallery/betts_10_1_opty_solution.npy b/examples-gallery/betts_10_1_opty_solution.npy new file mode 100644 index 0000000000000000000000000000000000000000..6bc6751af7d9ae1975b6cfb555619b022d9ab243 GIT binary patch literal 33784 zcmbTd`8!qr7dDOvA(1i_8jz^Bl9Xg2q9{a@2AR@i%sgZ!^E}VFzUiW>kSLn^#H`MM^Qdm*g3+fxZ)-xBBG#~zCR z9Wq`cdr#^!4Z@XC_kHZgyTQME?`|H!4#+>)c~a^|6VwGh819m(h2h8#dd(JPQ0{lf zD=nr7?r_U;#am@TivFV~mK(`%Z;$8pmLdR>SM7u9&DPrEqM@9ZS2;%7#^YB75pFP6e~KSp7JSILf}=m~`OvhYCKZ6$zhy@hTLN%;ipb>18;E;Idw#B~1R`64UX!S1Ahuk1 zR+J?fh!Rnvot`v-SpPJH&8i{*C3!P11!@K05$W?f8VmmTr035>ftf!l>WxZJEczjn zl*BjQCw>^U{)6_}dtX$RA4sb_>x;{h-{m|)eekL4EPKO}H*%cZpS6C)8^4vG&^!_H z7Sk3kPz?=up@YWM>^iF#9=IGCQzYYw%{K$R%pZB+#M9|xdHU{{axaG|Tg?qUt=J^A zxLmPoCux+S!WqxZd^bp-g^?K-;a_|T{UKz=V2+9Hf9%o zeO-bT)6Q?gatrWUc{&X(;qOr4sRK}h# z(qu~3S`zs<)t$cQ9<#-r|(XmQKUk$;9(G%Z1vV+ieQTngtr2u?%?x__G zr!Ou%=*rCN@j}j)kBsS$-I3x>dw4^MGhU0iXdM}Di}EIEwCb{E*r?z0wJGZf`enrg zftNTeo@c0W6?z6f3bGP=?QI~Jf2x_7?*{GPFYb@v^8?JXB(UuZ0nxqNzeTQx!+eXX zvP60$ForwnT8Kx(w^Kg0oE$OWWAKJC;aLpSZ_KqDNXLLs(c!U*yl4=vinG7f7zLdx z+~dmuk+9G8q1iNN1gLZ$L9bU~u*--^!P-9rrfzsl&btIb=LXeOa;zV`;OqV`bJ`2q zr3`a3=Uibf+gxx|$sW$VaR_NRY6^Q~g|c)s)M0N$!iU_)Pf#@7=XU&EOQca_Ea`anJ1X!oPG+2)74(Q0RMt4968e^y~O0;r-kfE%Q$is^Z_`$NAd|FZtb(XMd?y zo3aDW77L`E&v=1zGO6$TMQ;Id4oY(W+Q1i!8J1_UUZ81pC_Bv|7&H}!&ln6vz&^Qc zVxe{{_;&RwshcFg*9kA@U1o_Ot;T3@dpsVv@*W+pcpVK+y>IS~GljrN14?K(I03z~ z!kx%wbKL!l&vKGE6fHUz-mEDk;|qc}#U8fW!iRFU#4 z@y)!)dPHP7vcErYb@EUde&g~V`|+57Mq%Dw`A$WsP&xni9ps|hfj+|{Rp}^fnD9X? zJqf9kN~8X%Mxo*!juYRTgHTqvi7PJ81Nr{dKOTH)hk-+Z`o%`3c({*3b8X`VX5VeP zqhM=>R0U_aoTBYfoIv2Avv5b_dp25I6W+LQyPC__&KHS$%nRn8zC~-HY-JG@SG1uO zS1tKqjrBbZfhxydB55sGXvz9M_Fs37cq?!Tf1M3@u)uHz-dz!8PYS#OXXsmfPT^H> zxaic?WG)Xf>qn#HTBTs)h?>}0c|l<0I6hHfcLevgOH)iS0G@eWB7E$$G9G>1J7wsi zfy{>#YtHN6$GEq}D>D?lZ;oFcn&bL4cQl9BIFajRmCR;=Has_1 zD2FZZ1*!b{*=S3&ZQE7cb{0Zf%Jswp_Gp*#pKyYN9f}QmB};JHV%eG9`|GpFzb}rhA4z6(abOXtmIw<+N-6nPc#CH+|L3-@3qk zYT&ZdDGM+W9yGC$c7{CWGhz(1UhvI|*R$_b0My)w-#Azl0+ox;9&PYMK$OgBR=rcv zFn)jVWO`yO+&HkyH9$Nb_KO^+D9DS4Bf%lB4>!fZpI36V3$ih=ZXCs#xD^g(<{t6z zrw4+HP}*~2F(=5bVe21zr;X29vS{^ke6S#@F;;Fh8p{mnE1p%SqAhF6x6k!CNFNuh zVlY>P-dSaKF;9r7c>uQ7b4d7#)xA@kP>Jv{SD5>04VE9(JvBo>91$8bzt>iWOjnva zWbNwFWj9Yn>0&*4jF-A*9c{pmaqYdsTn%`Q;oo@GmwNoudG_pn{d$yVyl6j+b;uad zCcHR}cz12WAdR^egI(4Zb$3^z=RTd(vY`qTXPCT8Odw$z{RI=(N2M5hKIqP`;{>#z z(9%<+C`O|T7blBW^6}Rl*4<}Ba`4ar?)=rt43zl2BC?j2ilWaB_IWHOp*D#lTU;m} zr99HD~k3OJ&Wa%RA7q~T< z+4CRxtr)s}-Mj$?7HRGOT(qE2@4^8;89k6Zn69yA^&I3gT!$VHn?RD0>e)^4*Py2; zGd<{M1#?$?<>Z%ap>evh-BQyL5|_^N**tUx<>n6gQ+r(Dt#X0bxk5MCEInF0an%FZ z4a`q+y!Ql(wkeaV7vF+N%ORC_Z@oeM6&$-!>I0_q77RL-zRqec zfg%93UW!z(68wRBW>mcHjz2^^N$}2H^aIIaHQiNPKe$UD6MA9W7nUB>SZ3My!uznj z8LJCEu#co*prYjimsJdDUe$TS=g|ZnHy&@u7K&OIfBhB=OVgzStGz(m%}{G;!xQGG zE?X=d@q`me-o|Ind%#=Pu((tIxr5R{m+mYYHyFLFV#ZnI0!5cZ)!rR(hW`!~i%^F- zz}4M!tkHrt;Jruo3|pr;9FCg(UY-98%3rcG))?M`fA$b0YXZ2cf7ZVdS`V7|#4&Iq{V_3Ny6bwL-mPQE8KS z@JCw$^8Vg^sBs_(r*qnTnwV0sY|@dY<7Wy6`S&@zuuR1(Ki|hyIHclVb@9aat10-g z<*JxXT{0GTo(*eBNJNRt2Y1>J#$a3`qsW!E5PbcphcqGShOZrK>n`uV3mx>Kd3#I! zz`#X0Nq$#6#Mn7@QYU7>*FUiK_-P?X|Kaj@e~$>~>G^unPga2IH<|EVaO^v#j^j zt{!L(_dE*v)(C7$?Aetk(MQtNsyQ@0K6P^Zf{UW|z+XPWcGz#tgwnqdr2|&7)sFzxoKv zmkge39{UIsVHeb_9lGFN+?zTsmrl4t`GaTwgAd?ZQ}y;7`Ev?qzk*#(JG54`(tkbm z4$hnw8xu%rh39)Ry1p>AK(VCQ8(Qlo$e2;#x;)+hj5qwtjPKUNawc_hcOe3ocyD9a z$7=Xn^wV?eU?oswP%K$^k)WU7y@a2?6b`H2k(ym80rLQ@_m)dV@YDGFht+ch5PZ03 z*?Kky>QU#u=dVoIz9l)oa3>wg4&4rOqfUiP9P;pHElt0wFNxHZQ*WNF|$3jE#xUz-=9Vs zSovo9>HJe0aMb!}(Y$UAe%Z4t7AIpi@^$N54+eEp0)(l+2CXi zF$>^fq@HzSe+?>KYC%+kX3#Vo_v3@DDU2U|sP=fl2=-mPW*Hj)0*DtLCEukr1ka~E z(={jcV9ZkR+t~icaG{oFPp;Z+P!V|1=yOsDj3$rrE}9&J7X5~*kh{|O{NASkmI-zI z!TEdqvf)Eqe|zM2*WV$qZQ)o$?xZEYRdQ%cn&r zE4=SqadW)d24_F@E#B?2!^9~1eS9Sjc!R2K$4lo8PKTZxh#PXk!orae;UH%u>>QKt z;&8$E@#`Yq-YytCCDuGJ;ez_^*99_8xgzNU-`U-UuDGCe;owTBD_*msnP(@t;)vs; z*|u(1ytnLeb!NsDtJIDRZ~S(}P`bQ1fjw^cswMibF1s5xjg5EeA9KTNjta`K+iZo~U!fLnk=a9t}8QXaK3tXTlEIqX)7bYiu!|Fr|?=;efURmb)So zTye{tTYz(~50?29w=jGTLB740-vk)OqUwP?TFLFHc(VI&dcHy)nlt^2^jjz3?qIua zm60lxdGv6!MX?dDzDo98HE&1zn&cp^UER2m$j8O1`xyts?5=!xF@ki-@0Nd2P9guF z!nR1=MPv>f`kCqW6N8xz^e-xG;{VmJ{uh@)F_0qoI9m_SF5*JIynF;{}=0fp>$qY3C^_^!!KET zijn5Yng2#7i!f%u=s{dtA=bQ1KYk!BA9Wv}rn;`1i@Oy5bgt55W9ISOPt}7n@M!z! zwc_1rIAV8>x#(FkGPbaP>b8i-ZLBF1*onmP*mCW|wjsDYlI)QB#T%RJrkl@2IO5Bf zXA>T-KgJI=%PgfgUqU)(?K0(g7l>2k5>Y$v3n!l!N4(by0o9V_*Xzp>up&--u~Q}% z#=ZsU@n|GM>g`G`gRo?%4*fk@p_2;gk?NIEtm*Ll^NoZvSsCCbA@Q{JSQbQX$h=JT z%!YqHQH2cMIZ#tjUoJ_T2fl_+yzdF-LqcSMvV~Xy98x)|`9`D=Zv4pDn@nE>2M6EH z{bx}GOy!eKNBfI_>J#I;;QxxjY%Aa5wR15%Vz&Gf(pC&gls%H(v?Xw@dxn!BRRTGx zZ8AcJB~UeU_-w6535dJOWbchEfm#L*zv}1`sOSrcqw+6-{y?c)D`q9|OMfll!;KQ) z$|qRtJzWASlh-&0Pn7@+-2TF3Uji{8b*I|xaPd6Y~8 z-$7uA?S>ooI8;|m1pNt|0=mO82?^)spjoui-+yHhvL~kbg!-4^oNdtc*J-P8<{4x8 zi`#2JaoxG~;?O!={8k#wBD4XQ>AUZ@ifuyu-PU`$tDEpwPjfin@fO^dP0?Po*aGT~ z>G|{GTX2D6{=HbrCj5KfCgoAT0k(Fn8CPR}!Jkeg#ej?J5E6V$Zy@t0r0qUgs7hqpzQR;6p;bb3avQ&%NiuSkCX z&8Pw*GcJFBaI_ps=NJB~d0z&ZacxC+o|HmffX(4weFO-KKmKU>L+UpEFACtTJeZC2R%$xGYV&Yquq3d&rx06XzD6pE zko7}Ka=hb)7$;W%`ZaxfuvN2{NxKxon1Np}OE+NkGE_@Y>Y)WW62-S4y2p}qlGY3!cb(wnfXaf<)p>?Ran z+9k()a|_DE4Z4<|Z-LL?($hkzEfDnULFqQ~`+llvI%VG`%zF6;*Ku!vdbgM}>uC((Pi)@qw^SSSK%hl!Z_E7V`%WmbcD#R^$`+7O} z(!~h`vz7z zWY<;Is??N$zrp8ej<7PgaJlVs)~hloxqr=k?rIrWg>YW%I9dkA-G+r1D9Rw#dBc(K zdnwSeXtZ$>1-9QIE|l)2P&{>Wo99_6L^@3QKe|EA>l2-F#Y*A&HgWd` zS1I^1o)$U1s}vZW|5ZpY5J75G_UOTWBE*JI>Lk|@K{ENv(ED8Sxwew9dn6H_>wTUr zb|V6f$N@tzB|?QYL51xR5$<>9o;-M+2tJuZlY$aN*t*O^DRi6&H$~gZ*Y**?#qjpj zlWhWUdRwPR%@QE_Om5uGUIGlxuN~E>CBSms`uD920vPZ4{dC-)0My^FwI`Sp;DHY_ zAybn8>LAu@AVvVyo)zK$_7WgCy3JH|q6CJ!6(4exmcaC2jTVJ786WIMj0faPAX~OU z`w(>r6w$?TIbt#FKewRlZC6azoAgEkFBOAOPS|@F#$uSZ$Vl<|L-x5Qga2-D6~nEX zVg2)D+@j9%Ja|&B1a`+hA566-Kq;j?6!;b00q>*aOrnhB9SJ z{m}4bVvFVD7btEPs1C0l20zlL+6n4WDCA@{@V+|^e!4r_(_NFW?@T4lH-j1YK=rx! zJ@Y(hR0|s|H7tT@pfdm0i_4JHINQSIy8?Uv-P!oBYZWMu9PPSDu?AJ&Y4sC1)?tU4 zW=eejFOXIIXQ?~?3mS&0rX}JxKvw2TKp_7neDU$WeaLkadOIjyS{H8uv(g%`0AUjZ zp1$aa^4f%ih@6h2qMJ}cwSUh-!3M0E9?kr^djlROJu7s(^a~#O@62qet%HiJ?p(d{ z8i?nn)nthNgp6dj+cDIu5P5x=Dg6C1TpK!Bw8ZrTxK1e_8ID|liHP4xRUC8Rmg!Iv zkvG}9`*HO>P*$2e_{aDTnAXatT}PTh{kx#rQ=LXg)p&ou*1jID zUs+;EJWvO&!mfcL0=4iqd^Na^SPg#~-go`WtAdK$N)5Lkm0(V3?7wbM38y}|IoDED z!q(@*Y&*mXsO)HXy2ri(g#QM=3K~ZNo%wpsj$6J7cP6XUpL{k>_}D zM>!08`tn5Nl|!ob#~Z~y851{`wkMUYu*mp zA@g9@T<4_(NfN}H-0-m(D}%?}Uw+}6GT7FpU#!_%1~K~aZl8infufEs%8a%YIF|N3 zK5aq-kzy_;LK^{Ih{dG0?k7P1%5=vqwGz0q&0YH7O)=bWOA}yBC<0eL(g6eQLMSse zzWmlOALivw*wx<7h5Q#iiLt`jz@5F=%RHX}^oB2(jUvc+bzLI+m1qj2ZoZ{?Zj%T} zCuR1P+QmTH54w&vlW;IG{*&eOBnTwMG7nbfc|z=^L3wH;TTobE4bus{3zWvLX}=_2 zV3{fmQeY<+57e_w!c<*>Sh^)B95E1uS~*_yQ6yMu4oj_Qa;kyB*mgsyP$T* z?nLB&zC69SHyNdSDWVmfQjluma0T?GV#{BphXpL@DEF__;)QqyJ_%v=1ldgdoWA*P zLOcumAAAm6-=B@v;|G`KhO&`sDnfugAO~64seDpT=3+_p{K=4Fa^F6$?$CTR4{uj6 z&Fy=Whu_)dd+s#mp^OmkM8#?z5@?=3pJmL)D{oh&g7@d6q^xn?8ht({$1CfXuIFKT zRQM+8dmfgZH%R;2n}@D1IMa?!=V4p+Q4ML4e7xZvs`9ZW9}i4Xyk#a8;IYvUy1|@9 z=sDm(|Bt^IE7cgdIs8k|UnQx!TZ@P{^)D;mjV;6E7hg*w9#r7o|2AkFv#W7+^>0b< z^*Xe&S{AX6YQWS7bB&RT%{WPSD2+z14WHAzwDtJafid^&G_`0yV$P-SdtW~6#-6J* zC9;%#m@O~-BH`}8YE9wY zIjRV5hFMHsx%&YZ=g}j~%IoXDB{cggwd8bY1%uL*?!>mQ;*`UgxM+_x3>j2(tylVm z9FKTpr}u7R*{}AxfbK248mw@4Z@?eq@Lp{#y0nev4n*wh`-fKW;lx{EiaLhdiUCK4 zDe78`UraxAq^$eXw);Ep9;!OwHTv{u7pgj4qm{AiV^nolhCTZjgm={q(JNT)d$6mH z>Hd7uoY}6rQ0rhNACp~mjMa6imo#_P?RjtPoyxbX?$f+JV`U#zU7BjNccTSWodb)p zw$ct|oj$d#AA>bzUGH@!JBn`KmP@OykN7sQ_D`g|tJ6BRXOCSB9{-7Xe&Z%Hs;d~kZ|aRt z>oSVWJhKvuS;G1b#i88hc|2rtclN@8S=72^A>Zdag}#4YN;e5lpgYI#v4mfvcu20S zuXgY|*4ncU{^a|Lf&6O!KHz7R-nDwotgH{^0&@dy?&?AHfhQfJX&;e9*xUVFCbd1LjqPOcmoj~N`&R4qptwf3BQ_shvX>P(92lX5)AE_$v?uN({XZOr+|k>jm@ z5B=kE{E}TlwMVlYTbe|vQNA2oPyI;j`>!17GL&{R9Vo{w@!2@ZH4^IC#0}l}K*Du1 zX%5XK5{{78*&U2X820C$D5p3HStmHk(l^Rb_RxzLO{HZRCh6%bV_JrW?&JK7C(7{r zri-l9w^DSg^b4a3E5$F>r<~HRmmNWl@uLJ37Q3snwMvej;-oPeB3hM~SzB|FeYb1ZI=_e#(bV{jo7goXK2KJrA?gt^ zrsIFlTbXZtzcG=B)M@+7O6rN|V*ARTW|oMr7+2>Hm6DL_$@4KquPXdP^V9g{={nrF z$TYz%-iUNJ=qm4SwczFNVW(8mJ1{`y;T8+&BT~t*HpwXW;)d{~*oBJEcyzGj*Qoe6 z>}1zCzu+>43u`W2nb|00FV&uW|BFU&^Z z(bt#keXU<0`Y4;Vab_Px`@727NOwd2+Hd!`woVucRnk4as~wgUej9P`ZUNunFWfH0 z4Y2>1(>nAa6uoBKkUd`o0)H2Kf}fBe#pqLb_E{n%z8)R0Eh>WdreREdyakYZnkq%C zHxHJ4%vUIl9@y)IciZ~p!|#f0)^k1uP@=qI6jE6T=A^wfOddrbx7H9H zAYKg5=UEW{YIH5~223kO&t8nRhIF z=rDFIh3-Kau7RmiI4UwAyI-OVnkR~M{mjdt`+@8y!{{>jPzKH=;Q1Tn1@n^_*AlltJ)DXpkzI_xL}Y$v;HCUW_(r^KNP>a7X_ZwS7Ayh2J~S^FHnic$$^UBAD+MNtCQ)AW0IYKq~~T9(W`$70~is#?^O zD~7EWGv1?2#b9^rctqy6A_zhL&+6qxP{%{PcP_RFME+PrkNFnCQ{(;AZRE&!^XpaC zRKIu$n3Bji5T`h>`qd{`^YN4dwmPJa(f?CbZU70CSmOx>VdSTAJ8 z?5w>M_za_58THmfLoiT&%+Qd`lk@MKxXa@+3cq~ancf^9hmzShAF13YfW`4d^NZ+7 zU~YbHl6`Lq`i$kdbt{858 zXW_gxP(q%|$^19rScIX!1s{B+EWk@Q+U97T@^JJ@E9aw=d6@iJd}ga55B0&}`7xFP z>>K#OY%y7gpN(4zS{RB^c{p-aldKncYrhp`@FE~fi_?jy4kDI3yV+htQHJuhmz%^` zNZ9vd>Ewq#5;_=)S_tni#~yiOt{a8q`LxdJXP+9&Fe$i1Rn@c+}E+psu2rERUVjoTk`dr+2Q({1QF$2F3@0TEoC8L!&g;U$!C`^o`SRA4DMcRrd zN-{%FQPi?4MCO7M?73GHN2?nG89qJMI&E=~w9{S2XORl)((wvu8d;!lr6uj+S{__d zkjMy1DFUm)8~;s4m4KZyaWzVn2$?vo|A<%$w|zsLiN$4L=VD0wbD0Fg&9r%*W+Zs! z9>^QwK!QeUNt%BFB>3^4W+inmS)VExyOL>L4%EYczhX`WaJi0}Z!0v3P6Vp8bik+5elZojk6^d>$?%~}4}1?*l5jEV z2l9ZCky-y2IG{25Kc_!Tu&~}-U=qf5vAol~ zIR(rEe?JDDnucS(j|Drbrpa|lCt31y8s2>I^5f~5h9sipHtU0F__z1Xr#P~%8%^oH z6ly&Qv$vKEE+3izb|tEZy=1*Cb71Jq&btvf__-*-bZ{7|HzczRCmxd3ocvAY_dC^? zh@G!egQ{8xsHwN#K>I)m{$~*DEKOO68m*O-@6xkSG-maQ?(by0W0&h<@+J~lROX7k zk9lK6JAZfT@fY|@IbHdplsV|Bzo2{$o^X#~_~Bz}5Qv&(=zdWShg*pSlhrcOAgF2L zwm&Tn!c$D@HiHu3DwE19pS{Vj=580;C6WSF+&=@3oJ@s4;q>eu3~3;##a0&7p9YFM zZ;hP&(}7yKIyh&dsxNnM_xTu^3(cc>m%a>3-|c~5A2bh-q@ z9)(%UZxUc*D2TA7R|=j{5tq+iAwi=3k+PYB3h;O0qcg{9Xs^*+Yf-2J#_S;3_0R^$ zs3z7u)o2E<<8R~IecNE7t0H$*?mgJWM4KFv{s=id_YZuH?FL@IzjlNReIQ#?Vt8=h zXYjJQ_{^Vc2vklL#yQ;h23fIhPey(ifhROSBcp7`AYC=b*X_szWGrYN^QxMJI|75j zb(g2%+dz;|X~GP|nn;IM_RYe>#Jk&yLv!%_bkbST?0IlBbLVZ9T!3Wx&V_Zq1=xHZ zYATbk0J3&}S!<0Kpl+Xg4K2k2q=ai6Idpp-@)_OyFPYDQi=}8{^}|{4rFE1T;+lc{ zduMOD5vIUJ=;a;G)05yVWlJF!G7e(zKJn*Hk3v+iV)TIc2>6TCOnjpqhS?PveUaiX z&_rKL2z=cS&R1vc?eF)1XUHNW>E%ZtKBb)v*6kqUn487;;}*E%vFQG?rydfcoGp92 zs-X3(Z^wCEB7_TWJ-nx%2VFA(Ij!H4p-SP`&qnr{5f73R+#jd?3xClp@aYkVt@JfBa^P|x4$kG|~o5(NPvIKc9^>x4!GGQN?l z@fL}}L0ui+lLr$pDfIQLXqsdk%ws>QyekzK!_Hsup-e}mbBs>$+ZmW>CYeWE%R<$z z+2ix$ImnbJlBe93he3)0RDB5r*fN^h2xMQEO^Ac`C>g&^1hSb!J{P0j4+aLV=n@RI zIX9%COu!bkq`fPP1hn(`kGWaNQW#}`NNvba zQ%VwGevdrwOiLz5x{#_xoqLV58;5xbu~)r%(VE4riT2U}Y8!o-VU-v})9=ryt~(54 zO*RL=h{p(etdyyjn~WiwB##%3!vy{^V%)fuIfb3nfrD*}GZ?o>=WTd>9-l^+ic^yo zacp;TlJLQ0%mP+vX7v?(%qkU@?X`+=A>(ovN`GP>L$B{x?;6s48m<~2UB?qWj$a4^ zzp!1jvRAld12@JS9aD`qagjgTn~ru1U7AY^mh87so!|Igbo&+#5pSjLtZZQd8;5oW zIqqE0YSkd;a|jp`=CFl71K!OX*|UW)H$OFST5h6DSGD;5x(&SZXPTFB=@+(QCE+&J zFMN_!bk~_;9jm*dZz(SPL~i<`GKY7o7&cJ0wH~{Ig%7h6Z4H;vgd#cPq4pB8kmRh` zm*>$_xFPFMpb-#|NVN zo)-4i_}iFevEx`7UMUG)^CT3ZhDDcj+_@}F<*i@17?6Z~6Cx++Qp51uXSUG%r+D=u$Sc0SD`Yi%VCWfwMJr#!$}%xOZ1(-sp3Jv#}~pMDn^Z zH+^fS#|wBh^7w8{h#-0&%ll!IVTz=KF`TYV&UpBkSjc{PZxnc1UZYqYgv2l2U&eyM zut++New!f*7c%a?nLZYaid4G|UHap3R_ZG4qy0%}wk5V9&ym@*fJiax2Rk^bbRAfzI%;CI(ElPP%M%fIG0hZNiq zQ!R6)NkLCFNtY?7WZWsds$VsqgcM5IS#wHBc*LGNwd8an)@iHr59P<>)fOY!6U=c~ z!kqQN=T;0d)(@*3u!%zA$;ULiog=)vU27hz~pQF>J9VhB3V)2y?74nz)S&M$8C zzL=LZx)p0}$Ab%uVE`O7@=o{&WwUl}m&3$?br@{xofcvrY6y`U2c+(Z1zr#8aj?odzr zH=`&>WHhMksEYxEBM)wf&Bwv}YDS)|nFQEoSW15JAPJ~e>K@4mB|{_dPrMpVhV~%p z_HFAF`1(pvI%O*bnl1j-rs}1F*WoqS9#Sggb`EhGZKgtcUaaA=U>YnWmY4ogN`vl! zyN$Gu(||VZ#9@o)X~3I#S5x6#pM|(nZz>cP>+QXKy(%$D}~^=PD(m`zb)%`QnQv zYYH5e_+2E|lniJ452>mdCxh^eYxZyYWWd5^vRW1gKcs#Q?x;inOHtua!yiA0;Ivrp zW;ci9BbLys?1DAk*MC#%#Gs!=8vCA$*|^I5P4r9NINW%(HLBzQayI zSGzyr)#I*mgp2_cIpSPx`E&&Te!g|)x62gDSbqziPgz2t!NiaHWo!71{hN90;2*TK z*`jGpqpEuuuyym&DO6#C3JlW2x8yu$mO z#G3*1Q5Z2az1D*kGta#Bqq@-W^W;Ido_BbV^`XEBT{9}UKDzszww~P2)a@EpufkOu zvuM@GQXDw4C&eqL2r~v=E?%n5!rx3hS_6S#>88WZf6pKjp zbIgHr2P-tRlCmNFMfVl{f-Lwu*q?ByAQP;2_=rK78L&b9A{mm?A->)wy(B3OxGO9g zrn6IFX?J+?9ZUfuhRB2G#*#tMS6-`wH5phg&AWb7OM)jGvU#!+i4gqx^}f-Gcu4RI zD(K0IgP4ZvKiynnp*lj0G0=$Yua#N-{_-*!0w))JQtYEZ_-sUfKwKoS6Zi#cIwPP{ zxAc1-YXrO^&2I3$2nXt+W$&c^F!((7PWSACFnBS1E4T7bDAXxGX4R?*1ty`q_Tr3C zkgE-hI$swG^?CP>1+#@gPx-Sjs^T!n{UOCcN)HEi^U=d!Xd)rQGfkRnIts>ZO64@4 z#KOy+h0A7NU?cbpx3Gw|6S-#{UX$lE&ST2-W(gi zY2bMGAc^@X}))x3$^dNUXKj+h_k4Sr;80i101D6(Sg+gVU z@#xdzCv|$OQA^X|ckH7AtmN7t#qM&!byD=FE%Lk*JwF$9%{3ytA6A-V>8^%SZ&m{d z_ja@B{TCDV>T_P4k>TN1Cd9*K!04T`$DBvd%+Jb%K$weDaJZqu9}t@6)8#@M7&nIF@! zMnt|xAtf16CO2QxCkl6SPkuD?@WDU!6eUmP^f80<>CR-9BgFY9uog!LK|t4qbY0eH zNOoiSQD>Y4WvvgT&r_s9c=r?(N@s$8l<4<8Cv!mV_>97#`aBT2TDjazSqK&hf%PtA zo=a;SO7r?+G4%D670G@mhRh6q=XsG5aE%iowrQ3?+{H@H7c%6z!7&jPvpMp*9nYez zN-?}*y&+_*T?F|DX6$!eDu7|$szK>1P4)LeEly$Oze9)S`qf$^Oek~sDSmM;FN!*WPMrjNCVH=3K)D&G6{5y6v= zI`z>9bUWMhH&&+&#dj^p5qBCeEjW3SX`~Vlt93DdOD)7nhcC^l(_u(0G|t6*DFo6# zsu;Ulw^b-ex9{DXRTvBtC%4E>MiAr%h^q+* z1cKX7DNoV80btP5pek(a5A?OPvz{mXA*;uJ^|pXN%;o>=p1$o5JbmxKq-Og=BI|rQ zFM9w?+fDQvMg+hlZj_oG2?Q(NniKhiK*)XcXLR6r5HKq1#M+t$!9VGC*-uqLa9B}7 z{o-!&KGWri1zFW#sMr_#H!UIcgZJ(v&2?ctF-T`{@{`#!hD=UY>L!t9HjlWY@DE#v7-Rnpn27Aq^Up>1L z1}ewi{1@*S215%EDP79LfT^vq)nhFTr~+@$yglp^uKw)Cz=2X zC(O8&DiYzn|F@;%O-b-vWa}&RCWA?O*}BO7REWBnxJ>;v4JsDggyJq`0Fpv|=p?h? z|5WytaalayN z^ZOMMUpKRSxlmLDIf6~=H?WifXL1vMmHrH`R{aqy++P7n_ILQF`zqn<=fq}xNn*Xr z&mQ%o|Gq$IKmD=2Y}KIZG1O%1`vo}Pbgi(yQw0lpzH$j~D_|2HS9U1fXHb*(5RvLC z2J@>I)A^1SK*PLqwZ3BxToqdJr-L~K9^K^Q;5CSZu`eON2aO!z@Gdd7Ne)NMW&h)8 z(H)K~=TEVn9#6v6w>swK0y6NSQkd?nP%dV1=pH?CCm*>coy}fuC`1J=_lD;?ituUR zVZQ0WB4pmbZphKH2!Cn!7Kk(!Vja`blxr~s*ggK{O#AwL>@_{Z623JTH*UsfieEF4 z=iFhAPWn`Al$DY`GLeA6AtT2dQzM8tptd3AvOoU(RZFLE;tj@Uj+O8G`51KRN*nuI zTtP_QN%_}YAn>sKJG$`aBRGql^3x-}pE?)3a^vKQ2bLR^syc@fK$zZCin$~K!owKE zxpEVrsQAAZA$t=b`trT5Q0aJhMHkcKRuBV>LCxa5zLCUu+KsBmdp^QuhP3W0@xf5M z&P&2%e*oNBSAMxL+Y>e?sXlNsas~zZvRDD7_h1ra`GeQn0$vxt5wCc12lx-zhz{&K zjW^crZ@h5D2zkG8%jYE${bXwvZ`b%YsBl7I)iSWhRm{Or(&K z7mLSEeK6nrE*9mJ`Tuqkr2l+$!|VI8czV6-8une*btVu$`+#YF% zU{wCeXO+>Rn3x`O>mOmKiPHzZiLZ&kz}(S;;zyz|+NtH|bVm&S&G)_j%sL+LUGp6* z(oe((ohtV@qd(#4<;p&)mr`+M-md{SlMJ*w3mT=Ha&RJHpnvt&Japn{*tq;f0oIM{ z@Og0;<8;O@ORI@e%s*RemCRj%!qUdKxgUMO^t`Oft;GLxMfWqoh--Byv%BJH*0BbB zd?m!t=0GD_NS*$aDb|F$N{ekGRyE^5PU~T#`WCFQ{LG~Pq!m~1I&1Z;_b1Yd9oGGR z--ZK!=^9VXw_)S}t9q$oJ9hJ5ym|S*c3iQenEf^(uX{^>w|lA$*_uZBPnopggKWFP z|9pR<%-z<>p1vR0{^m=&aaS{b=TzVT?n?E*nbBXI37@jXJL<>lVFqYpsfiCZ57lQj})cs#n+d0;n=$2?Ur0(U4hLpdl%^h;`>rV^S!%|194@1 z3*SPMAKux)yPsvJH~zPs*5Taijtwf>C8Bh$xbxi&1s6jnthuyrXr+_`M)Hi`6}kNZ z-#0ny>UX@sp9zkcHEGt!88FLJ^W74Eb{XFBmnOdVB)<)fwRnL7L$cS~zdXjQ%uCNt zDj6cnF|be#yn{bXO!!x2YvEzJsXGtfoI@|VI$oCvPHgv`S6*dv2KX4*>cxz;U}}Y% zWBkjzpmPAs^qxHeIXjqZ?tKc?e~iMt?aV+!^?T7XX-nwkq-VERum&`_leEp?4R~cn zTb)gM4<(xpJq+it2Q~)X6n;fV@YuJ%x_p&0+}6Di`n%Z$#xy>)mnISKwclD9j2%6| z>2H*4nV}bi-|RHMd)@~|91Q&%_WMDY;=9)ED+A!dgrJ4}zW}h$Z+<4y5(xEXJJ*^L z{+A;q<2|e_yAnym9A8C>Z$G zuA1`^47m0UK7er8n)4Ku62jqLPfAxAX9Og?6wa2_jR1z@?^U>>Bfzlm$XUM52%w7| zSCM6ngd^wUTtDoHgkParbLSuu{>*y*>`{q?7nA+BqE18-{ZRK8s|i_*-(lW>771Ba zgDKngMgse2VJA1YDiDQN_s{0pST7R-IS~hp9X%_V4ZCa34#EAK1Y8g3D&CR4zv^}!^&NcUrZ0A0$cW&<|&VKV!g0idZ$SSeEg4pZRpMn zC|o@#$IG1tRhs5@f5{d6ip~PQVgN^iDDkb3W{GTwo>e(_>1BY#W9Lg6z3=|Yb5K*NRf#m z%m2u;)b-yfno;zm^!I|&Ye*G&zkvdhf1*ilGbMR?BXwN6Y*D|**Oe|(sj+C0ir4cO z*&Lj`$PMgii#+r!agoadVirl?5wS>R>5q#%pc}f#WqU#v8KN7!$l11_MgDtDaz=&Z zx-BGs(36xN3|zd<_$SHQdXhqwB!f#xel8$M%OhEZBri~0l9$deUBC4ErTZ=2f9ZKk z&u^j1!$FmAoGNcWRsIgDKKH5m#YvI%6{YGQBSY##wwKh;zZg>Y! z*YHyQV!at_N&QE!CGBCMMcT)YkF;0w3zAJyBnR?Ix>S-}S4&ci@&nHrQvON#gX$Mj zo~MpCQrBhGlgA|}KhmT4j-o%sD2fRblPJbhjG*X6(SqVxigVO+E>dKto^M`FK3}7p zq*xxw#2At>Z%I}tk<2P1^PhN$$)f(;R+IXx^37Y68-z*y&a#mD=Kn?NKSncl z@%*8z$VJjtgfFs`Cv1_O`J_Dqzmfi{CPR|(GwG+srXB zej{mKP4d-ek^;pf87P1LlS|4O|Kra~^3wUG>jkObTT%C0y8qJimY%;E_9Rjds|A{)9x5)-yvI)R})G8DO_mJ&SM=^a-o;j4f$U41*MFugF{&Rb0 z)S}GQPx_aTFd0X!zJx9wzg|VQ%QtICx|x&hvSl5~uZ5)ln!X_!eVJs4G|6y&l4(pN z|C=D&+pAwBue6aoNwwGB8d9#SAo-<~WONZpMXKFeQ2#^6sP^kbwOrUq${y~NeVCS$_R`QN?RQ>|w5K9v-wY|z-ZmMe{f|bH{-8Ka`o|H6`qM4Szur;)7DM?TC*_a3DgUgZ{B@9u2l|vhU!nZ_F6Hm@6pvHw;66nQimxf! zP_+7=^Gp4}C)9C6iuWk$QM^g>xKs{BNL<&+h0^CUO26|$ zq@HP%{?AbM7^Un}LD?&gvfl`0&(D;711NhdQua@v{9&Pq^ba1&U(zW5VK_kglO*L| zTPc6*aUlKAb%gXsZpuGvUz7f7M)|Kf<4(l}@#{r&RmXrQ+dUs(miC*D|*MGcCk6k{p=rYQZ5j3aIoS(`{Xh~l;Hq^#Ub@=goM z?>|U}{v_GePBOWRq#QBt?f<@CL=KSLKSEM$isTQv^2PIot4ZGBC26&TA_6RR^%*=zeG`S4R_?}rVTByLXaj^?qXdYZIy?G)4}^H9@SfmSG+XR4=un)Y=6Z#AtamuaYaayY-?7Om)XGMj+1 zAK<1RXC`fN z9@Z8gl7wJBGAm&xrkLy6c!ze|sRN6akGtCLuLZSU+0o}>RS?YSAN=3W61eZ4R{tWA z$b)e11or@u-zWU5*Lhc#PhcUlg72kl47{66v?w6v$(=pn?jaOQjI&$LnF&X}fu`i= z0sc9}eDVAOO~d1lQ2(9Nn)(5ITzhiuE$K=Zd~^HjpSCi0>``)bFc|hkStkqnhtb~n zyNKi2@?E|-vaf2px|<*R5s5!>)E`|ko>xtV2jC4=WBHUTfjAVqF#U@&2%WzL-%zRz z!rjq_PM>!R#x18KEDjUnPaTc)eTRiZQFZNubBukV=oaj!qTm&V^vs8Z_eg%k4CODk z^#eZQeQS%=zkYngH69b|>$*PTrHu0@g;PGF!3S;|BaM&9V)HxJr8*3gSMu3MiiP3h zwX*6Aw?mQpwdr(}dkE^>7WXzx4MxkxzXlwYL3n)LTSK=w5dZ6paPsO3z?b`rOW)S} zW6a-$V{4jyF_vXb&AwJ|WInd}WM-=;Hm^wJh--1jpAXNqnl-xOv<%aQsY+*Ly%MB) zD#sB=D)08LBl63$e7ZNWj#$SvNX#%v@Oh2qU+y|5t+K=&-)*>AUO&g^Lfz{smmlD6 zk#&a{cbrGpJA?V7MOR^*c^9{}$YaP~{o~udCMytdXlQ?K;RqwQc8%2+`@o;{8k_9y zk6>PLcWvWb66CZEdH6lf1- zx3({X58JHNzY%$Bgw^CagKhm#Dav?rk8&Q~{r2# z^O0Ey($))amm%_6lz4;}vP%b`M5Ngn>x*6(CgP?c!sUP`E)93q&>ldBaP#%f4MA{} z?x+~Y?s#ac+wO5}BndvYd^L_YO@{2ycbCe>lOZ!b_O6z&YGo%sw0i!dKndr-^Zgs>riZ(Fz}-D(ijBrne!WP$#Ch^?V@grypLOpXCK5 z+&pxZ*BoHgqoH*dLTwWPXV{Aa5rawLj=bU=}$r+hn9!q00 z`VSKHPjlZi_XiFm0Bg_M_(J9y@OdDe-yq)z=Z}o_ zik)dC@>k6UeVgfkUU9a!%cpx_&z^0kU8a9Qj<(7UyUIcM^yz!G>*irFk2RiXvLAu9 zeJg|4hmAt49!rn&r7?J}WZTvsF$SsT1JB<0j6scpx$?>bV^Hz;vCEMUqmU)8bY_R| z2*|RS1Xw&Ch84~S89W4r;IxZeLh-3T;2PL3-nr!$yvVe@(OA|4MRT=W+?P7xxti$E zp^Tq!)rh+ys`wi)FONO-_iqif^Cb-RxfeszI?p(Uv?RD5Hh7!Y%@K`%9||ZK&cJ1> z{oD=OO7O?$7m4=K)o6QeCTaUXJsOR@y6h4C4cDwTTzO%z84IJ6JsLl?Vo>RZJgHZ0 zxTWlHNtkgvZhBlRc4c1&GVY><^gXlp%yuB3O|e$HR|npazY=hgu^q2?-uZDq zr5W#honVw&RfltXi+oQz6rjF`-hIaU06g1g>v{WDFwlCn!lE`Oz;*|oSF{~TaIK3q zrIaTLwD+!!ha<##>A~W{2Vx1pE~M(hz9JrS7&ycOzQsV2g?QA8(kS3DDY|!{H39^0 z^Rdlu2nTa9{sa6^!@&G~RL!rk5TF2GCQKOuLAPS-PsU7t z;Hr8b^~~8997;2}LT$alKYjN;nU9`O*UXhG-s%qSr^HL*cDupxVrM!|Ul-V^9-v-dyt!B{#8as9}+$bNz#vK180x8*7BO` zu+u(qjZ)Ah*fqH)Nb1KK;J1q3_sBy5EDR4H{(M>lbgudJ8cRx|DBm^bk{BhNf6mRT zy6Y@TNU?1?RCEz}YMN6g$N$6L-DN?3N3}59+^F2wTpOK!^<1W}xrN^^?q0^@bPFHp z6oz^yXkp`F7oY207f|YP4?UBoG?emL9ImOn3tx*AE=mhpLIS6zZFHw2bb6WRM40(Q zvXABF?-|5+p@)v4Kt&SNo_A>*YR>{zfq?Y>wS~ZCah>i{+-I=fZMpK@x*BL~vKO@9 z-T-@UbDd(cX@Y@!Gl$iAT7hlC(dnd4I|v%}9d9e?1n=V+Hg!Thup;2|D`x{Dzm>&3 z-Zd`2U_XO3&mGGFSOYiBUI`C^!m7LBzYPAuW79LU%pyZDxpE(lTpI$7?@#L_uMqjK zM4n%_-7*AuYt968-24mI)aaQmhz&y2!TRyXz5}pu_SBvDieIpA^|hPQ1$}TxF~%)PXi z@0iW>G9lHf8T)@&+)uw( zH2e*Zo*R~9vung$x{8+DZ4K!2?e&(?jym+~+7fuxp$7XU|0q9bsKoh{`qpQwOOapw zXzPZp`S{XZ<=W7pbbJ%Q^HsSg4&BbFIV$7@qQ;W2P%7Rz_ctFaLo%VYF5&LhP&sVpu@PwJo zVgI$t67vjGwhq$5y}-`38MgC5Tn~==Pv69j2(xr1hx?K63qh)72o*-|Mw4VjsgK-%r3sZ>s0`7@~671 z4+KEalS?9;R|4Th&(o5l)q&u%KC)&)H3;I1=aTnU1;L4t_LJQgi1p9CEW(Z5!EkWj z4uzQKA;7KB$y>RO$gkQND0L($6qqEs({0tmfI+Fw#b`VXMm}_#M!S9lGouO@ZwX?2 zz0IfdAhG^AZll@9n>!;wM)z%Fz#C%yb4YEaOKt=h&kOCiNR0POWk?7f&5i(3MN#gI zR}pZi_TW#I?Gdmp$4l&KN;vSGE_v2W+;3}{N{20R|98%1Iitg2z}n~hFj*-K^hG|K zRJ(-&m;EuD*sc(`A}sz=Rw@K$_YAD)w+M!W4T2&iH9;VGmP`81wjdZuR?+HvL9APq zp)-jl_V=;dcw~dh&H!K?7F#8^+8;(drUaAYe2I3Jd-#E@58Sbu+Fs(~1((#M4&(%R zKqIHV=F@#{FyhMNz`oKMCX9C+S#D_y^Rn3<&!S&~)6T{MPb+z_TstvX_0k$6+G{L- ziFo13tt&!9q(9;^D}F1E)rqKVEXeVwKOH~&wv@5D=3#bh$GP22#hB;ycfWO91@_5h z)LXo&!Hjj&cjEi%vEvi-vB;Oi|Kil{aa3u+f~nP2u&E7;-+GR}*xrHJir*i+RqDe1 zbqOxk20bVfr^IdQODvpusStO#@)z!F@h&?sK7i6!?uHt%{l#WA`mIVfL&)CRY4`8Z zFn*{G`)Q{Gz|!?gK|G=f4rW_3P)>8={CdFLQKyEQ2`N z9vo-&^f!t#`jq~|KIBXlI`OZ+8&BWU8rr_B6B~Y8JoQ-BhPPvVr#Fjz$9_)1x$BDc z7;{|L!hN&?_m1CO8x@*|t$fFSWVb|O6Z?(z9|Zlt>A>yBpPyuc^3O7Eet{B@U3HFY ztoIA7v0Bb^ouM9PzY6!P$ZUiyE7C{gd78nh#{BTz4Xq&RSu)uEUmIL3%e|Rm*$#}% zTF#MAJAn7GgN@UXPPo5|rNgzg6Y}_f-WER71t`JkJMPs56DDn|hqAkXC0?jbEUOD@ zPX1o?)4hwxXREpG*10aI7y7teaJUn?j`22IKI(+)5yjd=Egi5YVr_iUt`2Y(N}iK9 zXoqX!wHdFy+CV#OhV#$epP-+fD)rRn2iOE|i_)`chV_q9nCcChz}b0D^8hr$@gjw7 zU#99o)A`@S?E%DoXZ1yyPW!5%&p{|1UzS5LpN?INPcg{7Shp{ic5snx7w-9m=&J~Pfc_0n94>lk0fXjczvq{q;>&F9qo2oY zu%^;KWthPM-@F=GuqEc1x<{86H<-Ji%VnKCBn{n&gHLB{X~l zdfkvYAt>hPtQ)>8usOVUnLCPX);={jBnEagox#9a`tlLzE-SAfhqtX|GkNDrs z4^FzFZ`*IP5+-+)WMr;bzS$i=S9Yyp1$WdoUw&?#o;%)oV9P7#<&K*b_*aECxTD9R z+YUy&9{79Td9L4j9%#XRc|kwV0}XFoJ!vN7i3;*D(fv-IctPLOC6d((kDCfdGPrr6 z$zE9lO9^j$M`!ep@vApp6@TlL zp#Wt2AeTMg7Jz2TF4ZcYf!I9iu;8j2gmQ^F-$YIY;~E(!2Sv3I9Jy=Dv-Msm&h6WD zuQY<#$D~&}Xm~|9Znj<1*J~Yt{Vz__sf$Em-@bOUmA|7=@%rW~^Cf z8V)AVF3PP*$G~FtcUuS2(NkbU$UL#%lCiGMvMOR-8cWp`JN?E?tkB)`;T0_l%gTEb zbRDv=<)N~yMR^wTwd@s78_&X+41c#4u52uq@yOh`I~$D^YYmDOv#~MZ{+gYavvFZ; zcK(||Haabw)=RX?#?gbNX(~S1#C+8hWrmb&JnJN(DOsD1S|_wc*9~W5n?>KhFB~~| zc9Y%vC(=2%+vk?$A?(279-hPKY+H@vHyLu} zuM6{9@v-Ep$<+0a9$A|!YW#A?AJ+A`d;T~|?@-{&mCvEuR`O9GH>&?w_x>NkxoL&| zV*k7)ato@z<>~C$lbb8T@Yg?8CimX^%UMsaAI$A>pA&pzsQv$a$NKrJ!uvFlKak&G z>$j6-wHQ08xw>Dg4$oWNFsZmuje{0`S9k6xM>F9SbcZ#IaX0(lGtt}g(Ooh^$h|8I zotu~($_}L9g=tqaB{5=uA0tinLu*3u%@+QhX)B%a-Ak2U1&vj3X~uiq86wa1>ptb- zTIC3s&Pc6b(eZ>8`*j!o23dim278kew;HaF-KEs0WR1qVGFx}MJLACM9Iy5v9~`V* zKgM&NSTDROL@>Z33{}LuE#DVM;w!pGts0%NX#Rw=*uXFW3(iQguxlovfsy^G-@B4A z?odMSzq4sg7OptQr|nspg;P;StU8mkFjMnhi-vs`vX-j5z|3u_DrT38kwlDedlF1)=d23{W$UuX5a^7dG<38GO*wq zPhI=Q46K+gI6qpFj=U=aqRsB5BL*DP@+!ZVMsu01CO)+>=VfhkM+1$0QQ- zq%FPID``QtEHbKdggLKJq66{RgRj6lm*ifM^j;rQ-Q=+^U$ktj3iW*TiChi}Bz zX1r?sgi77wpLrXzP*-)^agP0kXxn`-6|Prcf{aIDZFwE)R@TpOLNog2$CTW^)`2BD zPpyV|`%se4_kvX1AYLlTO?}rlf_bf6neC3_xQByTnoV#TUEeL(%;?Qw@81XUiW`^Z z{$MLzwLyk1w`WqLL*gG@ZV{dK(BLupT%!qQw_tJl+`Hkc^7{Pfa!YRx1hfS&%T3LC z6v_T>9uL^hI%x~cVAPD@^P^G|X!lyJrMzbZr3JNRHhdevn%4?zFSd5#LDN3RmfK&k zp#BK!uzWgNeA(>;mbq$=dlt|pqmz?XVXFC6REo@DS0aj4K_>v)g$tr8=Uerw-gV>i_truXxsUs z%Egmc+2n1J<9c^S1HUF5vyXLL@y{7TE`=A?#RWq^ULx_==T-$%4QPgW~HPm|V! zw6jIXEn+)9WmAMsDLj+QnF?|4BbU?NA9*-tmCkbjagMOz+kH^d()cq^bN<-(Pm7mc>j;TpVa6CUd$K!Z(3sv zT|0K_2!9~*YUJe|wY<@fBI=rk=gQkrLZkZT{w)pIJ(B)#D+gj7eMo;rbUe%zzCUW% zTLH>+<=$t9zCzD@_LpYIA7IO-`<8QOJIDq9yHy+C0kh&eQ||oN3G7o|eLLklL82&O zb*okf=n8EO?~Q8%yOrM?qJ&%E)^6WeE8PYNRy^pmz)}e_O=;Cvh&=fZ`dyy5ihTl( z<3jFYtpV`m_O=ecel09!I>BWW?uomKGQ1q?hjIx6_P91yOxTrD-BR@WaAcQTrIIwS9EPaC0mQULEKf zkBr4XY!fcLl5rR+ClVu{5QkFcUv*!tkH>D+kJ{X7@mTTxt<7a(y~$_S*qF;{@yJ!S zYiM_0JZgwFew$sFfP7bVd)f{spxlOXxn(*D=n=Vdhv~Zntl~d&=4o^SeksW3UR#lX zl0yDZulFZl%Ql%GmW+wGUu9~Q-Nr<`q7lyJwl@)*YKo6Vt0ZEi;_-#ntHk-A5l>|e z5;4!B!nVLP5mk7T)nC3(MBhE%oDAI)aU=iL4NifHxS*ikSs#&z-cK9mMdK6E@uJ2? zqfd!Mo}%;Un39P1wbZ6h6I6deb9|VTh$(SDNoE2us`_Y z=jCDv_>kqpkYs*5&XtKs2_B0_cd@g|d5_}I0OzbE(m_AMjtL3|v4=V6EBd`6rZgWJE@(0cUns%Nd}Z#b z9F_Qn(K7WSPc3@f)9(qs{}nHc_858uHDk2I*1fk|+i+9rE2J1kCz=o5y|lU(^=oIqc!!I{9~Y26f_&SDS5_#>s2?2UA!k@zF1tyH-X2&@bos zg$FlA(MRBpqUi4-v=R@xDb_ZCEJa57-BZ1o=y=g!{6Hs$<%m4K^xy{`6;|D2x1s?J z)^X&$kNk`b+d(z7G#&e8nia!Gi1DlB6?<-dO9a-MCqMRd5&M2Q?0x^uqZGJxekCsR zsf4sIPt?s`RfF~Dn9DAc8km;bvZJ%K2G$+!^XPtF13#ar@=FI-gL{P5nt+fhIN=v8 z-6~lQTrQnAFvk_2?~8{B`Mnb>b+PdBg6?s#I~ z%t@{r7O|N2>glfQ2GLm0Q0J7mBNA7y!y3L6@G*Jq*GueYF?VB^^_W~5s9LVpJ-j>} zM3^M((!JBcn(wXsAR{r(oOXHas8R;xJB#!T6M6MrRakY_7O{7c-!2pOIkR z-VETjrg2_m$bkJL5_j{8)8S}Gr5MM{bXakecYWEObeM_x&$;tw8bn1HZX2*B_T6~6 zN{CT74UW9czxE4L;mty*&9ZZ;aBR%?+?m!Cm^>C*>3KZ`2F?9AE`3V|R{PaQf6FI> zh__DlCnDeK@OQ`JTSG};QPR|KbypJb<%e|c)=Y%f^gma%EEB+E)8!Mwe(|93SihF` zI1VKK(YM?yj{z1@?*3I*qT!ZfS)QPB6j=34j&{Tn`2uJ$?(ZWap@9FSd4ggTEZh4* zt=cXc{+`Wa@^OoW6w9PZtA$p{ECJa36>6BQR7Ku@ zC-xT!dE%Ay6U?2w;D^^jww`Pcd@bZXHqZJC%f#RYhCO=OcYKe4{$Vb9>SIX!Jp zLS1sw7Y5oa&Fh>t<8oT{h1$Y+?iI8J*J(~}Mkd6+ z(P|c2rUg!H_nm@m&sX+s%@_yqY;om#U8B(d<*j{~&@i~w%TM(({s9;6yJO|deV|q8 z`M6HA6PCS_uVdcU3JpwJWz%IX<~^KWpL9XVkqTyes!^DX|}0Xk6%A5uQ4+^83)C0(`zB zeA4DfE?PfmY>YUTfl&o7tdraWO=hq^2Eb8NGD;c6B6+l zcr2`gBQsne+%)9%jrRdW9KNhSWfclTP8=8HpGUy**xQcgKGAU5*JS?{o;YZH$E#fM zHXhz9|Gl&O0tB%IAA+f$OkWc`Q8&0mvYHZa;R=r6Hu%-*wu*&15Cfh(V40u->{a@Y`KjBBVEk4QMdbhWgo5+XhMN8y4 z#-9wTWjB=?iT9KHRXm&2ZV~yG+a!apk0rtKOY;XZ%b61mZaV~JiQlOeH)CCF-L=B|H5P5(srUDlVTp)g0^PWMc z3n&zbwy$Avg{$jRfAQ^cg=D9lM^mr6!giLb!2m~MyvEjV`*e{j6o^h*^AO`V=I+Yn z`Ul+LL5gP1tru=^!0<}Wn>-@#kV=Wg8CG|&<1JQ@({zU$uk6<6r@F&;wT5n0ZV%Y@ zNMrJdnFn~hA93I$#&@zx)W$*E6BJKG>nL@3!tCFvbI*x9M>Cv?RRarN@c70yj(l%# zU`uHhsXpohYkpmR5jIf>hs>^xMXx9)-Y-{e^dj~#-N1SzVawY(Sf2b#`zZ0= 1.0 always. +print(f'minimum value of fact = {np.min(solution[4*num_nodes: 5*num_nodes]):.4e}') + +# %% +# sphinx_gallery_thumbnail_number = 2 + From 3d631eda17bd9f43b93a6a6e0c94e6a12bb414b7 Mon Sep 17 00:00:00 2001 From: Peter Stahlecker <83544146+Peter230655@users.noreply.github.com> Date: Sun, 24 Nov 2024 09:49:19 +0100 Subject: [PATCH 12/12] Delete examples-gallery/plot_betts_10_144_145.py --- examples-gallery/plot_betts_10_144_145.py | 220 ---------------------- 1 file changed, 220 deletions(-) delete mode 100644 examples-gallery/plot_betts_10_144_145.py diff --git a/examples-gallery/plot_betts_10_144_145.py b/examples-gallery/plot_betts_10_144_145.py deleted file mode 100644 index 6c626b63..00000000 --- a/examples-gallery/plot_betts_10_144_145.py +++ /dev/null @@ -1,220 +0,0 @@ -""" -Van der Pol Oscillator -====================== - -https://en.wikipedia.org/wiki/Van_der_Pol_oscillator - -These are example 10.144 / 145 from Betts' book "Practical Methods for Optimal Control -Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems. -It is described in more detail in section 4.14. example 4.11 of the book. - -""" -import numpy as np -import sympy as sm -import sympy.physics.mechanics as me -from opty.direct_collocation import Problem -from opty.utils import create_objective_function -import matplotlib.pyplot as plt - -# %% -# First version of the problem (10.144). -# -------------------------------------- -# **States** -# :math:`y_1, y_2` : state variables -# -# **Controls** -# :math:`u` : control variables -# -# Equations of Motion. -# -------------------- -t = me.dynamicsymbols._t - -y1, y2 = me.dynamicsymbols('y1, y2') -u = me.dynamicsymbols('u') - -eom = sm.Matrix([ - -y1.diff() + y2, - -y2.diff(t) + (1 - y1**2)*y2 - y1 + u, -]) -sm.pprint(eom) - -# %% -# Define and Solve the Optimization Problem. -# ------------------------------------------ -num_nodes = 2001 -iterations = 1000 - -t0, tf = 0.0, 5.0 -interval_value = (tf - t0)/(num_nodes - 1) - -state_symbols = (y1, y2) -unkonwn_input_trajectories = (u, ) - -objective = sm.integrate(u**2 + y1**2 + y2**2, t) -obj, obj_grad = create_objective_function(objective, - state_symbols, - unkonwn_input_trajectories, - tuple(), - num_nodes, - interval_value -) - -instance_constraints = ( - y1.func(t0) - 1.0, - y2.func(t0), -) - -bounds = { - y2: (-0.4, np.inf), -} - -prob = Problem( - obj, - obj_grad, - eom, - state_symbols, - num_nodes, - interval_value, - instance_constraints= instance_constraints, - bounds=bounds, -) - -# %% -# Solve the optimization problem. -# Give some rough estimates for the trajectories. - -initial_guess = np.ones(prob.num_free) - -# Find the optimal solution. -for i in range(1): - solution, info = prob.solve(initial_guess) - initial_guess = solution - print(info['status_msg']) - print(f'Objectve is: {info['obj_val']:.8f}, ' + - f'as per the book it is {2.95369916}, so the deviation is: ' + - f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %') -solution1 = solution -# %% -# Plot the optimal state and input trajectories. -prob.plot_trajectories(solution) - -# %% -# Plot the constraint violations. -prob.plot_constraint_violations(solution) - -# %% -# Plot the objective function as a function of optimizer iteration. -prob.plot_objective_value() - -# %% -# Second version of the problem (10.145). -# --------------------------------------- -# This is example 10.145 from Betts' book "Practical Methods for Optimal Control -# Using NonlinearProgramming", 3rd edition, Chapter 10: Test Problems. -# It is same as problem 10.144 but a bit reformulated. -# It has two control variables and one additional algebraic equation of motion. -# As opty needs as many state variables as equations of motion, I simply call -# the additional control variable v a state variable. -# It is described in more detail in section 4.14, example 4.11 of the book. -# -# This formulation seems to be more accurate compared to the above, when -# considering the violations of the constraints. -# -# **States** -# :math:`y_1, y_2, v` : state variables -# -# **Controls** -# :math:`u` : control variables -# -# Equations of Motion. -# -------------------- -y1, y2, v = me.dynamicsymbols('y1, y2, v') -u = me.dynamicsymbols('u') - -eom = sm.Matrix([ - -y1.diff() + y2, - -y2.diff(t) + v - y1 + u, - v - (1-y1**2)*y2, -]) -sm.pprint(eom) - -# %% -# Define and Solve the Optimization Problem. -# ------------------------------------------ -state_symbols = (y1, y2, v) -unkonwn_input_trajectories = (u,) - -objective = sm.integrate(u**2 + y1**2 + y2**2, t) -obj, obj_grad = create_objective_function(objective, - state_symbols, - unkonwn_input_trajectories, - tuple(), - num_nodes, - interval_value -) - -instance_constraints = ( - y1.func(t0) - 1.0, - y2.func(t0), -) - -bounds = { - y2: (-0.4, np.inf), -} - -prob = Problem( - obj, - obj_grad, - eom, - state_symbols, - num_nodes, - interval_value, - instance_constraints= instance_constraints, - bounds=bounds, -) - -# %% -# Solve the optimization problem. -# Give some rough estimates for the trajectories. - -initial_guess = np.ones(prob.num_free) - -# Find the optimal solution. -for i in range(1): - solution, info = prob.solve(initial_guess) - initial_guess = solution - print(info['status_msg']) - print(f'Objectve is: {info['obj_val']:.8f}, ' + - f'as per the book it is {2.95369916}, so the deviation is: ' + - f'{(info['obj_val'] -2.95369916) /2.95369916*100 :.5e} %') -# %% -# Plot the optimal state and input trajectories. -prob.plot_trajectories(solution) - -# %% -# Plot the constraint violations. -prob.plot_constraint_violations(solution) - -# %% -# Plot the objective function as a function of optimizer iteration. -prob.plot_objective_value() - -# %% -# Plot the Difference between the two Solutions. -# ---------------------------------------------- -diffy1 = solution1[: num_nodes] - solution[: num_nodes] -diffy2 = solution1[num_nodes: 2*num_nodes] - solution[num_nodes: 2*num_nodes] -diffu = solution1[2*num_nodes:] - solution[3*num_nodes:] -times = np.linspace(t0, tf, num_nodes) - -fig, ax = plt.subplots(3, 1, figsize=(7, 4), sharex=True, constrained_layout=True) -ax[0].plot(times, diffy1, label='Delta y1') -ax[0].legend() -ax[1].plot(times, diffy2, label='Delta y2') -ax[1].legend() -ax[2].plot(times, diffu, label='Delta u') -ax[2].legend() -ax[2].set_xlabel('Time') -ax[0].set_title('Differences between the two solutions') -avoid_printing = True -# sphinx_gallery_thumbnail_number = 2