Try out CasADi in your browser.
To run, click in the code cell below, and hit Shift+Enter. You may need to do 'Kernel->Restart' if you paused a while...
import casadi.*
x = SX.sym('x');
jacobian(sin(x),x)
Let us start out by finding the minimum of Rosenbrock's banana-valley function:
\begin{equation} \underset{\begin{array}{c}x, y\end{array}}{\text{minimize}} \quad x^2 + 100 \, (y-(1-x)^2)^2 \label{eq:rosenbrock} \end{equation}By inspection, we can see that its unique solution is $(0,1)$. The problem can be formulated and solved with CasADi as follows:
import casadi.*
% Symbolic representation
x=SX.sym('x');
y=SX.sym('y');
z=y-(1-x)^2;
f=x^2+100*z^2;
P=struct('x',[x;y],'f',f);
% Create solver instance
F=nlpsol('F','ipopt',P);
% Solve the problem
r=F('x0',[2.5 3.0]);
disp(r.x)
Let us reformulate the above as a constrained optimization problem, introducing a decision variable corresponding to z
above:
The problem can be formulated and solved with CasADi as follows:
import casadi.*
% Formulate the NLP
x=SX.sym('x');
y=SX.sym('y');
z=SX.sym('z');
f=x^2+100*z^2;
g=z+(1-x)^2-y;
P=struct('x',[x;y;z],'f',f,'g',g);
% Create solver instance
F=nlpsol('F','ipopt',P);
% Solve the problem
r=F('x0',[2.5 3.0 0.75],'ubg',0,'lbg',0);
disp(r.x)
Since the paper came out, a higher-level abstraction for modeling NLPs was introduced. Here's the above example remodeled:
import casadi.*
opti = Opti();
x = opti.variable();
y = opti.variable();
z = opti.variable();
opti.minimize(x**2+100*z**2);
opti.subject_to(z+(1-x)**2-y==0);
opti.set_initial(x,2.5);
opti.set_initial(y,3.0);
opti.set_initial(z,0.75);
opti.solver('ipopt');
sol = opti.solve();
[sol.value(x),sol.value(y)]
We now shift our attention to simulation and sensitivity analysis using
the CasADi's integrator
objects.
Consider the following initial-value problem in ODE corresponding to a Van der Pol oscillator:
With $p$ fixed to 0.1, we wish to solve for $x_{\text{f}} := x(1)$. This can be solved as follows:
import casadi.*
% Formulate the ODE
x=SX.sym('x',2);
p=SX.sym('p');
z=1-x(2)^2;
f=[z*x(1)-x(2)+p;x(1)];
dae=struct('x',x,'p',p,'ode',f);
% Create solver instance
op=struct('t0',0,'tf',1);
F=integrator('F','cvodes',dae,op);
% Solve the problem
r=F('x0',[0,1],'p',0.1);
disp(r.xf)
% Create Jacobian function
D=F.factory('D',{'x0','p'},{'jac:xf:x0'});
% Solve the problem
r=D('x0',[0,1],'p',0.1);
disp(r.jac_xf_x0)
By combining the nonlinear programing example with the embeddable integrator, we can implement the direct multiple shooting method by Bock and Plitt. We will consider a simple OCP reformulated as a DAE and with $p$ replaced by a time-varying control $u$:
\begin{align} \displaystyle \underset{\begin{array}{c}x(\cdot), z(\cdot), u(\cdot)\end{array}} {\text{minimize}}\quad &\displaystyle \int_{0}^{T}{ \left( x_1(t)^2 + x_2(t)^2 + u(t)^2 \right) \, dt} \\ \text{subject to} \, \quad & \left\{ \begin{array}{l} \dot{x}_1(t) = z(t) \, x_1(t) - x_2(t) + u(t) \\ \dot{x}_2(t) = x_1(t) \\ 0 = x_2(t)^2 + z(t) - 1\\ -1.0 \le u(t) \le 1.0, \quad x_1(t) \ge -0.25 \end{array} \right. \quad t \in [0,T] \\ & x_1(0)=0, \quad x_2(0)=1 \label{eq:vdp} \end{align} where $x(\cdot) \in \mathbb{R}^2$ is the (differential) state, $z(\cdot) \in \mathbb{R}$ is the algebraic variable and $u(\cdot) \in \mathbb{R}$ is the control. We let $T=10$.
Our goal is to transcribe the OCP to NLP form. In the direct approach, the first step in this process is a parameterization of the control trajectory. For simplicity, we assume a uniformly spaced, piecewise constant control trajectory: $$ u(t) := u_k \quad \text{for $t \in [t_k, t_{k+1}), \quad k=0, \ldots, N-1$} \quad \text{with $t_k := k \, T / N$.} $$
With the control fixed over one interval, we can use an integrator to reformulate the problem from continuous time to discrete time. Since we now have a DAE, we introduce an algebraic variable $z$ and the corresponding algebraic equation $g$ and use IDAS instead of CVODES. We also introduce a quadrature for calculating the contributions to the cost function.
import casadi.*
% Formulate the DAE
x=SX.sym('x',2);
z=SX.sym('z');
u=SX.sym('u');
f=[z*x(1)-x(2)+u;x(1)];
g=x(2)^2+z-1;
h=x(1)^2+x(2)^2+u^2;
dae=struct('x',x,'p',u,'ode',f,'z',z,'alg',g,'quad',h);
% Create solver instance
T = 10; % end time
N = 20; % discretization
op=struct('t0',0,'tf',T/N);
F=integrator('F','idas',dae,op);
% Empty NLP
w={}; lbw=[]; ubw=[];
G={}; J=0;
% Initial conditions
Xk=MX.sym('X0',2);
w{end+1}=Xk;
lbw=[lbw;0;1];
ubw=[ubw;0;1];
for k=1:N
% Local control
name=['U' num2str(k-1)];
Uk=MX.sym(name);
w{end+1}=Uk;
lbw=[lbw;-1];
ubw=[ubw; 1];
% Call integrator
Fk=F('x0',Xk,'p',Uk);
J=J+Fk.qf;
% Local state
name=['X' num2str(k)];
Xk=MX.sym(name,2);
w{end+1}=Xk;
lbw=[lbw;-.25;-inf];
ubw=[ubw; inf; inf];
G{end+1}=Fk.xf-Xk;
end
% Create NLP solver
nlp=struct('f',J,'g',vertcat(G{:}),'x',vertcat(w{:}));
S=nlpsol('S','ipopt',nlp,struct('ipopt',struct('tol',1e-6)));
% Solve NLP
r=S('lbx',lbw,'ubx',ubw,'x0',0,'lbg',0,'ubg',0);
disp(r.x);
Again, we note that the same can be achieved using a higher-le}vel abstraction for modeling NLPs:
import casadi.*
% Formulate the DAE
x=SX.sym('x',2);
z=SX.sym('z');
u=SX.sym('u');
f=[z*x(1)-x(2)+u;x(1)];
g=x(2)^2+z-1;
h=x(1)^2+x(2)^2+u^2;
dae=struct('x',x,'p',u,'ode',f,'z',z,'alg',g,'quad',h);
% Create solver instance
T = 10; % end time
N = 20; % discretization
op=struct('t0',0,'tf',T/N);
F=integrator('F','idas',dae,op);
opti = Opti();
% Build up objective
J = 0;
% Initial conditions
Xk=opti.variable(2);
opti.subject_to(Xk==[0;1]);
for k=1:N
% Local control
Uk=opti.variable();
opti.subject_to(-1<=Uk<=1);
% Call integrator
Fk=F('x0',Xk,'p',Uk);
J=J+Fk.qf;
% Local state
Xk=opti.variable(2);
opti.subject_to(Xk(1)>=-0.25);
opti.subject_to(Fk.xf==Xk);
end
opti.minimize(J);
% Create NLP solver
opti.solver('ipopt',struct,struct('tol',1e-6));
% Solve NLP
sol = opti.solve();
sol.value(opti.x')