.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\plot_betts_10_133_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_133_opty.py: Two-Strain Tuberculosis Model (betts_10_133) ============================================ This is example 10.133 from *John T.Betts, Practical Methods for Optimal Control Using NonlinearProgramming*, 3rd edition, Chapter 10: Test Problems. More details may be found in chapter 8.17 of this book. **States** - :math:`S, T, L_1, I_1, L_", I_2` : state variables **Controls** - :math:`u_1, u_2` : control variables .. GENERATED FROM PYTHON SOURCE LINES 20-28 .. 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 create_objective_function, MathJaxRepr .. GENERATED FROM PYTHON SOURCE LINES 29-31 Equations of Motion. -------------------- .. GENERATED FROM PYTHON SOURCE LINES 31-55 .. code-block:: Python start = time.time() t = me.dynamicsymbols._t S, T, L1, I1, L2, I2 = me.dynamicsymbols('S T L1 I1 L2 I2') u1, u2 = me.dynamicsymbols('u1, u2') beta1, d2, r2, betastar, beta2, k1, p, B1, nu, k2, q, B2, d1, r1, N = \ sm.symbols('beta1 d2 r2 betastar beta2 k1 p B1 nu k2 q B2 d1 r1 N') Lambda = nu * N eom = sm.Matrix([ -S.diff(t) + Lambda - beta1*S*I1/N - betastar*S*I2/N - nu*S, -T.diff(t) + (u1*r1*L1 - nu*T + (1-(1-u2)*(p+q))*r2*I1 - beta2*T*I1/N - betastar*T*I2/N), -L1.diff(t) + (beta1*S*I1/N - (nu+k1)*L1 - u1*r1*L1 + (1-u2)*p*r2*I1 + beta2*T*I1/N - betastar*L1*I2/N), -L2.diff(t) + (1-u2)*q*r2*I1 - (nu+k2)*L2 + betastar*(S+L1+T)*I2/N, -I1.diff(t) + k1*L1 - (nu+d1)*I1 - r2*I1, -I2.diff(t) + k2*L2 - (nu+d2)*I2, ]) MathJaxRepr(eom) .. raw:: html
$$\left[\begin{matrix}N \nu - \nu S - \dot{S} - \frac{\beta_{1} I_{1} S}{N} - \frac{betastar I_{2} S}{N}\\- \nu T + r_{1} L_{1} u_{1} + r_{2} \left(- \left(1 - u_{2}\right) \left(p + q\right) + 1\right) I_{1} - \dot{T} - \frac{\beta_{2} I_{1} T}{N} - \frac{betastar I_{2} T}{N}\\p r_{2} \left(1 - u_{2}\right) I_{1} - r_{1} L_{1} u_{1} - \left(k_{1} + \nu\right) L_{1} - \dot{L}_{1} + \frac{\beta_{1} I_{1} S}{N} + \frac{\beta_{2} I_{1} T}{N} - \frac{betastar I_{2} L_{1}}{N}\\q r_{2} \left(1 - u_{2}\right) I_{1} - \left(k_{2} + \nu\right) L_{2} - \dot{L}_{2} + \frac{betastar \left(L_{1} + S + T\right) I_{2}}{N}\\k_{1} L_{1} - r_{2} I_{1} - \left(d_{1} + \nu\right) I_{1} - \dot{I}_{1}\\k_{2} L_{2} - \left(d_{2} + \nu\right) I_{2} - \dot{I}_{2}\end{matrix}\right]$$


