Consider the initial value problem given below where the behaviour of the system is studied with respect to given set of inputs and initial condition. Lets make it parameteric by adding a scalar parameter to the ode and setting its value to zero.

ODE 1

\begin{aligned} (t_0,t_f,N)&=(0,10,25)\\ x(t_0)&=\begin{bmatrix} 0\\ 0 \end{bmatrix} \\ u(t)&=0.2t-1\\ p&=0\\ \begin{bmatrix} \dot{x}_0\\ \dot{x}_1 \end{bmatrix}&= \begin{bmatrix} (1-x_1^2)x_0-x_1+u\\ x_0 \end{bmatrix}+p \end{aligned}

The problem can be solved in 2 ways using Simulate.

1. eliminate controls from the ode

2. supply controls evaluated on the grid points

Elimination

• Firstly, the ode needs to modelled for which we call the get_symbols() function. The arguments are the number of states, controls, parameters and intervals for the grid. This also has an additional option to model the ode using symbolics of the SX or MX class. MX is usually preferred due to its generality over SX. The get_symbols() return a dictionary with required symbolics to construct the ode. This dictionary will keep updating as we pass it through each function.

• Using the symbolics stored in the dictionary, the ode is constructed. The ode will be passed to the scale_ode() method where, a time domain transformation takes place. In the scaled ode, the initial and final time are simply parameters.

• Post domain transformation, the numerical solver for integration and its options must be set. Currently, CVODES, rk4, IRK methods are supported. For detailed options list, refer to integrator section of CasADi API. Calling the integrator() creates an integrator object for the scaled ode.

• Lastly, simulate() is called where the initial state, initial and final time, control and parameter if any is passed. The system is integrated and results are stored. This can be accessed using keys ‘x_res’,’u_res’ and ‘t_grid’ if required.

import simulatetraj as a
import casadi as cs
sym_dict=a.get_symbols(n_x=2,n_u=0,n_p=1,N=25)
t=sym_dict['t']
x=sym_dict['x']
u=sym_dict['u']
p=sym_dict['p']
f=cs.vertcat((1-x[1]**2)*x[0]-x[1]+0.2*t-1,x[0])+p
_=a.scale_ode(sym_dict=sym_dict,f=f)
_=a.integrator(sym_dict=sym_dict,options=None)
a.simulate(t0=0,tf=10,x0=cs.DM([0,0]),p=cs.DM(0),sym_dict=sym_dict)
a.plot_sol(sym_dict=sym_dict)

without control elimination

The control input is approximated as a piecewise constant function.

import simulatetraj as a
import casadi as cs
sym_dict=a.get_symbols(n_x=2,n_u=1,n_p=1,N=25)
t=sym_dict['t']
x=sym_dict['x']
u=sym_dict['u']
p=sym_dict['p']
f=cs.vertcat((1-x[1]**2)*x[0]-x[1]+u,x[0])+p
_=a.scale_ode(sym_dict=sym_dict,f=f)
_=a.integrator(sym_dict=sym_dict,options=None)
a.simulate(t0=0,tf=10,x0=cs.DM([0,0]),u=cs.linspace(-1,1,25).T,p=cs.DM(0),sym_dict=sym_dict)
a.plot_sol(sym_dict=sym_dict)

ODE 2

In this example, the derivative of the square function is integrated.

\begin{aligned} (t_0,t_f,N)&=(0,10,100000)\\ x(t_0)&=0\\ \dot{x}&=2t \end{aligned}

sym_dict=get_symbols(n_x=1,n_u=0,n_p=0,N=100000)
t=sym_dict['t']
f=2*t
_=scale_ode(sym_dict=sym_dict,f=f)
_=integrator(sym_dict=sym_dict,options=None)
x0=cs.DM([0])
t0,tf=0,10
simulate(t0=t0,tf=tf,x0=x0,sym_dict=sym_dict)
plot_sol(sym_dict=sym_dict)
print('Global error x(N+1) for xdot=2t:',sym_dict['x_res'][:,-1]-100)

Global error x(N+1) for xdot=2t: DM(1.49132e-08)

Lotka voltera/prey-predator model

\begin{aligned} (t_0,t_f,N)&=(0,15,1000)\\ (\alpha,\beta)&=(0.01,0.02)\\ x(t_0)&= \begin{bmatrix} 20\\ 20 \end{bmatrix}\\ \begin{bmatrix} \dot{x}_0\\ \dot{x}_1 \end{bmatrix}&= \begin{bmatrix} x_0-\alpha x_0 x_1\\ -x_1+\beta x_0 x_1 \end{bmatrix} \end{aligned}

sym_dict=get_symbols(n_x=2,n_u=0,n_p=2,N=1000)
x=sym_dict['x']
p=sym_dict['p']
f=cs.vertcat(x[0]-p[0]*x[0]*x[1],-x[1]+p[1]*x[0]*x[1])
_=scale_ode(sym_dict=sym_dict,f=f)
_=integrator(sym_dict=sym_dict,options=None)
x0=cs.DM([20,20])
t0,tf=0,15
p=cs.DM([0.01,0.02])
simulate(t0=t0,tf=tf,x0=x0,p=p,sym_dict=sym_dict)
plot_sol(sym_dict=sym_dict)
plt.plot(sym_dict['x_res'][0,:].full().flatten(),sym_dict['x_res'][1,:].full().flatten(),'o')
plt.show()