CasADi demo

What is CasADi?

A tool for quick & efficient implementation of algorithms for dynamic optimization
Open source, LGPL-licensed, casadi.org
C++ / C++11
Interfaces to Python, Haskell, (Matlab?)
Numerical backends: IPOPT, Sundials, ...
Developers in group of Moritz Diehl:
- Joel Andersson
- Joris Gillis
- Greg Horn

Outline of demo

Scalar expression (SX) graphs
Functions of SX graphs
Matrices of scalar expressions
Automatic differentiation (AD)
Integrators
Matrix expression (MX) graphs
Functions of MX graphs
Solving an optimal control problem

Scalar expression (SX) graphs



In [ ]:

    
from __future__ import print_function
from pylab import *
from casadi import *
from casadi.tools import *  # for dotdraw
from IPython.display import Image, display
%matplotlib inline

def view_dot(graph):
   plt = Image(graph.create_png())
   display(plt)



In [ ]:

    
x = SX.sym("x")  # scalar symbolic primitives
y = SX.sym("y")

z = x*sin(x+y)   # common mathematical operators



In [ ]:

    
print(z)



In [ ]:

    
view_dot(dotgraph(z,direction="BT"))



In [ ]:

    
J = jacobian(z,x)
print(J)



In [ ]:

    
view_dot(dotgraph(J,direction="BT"))

Note 1: subexpressions are shared.

Graph $\leftrightarrow$ Tree
Different from Maple, Matlab symbolic, sympy, ...

A (very) little bit of Computer Algebra



In [ ]:

    
print(x*y/x-y)



In [ ]:

    
H = hessian(z,x)
print(H)

Functions of SX graphs

Sort graph into algorithm



In [ ]:

    
f = Function("f",[x,y],[z])

f.disp(True)

Note 2: re-use of tape variables: live-variables



In [ ]:

    
print(f(1.2,3.4))



In [ ]:

    
print(f(1.2,x+y))



In [ ]:

    
f.generate("f")
print(open("f.c").read())

Matrices of scalar expressions



In [ ]:

    
A = SX.sym("A",3,3)
B = SX.sym("B",3)
print(A)



In [ ]:

    
print(solve(A,B))



In [ ]:

    
print(trace(A))   # Trace



In [ ]:

    
print(mtimes(A,B))   # Matrix multiplication
#print(A @ B)        # Matrix multiplication in Python3



In [ ]:

    
print(norm_fro(A))  # Frobenius norm



In [ ]:

    
print(A[2,:])     # Slicing

Rule 1: Everything is a matrix



In [ ]:

    
print(A.shape, z.shape)



In [ ]:

    
I = SX.eye(3)
print(I)



In [ ]:

    
Ak = kron(I,A)
print(Ak)

Rule 1: Everything is a sparse matrix



In [ ]:

    
Ak.sparsity().spy()



In [ ]:

    
A.sparsity().spy()



In [ ]:

    
z.sparsity().spy()

Automatic differentiation (AD)

Consider an ode:

\begin{equation} \dot{p} = (1 - q^2)p-q+u \end{equation}\begin{equation} \dot{q} = p \end{equation}\begin{equation} \dot{c} = p^2+q^2+u^2 \end{equation}



In [ ]:

    
t = SX.sym("t")    # time
u = SX.sym("u")    # control
p = SX.sym("p")
q = SX.sym("q")
c = SX.sym("c")
x = vertcat(p,q,c) # state



In [ ]:

    
ode = vertcat((1 - q**2)*p - q + u, p, p**2+q**2+u**2)
print(ode, ode.shape)



In [ ]:

    
J = jacobian(ode,x)
print(J)



In [ ]:

    
f = Function("f",[t,u,x],[ode])

ffwd = f.forward(1)

fadj = f.reverse(1)

# side-by-side printing
print('{:*^24} || {:*^28} || {:*^28}'.format("f","ffwd","fadj"))
def short(f):
    import re
    return re.sub(r", a\.k\.a\. \"(\w+)\"",r". \1",f.str(True).replace(", No description available","").replace("Input ","").replace("Output ",""))
for l in zip(short(f).split("\n"),short(ffwd).split("\n"),short(fadj).split("\n")):
  print('{:<24} || {:<28} || {:<28}'.format(*l))

Performing forward sweeps gives the columns of J



In [ ]:

    
print(I)



In [ ]:

    
for i in range(3):
  print(ffwd(t,u,x, ode, 0,0,I[:,i]))



In [ ]:

    
print(J)

Performing adjoint sweeps gives the rows of J



In [ ]:

    
for i in range(3):
  print(fadj(t,u,x, ode, I[:,i])[2])

Often, you can do better than slicing with unit vectors

Note 3: CasADi does graph coloring for efficient sparse jacobians

Integrators

$\dot{x}=f(x,u,t)$ with $x = [p,q,c]^T$



In [ ]:

    
f = {'x':x,'t':t,'p':u,'ode':ode}

Construct an integrating block $x_{k+1} = \Phi(f;\Delta t;x_k,u_k)$



In [ ]:

    
tf = 10.0
N = 20
dt = tf/N



In [ ]:

    
Phi = integrator("Phi","cvodes",f,{"tf":dt})

x0 = DM([0,1,0])

print(Phi(x0=x0))



In [ ]:

    
x = x0
xs = [x]

for i in range(N):
  x = Phi(x0=x)["xf"]
  
  xs.append(x)



In [ ]:

    
plot(horzcat(*xs).T)
legend(["p","q","c"])

