.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/intermediate/plot_non_contiguous_parameter_identification.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_intermediate_plot_non_contiguous_parameter_identification.py: Parameter Identification from Noncontiguous Measurements ======================================================== 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 Example ========================== The position of a simple system consisting of a mass connected to a fixed point by a spring and a damper is simulated and recorded as noisy measurements. The spring constant and the damping coefficient will be identified. **State Variables** - :math:`x_1`: position of the mass of the first measurement trial [m] - :math:`x_2`: position of the mass of the second measurement trial [m] - :math:`x_3`: position of the mass of the third measurement trial [m] - :math:`x_4`: position of the mass of the fourth measurement trial [m] - :math:`u_1`: speed of the mass of the first measurement trial [m/s] - :math:`u_2`: speed of the mass of the second measurement trial [m/s] - :math:`u_3`: speed of the mass of the third measurement trial [m/s] - :math:`u_4`: speed of the mass of the fourth 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] .. GENERATED FROM PYTHON SOURCE LINES 44-51 .. code-block:: Python import sympy as sm import sympy.physics.mechanics as me import numpy as np from scipy.integrate import solve_ivp from opty import Problem import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 52-58 Equations of Motion ------------------- Set up the four sets of equations of motion, one for each of the four measurements. .. GENERATED FROM PYTHON SOURCE LINES 58-76 .. code-block:: Python x1, x2, x3, x4 = me.dynamicsymbols('x1, x2, x3, x4') u1, u2, u3, u4 = me.dynamicsymbols('u1, u2, u3, u4') m, c, k, l0 = sm.symbols('m, c, k, l0') t = me.dynamicsymbols._t eom = sm.Matrix([ x1.diff(t) - u1, x2.diff(t) - u2, x3.diff(t) - u3, x4.diff(t) - u4, m*u1.diff(t) + c*u1 + k*(x1 - l0), m*u2.diff(t) + c*u2 + k*(x2 - l0), m*u3.diff(t) + c*u3 + k*(x3 - l0), m*u4.diff(t) + c*u4 + k*(x4 - l0), ]) 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 ⎥ ⎢c⋅u₁(t) + k⋅(-l₀ + x₁(t)) + m⋅──(u₁(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢c⋅u₂(t) + k⋅(-l₀ + x₂(t)) + m⋅──(u₂(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢c⋅u₃(t) + k⋅(-l₀ + x₃(t)) + m⋅──(u₃(t))⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢c⋅u₄(t) + k⋅(-l₀ + x₄(t)) + m⋅──(u₄(t))⎥ ⎣ dt ⎦ .. GENERATED FROM PYTHON SOURCE LINES 77-84 Generate Noisy Measurement Data ------------------------------- Create four 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. .. GENERATED FROM PYTHON SOURCE LINES 84-116 .. code-block:: Python rhs = sm.Matrix([ u1, u2, u3, u4, 1/m*(-c*u1 - k*(x1 - l0)), 1/m*(-c*u2 - k*(x2 - l0)), 1/m*(-c*u3 - k*(x3 - l0)), 1/m*(-c*u4 - k*(x4 - l0)), ]) states = [x1, x2, x3, x4, u1, u2, u3, u4] parameters = [m, c, k, l0] par_vals = [1.0, 0.25, 1.0, 1.0] eval_rhs = sm.lambdify(states + parameters, rhs) t0, tf = 0.0, 20.0 num_nodes = 500 times = np.linspace(t0, tf, num=num_nodes) measurements = [] np.random.seed(123) for i in range(4): x0 = 4.0*np.random.randn(8) sol = solve_ivp(lambda t, x, p: eval_rhs(*x, *p).squeeze(), (t0, tf), x0, t_eval=times, args=(par_vals,)) measurements.append(sol.y[0, :] + 2.0*np.random.randn(len(sol.t))) measurements = np.array(measurements) print(measurements.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (4, 500) .. GENERATED FROM PYTHON SOURCE LINES 117-135 Setup the Identification Problem -------------------------------- The goal is to identify the damping coefficient :math:`c` and the spring constant :math:`k`. 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. .. math:: J(x_1, x_2, x_3, x_4) = \\ \int_{t_0}^{t_f} \left( w_1 (x_1 - x_1^m)^2 + w_2 (x_2 - x_2^m)^2 + w_3 (x_3 - x_3^m)^2 + w_4 (x_4 - x_4^m)^2 \right) dt .. GENERATED FROM PYTHON SOURCE LINES 135-161 .. code-block:: Python interval_value = (tf - t0) / (num_nodes - 1) w = [1.0, 1.0, 1.0, 1.0] def obj(free): return interval_value*np.sum( w[0]*(free[0*num_nodes:1*num_nodes] - measurements[0])**2 + w[1]*(free[1*num_nodes:2*num_nodes] - measurements[1])**2 + w[2]*(free[2*num_nodes:3*num_nodes] - measurements[2])**2 + w[3]*(free[3*num_nodes:4*num_nodes] - measurements[3])**2) def obj_grad(free): grad = np.zeros_like(free) grad[:num_nodes] = 2*w[0]*interval_value*( free[0*num_nodes:1*num_nodes] - measurements[0]) grad[num_nodes:2*num_nodes] = 2*w[1]*interval_value*( free[1*num_nodes:2*num_nodes] - measurements[1]) grad[2*num_nodes:3*num_nodes] = 2*w[2]*interval_value*( free[2*num_nodes:3*num_nodes] - measurements[2]) grad[3*num_nodes:4*num_nodes] = 2*w[3]*interval_value*( free[3*num_nodes:4*num_nodes] - measurements[3]) return grad .. GENERATED FROM PYTHON SOURCE LINES 162-164 By not including :math:`c` and :math:`k` in the parameter map, they will be treated as unknown parameters. .. GENERATED FROM PYTHON SOURCE LINES 164-167 .. code-block:: Python par_map = {m: par_vals[0], l0: par_vals[3]} print(par_map) .. rst-class:: sphx-glr-script-out .. code-block:: none {m: 1.0, l0: 1.0} .. GENERATED FROM PYTHON SOURCE LINES 168-185 .. code-block:: Python bounds = { c: (0.01, 2.0), k: (0.1, 10.0), } 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, ) .. GENERATED FROM PYTHON SOURCE LINES 186-194 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 two values are the guesses for :math:`c` and :math:`k`, respectively. .. GENERATED FROM PYTHON SOURCE LINES 194-198 .. code-block:: Python initial_guess = np.hstack((np.array(measurements).flatten(), # x np.zeros(4*num_nodes), # u [0.1, 3.0])) # c, k .. GENERATED FROM PYTHON SOURCE LINES 199-203 Solve the Optimization Problem ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 203-206 .. code-block:: Python solution, info = problem.solve(initial_guess) print(info['status_msg']) .. 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 207-209 .. code-block:: Python _ = problem.plot_objective_value() .. image-sg:: /examples/intermediate/images/sphx_glr_plot_non_contiguous_parameter_identification_001.png :alt: Objective Value :srcset: /examples/intermediate/images/sphx_glr_plot_non_contiguous_parameter_identification_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 210-212 .. code-block:: Python _ = problem.plot_constraint_violations(solution) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_non_contiguous_parameter_identification_002.png :alt: Constraint violations :srcset: /examples/intermediate/images/sphx_glr_plot_non_contiguous_parameter_identification_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 213-215 The identified parameters are: .. GENERATED FROM PYTHON SOURCE LINES 215-218 .. code-block:: Python print(f'Estimate of damping coefficient is {solution[-2]: 1.2f}') print(f'Estimate of the spring constant is {solution[-1]: 1.2f}') .. rst-class:: sphx-glr-script-out .. code-block:: none Estimate of damping coefficient is 0.25 Estimate of the spring constant is 1.01 .. GENERATED FROM PYTHON SOURCE LINES 219-222 Plot the Measurements and the Estimated Trajectories ---------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 222-228 .. code-block:: Python fig, ax = plt.subplots(8, 1, figsize=(6, 8), sharex=True) for i in range(4): ax[i].plot(times, measurements[i]) problem.plot_trajectories(solution, axes=ax) plt.show() .. image-sg:: /examples/intermediate/images/sphx_glr_plot_non_contiguous_parameter_identification_003.png :alt: State Trajectories :srcset: /examples/intermediate/images/sphx_glr_plot_non_contiguous_parameter_identification_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 33.160 seconds) .. _sphx_glr_download_examples_intermediate_plot_non_contiguous_parameter_identification.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_non_contiguous_parameter_identification.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_non_contiguous_parameter_identification.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_non_contiguous_parameter_identification.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_