In [69]:
from casadi import *
from casadi.tools import * # for dotdraw
In [70]:
x = SX.sym("x") # scalar symbolic primitives
y = SX.sym("y")
z = x*sin(x+y) # common mathematical operators
In [71]:
print z
In [72]:
dotdraw(z,direction="BT")
In [73]:
J = jacobian(z,x)
print J
In [74]:
dotdraw(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 [13]:
print x*y/x-y
In [9]:
H = hessian(z,x)
print H
Sort graph into algorithm
In [14]:
f = SXFunction([x,y],[z])
f.init()
print f
Note 2: re-use of tape variables: live-variables
In [15]:
f.setInput(1.2,0)
f.setInput(3.4,1)
f.evaluate()
print f.getOutput(0)
In [16]:
print f([1.2,3.4])
In [18]:
print f([1.2,x+y])
In [19]:
f.generateCode("f.c")
print file("f.c").read()
In [20]:
A = SX.sym("A",3,3)
B = SX.sym("B",3)
print A
In [21]:
print solve(A,B)
In [22]:
print trace(A) # Trace
In [23]:
print mul(A,B) # Matrix multiplication
In [24]:
print norm_F(A) # Frobenius norm
In [25]:
print A[2,:] # Slicing
Rule 1: Everything is a matrix
In [26]:
print A.shape, z.shape
In [27]:
I = SX.eye(3)
print I
In [28]:
Ak = kron(I,A)
print Ak
Rule 1: Everything is a sparse matrix
In [29]:
with nice_stdout():
Ak.sparsity().spy()
In [30]:
with nice_stdout():
A.sparsity().spy()
In [31]:
with nice_stdout():
z.sparsity().spy()
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 [32]:
t = SX.sym("t") # time
u = SX.sym("u") # control
x = SX.sym("[p,q,c]") # state
p,q,c = x
In [34]:
ode = vertcat([(1 - q**2)*p - q + u, p, p**2+q**2+u**2])
print ode, ode.shape
In [35]:
J = jacobian(ode,x)
print J
In [36]:
f = SXFunction([t,u,x],[ode])
f.init()
ffwd = f.derivative(1,0)
ffwd.init()
fadj = f.derivative(0,1)
fadj.init()
# side-by-side printing
print '{:*^24} || {:*^30} || {:*^30}'.format("f","ffwd","fadj")
for l in zip(str(f).split("\n"),str(ffwd).split("\n"),str(fadj).split("\n")):
print '{:<24} || {:<30} || {:<30}'.format(*l)
Performing forward sweeps gives the columns of J
In [37]:
print I
In [38]:
for i in range(3):
print ffwd([ t,u,x, 0,0,I[:,i] ]) [1]
In [39]:
print J
Performing adjoint sweeps gives the rows of J
In [40]:
for i in range(3):
print fadj([ t,u,x, I[:,i] ]) [3]
Often, you can do better than slicing with unit vectors
Note 3: CasADi does graph coloring for efficient sparse jacobians
$\dot{x}=f(x,u,t)$ with $x = [p,q,c]^T$
In [41]:
f = SXFunction(daeIn(x=x,t=t,p=u),daeOut(ode=ode))
f.init()
print str(f)[:327]
Construct an integrating block $x_{k+1} = \Phi(f;\Delta t;x_k,u_k)$
In [43]:
tf = 10.0
N = 20
dt = tf/N
In [44]:
Phi = Integrator("cvodes",f)
Phi.setOption("name","Phi")
Phi.setOption("tf",dt)
Phi.init()
x0 = vertcat([0,1,0])
Phi.setInput(x0,"x0")
Phi.evaluate()
print Phi.getOutput()
In [45]:
x = x0
xs = [x]
for i in range(N):
Phi.setInput(x,"x0")
Phi.evaluate()
x = Phi.getOutput()
xs.append(x)
In [46]:
plot(horzcat(xs).T)
legend(["p","q","c"])
Out[46]:
Rule 2: Everything is a Function (see http://docs.casadi.org)
Note 4: this is what makes CasADi stand out among AD tools
Recall
In [47]:
n = 3
A = SX.sym("A",n,n)
B = SX.sym("B",n,n)
C = mul(A,B)
print C
dotdraw(C,direction='BT')
What if you don't want to expand into scalar operations? ( avoid $O(n^3)$ storage)
In [48]:
A = MX.sym("A",n,n)
B = MX.sym("B",n,n)
C = mul(A,B)
print C
dotdraw(C,direction='BT')
What if you cannot expand into matrix operations? ( numerical algorithm )
In [49]:
C = solve(A,B)
print C
dotdraw(C,direction='BT')
In [50]:
X0 = MX.sym("x",3)
res = Phi(integratorIn(x0=X0))
XF, = integratorOut(res,"xf")
print XF
In [51]:
expr = sin(XF)+X0
dotdraw(expr,direction='BT')
In [52]:
F = MXFunction([X0],[ expr ])
F.init()
print F
In [53]:
F.setInput(x0)
F.evaluate()
print F.getOutput()
In [54]:
J = F.jacobian()
J.init()
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
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 [55]:
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 [56]:
print X.shape
print (N+1)*3+N
Demo: $\Phi(x_0,u_0)$
In [57]:
Xf = Phi( integratorIn(x0=X["x",0],p=X["u",0]) )[0]
print Xf
$ x_{k+1} - \Phi(x_k,u_k) = 0 , \quad \quad k = 0,1,\ldots, (N-1)$
In [58]:
g = [] # List of constraint expressions
for k in range(N):
Xf = Phi( integratorIn(x0=X["x",k],p=X["u",k]) )[0]
g.append( X["x",k+1]-Xf )
In [59]:
obj = X["x",N,"c"] # c_N
nlp = MXFunction( nlpIn(x=X), nlpOut(g=vertcat(g),f=obj) )
nlp.init()
print nlp
Block structure in the constraint Jacobian
In [60]:
jacG = nlp.jacobian("x","g")
jacG.init()
print jacG.getOutput().shape
with nice_stdout():
jacG.getOutput()[:20,:20].sparsity().spy()
Recall
\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 [61]:
solver = NlpSolver("ipopt",nlp)
solver.init()
solver.setInput(0,"lbg") # Equality constraints for shooting constraints
solver.setInput(0,"ubg") # 0 <= g <= 0
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
solver.setInput(lbx,"lbx")
solver.setInput(ubx,"ubx")
In [62]:
with nice_stdout():
solver.evaluate()
In [ ]:
print solver.getOutput()
In [63]:
sol = X(solver.getOutput())
In [64]:
plot(horzcat(sol["x",:]).T)
Out[64]:
In [65]:
step(range(N),sol["u",:])
Out[65]:
Showcase: kite-power optimization by Greg Horn, using CasADi backend
In [67]:
from IPython.display import YouTubeVideo
YouTubeVideo('tmjIBpb43j0')
Out[67]:
In [68]:
YouTubeVideo('SW6ZJzcMWAk')
Out[68]:
Distinction with other software:
| ACADOtoolkit | CasADi |
|---|---|
|
|
| Other operator-overloading AD tools | CasADi |
|---|---|
|
|
Closest similarity: AMPL