.. GENERATED FROM PYTHON SOURCE LINES 56-58 Define and Solve the Optimization Problem ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 58-84 .. code-block:: Python num_nodes = 801 t0, tf = 0.0, 5.0 interval_value = tf / (num_nodes - 1) state_symbols = (S, T, L1, L2, I1, I2) unkonwn_input_trajectories = (u1, u2) par_map = { beta1: 13, d2: 0.0, r2: 1.0, betastar: 0.029, beta2: 13, k1: 0.5, p: 0.4, B1: 50, nu: 0.0143, k2: 1.0, q: 0.1, B2: 500, d1: 0.0, r1: 2.0, N: 30000, } .. 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 objective = sm.Integral((L2 + I2 + 0.5*B1*u1**2 + 0.5*B2*u2**2).subs(par_map), t) obj, obj_grad = create_objective_function( objective, state_symbols, unkonwn_input_trajectories, tuple(), num_nodes, interval_value, ) .. GENERATED FROM PYTHON SOURCE LINES 99-100 Instance constraints and bounds. .. GENERATED FROM PYTHON SOURCE LINES 100-114 .. code-block:: Python instance_constraints = ( S.func(t0) - 76*par_map[N]/120, T.func(t0) - par_map[N]/120, L1.func(t0) - 36*par_map[N]/120, I1.func(t0) - 4*par_map[N]/120, L2.func(t0) - 2*par_map[N]/120, I2.func(t0) - par_map[N]/120, ) bounds = { u1: (0.05, 0.95), u2: (0.05, 0.95), } .. GENERATED FROM PYTHON SOURCE LINES 115-116 Create the optimization problem. .. GENERATED FROM PYTHON SOURCE LINES 116-129 .. code-block:: Python prob = Problem( obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, known_parameter_map=par_map, bounds=bounds, time_symbol=t, ) .. GENERATED FROM PYTHON SOURCE LINES 130-132 Acceptable tolerance and iteration settings are lowered, to get convergence with a reasonable number of iterations. .. GENERATED FROM PYTHON SOURCE LINES 132-136 .. code-block:: Python prob.add_option('max_iter', 4000) prob.add_option('acceptable_tol', 1e-2) prob.add_option('acceptable_iter', 5) .. GENERATED FROM PYTHON SOURCE LINES 137-138 Rough initial guess. .. GENERATED FROM PYTHON SOURCE LINES 138-140 .. code-block:: Python initial_guess = np.ones(prob.num_free) .. GENERATED FROM PYTHON SOURCE LINES 141-142 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 142-160 .. code-block:: Python start_solve = time.time() solution, info = prob.solve(initial_guess) end_solve = time.time() print(f"Solve time: {end_solve - start_solve:.3f} seconds") print(info['status_msg']) Jstar = 5152.07310 print(f"Objective value achieved: {info['obj_val']:.4f}, ", f"as per the book it is {Jstar:.4f}, so the deviation is: ", f"{(info['obj_val'] - Jstar) / Jstar*100:.3f} %") Tbstar = 1123 print(f"Individuals infected with resistant Tb = ", f"{solution[4*num_nodes-1] + solution[6*num_nodes-1]:.3f}, ", f"vs. {Tbstar} from the book, hence deviation: ", f"{(solution[4*num_nodes-1] + solution[6*num_nodes-1] - Tbstar) / Tbstar*100:.3f} %") .. rst-class:: sphx-glr-script-out .. code-block:: none Solve time: 65.949 seconds b'Algorithm stopped at a point that was converged, not to "desired" tolerances, but to "acceptable" tolerances (see the acceptable-... options).' Objective value achieved: 5154.9859, as per the book it is 5152.0731, so the deviation is: 0.057 % Individuals infected with resistant Tb = 1121.198, vs. 1123 from the book, hence deviation: -0.160 % .. GENERATED FROM PYTHON SOURCE LINES 161-162 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 162-164 .. code-block:: Python _ = prob.plot_trajectories(solution, show_bounds=True) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_133_opty_001.png :alt: State Trajectories, Input Trajectories :srcset: /examples/images/sphx_glr_plot_betts_10_133_opty_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 165-166 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 166-168 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_133_opty_002.png :alt: Constraint violations :srcset: /examples/images/sphx_glr_plot_betts_10_133_opty_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 169-170 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 170-171 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/images/sphx_glr_plot_betts_10_133_opty_003.png :alt: Objective Value :srcset: /examples/images/sphx_glr_plot_betts_10_133_opty_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 14.299 seconds) .. _sphx_glr_download_examples_plot_betts_10_133_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_133_opty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_133_opty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_133_opty.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_