.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/intermediate/plot_betts_10_47.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_intermediate_plot_betts_10_47.py: Singular Arc Problem ==================== Objective --------- - Show how to optimize a multiphase problem phasewise. - Show how to use a state variable to enforce a boundary condition. Introduction ------------ This is example 10.47 from [Betts2010]_, Chapter 10: Test Problems. John T. Betts made available his detailed results of this example, and 'error' is the relative error in the value compared to the values given. This is a three stage problem, where the final state of one stage is the initial state of the next stage. At the start of phase 2 there is a boundary condition: :math:`m \cdot g - (1 + \dfrac{v}{c}) \cdot \sigma v^2 e^{-h/h_0} = 0` It is enforced at the end of phase 1, by introducing a new state variable :math:`h_\textrm{end}`, and set :math:`h_{\textrm{end}} = m \cdot g - (1 + \dfrac{v}{c}) \cdot \sigma v^2 e^{-h/h_0}` and set the instance constraint :math:`h_{end} = 0` at the end of phase 1. Notes ----- ``opty`` presently does not support simultaneous optimization of multiple phases, so they have to be optimized phasewise. Dr. Betts confirmed that his solution was achieved by simultaneous optimization. So, at least in this case, phasewise optimization gives pretty close results. Interesting: The last phase, which seems to be the simplest, just free gliding, needs by far the most iterations to converge. Phase 2, the sigular arc, seems to pose no problems at all. **Phase 1 Maximum Thrust** **States** - :math:`h, v, m` : state variables - :math:`h_{\textrm{end}}` : state variable to enforce a boundary condition .. GENERATED FROM PYTHON SOURCE LINES 54-63 .. code-block:: Python import numpy as np import sympy as sm import sympy.physics.mechanics as me from opty import Problem import matplotlib.pyplot as plt import os .. GENERATED FROM PYTHON SOURCE LINES 64-65 Equations of motion. .. GENERATED FROM PYTHON SOURCE LINES 65-85 .. code-block:: Python t = me.dynamicsymbols._t h, v, m, h_end = me.dynamicsymbols('h v m, h_end') h_fast = sm.symbols('h_fast') # Parameters, the same for all three phases Tm = 193.044 g = 32.174 # Imperial units sigma = 5.4915348492 * 1.e-5 c = 1580.942579 h0 = 23800 eom = sm.Matrix([ -h.diff(t) + v, -v.diff(t) + 1/m * (Tm - sigma*v**2*sm.exp(-h/h0)) - g, -m.diff(t) - Tm/c, -h_end + m*g - (1 + v/c) * sigma*v**2*sm.exp(-h/h0), ]) sm.pprint(eom) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ↪ ⎢ v(t) - ──(h(t)) ↪ ⎢ dt ↪ ⎢ ↪ ⎢ -h(t) ↪ ⎢ ────── ↪ ⎢ 2 23800 ↪ ⎢ - 5.4915348492e-5⋅v (t)⋅ℯ + 193.044 d ↪ ⎢ ───────────────────────────────────────── - ──(v(t)) - 32.174 ↪ ⎢ m(t) dt ↪ ⎢ ↪ ⎢ d ↪ ⎢ - ──(m(t)) - 0.122106901644781 ↪ ⎢ dt ↪ ⎢ ↪ ⎢ -h(t) ↪ ⎢ ────── ↪ ⎢ 2 23800 ↪ ⎣- (3.47358273611277e-8⋅v(t) + 5.4915348492e-5)⋅v (t)⋅ℯ - h_end(t) + 32. ↪ ↪ ⎤ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ 174⋅m(t)⎦ .. GENERATED FROM PYTHON SOURCE LINES 86-92 .. code-block:: Python num_nodes = 101 t0, tf = 0*h_fast, (num_nodes-1) * h_fast interval_value = h_fast state_symbols = (h, v, m, h_end) .. GENERATED FROM PYTHON SOURCE LINES 93-96 Specify the objective function and form the gradient. h(:math:`t_f`) is to be maximized. Note that :math:`h_\textrm{fast}` is the last entry of free. .. GENERATED FROM PYTHON SOURCE LINES 96-122 .. code-block:: Python def obj(free): return -free[num_nodes-1] * free[-1] def obj_grad(free): grad = np.zeros_like(free) grad[num_nodes-1] = -free[-1] grad[-1] = -free[num_nodes-1] return grad # Specify the symbolic instance constraints, as per the example instance_constraints = ( h.func(0*h_fast), v.func(0*h_fast), m.func(0*h_fast) - 3.0, h_end.func(tf), ) bounds = { h_fast: (0.0, 0.5), m: (1.0, 3.0) } .. GENERATED FROM PYTHON SOURCE LINES 123-124 Create the optimization problem and set any options. .. GENERATED FROM PYTHON SOURCE LINES 124-141 .. code-block:: Python prob = Problem(obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, time_symbol=t, backend='numpy', ) prob.add_option('max_iter', 1000) # Give some rough estimates for the trajectories. initial_guess = np.ones(prob.num_free) * 0.1 .. GENERATED FROM PYTHON SOURCE LINES 142-143 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 143-148 .. code-block:: Python solution, info = prob.solve(initial_guess) print(info['status_msg']) t_phase1 = solution[-1] solution1 = solution .. 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).' .. GENERATED FROM PYTHON SOURCE LINES 149-150 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 150-152 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_001.png :alt: State Trajectories :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 153-154 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 154-156 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_002.png :alt: Constraint violations :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 157-158 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 158-170 .. code-block:: Python prob.plot_objective_value() error_t = (t_phase1*(num_nodes-1) - 13.726485)/13.726485 * 100 error_v = (solution[2*num_nodes-1] - 791.2744)/791.2744 * 100 error_h = (solution[num_nodes-1] - 4560.8912)/4560.8912 * 100 error_m = (solution[3*num_nodes-1] - 1.323901)/1.323901 * 100 print((f"duration of phase 1 is {(num_nodes-1) * solution[-1]:.3f}, " f"error is {error_t:.3f} %")) print((f"Height achieved is {- info['obj_val']/solution[-1]:.3f}, " f"error is {error_h:.3f} %")) print((f"Velocity achieved is {solution[2*num_nodes-1]:.3f}, " f"error is {error_v:.3f} %")) print(f"Final mass is {solution[3*num_nodes-1]:.3f}, ") .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_003.png :alt: Objective Value :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none duration of phase 1 is 13.729, error is 0.015 % Height achieved is 4557.723, error is -0.069 % Velocity achieved is 791.166, error is -0.014 % Final mass is 1.324, .. GENERATED FROM PYTHON SOURCE LINES 171-179 **Phase 2 Singular Arc** Now the thrust T(t) becomes a state variable. **States** - :math:`h, v, m, T`: state variables .. GENERATED FROM PYTHON SOURCE LINES 179-197 .. code-block:: Python # Set up the eoms T = me.dynamicsymbols('T') eom = sm.Matrix([ -h.diff(t) + v, -v.diff(t) + 1/m * (T - sigma*v**2*sm.exp(-h/h0)) - g, -m.diff(t) - T/c, (T - sigma*v**2*sm.exp(-h/h0) - m*g -m*g / (1 + 4*c/v + 2*c**2/v**2) * (c**2 / (h0 * g) * (1 + v/c) - 1 - 2 * c/v)), ]) sm.pprint(eom) state_symbols = (h, v, m, T) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ↪ ⎢ v(t) - ──(h(t)) ↪ ⎢ dt ↪ ⎢ ↪ ⎢ -h(t) ↪ ⎢ ────── ↪ ⎢ 2 23800 ↪ ⎢ T(t) - 5.4915348492e-5⋅v (t)⋅ℯ d ↪ ⎢ ──────────────────────────────────── - ──(v ↪ ⎢ m(t) dt ↪ ⎢ ↪ ⎢ d ↪ ⎢ -0.000632534042212042⋅T(t) - ──(m( ↪ ⎢ dt ↪ ⎢ ↪ ⎢ -h(t) ⎛ ↪ ⎢ ────── 32.174⋅⎜0.00206459124701 ↪ ⎢ 2 23800 ⎝ ↪ ⎢T(t) - 32.174⋅m(t) - 5.4915348492e-5⋅v (t)⋅ℯ - ──────────────────────── ↪ ⎢ 63 ↪ ⎢ 1 + ── ↪ ⎢ ↪ ⎣ ↪ ↪ ⎤ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ (t)) - 32.174 ⎥ ↪ ⎥ ↪ ⎥ ↪ ⎥ ↪ t)) ⎥ ↪ ⎥ ↪ ⎥ ↪ 3161.885158⎞ ⎥ ↪ 662⋅v(t) + 2.26400021063928 - ───────────⎟⋅m(t)⎥ ↪ v(t) ⎠ ⎥ ↪ ───────────────────────────────────────────────⎥ ↪ 23.770316 4998758.87619034 ⎥ ↪ ───────── + ──────────────── ⎥ ↪ v(t) 2 ⎥ ↪ v (t) ⎦ .. GENERATED FROM PYTHON SOURCE LINES 198-200 The instance constraints (of course) are the values achieved at the end of the first phase. .. GENERATED FROM PYTHON SOURCE LINES 200-216 .. code-block:: Python instance_constraints = ( h.func(0*h_fast) - solution[num_nodes-1], v.func(0*h_fast) - solution[2*num_nodes-1], m.func(0*h_fast) - solution[3*num_nodes-1], T.func(0*h_fast) - Tm, m.func((num_nodes-1)*h_fast) - 1.0, ) # As per the example, the mass m >= 1.0 bounds = { h_fast: (0.0, 0.5), T: (0.0, Tm), m: (1.0, solution[3*num_nodes-1]) } .. GENERATED FROM PYTHON SOURCE LINES 217-218 Create the optimization problem and set any options. .. GENERATED FROM PYTHON SOURCE LINES 218-239 .. code-block:: Python prob = Problem(obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, time_symbol=t, backend='numpy', ) prob.add_option('max_iter', 3000) # Give some rough estimates for the trajectories. # The result of a previous run are used to save time. initial_guess = np.array((*[solution[num_nodes-1] for _ in range(num_nodes)], *[solution[2*num_nodes-1] for _ in range(num_nodes)], *[solution[3*num_nodes-1] for _ in range(num_nodes)], *[Tm/c for _ in range(num_nodes)], solution[-1])) .. GENERATED FROM PYTHON SOURCE LINES 240-241 Find the optimal solution. .. GENERATED FROM PYTHON SOURCE LINES 241-246 .. code-block:: Python solution, info = prob.solve(initial_guess) initial_guess = solution print(info['status_msg']) t_phase2 = solution[-1] solution2 = solution .. 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).' .. GENERATED FROM PYTHON SOURCE LINES 247-248 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 248-250 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_004.png :alt: State Trajectories :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 251-252 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 252-254 .. code-block:: Python _ = prob.plot_constraint_violations(solution) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_005.png :alt: Constraint violations :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 255-256 Plot the objective function as a function of optimizer iteration. .. GENERATED FROM PYTHON SOURCE LINES 256-271 .. code-block:: Python _ = prob.plot_objective_value() error_t = ((t_phase2*(num_nodes-1) - (22.025604-13.726485)) / (22.025604-13.726485) * 100) error_v = (solution[2*num_nodes-1] - 789.64098)/789.64098 * 100 error_h = (solution[num_nodes-1] - 11121.110)/11121.110 * 100 error_m = (solution[3*num_nodes-1] - 1.0)/1.0 * 100 print((f'duration of phase 2 is {(num_nodes-1) * solution[-1]:.3f}, ' f'error is {error_t:.3f} %')) print((f"Height achieved is {- info['obj_val']/solution[-1]:.3f}, " f'error is {error_h:.3f} %')) print((f'Velocity achieved is {solution[2*num_nodes-1]:.3f}, ' f'error is {error_v:.3f} %')) print(f'Final mass is {solution[3*num_nodes-1]:.3f}, ') .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_006.png :alt: Objective Value :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none duration of phase 2 is 8.310, error is 0.126 % Height achieved is 11125.188, error is 0.037 % Velocity achieved is 789.517, error is -0.016 % Final mass is 1.000, .. GENERATED FROM PYTHON SOURCE LINES 272-281 **Phase 3 No Thrust** Now it is free gliding without thrust. **States** - :math:`h, v, m`: state variables. Set up the eoms. .. GENERATED FROM PYTHON SOURCE LINES 281-290 .. code-block:: Python eom = sm.Matrix([ -h.diff(t) + v, -v.diff(t) - sigma*v**2*sm.exp(-h/h0)/m - g, -m.diff(t) - 0, ]) sm.pprint(eom) state_symbols = (h, v, m) .. rst-class:: sphx-glr-script-out .. code-block:: none ⎡ d ⎤ ⎢ v(t) - ──(h(t)) ⎥ ⎢ dt ⎥ ⎢ ⎥ ⎢ -h(t) ⎥ ⎢ ──────⎥ ⎢ 2 23800 ⎥ ⎢ d 5.4915348492e-5⋅v (t)⋅ℯ ⎥ ⎢- ──(v(t)) - 32.174 - ─────────────────────────────⎥ ⎢ dt m(t) ⎥ ⎢ ⎥ ⎢ d ⎥ ⎢ -──(m(t)) ⎥ ⎣ dt ⎦ .. GENERATED FROM PYTHON SOURCE LINES 291-297 The maximum height is achieved when the velocity is zero. The objective function to be minimized is the square of the final speed. There is really nothing here to be optimized, just stopp when the velocity is zero. However, the 'obvious' choice of setting an instance constraint v(:math:`t_f`) = 0.0 creates an error message: *Problem has too few degrees of freedom*. .. GENERATED FROM PYTHON SOURCE LINES 297-325 .. code-block:: Python def obj(free): return free[2*num_nodes-1]**2 def obj_grad(free): grad = np.zeros_like(free) grad[2*num_nodes-1] = 2*free[2*num_nodes-1] grad[-1] = 0.0 return grad # Specify the symbolic instance constraints. The starting values are the final # values of the previous phase. The mass must be 1.0 at the end. instance_constraints = ( h.func(0*h_fast) - solution[num_nodes-1], v.func(0*h_fast) - solution[2*num_nodes-1], m.func(0*h_fast) - solution[3*num_nodes-1], m.func((num_nodes-1)*h_fast) - 1.0, ) bounds = { h_fast: (0.0, 0.5), m: (solution[3*num_nodes-1], 1.0), v: (0.0, np.inf), } .. GENERATED FROM PYTHON SOURCE LINES 326-327 Create the optimization problem and set any options. .. GENERATED FROM PYTHON SOURCE LINES 327-366 .. code-block:: Python prob = Problem(obj, obj_grad, eom, state_symbols, num_nodes, interval_value, instance_constraints=instance_constraints, bounds=bounds, time_symbol=t, backend='numpy', ) prob.add_option('max_iter', 20000) fname = f'betts_10_47_phase3_{num_nodes}_nodes_solution.csv' # Use a given solution if available, else give an initial guess and solve the # problem. if os.path.exists(fname): solution = np.loadtxt(fname) # Use solution given. solution = np.loadtxt(fname) else: # Give some rough estimates for the trajectories. initial_guess = np.array((*[solution[num_nodes-1] for _ in range(num_nodes)], *[solution[2*num_nodes-1] for _ in range(num_nodes)], *[solution[3*num_nodes-1] for _ in range(num_nodes)], solution[-1]) ) # Find the optimal solution. solution, info = prob.solve(initial_guess) print(info['status_msg']) # Plot the objective function as a function of optimizer iteration. _ = prob.plot_objective_value() t_phase3 = solution[-1] solution3 = solution # np.savetxt(fname, solution, fmt='%.12f') .. GENERATED FROM PYTHON SOURCE LINES 367-368 Plot the optimal state and input trajectories. .. GENERATED FROM PYTHON SOURCE LINES 368-370 .. code-block:: Python _ = prob.plot_trajectories(solution) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_007.png :alt: State Trajectories :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 371-372 Plot the constraint violations. .. GENERATED FROM PYTHON SOURCE LINES 372-390 .. code-block:: Python _ = prob.plot_constraint_violations(solution) print(f'final speed is {solution[2*num_nodes-1]:.3e}, ideally it should be 0') print('\n') error_h = (solution[num_nodes-1] - 18664.187)/18664.187 * 100 print((f'final height achieved is is {solution[num_nodes-1]:.3f}, error is ' f' {error_h:.3f} %')) error_t1 = (t_phase1*(num_nodes-1) - 13.751270)/13.751270 * 100 error_t2 = ((t_phase1+t_phase2)*(num_nodes-1) - 21.987)/21.987 * 100 print(('duration of first and second phase is ' f'{(t_phase1+t_phase2)*(num_nodes-1):.3f}, error is ' f'{error_t2:.3f} %')) error_t3 = (((t_phase1+t_phase2+t_phase3)*(num_nodes-1) - 42.641355)/42.641355 * 100) print((f'total duration is {(t_phase1+t_phase2+t_phase3) * (num_nodes-1):.3f}' f' sec, error is {error_t3:.3f} %')) .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_008.png :alt: Constraint violations :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none final speed is 1.143e-02, ideally it should be 0 final height achieved is is 18661.338, error is -0.015 % duration of first and second phase is 22.038, error is 0.232 % total duration is 43.041 sec, error is 0.937 % .. GENERATED FROM PYTHON SOURCE LINES 391-392 Plot complete journey. .. GENERATED FROM PYTHON SOURCE LINES 392-415 .. code-block:: Python h_list = np.concatenate((solution1[:num_nodes], solution2[:num_nodes], solution3[:num_nodes])) v_list = np.concatenate((solution1[num_nodes:2*num_nodes], solution2[num_nodes:2*num_nodes], solution3[num_nodes:2*num_nodes])) m_list = np.concatenate((solution1[2*num_nodes:3*num_nodes], solution2[2*num_nodes:3*num_nodes], solution3[2*num_nodes:3*num_nodes])) times = np.linspace(0, (num_nodes-1)*(t_phase1+t_phase2+t_phase3), num_nodes*3) fig, ax = plt.subplots(3, 1, figsize=(6.8, 4), sharex=True, layout='constrained') ax[0].plot(times, h_list) ax[0].set_ylabel('height [ft]') ax[0].set_title('Trajectory of total journey') ax[1].plot(times, v_list) ax[1].set_ylabel('velocity [ft/s]') ax[1].set_title('Velocity') ax[2].plot(times, m_list) ax[2].set_ylabel('mass [lbs]') ax[2].set_xlabel('time [s]') _ = ax[2].set_title('Mass') .. image-sg:: /examples/intermediate/images/sphx_glr_plot_betts_10_47_009.png :alt: Trajectory of total journey, Velocity, Mass :srcset: /examples/intermediate/images/sphx_glr_plot_betts_10_47_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 416-417 sphinx_gallery_thumbnail_number = 9 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 7.254 seconds) .. _sphx_glr_download_examples_intermediate_plot_betts_10_47.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_47.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_betts_10_47.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_betts_10_47.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_