.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\plot_betts_10_90_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_90_opty.py: Nonconvex Delay (betts_10_90) ============================= This is example 10.90 from *John T. Betts, Practical Methods for Optimal Control Using NonlinearProgramming*, 3rd edition, Chapter 10: Test Problems. As reference, Betts gives: *H. Mauer, Theory and applications of optimal control problems with control and state delays*, Proceedings of the 53rd Australian Mathematical Confewrence, Adelaide, 2009. This simulation shows the use of instance constrains of the form: :math:`x_i(t_0) = x_j(t_f)`. **States** - :math:`x_1,....., x_N` : state variables **Controls** - :math:`u_1,....., u_N` : control variables .. GENERATED FROM PYTHON SOURCE LINES 28-35 .. code-block:: Python import numpy as np import sympy as sm import sympy.physics.mechanics as me from opty import Problem from opty.utils import create_objective_function, MathJaxRepr .. GENERATED FROM PYTHON SOURCE LINES 36-38 Equations of Motion ------------------- .. GENERATED FROM PYTHON SOURCE LINES 38-54 .. code-block:: Python t = me.dynamicsymbols._t N = 20 tf = 0.1 sigma = 5 x = me.dynamicsymbols(f'x:{N}') u = me.dynamicsymbols(f'u:{N}') eom = sm.Matrix([ *[-x[i].diff(t) + 1**2 - u[i] for i in range(sigma)], *[-x[i].diff(t) + x[i-sigma]**2 - u[i] for i in range(sigma, N)], ]) MathJaxRepr(eom) .. raw:: html
$$\left[\begin{matrix}- u_{0} - \dot{x}_{0} + 1\\- u_{1} - \dot{x}_{1} + 1\\- u_{2} - \dot{x}_{2} + 1\\- u_{3} - \dot{x}_{3} + 1\\- u_{4} - \dot{x}_{4} + 1\\- u_{5} + x_{0}^{2} - \dot{x}_{5}\\- u_{6} + x_{1}^{2} - \dot{x}_{6}\\- u_{7} + x_{2}^{2} - \dot{x}_{7}\\- u_{8} + x_{3}^{2} - \dot{x}_{8}\\- u_{9} + x_{4}^{2} - \dot{x}_{9}\\- u_{10} + x_{5}^{2} - \dot{x}_{10}\\- u_{11} + x_{6}^{2} - \dot{x}_{11}\\- u_{12} + x_{7}^{2} - \dot{x}_{12}\\- u_{13} + x_{8}^{2} - \dot{x}_{13}\\- u_{14} + x_{9}^{2} - \dot{x}_{14}\\- u_{15} + x_{10}^{2} - \dot{x}_{15}\\- u_{16} + x_{11}^{2} - \dot{x}_{16}\\- u_{17} + x_{12}^{2} - \dot{x}_{17}\\- u_{18} + x_{13}^{2} - \dot{x}_{18}\\- u_{19} + x_{14}^{2} - \dot{x}_{19}\end{matrix}\right]$$


.. GENERATED FROM PYTHON SOURCE LINES 55-57 Define and Solve the Optimization Problem ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 57-64 .. code-block:: Python num_nodes = 501 t0 = 0.0 interval_value = (tf - t0) / (num_nodes - 1) state_symbols = x unkonwn_input_trajectories = u .. GENERATED FROM PYTHON SOURCE LINES 65-66 Specify the objective function and form the gradient. .. GENERATED FROM PYTHON SOURCE LINES 66-77 .. code-block:: Python objective = sm.Integral(sum([x[i]**2 + u[i]**2 for i in range(N)]), t) obj, obj_grad = create_objective_function( objective, state_symbols, unkonwn_input_trajectories, tuple(), num_nodes, interval_value ) .. GENERATED FROM PYTHON SOURCE LINES 78-79 Instance Constraints and Bounds. .. GENERATED FROM PYTHON SOURCE LINES 79-89 .. code-block:: Python instance_constraints = ( x[0].func(t0) - 1.0, *[x[i].func(t0) - x[i-1].func(tf) for i in range(1, N)], ) limit_value = np.inf bounds = { x[i]: (0.7, limit_value) for i in range(0, N) } .. GENERATED FROM PYTHON SOURCE LINES 90-92 Set Up the Control Problem and Solve it --------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 92-105 .. code-block:: Python prob = Problem( obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, time_symbol=t, ) .. GENERATED FROM PYTHON SOURCE LINES 106-107 Rough initial guess. .. GENERATED FROM PYTHON SOURCE LINES 107-109 .. code-block:: Python initial_guess = np.ones(prob.num_free) * 0.1 .. GENERATED FROM PYTHON SOURCE LINES 110-111 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 111-119 .. code-block:: Python solution, info = prob.solve(initial_guess) print(info['status_msg']) Jstr = 2.79685770 print(f"Objective value achieved: {info['obj_val']:.4f}, as per the book ", f"it is {Jstr:.4f}, so the error is: ", f"{(info['obj_val'] - Jstr) / Jstr * 100:.4f} % ") print('\n') .. 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).' Objective value achieved: 2.7967, as per the book it is 2.7969, so the error is: -0.0040 % .. GENERATED FROM PYTHON SOURCE LINES 120-121 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 121-123 .. code-block:: Python _ = prob.plot_trajectories(solution, show_bounds=True) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_90_opty_001.png :alt: State Trajectories, Input Trajectories :srcset: /examples/images/sphx_glr_plot_betts_10_90_opty_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 124-125 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 125-127 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_90_opty_002.png :alt: Constraint violations :srcset: /examples/images/sphx_glr_plot_betts_10_90_opty_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 128-129 Plot the objective function. .. GENERATED FROM PYTHON SOURCE LINES 129-132 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/images/sphx_glr_plot_betts_10_90_opty_003.png :alt: Objective Value :srcset: /examples/images/sphx_glr_plot_betts_10_90_opty_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 133-134 Are the inequality constraints strictly satisfied? .. GENERATED FROM PYTHON SOURCE LINES 134-141 .. code-block:: Python min_x = np.min(solution[0: N*num_nodes-1]) if min_x >= 0.7: print(f"Minimal value of the x\u1D62 is: {min_x:.12f} >= 0.7", "hence satisfied") else: print(f"Minimal value of the x\u1D62 is: {min_x:.12f} < 0.7", "hence not satisfied") .. rst-class:: sphx-glr-script-out .. code-block:: none Minimal value of the xᵢ is: 0.700000012530 >= 0.7 hence satisfied .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 14.750 seconds) .. _sphx_glr_download_examples_plot_betts_10_90_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_90_opty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_90_opty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_90_opty.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_