.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\plot_betts_10_141_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_141_opty.py: Tumor Antigionesis (betts_10_141) ================================= This is example 10.141 from *John T. Betts, Practical Methods for Optimal Control Using NonlinearProgramming*, 3rd edition, Chapter 10: Test Problems. It is described in more detail in section 8.17 of the book. Notes ----- - The variable :math:`\textrm{hilfs}` ist itroduced to eneforce :math:`y(tf) \leq A` as an instance inequality constraint, which opty presently does not support. (of course, looking at the equations of motion, it is clear that :math:`\dfrac{d}{dt}y(t) \geq 0`, so bounding y(t) would have been sufficient - which opty does support) - Unless :math:`\left[ \textrm{num}_{\textrm{nodes}} - 1 \right] \cdot h_{\textrm{fast}}` is chosen only a bit larger than the result, convergence becomes difficult. (trial and error approach) **States** - :math:`p, q, y` : state variables **Controls** - :math:`u` : control variable .. GENERATED FROM PYTHON SOURCE LINES 30-37 .. code-block:: Python import numpy as np import sympy as sm import sympy.physics.mechanics as me import time from opty import Problem from opty.utils import MathJaxRepr .. GENERATED FROM PYTHON SOURCE LINES 38-40 Equations of Motion -------------------- .. GENERATED FROM PYTHON SOURCE LINES 40-48 .. code-block:: Python t = me.dynamicsymbols._t p, q, y = me.dynamicsymbols('p, q, y') u = me.dynamicsymbols('u') hilfs = me.dynamicsymbols('hilfs') h_fast = sm.symbols('h_fast') .. GENERATED FROM PYTHON SOURCE LINES 49-50 Parameters .. GENERATED FROM PYTHON SOURCE LINES 50-63 .. code-block:: Python chi = 0.084 G = 0.15 b = 5.85 nu = 0.02 d = 0.00873 a = 75.0 A = 15.0 pbar = ((b - nu)/d)**(3/2) qbar = pbar p0 = pbar/2.0 q0 = qbar/4.0 .. GENERATED FROM PYTHON SOURCE LINES 64-65 Equations of motion. .. GENERATED FROM PYTHON SOURCE LINES 65-74 .. code-block:: Python eom = sm.Matrix([ -p.diff(t) - chi*p*sm.ln(p/q), -q.diff(t) + q*(b - (nu + d*p**(2/3) + G*u)), -y.diff(t) + u, hilfs * A - y ]) MathJaxRepr(eom) .. raw:: html
$$\left[\begin{matrix}- 0.084 p \log{\left(\frac{p}{q} \right)} - \dot{p}\\\left(- 0.00873 p^{0.666666666666667} - 0.15 u + 5.83\right) q - \dot{q}\\u - \dot{y}\\15.0 hilfs - y\end{matrix}\right]$$


.. GENERATED FROM PYTHON SOURCE LINES 75-77 Define and Solve the Optimization Problem. ------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 77-84 .. code-block:: Python num_nodes = 501 interval_value = h_fast t0, tf = 0*h_fast, (num_nodes-1)*h_fast state_symbols = (p, q, y) .. GENERATED FROM PYTHON SOURCE LINES 85-86 Specify the objective function and form the gradient. .. GENERATED FROM PYTHON SOURCE LINES 86-98 .. code-block:: Python def obj(free): return free[num_nodes-1] def obj_grad(free): grad = np.zeros_like(free) grad[num_nodes-1] = 1.0 return grad .. GENERATED FROM PYTHON SOURCE LINES 99-100 Define the instance constraints, bounds, and eom bounds. .. GENERATED FROM PYTHON SOURCE LINES 100-117 .. code-block:: Python instance_constraints = ( p.func(t0) - p0, q.func(t0) - q0, y.func(t0), hilfs.func(tf) - 1.0, ) bounds = { h_fast: (0.0, 2.5 / (num_nodes-1)), p: (0.01, pbar), q: (0.01, qbar), y: (0.0, np.inf), u: (0, a), } eom_bounds = {3: (0.0, np.inf)} .. GENERATED FROM PYTHON SOURCE LINES 118-119 Set up the Problem. .. GENERATED FROM PYTHON SOURCE LINES 119-131 .. code-block:: Python prob = Problem( obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, eom_bounds=eom_bounds, time_symbol=t ) .. GENERATED FROM PYTHON SOURCE LINES 132-133 Give some rough estimates for the trajectories. .. GENERATED FROM PYTHON SOURCE LINES 133-135 .. code-block:: Python initial_guess = np.ones(prob.num_free) .. GENERATED FROM PYTHON SOURCE LINES 136-137 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 137-141 .. code-block:: Python start = time.time() solution, info = prob.solve(initial_guess) print(f"Solved in {time.time() - start:.2f} seconds.") .. rst-class:: sphx-glr-script-out .. code-block:: none Solved in 29.24 seconds. .. GENERATED FROM PYTHON SOURCE LINES 142-143 Print some information about the solution. .. GENERATED FROM PYTHON SOURCE LINES 143-152 .. code-block:: Python print(info['status_msg']) tfstar = 1.1961336 print(f"Duration is: {solution[-1]*(num_nodes-1):.4f} " f"as per the book it is {tfstar:.3f}, so the deviation is: " f"{(solution[-1]*(num_nodes-1) - tfstar)/tfstar*100:.3f} %") Jstar = 7571.6712 print(f"p(tf) = {solution[num_nodes-1]:.4f}" + f"as per the book it is {Jstar:.4f}, so the deviation is: " + f"{(solution[num_nodes-1] - Jstar)/Jstar*100:.3f} %") .. 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).' Duration is: 1.1788 as per the book it is 1.196, so the deviation is: -1.450 % p(tf) = 7586.6171as per the book it is 7571.6712, so the deviation is: 0.197 % .. GENERATED FROM PYTHON SOURCE LINES 153-154 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 154-156 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_141_opty_001.png :alt: State Trajectories, Input Trajectories :srcset: /examples/images/sphx_glr_plot_betts_10_141_opty_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 157-158 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 158-160 .. code-block:: Python _ = prob.plot_constraint_violations(solution, subplots=True, show_bounds=True) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_141_opty_002.png :alt: Constraint violations Values of bounded EoMs :srcset: /examples/images/sphx_glr_plot_betts_10_141_opty_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 161-162 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 162-163 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/images/sphx_glr_plot_betts_10_141_opty_003.png :alt: Objective Value :srcset: /examples/images/sphx_glr_plot_betts_10_141_opty_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 164-165 sphinx_gallery_thumbnail_number = 2 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 37.711 seconds) .. _sphx_glr_download_examples_plot_betts_10_141_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_141_opty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_141_opty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_141_opty.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_