Question(s) about code generation & usage in Matlab/Simulink

[ottarg]

Hello everyone!

I will start by doing my best to avoid the x-y problem thing, by first stating what it is that I actually want as an end result: My ultimate goal is to generate C/C++ code from https://pybamm.org/ models and run them in Simulink and on embedded systems. I am approaching this by combining some built in methods of PyBamm/CasADi, but if I am going at this from the wrong way, I would be more than happy to do something else.

What I understand so far:

• PyBamm models have a model.generate() method which produces generated CasADi code. Some models are ODE models, some are DAE, but all produce functions for the right-hand-side, algebraic equations & Jacobians upon generation. Example
  • To use these methods, I need to pair them with a solver of some sort, since only the constituent system equations are exported, not the solver
  • I have been successful in running these models in Matlab/Simulink with ode15s (matlab) and in Simulink using interpreted matlab functions + the algebraic constraint block
    • The interpreted matlab blocks are quite slow, and I know that I can/should use S-functions instead, but I don't yet know how to do that.

My question:

• Do I understand correctly that if I were instead to construct a system model + solver in CasADi, I could generate code from that entire thing and use it in an s-function/on embedded hardware as a standalone block?
  • If so, I would very much appreciate some help with that! I have seen some examples and discussion about which solvers are compatible with code generation, but when I try to construct a CasADi system from my C++ methods, I get errors which I can't reconcile.
  • I would then also wonder if I should not be invoking the code generation method at a different point in PyBamm? Surely it must somewhere construct a complete CasADi system model with solver, which presumably would then be a "Function" class object, which could generate code for me?

I have tried my best to follow examples/tutorials to build the system in CasADi, but the error I get is that my rhs function does not have a derivative. Below is my test code and error:

import casadi.*
clearvars 
% Load external functions
x0_fn = external('x0', './dfn.so');
z0_fn = external('z0', './dfn.so');
rhs_fn = external('rhs_', './dfn.so');
alg_fn = external('alg_', './dfn.so');
jac_rhs_fn = external('jac_rhs', './dfn.so');
jac_alg_fn = external('jac_alg', './dfn.so');
var_fn = external('variables', './dfn.so');

% Define state and algebraic variable initial conditions
x0 = x0_fn();
z0 = z0_fn();
% Define number of states and algebraic variables
nx = length(x0.o0); % Number of ODE states
nz = length(z0.o0); % Number of algebraic variables

% Define symbolic variables
t = MX.sym('t');
x = MX.sym('x', nx);
z = MX.sym('z', nz);
p = MX.sym('p', 1); % My model has just one parameter, the input cell current

odes = rhs_fn(t,x, z, p);
algs= alg_fn(t,x, z, p);

opts = struct;
% Create CasADi integrator for DAE system
dae = struct('t',t,'x',x,'z',z,'p',p,'ode',odes,'alg',algs);
dae_integrator = integrator('F', 'idas', dae,opts);

I should note here, that when I tried myself to make a simple newton method solver with these generated functions, I found that the Jacobian of the rhs (from jac_rhs)was not full rank, which was a stopping point for me. I wonder if this issue is related? I know I am not providing that jac_rhs function at all here, but maybe the system is truly somehow deficient? Very confusing to me though, is that the system does definitely still work both in ode15s and Simulink. In ode15s, if I provide no Jacobian function, it works but when I do supply it with the (rank deficient) Jacobian it also still works -- and in fact is more stable and faster, which perplexes me.

I am also not sure if I should be trying to make a solver object or an integrator object? I am not trying to do any NLP/optimization, like I have seen in CasADi-driven MPC in Simulink.

Any help or pointers are very welcome!

[ernestds]

You have to provide to the rhs function the jacobian as a cached function.
In Python you'd do something like this:

ca.external('name',folder,{'cache': {'jac_name': jac_rhs_fn}})

the name, input and output sparsities must match exactly for the jacobian, so for example the jac_rhs_fn.name() should be 'jac_name'

[ottarg]

Thanks for the response! It seemed like there must be something like that. I will give it a try!

Since writing my original question, I was able to find an alternative way to construct the integrator not by using the compiled C++ methods, but instead by digging into the PyBaMM source a bit and finding where the integrator is constructed in their process and then applying the age old trick of "monkey see, monkey do"

Unfortunately, I also learned that the model I want to compile to C/C++ is a DAE system and they are using the 'idas' integrator, which apparently can not be used with code generation. From what I can see, the only other integrator type that can handle DAEs is the collocation one. I don't know much about collocation, but from what I do know it seems like it always solves the entire timeframe simultaneously and so would never be suitable for an embedded solution. In general I am very new to DAEs so I'm not really sure what my options are. Do you know of any integrator plug-ins which are perhaps not shipped with the default package that are compatible with code-gen and can do DAEs? Is there any way for me to convert a DAE system into an ODE one?

Sorry for bombarding you with questions, I hate to punish you for the good deed of being helpful, but the help is very much appreciated!