.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\plot_betts_10_1_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_betts_10_1_opty.py: Stiff Set of DAEs (betts-10-1) ============================== This is example 10.1 from *John T. Betts Practical Methods for Optimal Control Using NonlinearProgramming*, 3rd edition, Chapter 10: Test Problems. More details are in section 4.11.7 of the book, where it says it was designed specifically to be **hard** to solve. Notes: ------ - The equations of motion do not seem to have any 'physical' meaning. - Only by relenting the bound on the inequality constraint from :math:`y^Ty \geq \textrm{function(time, parameters)}` to :math:`y^Ty \geq \textrm{function(time, parameters)} - 0.01`, did the problem reach convergence. - seems very sensitive to the number of nodes used. **States** - :math:`y_1, .., y_4` : state variables as per the book - :math:`T` : variable to keep track of the time. **Specifieds** - :math:`u_1, u_2` : control variables .. GENERATED FROM PYTHON SOURCE LINES 32-41 .. code-block:: Python import time 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, MathJaxRepr .. GENERATED FROM PYTHON SOURCE LINES 43-46 Equations of Motion ------------------- .. GENERATED FROM PYTHON SOURCE LINES 46-50 .. code-block:: Python t = me.dynamicsymbols._t y1, y2, y3, y4 = me.dynamicsymbols('y1 y2 y3 y4') u1, u2 = me.dynamicsymbols('u1 u2') .. GENERATED FROM PYTHON SOURCE LINES 51-53 As the time occurs explicitly in the equations of motion, and presently opty does not support that, a new state T(t) is introduced to keep track of time. .. GENERATED FROM PYTHON SOURCE LINES 53-56 .. code-block:: Python T = sm.symbols('T', cls=sm.Function) .. GENERATED FROM PYTHON SOURCE LINES 57-58 Get the rhs of the algebraic inequality. .. GENERATED FROM PYTHON SOURCE LINES 58-68 .. code-block:: Python def p(t, a, b): return sm.exp(-b*(t - a)**2) terrain = (3.0*(p(T(t), 3, 12) + p(T(t), 6, 10) + p(T(t), 10, 6)) + 8.0*p(T(t), 15, 4) + 0.01) .. GENERATED FROM PYTHON SOURCE LINES 69-70 Formulate the equations of motion as per the book. .. GENERATED FROM PYTHON SOURCE LINES 70-82 .. code-block:: Python eom = sm.Matrix([ -y1.diff(t) - 10*y1 + u1 + u2, -y2.diff(t) - 2*y2 + u1 - 2*u2, -y3.diff(t) - 3*y3 + 5*y4 + u1 - u2, -y4.diff(t) + 5*y3 - 3*y4 + u1 + 3*u2, y1**2 + y2**2 + y3**2 + y4**2 - terrain, T(t).diff(t) - 1, ]) MathJaxRepr(eom) .. raw:: html
$$\left[\begin{matrix}u_{1} + u_{2} - 10 y_{1} - \dot{y}_{1}\\u_{1} - 2 u_{2} - 2 y_{2} - \dot{y}_{2}\\u_{1} - u_{2} - 3 y_{3} + 5 y_{4} - \dot{y}_{3}\\u_{1} + 3 u_{2} + 5 y_{3} - 3 y_{4} - \dot{y}_{4}\\y_{1}^{2} + y_{2}^{2} + y_{3}^{2} + y_{4}^{2} - 0.01 - 3.0 e^{- 12 \left(T - 3\right)^{2}} - 3.0 e^{- 10 \left(T - 6\right)^{2}} - 3.0 e^{- 6 \left(T - 10\right)^{2}} - 8.0 e^{- 4 \left(T - 15\right)^{2}}\\\dot{T} - 1\end{matrix}\right]$$


