4.1.1. Van der pol Oscillator#
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: vanderpol.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.1.1. OCP definition#
Van der Pol OCP: https://web.casadi.org/docs/#a-simple-test-problem
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)
[3]:
def dynamics(x, u, t):
return [(1 - x[1] * x[1]) * x[0] - x[1] + u[0], x[0]]
[4]:
def running_cost(x, u, t):
return x[0] * x[0] + x[1] * x[1] + u[0] * u[0]
[5]:
ocp.dynamics[0] = dynamics
ocp.running_costs[0] = running_cost
Initial state
[6]:
ocp.x00[0] = [0, 1]
Box constraints
[7]:
ocp.lbu[0] = -1.0
ocp.ubu[0] = 1.0
ocp.lbx[0][1] = -0.25
ocp.lbtf[0] = 10.0
ocp.ubtf[0] = 10.0
[8]:
ocp.validate()
4.1.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]:
mpo, post = mp.solve(ocp, n_segments=1, poly_orders=25, 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............................: 76
variables with only lower bounds: 25
variables with lower and upper bounds: 26
variables with only upper bounds: 0
Total number of equality constraints.................: 52
Total number of inequality constraints...............: 25
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 25
inequality constraints with only upper bounds: 0
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 2.8734932991287625e+00 2.8734932991287625e+00
Dual infeasibility......: 3.4605651677566129e-13 3.4605651677566129e-13
Constraint violation....: 4.6562753652779065e-13 4.6562753652779065e-13
Complementarity.........: 2.5153452656320664e-09 2.5153452656320664e-09
Overall NLP error.......: 2.5153452656320664e-09 2.5153452656320664e-09
Number of objective function evaluations = 12
Number of objective gradient evaluations = 12
Number of equality constraint evaluations = 12
Number of inequality constraint evaluations = 12
Number of equality constraint Jacobian evaluations = 12
Number of inequality constraint Jacobian evaluations = 12
Number of Lagrangian Hessian evaluations = 11
Total CPU secs in IPOPT (w/o function evaluations) = 0.016
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 53.00us ( 4.42us) 52.65us ( 4.39us) 12
nlp_g | 327.00us ( 27.25us) 323.85us ( 26.99us) 12
nlp_grad | 58.00us ( 58.00us) 55.45us ( 55.45us) 1
nlp_grad_f | 77.00us ( 5.92us) 77.00us ( 5.92us) 13
nlp_hess_l | 112.00us ( 10.18us) 115.53us ( 10.50us) 11
nlp_jac_g | 515.00us ( 39.62us) 515.59us ( 39.66us) 13
total | 19.01ms ( 19.01ms) 19.51ms ( 19.51ms) 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 10 (Exact), {x[-1]}")
Terminal time, state : 10.0000 vs 10 (Exact), [-0.00178384 -0.00014044]
4.1.1.3. Solve again with Chebyshev-Gauss-Lobatto (CGL) roots#
[11]:
mpo, post = mp.solve(ocp, n_segments=1, poly_orders=25, scheme="CGL", plot=True)
Total number of variables............................: 76
variables with only lower bounds: 25
variables with lower and upper bounds: 26
variables with only upper bounds: 0
Total number of equality constraints.................: 52
Total number of inequality constraints...............: 25
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 25
inequality constraints with only upper bounds: 0
Number of Iterations....: 12
(scaled) (unscaled)
Objective...............: 2.8734060736207896e+00 2.8734060736207896e+00
Dual infeasibility......: 2.6383995374325467e-12 2.6383995374325467e-12
Constraint violation....: 2.1178614417749486e-12 2.1178614417749486e-12
Complementarity.........: 2.5713635670826573e-09 2.5713635670826573e-09
Overall NLP error.......: 2.5713635670826573e-09 2.5713635670826573e-09
Number of objective function evaluations = 13
Number of objective gradient evaluations = 13
Number of equality constraint evaluations = 13
Number of inequality constraint evaluations = 13
Number of equality constraint Jacobian evaluations = 13
Number of inequality constraint Jacobian evaluations = 13
Number of Lagrangian Hessian evaluations = 12
Total CPU secs in IPOPT (w/o function evaluations) = 0.013
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 58.00us ( 4.46us) 57.81us ( 4.45us) 13
nlp_g | 358.00us ( 27.54us) 355.20us ( 27.32us) 13
nlp_grad | 68.00us ( 68.00us) 67.78us ( 67.78us) 1
nlp_grad_f | 86.00us ( 6.14us) 83.08us ( 5.93us) 14
nlp_hess_l | 127.00us ( 10.58us) 126.75us ( 10.56us) 12
nlp_jac_g | 536.00us ( 38.29us) 610.09us ( 43.58us) 14
total | 20.51ms ( 20.51ms) 20.98ms ( 20.98ms) 1
[12]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 10s (Exact), {x[-1]}")
Terminal time, state : 10.0000 vs 10s (Exact), [-0.00178269 -0.00013973]
4.1.1.3.1. Solve again with Legendre-Gauss-Lobatto (LGL) roots#
[13]:
mpo, post = mp.solve(ocp, n_segments=1, poly_orders=25, scheme="LGL", plot=True)
Total number of variables............................: 76
variables with only lower bounds: 25
variables with lower and upper bounds: 26
variables with only upper bounds: 0
Total number of equality constraints.................: 52
Total number of inequality constraints...............: 25
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 25
inequality constraints with only upper bounds: 0
Number of Iterations....: 11
(scaled) (unscaled)
Objective...............: 2.8734849959084205e+00 2.8734849959084205e+00
Dual infeasibility......: 8.2869030765986375e-12 8.2869030765986375e-12
Constraint violation....: 6.0937921375625592e-12 6.0937921375625592e-12
Complementarity.........: 2.6832881949699023e-09 2.6832881949699023e-09
Overall NLP error.......: 2.6832881949699023e-09 2.6832881949699023e-09
Number of objective function evaluations = 12
Number of objective gradient evaluations = 12
Number of equality constraint evaluations = 12
Number of inequality constraint evaluations = 12
Number of equality constraint Jacobian evaluations = 12
Number of inequality constraint Jacobian evaluations = 12
Number of Lagrangian Hessian evaluations = 11
Total CPU secs in IPOPT (w/o function evaluations) = 0.015
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
solver : t_proc (avg) t_wall (avg) n_eval
nlp_f | 57.00us ( 4.75us) 54.18us ( 4.51us) 12
nlp_g | 330.00us ( 27.50us) 326.96us ( 27.25us) 12
nlp_grad | 53.00us ( 53.00us) 52.09us ( 52.09us) 1
nlp_grad_f | 77.00us ( 5.92us) 76.10us ( 5.85us) 13
nlp_hess_l | 119.00us ( 10.82us) 116.95us ( 10.63us) 11
nlp_jac_g | 510.00us ( 39.23us) 516.94us ( 39.76us) 13
total | 18.17ms ( 18.17ms) 17.80ms ( 17.80ms) 1
[14]:
x, u, t, a = post.get_data()
print(f"Terminal time, state : {t[-1][0]:.4f} vs 10s (Exact), {x[-1]}")
Terminal time, state : 10.0000 vs 10s (Exact), [-0.00178275 -0.00014031]
[ ]: