.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples\plot_betts_10_58_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_58_opty.py: Heat Equation (betts_10_58) =========================== This is example 10.58 from *John T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming*, 3rd edition, chapter 10: Test Problems. It deals with the 'discretization' of a PDE. Note: ----- - While it converges rapidly to a solution, the objective value achieved is substantially different from that reported in the book. **States** - :math:`y_0, .....y_9, w` : state variables **Specifieds** - :math:`v, q_{00}, q_{11}` : control variables .. GENERATED FROM PYTHON SOURCE LINES 24-31 .. code-block:: Python import numpy as np import sympy as sm import sympy.physics.mechanics as me from opty.direct_collocation import Problem .. GENERATED FROM PYTHON SOURCE LINES 33-35 Equations of Motion. -------------------- .. GENERATED FROM PYTHON SOURCE LINES 35-42 .. code-block:: Python t = me.dynamicsymbols._t q = list(me.dynamicsymbols(f'q:{10}')) uq = [q[i].diff(t) for i in range(10)] w = me.dynamicsymbols('w') v, q00, q11 = me.dynamicsymbols('v q00 q11') .. GENERATED FROM PYTHON SOURCE LINES 43-44 Parameters fom the example. .. GENERATED FROM PYTHON SOURCE LINES 44-49 .. code-block:: Python qa = 0.2 gamma = 0.04 h = 10.0 delta = 1.0 / 9.0 .. GENERATED FROM PYTHON SOURCE LINES 50-51 Equations of motion as per the book. .. GENERATED FROM PYTHON SOURCE LINES 51-61 .. code-block:: Python eom = sm.Matrix([ -uq[0] + 1/delta**2 * (q[1] - 2*q[0] + q00), *[-uq[i] + 1/delta**2 * (q[i+1] - 2*q[i] + q[i-1]) for i in range(1, 9)], -uq[-1] + 1/delta**2 * (q11 - 2*q[-1] + q[-2]), -w.diff(t) + 1/gamma * (v - w), h*(q[0] - w) - 1/(2*delta)*(q[1] - q00), 1/(2*delta)*(q11 - q[-2]), ]) .. GENERATED FROM PYTHON SOURCE LINES 62-64 Optimization ------------ .. GENERATED FROM PYTHON SOURCE LINES 64-98 .. code-block:: Python t0, tf = 0.0, 0.2 num_nodes = 501 interval_value = (tf - t0) / (num_nodes - 1) state_symbols = q + [w] def obj(free): value1 = 1 / (2 * delta) * (qa - free[num_nodes-1])**2 value2 = 1 / (2 * delta) * (qa - free[10 * num_nodes - 1])**2 value3 = 1 / delta * np.sum([(qa - free[(i + 1) * num_nodes - 1])**2 for i in range(1, 9)]) return value1 + value2 + value3 def obj_grad(free): grad = np.zeros_like(free) grad[num_nodes - 1] = -1 / delta * (qa - free[num_nodes - 1]) grad[10 * num_nodes - 1] = -1 / delta * (qa - free[10 * num_nodes - 1]) for i in range(1, 9): grad[(i + 1) * num_nodes - 1] = (-2 / delta * (qa - free[(i + 1) * num_nodes - 1])) return grad # Specify the symbolic instance constraints and the bound, as per the example. instance_constraints = ( *[q[i].func(t0) - 0 for i in range(10)], w.func(t0) - 0, ) bounds = {v: (0.0, 1.0)} .. GENERATED FROM PYTHON SOURCE LINES 99-100 Create the optimization problem. .. GENERATED FROM PYTHON SOURCE LINES 100-113 .. 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 114-115 Give some rough estimates for the trajectories. .. GENERATED FROM PYTHON SOURCE LINES 115-117 .. code-block:: Python initial_guess = np.zeros(prob.num_free) .. GENERATED FROM PYTHON SOURCE LINES 118-119 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 119-126 .. code-block:: Python solution, info = prob.solve(initial_guess) print(info['status_msg']) print(f"Objective value achieved: {info['obj_val']:.4e}, ", f"as per the book it is {2.45476113*1.e-3:.4e} \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: 1.9996e-01, as per the book it is 2.4548e-03 .. GENERATED FROM PYTHON SOURCE LINES 127-128 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 128-130 .. code-block:: Python _ = prob.plot_trajectories(solution, show_bounds=True) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_58_opty_001.png :alt: State Trajectories, Input Trajectories :srcset: /examples/images/sphx_glr_plot_betts_10_58_opty_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 131-132 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 132-134 .. code-block:: Python _ = prob.plot_constraint_violations(solution, subplots=True) .. image-sg:: /examples/images/sphx_glr_plot_betts_10_58_opty_002.png :alt: Constraint violations :srcset: /examples/images/sphx_glr_plot_betts_10_58_opty_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 135-136 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 136-137 .. code-block:: Python _ = prob.plot_objective_value() .. image-sg:: /examples/images/sphx_glr_plot_betts_10_58_opty_003.png :alt: Objective Value :srcset: /examples/images/sphx_glr_plot_betts_10_58_opty_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.376 seconds) .. _sphx_glr_download_examples_plot_betts_10_58_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_58_opty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_58_opty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_58_opty.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_