4.1.2. Moon lander (2D) Example#
2023 Devakumar Thammisetty
MPOPT is an open-source Multi-phase Optimal Control Problem (OCP) solver based on pseudo-spectral collocation with customized adaptive grid refinement techniques.
Download this notebook: moon_lander.ipynb
Install mpopt from pypi using the following. Disable after first usage
Import mpopt (Contains main solver modules)
[1]:
#!pip install mpopt
from mpopt import mp
4.1.2.1. Defining OCP#
The Fuel optimal solution to the moon-lander OCP is known to have bang-bang thrust profile. The selected OCP has one discontinuity.
We first create an OCP object and then polulate the object with dynamics, path_constraints, terminal_constraints and objective (running_costs, terminal_costs)
[2]:
ocp = mp.OCP(n_states=2, n_controls=1, n_phases=1)
[3]:
ocp.dynamics[0] = lambda x, u, t: [x[1], u[0] - 1.5]
[4]:
ocp.running_costs[0] = lambda x, u, t: u[0]
[5]:
ocp.terminal_constraints[0] = lambda xf, tf, x0, t0: [xf[0], xf[1]]
Initial state
[6]:
ocp.x00[0] = [10.0, -2.0]
Box constraints
[7]:
ocp.lbx[0][0] = 0.0
ocp.lbu[0], ocp.ubu[0] = 0, 3
[8]:
ocp.validate()
4.1.2.2. Solve and plot the results in one line#
Lets solve the OCP using following pseudo-spectral approximation * Collocation using Legendre-Gauss-Radau roots * Let’s plot the position and velocity evolution with time starting from 0.
The OCP is a free final time formulation,
[9]:
mpo, post = mp.solve(ocp, n_segments=10, poly_orders=6, scheme="LGR", plot=True)
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
Total number of variables............................: 182
variables with only lower bounds: 61
variables with lower and upper bounds: 61
variables with only upper bounds: 0
Total number of equality constraints.................: 124
Total number of inequality constraints...............: 60
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 60
inequality constraints with only upper bounds: 0
Number of Iterations....: 32
(scaled) (unscaled)
Objective...............: 8.2477255075783038e+00 8.2477255075783038e+00
Dual infeasibility......: 1.9684090468707893e-10 1.9684090468707893e-10
Constraint violation....: 4.0178316229599886e-10 4.0178316229599886e-10
Complementarity.........: 4.4108642187588750e-09 4.4108642187588750e-09
Overall NLP error.......: 4.4108642187588750e-09 4.4108642187588750e-09
Number of objective function evaluations = 33
Number of objective gradient evaluations = 33
Number of equality constraint evaluations = 33
Number of inequality constraint evaluations = 33
Number of equality constraint Jacobian evaluations = 33
Number of inequality constraint Jacobian evaluations = 33
Number of Lagrangian Hessian evaluations = 32
Total CPU secs in IPOPT (w/o function evaluations) = 0.046
Total CPU secs in NLP function evaluations = 0.004
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 145.00us ( 4.39us) 144.99us ( 4.39us) 33
nlp_g | 781.00us ( 23.67us) 774.17us ( 23.46us) 33
nlp_grad | 47.00us ( 47.00us) 44.91us ( 44.91us) 1
nlp_grad_f | 213.00us ( 6.26us) 212.28us ( 6.24us) 34
nlp_hess_l | 274.00us ( 8.56us) 269.97us ( 8.44us) 32
nlp_jac_g | 1.03ms ( 30.24us) 1.04ms ( 30.44us) 34
total | 52.43ms ( 52.43ms) 51.88ms ( 51.88ms) 1
Lets retrive the solution to see the terminal time.
x: states, u: Controls, t:time, a:Algebraic variables in case OCP has differential algebraic equations (DAEs)
Last element of t and x gives the terminal values. Exact terminal time from the analytical solution is 4.1641s.
[10]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 4.1641s (Exact), {x[-1]}")
Terminal time, state : 4.1652 vs 4.1641s (Exact), [9.85970078e-37 0.00000000e+00]
4.1.2.3. Solve again with Chebyshev-Gauss-Lobatto (CGL) roots#
[11]:
mpo, post = mp.solve(ocp, n_segments=2, poly_orders=30, scheme="CGL", plot=True)
Total number of variables............................: 182
variables with only lower bounds: 61
variables with lower and upper bounds: 61
variables with only upper bounds: 0
Total number of equality constraints.................: 124
Total number of inequality constraints...............: 60
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 60
inequality constraints with only upper bounds: 0
Number of Iterations....: 41
(scaled) (unscaled)
Objective...............: 8.2457172048588543e+00 8.2457172048588543e+00
Dual infeasibility......: 1.3972013274602247e-12 1.3972013274602247e-12
Constraint violation....: 3.4862335240859466e-11 3.4862335240859466e-11
Complementarity.........: 3.9978925923296842e-09 3.9978925923296842e-09
Overall NLP error.......: 3.9978925923296842e-09 3.9978925923296842e-09
Number of objective function evaluations = 42
Number of objective gradient evaluations = 42
Number of equality constraint evaluations = 42
Number of inequality constraint evaluations = 42
Number of equality constraint Jacobian evaluations = 42
Number of inequality constraint Jacobian evaluations = 42
Number of Lagrangian Hessian evaluations = 41
Total CPU secs in IPOPT (w/o function evaluations) = 0.150
Total CPU secs in NLP function evaluations = 0.009
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 188.00us ( 4.48us) 183.99us ( 4.38us) 42
nlp_g | 2.48ms ( 58.98us) 2.49ms ( 59.19us) 42
nlp_grad | 114.00us (114.00us) 113.29us (113.29us) 1
nlp_grad_f | 273.00us ( 6.35us) 266.81us ( 6.20us) 43
nlp_hess_l | 336.00us ( 8.20us) 332.45us ( 8.11us) 41
nlp_jac_g | 3.81ms ( 88.63us) 3.80ms ( 88.47us) 43
total | 167.97ms (167.97ms) 167.12ms (167.12ms) 1
[12]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 4.1641s (Exact), {x[-1]}")
Terminal time, state : 4.1661 vs 4.1641s (Exact), [0. 0.]
4.1.2.4. Solve again with Legendre-Gauss-Lobatto (LGL) roots#
[13]:
mpo, post = mp.solve(ocp, n_segments=2, poly_orders=30, scheme="LGL", plot=True)
Total number of variables............................: 182
variables with only lower bounds: 61
variables with lower and upper bounds: 61
variables with only upper bounds: 0
Total number of equality constraints.................: 124
Total number of inequality constraints...............: 60
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 60
inequality constraints with only upper bounds: 0
Number of Iterations....: 42
(scaled) (unscaled)
Objective...............: 8.2425586640613506e+00 8.2425586640613506e+00
Dual infeasibility......: 6.8055248936271795e-11 6.8055248936271795e-11
Constraint violation....: 6.9194650009762881e-12 8.8133944586843427e-12
Complementarity.........: 4.7886457777884481e-09 4.7886457777884481e-09
Overall NLP error.......: 4.7886457777884481e-09 4.7886457777884481e-09
Number of objective function evaluations = 43
Number of objective gradient evaluations = 43
Number of equality constraint evaluations = 43
Number of inequality constraint evaluations = 43
Number of equality constraint Jacobian evaluations = 43
Number of inequality constraint Jacobian evaluations = 43
Number of Lagrangian Hessian evaluations = 42
Total CPU secs in IPOPT (w/o function evaluations) = 0.139
Total CPU secs in NLP function evaluations = 0.010
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 184.00us ( 4.28us) 181.26us ( 4.22us) 43
nlp_g | 2.61ms ( 60.79us) 2.62ms ( 60.86us) 43
nlp_grad | 113.00us (113.00us) 112.84us (112.84us) 1
nlp_grad_f | 283.00us ( 6.43us) 276.32us ( 6.28us) 44
nlp_hess_l | 350.00us ( 8.33us) 342.76us ( 8.16us) 42
nlp_jac_g | 3.91ms ( 88.80us) 3.93ms ( 89.37us) 44
total | 149.56ms (149.56ms) 149.21ms (149.21ms) 1
[14]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 4.1641s (Exact), {x[-1]}")
Terminal time, state : 4.1662 vs 4.1641s (Exact), [3.68621737e-35 0.00000000e+00]