In [13]:
%matplotlib inline
import numpy as np
import casadi as ca
import matplotlib.pyplot as plt
from operator import itemgetter
In [16]:
nk = 100 # Control discretization
tf = 10.0 # End time
# Declare variables (use scalar graph)
u = ca.SX.sym( 'u' ) # control
x = ca.SX.sym( 'x', 3 ) # states
This notebook will try to solve third order system in form of: $$ \dddot x + \ddot x + \dot x + x - u = 0 $$ Or alternatively in the 1st order form: $$\begin{cases} \dot x_0 = -x_0 - x_1 - x_2 + u \\ \dot x_1 = x_0 \\ \dot x_2 = x_1 \end{cases} $$ Subject to: $$ \min \int_{0}^{T}x_2^2 + u^2 \,dx $$
In [17]:
# ODE right hand side and quadratures
xdot = ca.vertcat( [-x[0] - x[1] - x[2] + u, x[0], x[1]] )
qdot = x[1]*x[1] + x[2]*x[2] + u*u
# DAE residual function
dae = ca.SXFunction("dae", ca.daeIn(x=x, p=u), ca.daeOut(ode=xdot, quad=qdot))
# Create an integrator
integrator = ca.Integrator("integrator", "cvodes", dae, {"tf":tf/nk})
# All controls (use matrix graph)
x = ca.MX.sym("x",nk) # nk-by-1 symbolic variable
U = ca.vertsplit(x) # cheaper than x[0], x[1], ...
# The initial state (x_0=0, x_1=1)
X = ca.MX([0,1, 1])
In [18]:
# Objective function
f = 0
# Build a graph of integrator calls
for k in range(nk):
X,QF = itemgetter('xf','qf')(integrator( {'x0':X, 'p':U[k] } ) )
f += QF
# Terminal constraints: x_0(T)=x_1(T)=0
g = X
# Allocate an NLP solver
opts = {'linear_solver': 'ma27'}
nlp = ca.MXFunction("nlp", ca.nlpIn(x=x), ca.nlpOut(f=f,g=g))
solver = ca.NlpSolver("solver", "ipopt", nlp, opts)
# Solve the problem
sol = solver({"lbx" : -0.75,
"ubx" : 1,
"x0" : 0,
"lbg" : 0,
"ubg" : 0})
In [21]:
# Retrieve the solution
u_opt = np.array(sol["x"])
# Time grid
tgrid_x = np.linspace(0,10,nk+1)
tgrid_u = np.linspace(0,10,nk)
# Plot the results
plt.figure(1)
plt.clf()
plt.plot(tgrid_u,u_opt,'b-')
plt.title("third order system optimisation - single shooting")
plt.xlabel('time')
plt.legend(['u trajectory'])
plt.grid()
In [20]:
sol
Out[20]:
In [ ]: