Discretize nonlinear continuos-time system using RK45 instead of Euler

[JOlucak]

In some situation we want or have to deal with discrete time problems. Depending which method is used, the accuracy is affected by the method. A good trade-off seems to be RK45 so I tried to set up the discretization. For simulation in works fine. Though, if we want to use it to parameterize the problem in becomes more difficult. I think I mess up with the different data types.

Euler is rather straight forward i.e. xk1 = xk +dt*f(x,u) which could be also easily implemented in some SOS constraints.

Below my trial to compute a discrete-time terminal LQR controller using RK45:

% define symbolic variables to setup dynamics as a function (otherwise CaSoS throws an error when the function f is setup)
x = casos.PS.sym('x',2);
u = casos.PS.sym('u',1);

% continous-time system dynamics
f0 = [x(2);
(1-x(1)^2)*x(2)-x(1)];

gx = [0;1];

xdot = f0+gx*u;

f = casos.Function('f', {x, u}, {xdot});

% RK45 scheme

% step-length and sub-steplength
dt0 = 0.1;
rk = 4;

dt = dt0/rk;

% symbolic variables needed
k1 = casos.PS.sym('k1',[length(x),rk]);
k2 = casos.PS.sym('k2',[length(x),rk]);
k3 = casos.PS.sym('k3',[length(x),rk]);
k4 = casos.PS.sym('k4',[length(x),rk]);
xk = [x casos.PS.sym('r',[length(x),rk] )];

for j=1:rk
k1(:,j) = f(xk(:,j), u);
k2(:,j) = f(xk(:,j) + dtk1(:,j)/2, u);
k3(:,j) = f(xk(:,j) + dtk2(:,j)/2, u);
k4(:,j) = f(xk(:,j) + dtk3(:,j), u);
xk(:,j+1) = xk(:,j) + dt(k1(:,j) + 2k2(:,j) + 2k3(:,j) + k4(:,j))/6;
end

% setup a function to evaluate
xk_1 = casos.Function('fk', {x,u}, {xk(:,end)}, {'x0' 'uk'}, {'xk'});

%% here comes the part I am stuck
x = casos.PS.Indeterminates('x',2);
u = casos.PS.Indeterminates('u',1);

% we want to linearize the system to obtain a discrete-time controller
[Ad,Bd] = plinearize(xk_1(x,u),x,u);

Ad = full(casadi.DM(Ad));
Bd = full(casadi.DM(Bd));

% --- Both Ad and Bd are full of zeros i.e. the derivative is wrongly calculated probably because the variables are probably messed up.

% example terminal LQR
[k,P] = dlqr(Ad, Bd, eye(2), 1);

The second issue occurs when we set up a constraint i.e.

s(V-gamma) - V( x_k1 ) - V(x_k) ...

During the build process the message with free variables appears.

My question:

• How could that be fixed? (is it even possible to use functions for that; do we have to use subs(.) for that?