c code generation for ODE/DAE

[Horea Alexandru Caramizaru]

Hi All,

I was wondering if there is someone here who can help me with some problems I got working with CasADi.

I've seen that in the current version is possible to generate c code for different symbolic expressions and then reuse by doing a function binding. ( in matlab )

I would also want to generate c code for some ODE/DAE and use it afterwords by calling them inside CasADi Integrator. But this doesn't seem possible since CasADi's Integrator doesn't accept Functors ( binding functions ) as a parameter but only symbolic/struct parameter.

Is there a workaround / different flow for this use case or current version doesn't support code generation for ODE/DAE?

In case of a solution, I would appreciate an example since currently, it seems that only simple expressions can be considered for c code generation ( and even if you are allowed to generate you can't use them because of parameter restriction inside "Integrator" )

Many thanks!

Horea Caramizaru

[Joris Gillis]

Dear Horea,

Indeed, the integrator variant accepting a Function argument seems to be deactivated in Matlab/Python.

You could always call the function binding with symbolic arguments to get expressions again.

Alternatively, use the 'jit' family of options. For example: import casadi.*

x = SX.sym('x');


sys = struct;
sys.x = x;
sys.ode = x^2;

intg_opts = struct;
intg_opts.expand = true;

helper_opts = struct;
helper_opts.jit = true;
helper_opts.compiler = 'shell';
helper_opts.jit_options.verbose = true;
intg_opts.common_options.final_options = helper_opts;

intg = integrator('intg','cvodes',sys,intg_opts);

intg('x0',0);

Best regards,

  Joris

[Horea Alexandru Caramizaru]

Dear Joris,

Thanks for your response.

Unfortunately, the approach using 'jit' has its own limitations. Let's say I want to change the interval time ( t0, tf ) for integration. For that, I'll have to recreate an integrator ( e.g. to regenerate/recompile ) which makes the process useless since the reason behind using this approach was to minimize the time needed. ( or just change some other options passed to the integrator... )

You were also mentioning "You could always call the function binding with symbolic arguments to get expressions again." but I can't say it's clear for me what I need to do. Can you give me an example?

Also, is there a technical reason to have Function arguments deactivated in Matlab/Python? Any chance this to be changed in the near future?

Best regards,

Horea C.

[Horea Alexandru Caramizaru]

Hi all,

I was trying to find out whether the disable constructor "integrator" actually works. For that I've recompiled CasADi and regenerate the Matlab bindings with that specific call activated:

//#ifndef SWIG
CASADI_EXPORT Function integrator(const std::string& name, const std::string& solver,
const Function& dae, const Dict& opts=Dict());

//#endif // SWIG

Unfortunately, there seems to be a bug related to this constructor, as I get an 'out of bound' error wrt to the number of the parameters of my compiled/uncompiled functor that is built around my DAE.

The code I am trying is this:

import casadi.*

dae = DaeBuilder();

% Add input expressions

a = dae.add_p('a');

b = dae.add_p('b');

u = dae.add_u('u');

h = dae.add_x('h');

v = dae.add_x('v');

m = dae.add_x('m');

% Constants

g = 9.81; % gravity

% Add output expressions

hdot = v;

vdot = (u-a*v^2)/m-g;

mdot = -b*u^2;

dae.add_ode('hdot', hdot);

dae.add_ode('vdot', vdot);

dae.add_ode('mdot', mdot);

% Specify initial conditions

dae.set_start('h', 0);

dae.set_start('v', 0);

dae.set_start('m', 1);

% Add meta information

dae.set_unit('h','m');

dae.set_unit('v','m/s');

dae.set_unit('m','kg');

% Print DAE

disp(dae, true);

f = dae.create('f',...

     {'x','u','p'},{'ode'});

opts = struct('mex', true);

C = CodeGenerator('gen.c',opts);

C.add(f);

C.generate();

C = Importer('gen.c','shell');

f = external('f',C);

disp(f);

intg_opts.t0=0;

intg_opts.tf=0.5;

intg_opts.verbose = true;

intg = integrator('intg','idas',f);

The output and error is:

ns = 0, nx = 3, nz = 0, nq = 0, ny = 0, np = 2, nd = 0, nu = 1

Variables

  t = t

  x = [h, v, m]

  p =  [a, b]

  u =  [u]

Differential equations

  der_h == v

  der_v == (((u-(a*sq(v)))/m)-9.81)

  der_m == (-(b*sq(u)))

f:(x[3],u,p[2])->(ode[3]) MXFunction

Error using casadi.integrator (line 517)

.../casadi/core/function_internal.cpp:145: Error calling IdasInterface::init for 'intg':

Error in Function::sparsity_in for 'f' [MXFunction] at .../casadi/core/function.cpp:897:

vector::_M_range_check: __n (which is 3) >= this->size() (which is 3)

Error in untitled (line 86)

intg = integrator('intg','idas',f);

Like I said in the beginning, there seems to be a bug here. I was trying to see whether there is an obvious error ( something that I can see without putting in place a development environment ) but I wasn't successful in finding a solution. Does anyone have any idea about this problem so that I can fix it faster? :)

The second approach will be to do what Joris has proposed: " You could always call the function binding with symbolic arguments to get expressions again." But it's unclear to me what he means...

Any help will be greatly appreciated :)

Best regards,

Horea Caramizaru

[Joris Gillis]

Dear Horea,

You can work with a time-transformation to keep t0 and tf fixed at, say, 0 and 1.

"call the function binding with symbolic arguments to get expressions again":

x = MX.sym('x',f.size_in('x'));
u = MX.sym('u',f.size_in('u'));
p = MX.sym('p',f.size_in('p'));

intg = integrator('intg','idas',struct('x',x,'p',[u;p],'ode', f(x,u,p) ));

IDAS needs derivatives, so beter code-generate them as well:

C = CodeGenerator('gen.c',opts);
C.add(f);

C.add(f.jacobian());
C.generate();

Best regards,

  Joris

[Matias Tofteby]

Is there any way to include the quad function into the integrator in a similar fashion?

Something like:

f = dae.create('f',...

     {'x','u','p'},{'ode','quad'});

x = MX.sym('x',f.size_in('x'));
u = MX.sym('u',f.size_in('u'));
p = MX.sym('p',f.size_in('p'));

intg = integrator('intg','idas',struct('x',x,'p',[u;p],'ode', f(x,u,p) ,'quad',f.quad));

Best regards,

Matias.

[Joel Andersson]

You can divide up that expressions into two rows:

 

ode, quad = f(x,u,p)

intg = integrator('intg','idas',struct('x',x,'p',[u;p],'ode', ode ,'quad',quad));

But it's a bit of a roundabout way of doing it: From DaeBuilder you already have the expressions, so creating a Function instance only to get back expressions is unnecessary.

I created an issue to make provide a simpler syntax for this: https://github.com/casadi/casadi/issues/2783

There is a lot of work going on right now to make DaeBuilder more mature and align with the FMI standard. There is also work to improve the interoperability with Modelica.

As of now, DaeBuilder is not very mature, so if you want to work with a mature code, you might want to avoid using DaeBuilder altogether for now.

Joel 

[Joel Andersson]

You can try the following by the way:

x = vertcat(dae.x);

u = vertcat(dae.u);

p = vertcat(dae.p);

ode = vertcat(dae.ode);

quad = vertcat(dae.quad);

dae = struct('x', x, 'p', [u;p], 'ode', ode, 'quad', quad);

intg = integrator('intg','idas',dae);

...

That way, you'll avoid the (unnecessary) Function creation.

Joel