.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/beginner/plot_betts2003.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_betts2003.py: Parameter Identification: Betts & Huffman 2003 ============================================== This is the single parameter identification problem presented in section 7 of [Betts2003]_. Objectives ---------- - Show how to execute a parameter identification. - Show how to handle time being explicit in the equations of motion using a known trajectory set to time. .. GENERATED FROM PYTHON SOURCE LINES 16-26 .. code-block:: Python import numpy as np import sympy as sym import matplotlib.pyplot as plt from opty import Problem duration = 1.0 # seconds num_nodes = 100 interval = duration / (num_nodes - 1) .. GENERATED FROM PYTHON SOURCE LINES 27-28 Define the variables. .. GENERATED FROM PYTHON SOURCE LINES 28-34 .. code-block:: Python mu, p, t = sym.symbols('mu, p, t') y1, y2, T = sym.symbols('y1, y2, T', cls=sym.Function) state_symbols = (y1(t), y2(t)) constant_symbols = (mu, p) .. GENERATED FROM PYTHON SOURCE LINES 35-39 Define the symbolic equations of motion Use a trick here because time was explicit in the eoms. Set T(t) to a function of time and then pass in time as a known trajectory in the problem. .. GENERATED FROM PYTHON SOURCE LINES 39-44 .. code-block:: Python eom = sym.Matrix([ y1(t).diff(t) - y2(t), y2(t).diff(t) - mu**2 * y1(t) + (mu**2 + p**2) * sym.sin(p * T(t)), ]) .. GENERATED FROM PYTHON SOURCE LINES 45-46 Specify the known system parameters. .. GENERATED FROM PYTHON SOURCE LINES 46-48 .. code-block:: Python par_map = {mu: 60.0} .. GENERATED FROM PYTHON SOURCE LINES 49-50 Generate data representing measurements of the system's motion. .. GENERATED FROM PYTHON SOURCE LINES 50-56 .. code-block:: Python time = np.linspace(0.0, 1.0, num=num_nodes) y1_m = np.sin(np.pi * time) + np.random.normal(scale=0.05, size=len(time)) y2_m = np.pi * np.cos(np.pi * time) + np.random.normal(scale=0.05, size=len(time)) .. GENERATED FROM PYTHON SOURCE LINES 57-59 Specify the objective function and its gradient to minimize the sum of the squares of ``y1``. .. GENERATED FROM PYTHON SOURCE LINES 59-69 .. code-block:: Python def obj(free): return interval * np.sum((y1_m - free[:num_nodes])**2) def obj_grad(free): grad = np.zeros_like(free) grad[:num_nodes] = 2.0 * interval * (free[:num_nodes] - y1_m) return grad .. GENERATED FROM PYTHON SOURCE LINES 70-71 Specify the symbolic instance constraints, i.e. initial and end conditions. .. GENERATED FROM PYTHON SOURCE LINES 71-73 .. code-block:: Python instance_constraints = (y1(0.0), y2(0.0) - np.pi) .. GENERATED FROM PYTHON SOURCE LINES 74-75 Create an optimization problem. .. GENERATED FROM PYTHON SOURCE LINES 75-84 .. code-block:: Python prob = Problem(obj, obj_grad, eom, state_symbols, num_nodes, interval, known_parameter_map=par_map, known_trajectory_map={T(t): time}, instance_constraints=instance_constraints, time_symbol=t, integration_method='midpoint') .. GENERATED FROM PYTHON SOURCE LINES 85-86 Use a random positive initial guess. .. GENERATED FROM PYTHON SOURCE LINES 86-88 .. code-block:: Python initial_guess = np.random.randn(prob.num_free) .. 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) 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 94-95 Print results. .. GENERATED FROM PYTHON SOURCE LINES 95-104 .. code-block:: Python known_msg = "Known value of p = {}".format(np.pi) identified_msg = "Identified value of p = {}".format(solution[-1]) divider = '=' * max(len(known_msg), len(identified_msg)) print(divider) print(known_msg) print(identified_msg) print(divider) .. rst-class:: sphx-glr-script-out .. code-block:: none ========================================== Known value of p = 3.141592653589793 Identified value of p = 3.1409353268742914 ========================================== .. GENERATED FROM PYTHON SOURCE LINES 105-106 Plot results. .. GENERATED FROM PYTHON SOURCE LINES 106-122 .. code-block:: Python fig_y1, axes_y1 = plt.subplots(3, 1, layout='constrained') legend = ['measured', 'initial guess', 'direct collocation solution'] axes_y1[0].plot(time, y1_m, '.k', time, initial_guess[:num_nodes], '.b', time, solution[:num_nodes], '.r') axes_y1[0].set_xlabel('Time [s]') axes_y1[0].set_ylabel('y1') 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_betts2003_001.png :alt: Initial Guess Constraint Violations, Solution Constraint Violations :srcset: /examples/beginner/images/sphx_glr_plot_betts2003_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 123-137 .. code-block:: Python fig_y2, axes_y2 = plt.subplots(3, 1, layout='constrained') axes_y2[0].plot(time, y2_m, '.k', time, initial_guess[num_nodes:-1], '.b', time, solution[num_nodes:-1], '.r') axes_y2[0].set_xlabel('Time [s]') axes_y2[0].set_ylabel('y2') 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:]) .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts2003_002.png :alt: Initial Guess Constraint Violations, Solution Constraint Violations :srcset: /examples/beginner/images/sphx_glr_plot_betts2003_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 138-140 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts2003_003.png :alt: Constraint violations :srcset: /examples/beginner/images/sphx_glr_plot_betts2003_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 141-144 .. code-block:: Python _ = prob.plot_trajectories(solution) plt.show() .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts2003_004.png :alt: State Trajectories, Input Trajectories :srcset: /examples/beginner/images/sphx_glr_plot_betts2003_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 31.614 seconds) .. _sphx_glr_download_examples_beginner_plot_betts2003.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_betts2003.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts2003.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts2003.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_