.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/beginner/plot_vyasarayani.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_beginner_plot_vyasarayani.py: Multi-minimum Parameter Identification ====================================== This example is taken from [Vyasarayani2011]_. In Section 3.1 there is a simple example of a single pendulum parameter identification that has many local minima. For the following differential equations that describe a single pendulum acting under the influence of gravity, the goals is to identify the parameter p given noisy measurements of the angle, y1. :: -- -- -- -- | y1' | | y2 | y' = f(y, t) = | | = | | | y2' | | -p*sin(y1) | -- -- -- -- .. GENERATED FROM PYTHON SOURCE LINES 22-29 .. code-block:: Python import numpy as np import sympy as sm from scipy.integrate import odeint import matplotlib.pyplot as plt from opty import Problem .. GENERATED FROM PYTHON SOURCE LINES 30-31 Specify the symbolic equations of motion. .. GENERATED FROM PYTHON SOURCE LINES 31-38 .. code-block:: Python p, t = sm.symbols('p, t') y1, y2 = [f(t) for f in sm.symbols('y1, y2', cls=sm.Function)] y = sm.Matrix([y1, y2]) f = sm.Matrix([y2, -p*sm.sin(y1)]) eom = y.diff(t) - f sm.pprint(eom) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ⎤ ⎢ -y₂(t) + ──(y₁(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢p⋅sin(y₁(t)) + ──(y₂(t))⎥ ⎣ dt ⎦ .. GENERATED FROM PYTHON SOURCE LINES 39-40 Generate some data by integrating the equations of motion. .. GENERATED FROM PYTHON SOURCE LINES 40-57 .. code-block:: Python duration = 50.0 num_nodes = 5000 interval = duration/(num_nodes - 1) time = np.linspace(0.0, duration, num=num_nodes) p_val = 10.0 y0 = [np.pi/6.0, 0.0] def eval_f(y, t, p): return np.array([y[1], -p*np.sin(y[0])]) y_meas = odeint(eval_f, y0, time, args=(p_val,)) y1_meas = y_meas[:, 0] y2_meas = y_meas[:, 1] .. GENERATED FROM PYTHON SOURCE LINES 58-59 Add measurement noise. .. GENERATED FROM PYTHON SOURCE LINES 59-63 .. code-block:: Python y1_meas += np.random.normal(scale=0.05, size=y1_meas.shape) y2_meas += np.random.normal(scale=0.1, size=y2_meas.shape) .. GENERATED FROM PYTHON SOURCE LINES 64-67 Setup the optimization problem to minimize the error in the simulated angle and the measured angle. The midpoint integration method is preferable to the backward Euler method because no artificial damping is introduced. .. GENERATED FROM PYTHON SOURCE LINES 67-83 .. code-block:: Python def obj(free): """Minimize the error in the angle, y1.""" return interval*np.sum((y1_meas - free[:num_nodes])**2) def obj_grad(free): grad = np.zeros_like(free) grad[:num_nodes] = 2.0*interval*(free[:num_nodes] - y1_meas) return grad prob = Problem(obj, obj_grad, eom, (y1, y2), num_nodes, interval, time_symbol=t, integration_method='midpoint') num_states = len(y) .. GENERATED FROM PYTHON SOURCE LINES 84-86 Give noisy measurements as the initial state guess and a random positive values as the parameter guess. .. GENERATED FROM PYTHON SOURCE LINES 86-88 .. code-block:: Python initial_guess = np.hstack((y1_meas, y2_meas, 100.0*np.random.random(1))) .. GENERATED FROM PYTHON SOURCE LINES 89-90 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 90-93 .. code-block:: Python solution, info = prob.solve(initial_guess) p_sol = solution[-1] .. GENERATED FROM PYTHON SOURCE LINES 94-95 Print the result. .. GENERATED FROM PYTHON SOURCE LINES 95-106 .. code-block:: Python known_msg = "Known value of p = {}".format(p_val) guess_msg = "Initial guess for p = {}".format(initial_guess[-1]) identified_msg = "Identified value of p = {}".format(p_sol) divider = '='*max(len(known_msg), len(identified_msg)) print(divider) print(known_msg) print(guess_msg) print(identified_msg) print(divider) .. rst-class:: sphx-glr-script-out .. code-block:: none ========================================== Known value of p = 10.0 Initial guess for p = 45.42404253743789 Identified value of p = 10.001366458774221 ========================================== .. GENERATED FROM PYTHON SOURCE LINES 107-108 Plot constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 108-110 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/beginner/images/sphx_glr_plot_vyasarayani_001.png :alt: Constraint violations :srcset: /examples/beginner/images/sphx_glr_plot_vyasarayani_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 111-112 Simulate with the identified parameter. .. GENERATED FROM PYTHON SOURCE LINES 112-116 .. code-block:: Python y_sim = odeint(eval_f, y0, time, args=(p_sol,)) y1_sim = y_sim[:, 0] y2_sim = y_sim[:, 1] .. GENERATED FROM PYTHON SOURCE LINES 117-118 Plot results .. GENERATED FROM PYTHON SOURCE LINES 118-137 .. code-block:: Python fig_y1, axes_y1 = plt.subplots(3, 1, layout='constrained') legend = ['measured', 'initial guess', 'direct collocation solution', 'identified simulated'] axes_y1[0].plot(time, y1_meas, '.k', time, initial_guess[:num_nodes], '.b', time, solution[:num_nodes], '.r', time, y1_sim, 'g') axes_y1[0].set_xlabel('Time [s]') axes_y1[0].set_ylabel('y1 [rad]') axes_y1[0].legend(legend) axes_y1[1].set_title('Initial Guess Constraint Violations') axes_y1[1].plot(prob.con(initial_guess)[:num_nodes - 1]) axes_y1[2].set_title('Solution Constraint Violations') axes_y1[2].plot(prob.con(solution)[:num_nodes - 1]) .. image-sg:: /examples/beginner/images/sphx_glr_plot_vyasarayani_002.png :alt: Initial Guess Constraint Violations, Solution Constraint Violations :srcset: /examples/beginner/images/sphx_glr_plot_vyasarayani_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 139-155 .. code-block:: Python fig_y2, axes_y2 = plt.subplots(3, 1, layout='constrained') axes_y2[0].plot(time, y2_meas, '.k', time, initial_guess[num_nodes:-1], '.b', time, solution[num_nodes:-1], '.r', time, y2_sim, 'g') axes_y2[0].set_xlabel('Time [s]') axes_y2[0].set_ylabel('y2 [rad]') axes_y2[0].legend(legend) axes_y2[1].set_title('Initial Guess Constraint Violations') axes_y2[1].plot(prob.con(initial_guess)[num_nodes - 1:]) axes_y2[2].set_title('Solution Constraint Violations') axes_y2[2].plot(prob.con(solution)[num_nodes - 1:]) plt.show() .. image-sg:: /examples/beginner/images/sphx_glr_plot_vyasarayani_003.png :alt: Initial Guess Constraint Violations, Solution Constraint Violations :srcset: /examples/beginner/images/sphx_glr_plot_vyasarayani_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 32.873 seconds) .. _sphx_glr_download_examples_beginner_plot_vyasarayani.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_vyasarayani.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_vyasarayani.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_vyasarayani.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_