4.2.1. Two-phase Schwartz OCP#
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: twophaseschwartz.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.2.1.1. Defining OCP#
OCP: https://tomopt.com/docs/propt/tomlab_propt123.php
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=2)
[3]:
# Step-1 : Define dynamics
def dynamics0(x, u, t):
return [x[1], u[0] - 0.1 * (1.0 + 2.0 * x[0] * x[0]) * x[1]]
ocp.dynamics = [dynamics0, dynamics0]
[4]:
# Step-2: Add path constraints
def path_constraints0(x, u, t):
return [
1.0 - 9.0 * (x[0] - 1) * (x[0] - 1) - (x[1] - 0.4) * (x[1] - 0.4) / (0.3 * 0.3)
]
ocp.path_constraints[0] = path_constraints0
[5]:
# Step-3: Add terminal cost
def terminal_cost1(xf, tf, x0, t0):
return 5 * (xf[0] * xf[0] + xf[1] * xf[1])
ocp.terminal_costs[1] = terminal_cost1
Initial state and Final guess
[6]:
ocp.x00[0] = [1, 1]
ocp.x00[1] = [1, 1]
ocp.xf0[0] = [1, 1]
ocp.xf0[1] = [0, 0]
Box constraints
[7]:
ocp.lbx[0][1] = -0.8
ocp.lbu[0], ocp.ubu[0] = -1, 1
ocp.lbt0[0], ocp.ubt0[0] = 0, 0
ocp.lbtf[0], ocp.ubtf[0] = 1, 1
ocp.lbtf[1], ocp.ubtf[1] = 2.9, 2.9
[8]:
ocp.validate()
4.2.1.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.
[9]:
# ocp.du_continuity[0] = 1
mpo, post = mp.solve(ocp, n_segments=1, poly_orders=20, 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............................: 125
variables with only lower bounds: 21
variables with lower and upper bounds: 21
variables with only upper bounds: 0
Total number of equality constraints.................: 88
Total number of inequality constraints...............: 41
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 20
inequality constraints with only upper bounds: 21
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 6.4094433643382967e-22 6.4094433643382967e-22
Dual infeasibility......: 4.9116413862415506e-11 4.9116413862415506e-11
Constraint violation....: 1.5082379789532752e-12 1.5181189638724391e-12
Complementarity.........: 2.5210268504311982e-09 2.5210268504311982e-09
Overall NLP error.......: 2.5210268504311982e-09 2.5210268504311982e-09
Number of objective function evaluations = 11
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 11
Number of inequality constraint evaluations = 11
Number of equality constraint Jacobian evaluations = 10
Number of inequality constraint Jacobian evaluations = 10
Number of Lagrangian Hessian evaluations = 9
Total CPU secs in IPOPT (w/o function evaluations) = 0.022
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 35.00us ( 3.18us) 35.40us ( 3.22us) 11
nlp_g | 335.00us ( 30.45us) 333.19us ( 30.29us) 11
nlp_grad | 57.00us ( 57.00us) 56.50us ( 56.50us) 1
nlp_grad_f | 47.00us ( 4.27us) 44.33us ( 4.03us) 11
nlp_hess_l | 122.00us ( 13.56us) 123.24us ( 13.69us) 9
nlp_jac_g | 523.00us ( 47.55us) 511.58us ( 46.51us) 11
total | 29.37ms ( 29.37ms) 28.73ms ( 28.73ms) 1
Retrive the solution
x: states, u: Controls, t:time, a:Algebraic variables
[10]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 2.9 (Exact), {x[-1]}")
Terminal time, state : 2.9000 vs 2.9 (Exact), [ 9.53605506e-12 -6.10348435e-12]
4.2.1.3. Solve again with Chebyshev-Gauss-Lobatto (CGL) roots#
[11]:
mpo, post = mp.solve(ocp, n_segments=1, poly_orders=20, scheme="CGL", plot=True)
Total number of variables............................: 125
variables with only lower bounds: 21
variables with lower and upper bounds: 21
variables with only upper bounds: 0
Total number of equality constraints.................: 88
Total number of inequality constraints...............: 41
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 20
inequality constraints with only upper bounds: 21
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 6.7496749138057129e-22 6.7496749138057129e-22
Dual infeasibility......: 5.0501516418909649e-11 5.0501516418909649e-11
Constraint violation....: 1.5804024755539103e-12 1.5804024755539103e-12
Complementarity.........: 2.5216728438899515e-09 2.5216728438899515e-09
Overall NLP error.......: 2.5216728438899515e-09 2.5216728438899515e-09
Number of objective function evaluations = 11
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 11
Number of inequality constraint evaluations = 11
Number of equality constraint Jacobian evaluations = 10
Number of inequality constraint Jacobian evaluations = 10
Number of Lagrangian Hessian evaluations = 9
Total CPU secs in IPOPT (w/o function evaluations) = 0.024
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 37.00us ( 3.36us) 35.08us ( 3.19us) 11
nlp_g | 338.00us ( 30.73us) 335.83us ( 30.53us) 11
nlp_grad | 58.00us ( 58.00us) 56.75us ( 56.75us) 1
nlp_grad_f | 50.00us ( 4.55us) 47.37us ( 4.31us) 11
nlp_hess_l | 127.00us ( 14.11us) 125.98us ( 14.00us) 9
nlp_jac_g | 515.00us ( 46.82us) 519.92us ( 47.27us) 11
total | 27.65ms ( 27.65ms) 27.21ms ( 27.21ms) 1
[12]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 2.9s (Exact), {x[-1]}")
Terminal time, state : 2.9000 vs 2.9s (Exact), [ 9.78414149e-12 -6.26610513e-12]
4.2.1.4. Solve again with Legendre-Gauss-Lobatto (LGL) roots#
[13]:
mpo, post = mp.solve(ocp, n_segments=1, poly_orders=20, scheme="LGL", plot=True)
Total number of variables............................: 125
variables with only lower bounds: 21
variables with lower and upper bounds: 21
variables with only upper bounds: 0
Total number of equality constraints.................: 88
Total number of inequality constraints...............: 41
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 20
inequality constraints with only upper bounds: 21
Number of Iterations....: 9
(scaled) (unscaled)
Objective...............: 6.7267963401089627e-22 6.7267963401089627e-22
Dual infeasibility......: 4.9171301510338480e-11 4.9171301510338480e-11
Constraint violation....: 1.5127898933542383e-12 1.5127898933542383e-12
Complementarity.........: 2.5210877471234159e-09 2.5210877471234159e-09
Overall NLP error.......: 2.5210877471234159e-09 2.5210877471234159e-09
Number of objective function evaluations = 11
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 11
Number of inequality constraint evaluations = 11
Number of equality constraint Jacobian evaluations = 10
Number of inequality constraint Jacobian evaluations = 10
Number of Lagrangian Hessian evaluations = 9
Total CPU secs in IPOPT (w/o function evaluations) = 0.021
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 35.00us ( 3.18us) 33.39us ( 3.04us) 11
nlp_g | 334.00us ( 30.36us) 332.27us ( 30.21us) 11
nlp_grad | 57.00us ( 57.00us) 56.27us ( 56.27us) 1
nlp_grad_f | 50.00us ( 4.55us) 46.60us ( 4.24us) 11
nlp_hess_l | 123.00us ( 13.67us) 123.49us ( 13.72us) 9
nlp_jac_g | 526.00us ( 47.82us) 532.02us ( 48.37us) 11
total | 27.35ms ( 27.35ms) 26.99ms ( 26.99ms) 1
[14]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 2.9s (Exact), {x[-1]}")
Terminal time, state : 2.9000 vs 2.9s (Exact), [ 9.76831449e-12 -6.25427524e-12]