.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/beginner/plot_betts_10_103_104.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_betts_10_103_104.py: Compare DAE vs. ODE Formulation =============================== Objective --------- - To compare the results of the DAE and ODE formulations of the same problem. In one case one equation of motion is an algebraic equation (case 10.103), in the other case :math:`\dfrac{d}{dt}(\textrm{algebraic equation})` is used. Introduction ------------ This is example 10.103 from [Betts2010]_'s test problems. It has four differential equations and one algebraic equation. This example compares two solutions, the first uses the algebraic equation directly and the second uses the derivative of the algebraic equation, i.e. converts the DAEs to ODEs. As expected, the DAE formulation gives a better result than the ODE formulation. The ODE formulation seems to run a bit faster. **States** - :math:`y_0,...,y_4` : state variables **Specifieds** - :math:`u` : control variable .. GENERATED FROM PYTHON SOURCE LINES 34-41 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 42-46 Information Common to Both Formulations --------------------------------------- Define the time varying variables. .. GENERATED FROM PYTHON SOURCE LINES 46-53 .. code-block:: Python t = me.dynamicsymbols._t y = me.dynamicsymbols('y0, y1, y2, y3, y4') u = me.dynamicsymbols('u') # Define the constant parameters. g = 9.81 .. GENERATED FROM PYTHON SOURCE LINES 54-55 Set up the time intervals. .. GENERATED FROM PYTHON SOURCE LINES 55-61 .. code-block:: Python tf = 3.0 num_nodes = 751 t0, tf = 0.0, tf interval_value = (tf - t0)/(num_nodes - 1) .. GENERATED FROM PYTHON SOURCE LINES 62-63 Specify the objective function and form the gradient. .. GENERATED FROM PYTHON SOURCE LINES 63-75 .. code-block:: Python obj_func = sm.Integral(u**2, t) sm.pprint(obj_func) obj, obj_grad = create_objective_function( obj_func, y, (u,), tuple(), num_nodes, node_time_interval=interval_value, time_symbol=t, ) .. rst-class:: sphx-glr-script-out .. code-block:: none ⌠ ⎮ 2 ⎮ u (t) dt ⌡ .. GENERATED FROM PYTHON SOURCE LINES 76-77 Specify the symbolic instance constraints. .. GENERATED FROM PYTHON SOURCE LINES 77-84 .. code-block:: Python instance_constraints = ( y[0].func(t0) - 1, *[y[i].func(t0) - 0 for i in range(1, 5)], y[0].func(tf) - 0, y[2].func(tf) - 0, ) .. GENERATED FROM PYTHON SOURCE LINES 85-86 Specify the bounds. .. GENERATED FROM PYTHON SOURCE LINES 86-94 .. code-block:: Python bounds = { y[0]: (-5, 5), y[1]: (-5, 5), y[2]: (-5, 5), y[3]: (-5, 5), y[4]: (-1, 15), } .. GENERATED FROM PYTHON SOURCE LINES 95-99 DAE Formulation --------------- Define the differential algebraic equations of motion. .. GENERATED FROM PYTHON SOURCE LINES 99-108 .. code-block:: Python eom = sm.Matrix([ -y[0].diff(t) + y[2], -y[1].diff(t) + y[3], -y[2].diff(t) - 2*y[4]*y[0] + u*y[1], -y[3].diff(t) - g - 2*y[4]*y[1] - u*y[0], y[2]**2 + y[3]**2 - 2*y[4] - g*y[1], ]) sm.pprint(eom) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ⎤ ⎢ y₂(t) - ──(y₀(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ y₃(t) - ──(y₁(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ u(t)⋅y₁(t) - 2⋅y₀(t)⋅y₄(t) - ──(y₂(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢-u(t)⋅y₀(t) - 2⋅y₁(t)⋅y₄(t) - ──(y₃(t)) - 9.81⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ 2 2 ⎥ ⎣ -9.81⋅y₁(t) + y₂ (t) + y₃ (t) - 2⋅y₄(t) ⎦ .. GENERATED FROM PYTHON SOURCE LINES 109-110 Create the optimization problem and set any options. .. GENERATED FROM PYTHON SOURCE LINES 110-138 .. code-block:: Python prob = Problem( obj, obj_grad, eom, y, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, ) prob.add_option('nlp_scaling_method', 'gradient-based') # Give rough estimates for the x and y trajectories. initial_guess = np.zeros(prob.num_free) # Find the optimal solution. solution, info = prob.solve(initial_guess) print(info['status_msg']) if info['obj_val'] < 12.8738850: msg = 'opty with DAE gets a better result' else: msg = 'opty with DAE gets a worse result' print(f'Minimal objective value achieved: {info["obj_val"]:.4f}, ' + f'as per the book it is {12.8738850}, {msg} ') obj_DAE = info["obj_val"] .. 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).' Minimal objective value achieved: 12.0434, as per the book it is 12.873885, opty with DAE gets a better result .. GENERATED FROM PYTHON SOURCE LINES 139-140 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 140-142 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_001.png :alt: State Trajectories, Input Trajectories :srcset: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 143-144 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 144-146 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_002.png :alt: Constraint violations :srcset: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 147-148 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 148-150 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_003.png :alt: Objective Value :srcset: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 151-156 ODE Formulation --------------- Differentiate the algebraic equation with respect to time and add it to the equations of motion instead of the algebraic equation, as before. .. GENERATED FROM PYTHON SOURCE LINES 156-166 .. code-block:: Python eom = sm.Matrix([ -y[0].diff(t) + y[2], -y[1].diff(t) + y[3], -y[2].diff(t) - 2*y[4]*y[0] + u*y[1], -y[3].diff(t) - g - 2*y[4]*y[1] - u*y[0], -y[4].diff(t) + y[2]*y[2].diff(t) + y[3]*y[3].diff(t) - g*y[1].diff(t)/2, ]) sm.pprint(eom) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ⎤ ⎢ y₂(t) - ──(y₀(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ y₃(t) - ──(y₁(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ u(t)⋅y₁(t) - 2⋅y₀(t)⋅y₄(t) - ──(y₂(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -u(t)⋅y₀(t) - 2⋅y₁(t)⋅y₄(t) - ──(y₃(t)) - 9.81 ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ d d d d ⎥ ⎢y₂(t)⋅──(y₂(t)) + y₃(t)⋅──(y₃(t)) - 4.905⋅──(y₁(t)) - ──(y₄(t))⎥ ⎣ dt dt dt dt ⎦ .. GENERATED FROM PYTHON SOURCE LINES 167-168 Create the optimization problem and set any options. .. GENERATED FROM PYTHON SOURCE LINES 168-179 .. code-block:: Python prob = Problem( obj, obj_grad, eom, y, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, ) .. GENERATED FROM PYTHON SOURCE LINES 180-181 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 181-192 .. code-block:: Python solution, info = prob.solve(initial_guess) print(info['status_msg']) if info['obj_val'] < 12.8738850: msg = 'opty with ODE gets a better result' else: msg = 'opty with ODE gets a worse result' print(f'Minimal objective value achieved: {info["obj_val"]:.4f}, ' + f'as per the book it is {12.8738850}, {msg} ') obj_ODE = info["obj_val"] .. 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).' Minimal objective value achieved: 12.5526, as per the book it is 12.873885, opty with ODE gets a better result .. GENERATED FROM PYTHON SOURCE LINES 193-194 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 194-196 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_004.png :alt: State Trajectories, Input Trajectories :srcset: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 197-198 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 198-200 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_005.png :alt: Constraint violations :srcset: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 201-202 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 202-204 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_006.png :alt: Objective Value :srcset: /examples/beginner/images/sphx_glr_plot_betts_10_103_104_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 205-207 Compare the results -------------------- .. GENERATED FROM PYTHON SOURCE LINES 207-214 .. code-block:: Python if obj_DAE < obj_ODE: value = (obj_ODE - obj_DAE)/obj_ODE*100 print(f'DAE formulation gives a better result by {value:.2f} %') else: value = (obj_DAE - obj_ODE)/obj_DAE*100 print(f'ODE formulation gives a better result by {value:.2f} %') .. rst-class:: sphx-glr-script-out .. code-block:: none DAE formulation gives a better result by 4.06 % .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 4.344 seconds) .. _sphx_glr_download_examples_beginner_plot_betts_10_103_104.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_betts_10_103_104.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_103_104.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_103_104.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_