In [1]:
from pylab import *
from casadi import *
from casadi.tools import * # for dotdraw
%matplotlib inline
In [2]:
x = SX.sym("x") # scalar symbolic primitives
y = SX.sym("y")
z = x*sin(x+y) # common mathematical operators
In [3]:
print z
In [4]:
dotdraw(z,direction="BT")
In [5]:
J = jacobian(z,x)
print J
In [6]:
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 [7]:
print x*y/x-y
In [8]:
H = hessian(z,x)
print H
Sort graph into algorithm
In [9]:
f = SXFunction("f",[x,y],[z])
print f
Note 2: re-use of tape variables: live-variables
In [10]:
print f([1.2,3.4])
In [11]:
print f([1.2,x+y])
In [12]:
f.generate("f")
print file("f.c").read()
In [13]:
A = SX.sym("A",3,3)
B = SX.sym("B",3)
print A
In [14]:
print solve(A,B)
In [15]:
print trace(A) # Trace
In [16]:
print mul(A,B) # Matrix multiplication
In [17]:
print norm_F(A) # Frobenius norm
In [18]:
print A[2,:] # Slicing
Rule 1: Everything is a matrix
In [19]:
print A.shape, z.shape
In [20]:
I = SX.eye(3)
print I
In [21]:
Ak = kron(I,A)
print Ak
Rule 1: Everything is a sparse matrix
In [22]:
Ak.sparsity().spy()
In [23]:
A.sparsity().spy()
In [24]:
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 [25]:
t = SX.sym("t") # time
u = SX.sym("u") # control
x = SX.sym("[p,q,c]") # state
p,q,c = x
In [26]:
ode = vertcat([(1 - q**2)*p - q + u, p, p**2+q**2+u**2])
print ode, ode.shape
In [27]:
J = jacobian(ode,x)
print J
In [28]:
f = SXFunction("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",str(f).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 [29]:
print I
In [30]:
for i in range(3):
print ffwd([ t,u,x, ode, 0,0,I[:,i] ])[0]
In [31]:
print J
Performing adjoint sweeps gives the rows of J
In [32]:
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
$\dot{x}=f(x,u,t)$ with $x = [p,q,c]^T$
In [33]:
f = SXFunction("f",daeIn(x=x,t=t,p=u),daeOut(ode=ode))
f.printDimensions()
Construct an integrating block $x_{k+1} = \Phi(f;\Delta t;x_k,u_k)$
In [34]:
tf = 10.0
N = 20
dt = tf/N
In [35]:
Phi = Integrator("Phi","cvodes",f,{"tf":dt})
x0 = vertcat([0,1,0])
print Phi({"x0":x0})
In [36]:
x = x0
xs = [x]
for i in range(N):
x = Phi({"x0":x})["xf"]
xs.append(x)
In [37]:
plot(horzcat(xs).T)
legend(["p","q","c"])
Out[37]:
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 [38]:
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 [39]:
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 [40]:
C = solve(A,B)
print C
dotdraw(C,direction='BT')
In [41]:
X0 = MX.sym("x",3)
XF = Phi({"x0":X0})["xf"]
print XF
In [42]:
expr = sin(XF)+X0
dotdraw(expr,direction='BT')
In [44]:
F = MXFunction("F",[X0],[ expr ])
print F
In [45]:
print F([x0])
In [46]:
J = F.jacobian()
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 [47]:
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 [48]:
print X.shape
print (N+1)*3+N
Demo: $\Phi(x_0,u_0)$
In [49]:
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 [50]:
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 [52]:
obj = X["x",N,"c"] # c_N
bugfix_140 = {"ad_weight_sp": 0} # This can be removed when 2.4.1 is released
nlp = MXFunction("nlp",nlpIn(x=X), nlpOut(g=vertcat(g),f=obj) ,bugfix_140)
print nlp
Block structure in the constraint Jacobian
In [53]:
jacG = nlp.jacobian("x","g")
S = jacG.sparsity_out(0)
print S.shape
DM.ones(S)[: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 [54]:
solver = NlpSolver("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 [55]:
sol_out = solver({
"lbg": 0, # Equality constraints for shooting constraints
"ubg": 0, # 0 <= g <= 0
"lbx": lbx,
"ubx": ubx
})
In [56]:
print sol_out["x"]
In [57]:
sol = X(sol_out["x"])
In [58]:
plot(horzcat(sol["x",:]).T)
Out[58]:
In [59]:
step(range(N),sol["u",:])
Out[59]:
Showcase: kite-power optimization by Greg Horn, using CasADi backend
In [60]:
from IPython.display import YouTubeVideo
YouTubeVideo('tmjIBpb43j0')
Out[60]:
In [61]:
YouTubeVideo('SW6ZJzcMWAk')
Out[61]:
Distinction with other software:
| ACADOtoolkit | CasADi |
|---|---|
|
|
| Other operator-overloading AD tools | CasADi |
|---|---|
|
|
Closest similarity: AMPL