.. GENERATED FROM PYTHON SOURCE LINES 83-85 Set up and Solve the Optimization Problem ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 85-95 .. code-block:: Python t0, tf = 0.0, 20.0 num_nodes = 501 interval_value = (tf - t0) / (num_nodes - 1) times = np.linspace(t0, tf, num_nodes) state_symbols = (y1, y2, y3, y4, T(t)) specified_symbols = (u1, u2) integration_method = 'backward euler' .. GENERATED FROM PYTHON SOURCE LINES 96-97 Specify the objective function and form the gradient. .. GENERATED FROM PYTHON SOURCE LINES 97-110 .. code-block:: Python obj_func = sm.Integral(100.0*(y1**2 + y2**2 + y3**2 + y4**2) + 0.01*(u1**2 + u2**2), t) obj, obj_grad = create_objective_function( obj_func, state_symbols, specified_symbols, tuple(), num_nodes, node_time_interval=interval_value, integration_method=integration_method, ) .. GENERATED FROM PYTHON SOURCE LINES 111-112 Specify the symbolic instance constraints, as per the example. .. GENERATED FROM PYTHON SOURCE LINES 112-124 .. code-block:: Python instance_constraints = ( T(t0) - t0, y1.func(t0) - 2.0, y2.func(t0) - 1.0, y3.func(t0) - 2.0, y4.func(t0) - 1.0, y1.func(tf) - 2.0, y2.func(tf) - 3.0, y3.func(tf) - 1.0, y4.func(tf) + 2.0, ) .. GENERATED FROM PYTHON SOURCE LINES 125-127 Specify the equation of motion bounds. Here the left side of the interval is changed from 0.0 to -0.01 to get convergence. .. GENERATED FROM PYTHON SOURCE LINES 127-129 .. code-block:: Python eom_bounds = {4: (-0.01, np.inf)} .. GENERATED FROM PYTHON SOURCE LINES 130-131 Create the optimization problem, set some options. .. GENERATED FROM PYTHON SOURCE LINES 131-148 .. code-block:: Python prob = Problem( obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, eom_bounds=eom_bounds, integration_method=integration_method, time_symbol=t, ) prob.add_option('nlp_scaling_method', 'gradient-based') prob.add_option('max_iter', 8000) .. GENERATED FROM PYTHON SOURCE LINES 149-150 Give some rough estimates for the trajectories. .. GENERATED FROM PYTHON SOURCE LINES 150-152 .. code-block:: Python initial_guess = np.zeros(prob.num_free) .. GENERATED FROM PYTHON SOURCE LINES 153-154 Solve the optimization problem. .. GENERATED FROM PYTHON SOURCE LINES 154-166 .. code-block:: Python start = time.time() solution, info = prob.solve(initial_guess) print(info['status_msg']) print(f"Objective value achieved: {info['obj_val']:.4f}, as per the book ", f"it is {2030.85609}, so the difference is: " f"{(info['obj_val'] - 2030.85609)/2030.85609*100:.3f} % \n ") end = time.time() print(f'Total time needed: {end - start:.2f} seconds') .. rst-class:: sphx-glr-script-out .. code-block:: none b'Algorithm stopped at a point that was converged, not to "desired" tolerances, but to "acceptable" tolerances (see the acceptable-... options).' Objective value achieved: 2009.5722, as per the book it is 2030.85609, so the difference is: -1.048 % Total time needed: 146.06 seconds .. GENERATED FROM PYTHON SOURCE LINES 167-168 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 168-170 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_1_opty_001.png :alt: State Trajectories, Input Trajectories :srcset: /examples/images/sphx_glr_plot_betts_10_1_opty_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 171-172 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 172-174 .. code-block:: Python _ =prob.plot_constraint_violations(solution, subplots=True) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_1_opty_002.png :alt: Constraint violations :srcset: /examples/images/sphx_glr_plot_betts_10_1_opty_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 175-176 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 176-177 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/images/sphx_glr_plot_betts_10_1_opty_003.png :alt: Objective Value :srcset: /examples/images/sphx_glr_plot_betts_10_1_opty_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 37.207 seconds) .. _sphx_glr_download_examples_plot_betts_10_1_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_betts_10_1_opty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_1_opty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_1_opty.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_