Basic continuous time simulation of a dae system

[Krishnakumar Gopalakrishnan]

I am really sorry if this question is too basic.  Before using all the advanced ideas in casADi, i.e. optimal control, sensitivity analysis etc, I just wish to make a simple time-domain simulation of my DAE model from t0 to tf. The idea is to simply plot the evolution of differential and algebraic variables over time (of the basic open loop system).

I have got started on the Robertson problem having formulated these as semi-explicit DAEs in MATLAB.

clear;close all;clc;format short g;
import casadi.*

y1 = SX.sym('y1');
y2 = SX.sym('y2');
y3 = SX.sym('y3');
dae = struct('x',[y1;y2],'z',y3,'ode',[-0.04*y1 + 1e4*y2*y3;0.04*y1 - 1e4*y2*y3 - (3e7)*y2.^2],'alg',y1+y2+y3-1);
F = integrator('F','idas',dae);
result = F('x0',[1;0],'z0',0);
disp(result.xf);

Any help on this is much appreciated.

[Adrian Bürger]

You are currently evaluating your integrator only once. For a simulation over a given time horizon you can put it into a loop, while initializing with the result of the previous step. Assuming you you have a time grid with equidistant time points, this could look something like this:

t0 = 0.0;

tf = 3;

N = 100;

dt = (tf - t0) / N;

y1 = SX.sym('y1');

y2 = SX.sym('y2');

y3 = SX.sym('y3');

dae = struct('x',[y1;y2],'z',y3,'ode',[-0.04*y1 + 1e4*y2*y3;0.04*y1 - 1e4*y2*y3 - (3e7)*y2.^2],'alg',y1+y2+y3-1);

F = integrator('F','idas',dae, struct('t0', 0.0, 'tf', dt));

x0 = [1;0];

z0 = 0;

x_sim = {x0};

z_sim = {z0};

for k = 1:N-1

    

    result = F('x0', x_sim{end}, 'z0', z_sim{end});

    x_sim = {x_sim{:} result.xf};

    z_sim = {z_sim{:} result.zf};

end

x_sim = full(horzcat(x_sim{:}));

z_sim = full(horzcat(z_sim{:}));

figure();

plot(x_sim');

figure();

plot(z_sim);

[Krishnakumar Gopalakrishnan]

Hi Adrien,

Thanks a lot for your reply and the excellent answer. I can confirm that it works

There is one major hiccup though. IDA is an adaptive time-step solver, isn't it, i.e. it can automatically vary the time-steps depending on the smoothness of the solution.  In my understanding, it will take short steps if it detects a large change in solution, and progressively longer steps when the solution dynamics are slow, as we progress in time.

By forcing the use of a fixed time-step, are we losing these benefits? 

Regards,

Krishna

PS:
It could also be that the integrator is internally using adaptive time-steps, and simply reports the value at the end of the interval.  Even in this case, using loops etc. in MATLAB could be slower. It would be nice if we can simply specifiy the simulation end-time and ask IDA to return the solution at some pre-specified grid of time-points.

[Adrian Bürger]

Hi again Krishnakumar,

the case is what you wrote in your PS: in the example above, we do not fix the step size of the integrator, but set the end time point for one single integration step. The whole simulation then is a sequence of such steps where each step starts at the end point of the previous one. So we basically determine at which points over whole horizon we want to access the values of the states, but not the steps that IDAS takes internally.

I do not know whether it is possible to access the state values at the steps that IDAS takes internally so that you could use one single integrator call for your whole horizon. However, there exist a CasADi functionality called mappacum that might help you to prevent looping in case this aspect would actually become time critical, but once this is not the case I would suggest to stick to the loop for now for the sake of simplicity.