.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\plot_identification_of_four_parameters_opty.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_plot_identification_of_four_parameters_opty.py: Parameter Identification with Non-Contiguous Measurements. ========================================================== Objective --------- - Show how to use opty to estimate parameters of a mechanical system from noisy measurements. - Show a way to handle explicit time in the equations of motion (very small point) Description ----------- For parameter estimation it is common to collect measurements of a system's trajectories for distinct experiments. For example, if you are identifying the parameters of a mass-spring-damper system, you will excite the system with different initial conditions multiple times. The data cannot simply be stacked and the identification run because the measurement data would be discontinuous between trials. A work around 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. For example: *number_of_measurements* measurements of the location of a simple system consisting of a mass connected to a fixed point by a spring and a damper, with a force = :math:`F_1 \sin(\omega \cdot t)` acting on the mass. The movement is in a horizontal direction. :math:`c, F_1, k, \omega` are to be identified. Notes ----- This is just a slightly more elaborate version of an example from opty: https://opty.readthedocs.io/stable/examples/intermediate/plot_non_contiguous_parameter_identification.html#sphx-glr-examples-intermediate-plot-non-contiguous-parameter-identification-py **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`: damping coefficient [Ns/m] - :math:`k`: spring constant [N/m] - :math:`l_0`: natural length of the spring [m] - :math:`F_1`: amplitude of the force [N] - :math:`\omega`: frequency of the force [rad/s] .. GENERATED FROM PYTHON SOURCE LINES 61-63 Set up the equations of motion and integrate them to get the measurements. .. GENERATED FROM PYTHON SOURCE LINES 63-90 .. code-block:: Python import sympy as sm import numpy as np import sympy.physics.mechanics as me import matplotlib.pyplot as plt from scipy.integrate import solve_ivp from opty import Problem number_of_measurements = 15 t0, tf = 0, 10 num_nodes = 400 times = np.linspace(t0, tf, num_nodes) t_span = (t0, tf) np.random.seed(1234) x = me.dynamicsymbols(f'x:{number_of_measurements}') u = me.dynamicsymbols(f'u:{number_of_measurements}') fx = me.dynamicsymbols('f_x') m, c, k, l0 = sm.symbols('m, c, k, l0') F1, omega = sm.symbols('F1, omega') t = me.dynamicsymbols._t T = me.dynamicsymbols('T') eom1 = sm.Matrix([x[i].diff(t) - u[i] for i in range(number_of_measurements)]) eom2 = sm.Matrix([m*u[i].diff(t) + c*u[i] + k*(x[i] - l0) - F1*sm.sin(omega*T) for i in range(number_of_measurements)]) eom = eom1.col_join(eom2) .. GENERATED FROM PYTHON SOURCE LINES 91-93 Print the equations of motion. .. GENERATED FROM PYTHON SOURCE LINES 93-94 .. code-block:: Python sm.pprint(eom) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ⎤ ⎢ -u₀(t) + ──(x₀(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₁(t) + ──(x₁(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₂(t) + ──(x₂(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₃(t) + ──(x₃(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₄(t) + ──(x₄(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₅(t) + ──(x₅(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₆(t) + ──(x₆(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₇(t) + ──(x₇(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₈(t) + ──(x₈(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₉(t) + ──(x₉(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₁₀(t) + ──(x₁₀(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₁₁(t) + ──(x₁₁(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₁₂(t) + ──(x₁₂(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₁₃(t) + ──(x₁₃(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u₁₄(t) + ──(x₁₄(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₀(t) + k⋅(-l₀ + x₀(t)) + m⋅──(u₀(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₁(t) + k⋅(-l₀ + x₁(t)) + m⋅──(u₁(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₂(t) + k⋅(-l₀ + x₂(t)) + m⋅──(u₂(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₃(t) + k⋅(-l₀ + x₃(t)) + m⋅──(u₃(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₄(t) + k⋅(-l₀ + x₄(t)) + m⋅──(u₄(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₅(t) + k⋅(-l₀ + x₅(t)) + m⋅──(u₅(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₆(t) + k⋅(-l₀ + x₆(t)) + m⋅──(u₆(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₇(t) + k⋅(-l₀ + x₇(t)) + m⋅──(u₇(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₈(t) + k⋅(-l₀ + x₈(t)) + m⋅──(u₈(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -F₁⋅sin(ω⋅T(t)) + c⋅u₉(t) + k⋅(-l₀ + x₉(t)) + m⋅──(u₉(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢-F₁⋅sin(ω⋅T(t)) + c⋅u₁₀(t) + k⋅(-l₀ + x₁₀(t)) + m⋅──(u₁₀(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢-F₁⋅sin(ω⋅T(t)) + c⋅u₁₁(t) + k⋅(-l₀ + x₁₁(t)) + m⋅──(u₁₁(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢-F₁⋅sin(ω⋅T(t)) + c⋅u₁₂(t) + k⋅(-l₀ + x₁₂(t)) + m⋅──(u₁₂(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢-F₁⋅sin(ω⋅T(t)) + c⋅u₁₃(t) + k⋅(-l₀ + x₁₃(t)) + m⋅──(u₁₃(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢-F₁⋅sin(ω⋅T(t)) + c⋅u₁₄(t) + k⋅(-l₀ + x₁₄(t)) + m⋅──(u₁₄(t))⎥ ⎣ dt ⎦ .. GENERATED FROM PYTHON SOURCE LINES 95-97 Equations of motion for the solve_ivp integration. .. GENERATED FROM PYTHON SOURCE LINES 97-121 .. code-block:: Python rhs1 = np.array([u[i] for i in range(number_of_measurements)]) rhs2 = np.array([1/m * (-c*u[i] - k*(x[i] - l0) + F1*sm.sin(omega*t)) for i in range(number_of_measurements)]) rhs = np.concatenate((rhs1, rhs2)) qL = x + u pL = [m, c, k, l0, F1, omega, t] rhs_lam = sm.lambdify(qL + pL, rhs) def gradient(t, x, args): args[-1] = t return rhs_lam(*x, *args).reshape(2 * number_of_measurements) # Different initial conditions for the different measurements. x0 = np.array([2 + (i / 4) * (-1)**i for i in range(number_of_measurements)] + [0 for _ in range(number_of_measurements)]) pL_vals = [1.0, 0.25, 2.0, 1.0, 6.0, 3.0, t0] resultat1 = solve_ivp(gradient, t_span, x0, t_eval=times, args=(pL_vals,)) resultat = resultat1.y.T .. GENERATED FROM PYTHON SOURCE LINES 122-125 Create the noisy measurements: simply add noise to the results of the integration.It is assumed, that only the locations are measured, not the speeds. .. GENERATED FROM PYTHON SOURCE LINES 125-132 .. code-block:: Python measurements = [] for i in range(number_of_measurements): measurements.append(resultat[:, i] + np.random.randn(resultat.shape[0]) * 1.0 / (i / 20 + 1) + np.random.randn(1) * 1.0 / (i / 20 + 1)) .. GENERATED FROM PYTHON SOURCE LINES 133-144 Set up the Estimation Problem. -------------------------------- The idea is Gauss' method of least squares: https://en.wikipedia.org/wiki/Least_squares If some measurement is considered more reliable, its weight w may be increased. objective = :math:`\int_{t_0}^{t_f} \left[ \sum_{i=1}^{ \textrm{number_of_measurements}} (w_i \cdot (x_i - x_i^m)^2 \right] dt` .. GENERATED FROM PYTHON SOURCE LINES 144-198 .. code-block:: Python # This is added to get the explicit time for opty. eom = eom.col_join(sm.Matrix([T.diff(t) - 1])) state_symbols = x + u + [T] interval_value = (tf - t0) / (num_nodes - 1) par_map = {m: pL_vals[0], l0: pL_vals[3], } # Weight vector. Here measurements with higher order number are considered # more reliable, so their weight is set to be larger. w = [1 + i for i in range(number_of_measurements)] def obj(free): return interval_value * np.sum([w[i] * np.sum( (free[i*num_nodes:(i+1)*num_nodes] - measurements[i])**2) for i in range(number_of_measurements)]) def obj_grad(free): grad = np.zeros_like(free) for i in range(number_of_measurements): grad[i*num_nodes: (i+1)*num_nodes] = 2*w[i]*interval_value*( free[i*num_nodes:(i+1)*num_nodes] - measurements[i] ) return grad bounds = { c: (0, 1), k: (1, 3), F1: (5, 10), omega: (2, 7), } instance_constraints = ( T.func(t0) - 0.0, ) problem = Problem( obj, obj_grad, eom, state_symbols, num_nodes, interval_value, known_parameter_map=par_map, instance_constraints=instance_constraints, bounds=bounds, time_symbol=me.dynamicsymbols._t, ) .. GENERATED FROM PYTHON SOURCE LINES 199-200 This gives the unknown parameters, and their sequence in the solution vector. .. GENERATED FROM PYTHON SOURCE LINES 200-201 .. code-block:: Python print('unknown parameters are:', problem.collocator.unknown_parameters) .. rst-class:: sphx-glr-script-out .. code-block:: none unknown parameters are: (F1, c, k, omega) .. GENERATED FROM PYTHON SOURCE LINES 202-203 Initial guess. .. GENERATED FROM PYTHON SOURCE LINES 203-213 .. code-block:: Python list1 = [list(measurements[i]) for i in range(number_of_measurements)] list1 = list(np.array(list1).flat) initial_guess = np.array( list1 + list(np.zeros(number_of_measurements*num_nodes)) + list(np.linspace(t0, tf, num_nodes)) + [0, 0, 0, 0] ) .. GENERATED FROM PYTHON SOURCE LINES 214-216 Solve the Optimization Problem. .. GENERATED FROM PYTHON SOURCE LINES 216-219 .. code-block:: Python solution, info = problem.solve(initial_guess) print(info['status_msg']) _ = problem.plot_objective_value() .. image-sg:: /examples/images/sphx_glr_plot_identification_of_four_parameters_opty_001.png :alt: Objective Value :srcset: /examples/images/sphx_glr_plot_identification_of_four_parameters_opty_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).' .. GENERATED FROM PYTHON SOURCE LINES 220-221 .. code-block:: Python _ = problem.plot_constraint_violations(solution) .. image-sg:: /examples/images/sphx_glr_plot_identification_of_four_parameters_opty_002.png :alt: Constraint violations :srcset: /examples/images/sphx_glr_plot_identification_of_four_parameters_opty_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 222-225 Results obtained ---------------- .. GENERATED FROM PYTHON SOURCE LINES 225-235 .. code-block:: Python print('As the true values are known in this example, an error can be given') print(f'Estimate of dampening constant is {solution[-3]:.2f} error is ' f'{(solution[-3] - pL_vals[1])/pL_vals[1] * 100:.2f} %') print(f'Estimate of spring constant is {solution[-2]:.2f} error is ' f'{(solution[-2] - pL_vals[2])/pL_vals[2] * 100:.2f} %') print(f'Estimate of force is {solution[-4]:.2f} error is ' f'{(solution[-4] - pL_vals[4])/pL_vals[4] * 100:.2f} %') print(f'Estimate of the frequency is {solution[-1]:.2f} error is ' f'{(solution[-1] - pL_vals[5])/pL_vals[5] * 100:.2f} %') .. rst-class:: sphx-glr-script-out .. code-block:: none As the true values are known in this example, an error can be given Estimate of dampening constant is 0.21 error is -14.30 % Estimate of spring constant is 1.98 error is -1.01 % Estimate of force is 5.93 error is -1.25 % Estimate of the frequency is 2.99 error is -0.48 % .. GENERATED FROM PYTHON SOURCE LINES 236-237 Plot the measurements. .. GENERATED FROM PYTHON SOURCE LINES 237-247 .. code-block:: Python fig, ax = plt.subplots(number_of_measurements, 1, figsize=(8, 1.25 * number_of_measurements), sharex=True, layout='constrained') for i in range(number_of_measurements): ax[i].plot(times, measurements[i]) ax[i].set_ylabel(f'Set - {i+1} [m]') ax[0].set_title('Measurements') _ = ax[-1].set_xlabel('Time [sec]') plt.show() .. image-sg:: /examples/images/sphx_glr_plot_identification_of_four_parameters_opty_003.png :alt: Measurements :srcset: /examples/images/sphx_glr_plot_identification_of_four_parameters_opty_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 22.450 seconds) .. _sphx_glr_download_examples_plot_identification_of_four_parameters_opty.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_identification_of_four_parameters_opty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_identification_of_four_parameters_opty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_identification_of_four_parameters_opty.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_