Rule 2: Everything is a Function (see http://docs.casadi.org)

Matrix expression (MX) graphs

Note 4: this is what makes CasADi stand out among AD tools

Recall



In [ ]:

    
n = 3

A = SX.sym("A",n,n)
B = SX.sym("B",n,n)
C = mtimes(A,B)
print(C)
view_dot(dotgraph(C,direction='BT'))

What if you don't want to expand into scalar operations? ( avoid $O(n^3)$ storage)



In [ ]:

    
A = MX.sym("A",n,n)
B = MX.sym("B",n,n)
C = mtimes(A,B)
print(C)
view_dot(dotgraph(C,direction='BT'))

What if you cannot expand into matrix operations? ( numerical algorithm )



In [ ]:

    
C = solve(A,B)
print(C)
view_dot(dotgraph(C,direction='BT'))



In [ ]:

    
X0 = MX.sym("x",3)

XF = Phi(x0=X0)["xf"]
print(XF)



In [ ]:

    
expr = sin(XF)+X0
view_dot(dotgraph(expr,direction='BT'))

Functions of MX graphs



In [ ]:

    
F = Function("F",[X0],[  expr  ])
print(F)



In [ ]:

    
print(F(x0))



In [ ]:

    
J = Function("J",[X0],[  jacobian(expr,X0)  ])

print(J(x0))

This shows how an integrator-call can be embedded in matrix graph.

More possibilities: external compiled library, a call to Matlab/Scipy

Solving an optimal control problem

\begin{equation} \begin{array}{cl} \underset{p(.),q(.),u(.)}{\text{minimize}} & \displaystyle \int_{0}^{T}{ p(t)^2 + q(t)^2 + u(t)^2 dt} \\\\ \text{subject to} & \dot{p} = (1 - q^2)p-q+u \\\\ & \dot{q} = p \\\\ & p(0) = 0, q(0) = 1 \\\\ &-1 \le u(t) \le 1 \end{array} \end{equation}

Remember, $\dot{x}=f(x,u,t)$ with $x = [p,q,c]^T$

\begin{equation} \begin{array}{cl} \underset{x(.),u(.)}{\text{minimize}} & c(T) \\\\ \text{subject to} & \dot{x} = f(x,u) \\\\ & p(0) = 0, q(0) = 1, c(0)= 0 \\\\ &-1 \le u(t) \le 1 \end{array} \end{equation}

Discretization with multiple shooting

\begin{equation} \begin{array}{cl} \underset{x_{\bullet},u_{\bullet}}{\text{minimize}} & c_N \\\\ \text{subject to} & x_{k+1} - \Phi(x_k,u_k) = 0 , \quad \quad k = 0,1,\ldots, (N-1) \\\\ & p_0 = 0, q_0 = 1, c_0 = 0 \\\\ &-1 \le u_k \le 1 , \quad \quad k = 0,1,\ldots, (N-1) \end{array} \end{equation}

Cast as NLP

\begin{equation} \begin{array}{cl} \underset{X}{\text{minimize}} & F(X,P) \\\\ \text{subject to} & \text{lbx} \le X \le \text{ubx} \\\\ & \text{lbg} \le G(X,P) \le \text{ubg} \\\\ \end{array} \end{equation}



In [ ]:

    
X = struct_symMX([
     (
      entry("x", repeat=N+1, struct=struct(["p","q","c"]) ),
      entry("u", repeat=N)
     )
    ])

X is a symbolic matrix primitive, but with fancier indexing



In [ ]:

    
print(X.shape)
print((N+1)*3+N)

Demo: $\Phi(x_0,u_0)$



In [ ]:

    
Xf = Phi( x0=X["x",0],p=X["u",0] )["xf"]

print(Xf)

$ x_{k+1} - \Phi(x_k,u_k) = 0 , \quad \quad k = 0,1,\ldots, (N-1)$



In [ ]:

    
g = [] # List of constraint expressions

for k in range(N):
  Xf = Phi( x0=X["x",k],p=X["u",k] )["xf"]
  g.append( X["x",k+1]-Xf )



In [ ]:

    
obj = X["x",N,"c"] # c_N

nlp = {"x":X, "g": vcat(g), "f": obj}

print(nlp)

Block structure in the constraint Jacobian



In [ ]:

    
jacG = jacobian(nlp["g"],nlp["x"])

S = jacG.sparsity()

print(S.shape)

DM.ones(S)[:20,:20].sparsity().spy()

Recall



In [ ]:

    
solver = nlpsol("solver","ipopt",nlp)

lbx = X(-inf)
ubx = X(inf)

lbx["u",:] = -1; ubx["u",:] = 1   #   -1 <= u(t) <= 1

lbx["x",0] = ubx["x",0] = x0      # Initial condition



In [ ]:

    
sol_out = solver(
           lbg = 0, # Equality constraints for shooting constraints
           ubg = 0, #    0 <= g <= 0
           lbx = lbx,
           ubx = ubx)



In [ ]:

    
print(sol_out["x"])



In [ ]:

    
sol = X(sol_out["x"])



In [ ]:

    
plot(horzcat(*sol["x",:]).T)



In [ ]:

    
step(range(N),sol["u",:])

Wrapping up

Showcase: kite-power optimization by Greg Horn, using CasADi backend



In [ ]:

    
from IPython.display import YouTubeVideo
YouTubeVideo('tmjIBpb43j0')



In [ ]:

    
YouTubeVideo('SW6ZJzcMWAk')

Distinction with other software:

ACADOtoolkit	CasADi
Black-box solver Standard-form OCP Good at small-scale real-time NMPC Easy to get started	Write your own solver using a pool of building-blocks No limitations on formulation Good at large-scale OCP Easy to extend

Other operator-overloading AD tools	CasADi
Scalar graphs only, checkpointing	Scalar and matrix graphs Batteries included: Ipopt, Sundials

Closest similarity: AMPL