Create a collection of DAE models

[hsorby]

We need a collection of DAE examples that we can use to build DAE code generation abilities in libCellML.

[agarny]

Parabola model using either an ODE or a DAE approach (using CellML 1.0):

• ODE model;
• DAE model; and
• Variant DAE model.

Another simple DAE model.

[agarny]

@nickerso mentioned https://github.com/SciML/CellMLToolkit.jl on Slack. It looks pretty cool, although I haven't checked it in details. However, I have had a quick look at SciML and DifferentialEquations.jl in particular. They have a DAE example (the Robertson biochemical reactions), which they solve using the IDA solver and, indeed, they need to provide both the y0 and yp0 values to IDA (see here). (Note that as we discussed this morning, we need to provide the initial value for the DAE state variable. From there, we can compute the yp0 values for both ODE and DAE state variables.)

Anyway, that's another very simple model, which can be expressed as either an ODE or DAE model. So, perfect for our testing.

I am therefore going to create a CellML version of both approaches and add them to tests/resources/generator, together with the ones in my previous comment.

[agarny]

FWIW, I am not adding the models mentioned in my original comment. They were more toy models for me to test things in OpenCOR.

[agarny]

Food for thought: the DAE form of the Robertson model is:

dy1/dt = -k1*y1 + k3*y2*y4
dy2/dt = k1*y1 - k2*y2^2 - k3*y2*y3
1 = y1+y2+y3

and in our CellML file, y1, y2 and y3 have an initial value.

Analysis of the model:

• By default, we would say that y1, y2 and y3 can be either a constant, an ODE state or a DAE state.
• Both y1 and y2 are used in an ODE, so they are an ODE state.
• y3 is not used in an ODE, so it has to be either a constant or a DAE state.
• y3 is used in 1 = y1+y2+y3 where we know that both y1 and y2 are ODE states. So, y3 cannot be a constant and must therefore be a DAE state.

Determine the intial value of dy1/dt, dy2/dt and dy3/dt:

• The first two equations are for dy1/dt and dy2/dt, so we can compute the initial value of dy1/dt and dy2/dt.
• From 1 = y1+y2+y3 and using the first two equations, we have dy3/dt = -(dy1/dt + dy2/dt) = k2*y2^2. From there, we can compute the initial value of dy3/dt.

If we can automate that last step (for dy3/dt), then we are all good. However, this may in some (most?) cases require resorting to some algorithmic differentiation.

Also, if we were able to differentiate automatically and obtain dy3/dt, then it means that we would have converted our DAE model into an ODE model. So, in that case, we ought to solve the model as an ODE system (since it would be faster than solving it as a DAE model). (In fact, with that model, y3 can be computed directly using y3 = 1-y1-y2, but this is really a "simple" model.)

Anyway, in a lot of cases, it's not going to be that "simple" to compute the inital value of all dy/dt values.

In fact, if we want to use the IDA solver then we will need the CellML file to contain the initial values of the different ODE/DAE states (already the case), as well as the initial value of their derivative (not the case). With that information, the model would be automatically considered as a DAE model and y1, y2 and y3 as ODE/DAE states, which would make things much easier for us.

So, at this stage, I would say that we need/want to be able to specify the initial value of an ODE/DAE state in a CellML file. If we can do that, then it's going to be relatively trivial to generate code for a DAE model for use with IDA.

[agarny]

FWIW, this is the Saucerman model I was talking about. However, trying it in OpenCOR, it opens as being underconstrained, which might not be too surprising since that model is "the original unchecked version of the model imported from the previous CellML model repository, 24-Jan-2006". So, probably not a good idea to add this model to our test suite.

[agarny]

Regarding my food for thought, I guess the solution would be to have a new attribute for the variable element: initial_pvalue. Thus:

• if we are dealing with an algebraic or ODE model, then the value of initial_pvalue is ignored. If initial_pvalue is set for a variable, then the analyser will generate a warning (or the validator an error?).
• if we are dealing with a DAE model, then all ODE and DAE states should have initial_pvalue set (in addition to initial_value, of course). If a variable has initial_value set, but not initial_pvalue set, then it will be automatically considered as a constant. (I currently feel that all of this is needed if we want to use IDA as our DAE solver.)

[anandijain]

Hello, I am helping work on CellMLToolkit.jl and interested in helping build a good collection of DAE CellML models.

One neat thing I'll mention for SciML is the automated index reduction work that has recently been released. So I'm interested in seeing how our system tears biological DAEs.

We are still looking to implement benchmarks and I think it'd be great to work together to find a common set of models to benchmark the libraries against.

Currently, we have just been testing on the Physiome https://www.cellml.org/models.

[agarny]

Hi @anandijain, as you can tell, we are still very much at an early stage when it comes to supporting DAE models in libCellML.

If I am not mistaken (@nickerso to correct me, if needed), we only intend to support index-1 DAE models in libCellML since we don't have access to a symbolic math library and therefore can't do automated index reduction. (We were thinking of using SymEngine, but if I recall correctly we needed a way to convert MathML to SymEngine's internal format.)

Anyway, this means that DAE models of a higher index will first have to be converted to an index-1 DAE model before they can be used with libCellML (and IDA). When it comes to the CellML version of a DAE model, I feel like we need to provide the initial value of the state variables, but also the initial value of their derivative. IOW, we need to provide both y0 and yp0 (if we are to solve DAE models using IDA). However, this needs to be confirmed.

In summary, our collection of DAE models will be index-1 DAE models. As for our ODE test models, all of our DAE test models will be available in our GitHub repository (somewhere here), but published DAE models will first be made available on our Physiome Model Repository (PMR), to which you are more than welcome to contribute (including for ODE/DAE test models).

[anandijain]

Yes the differential states must be given u0s. https://github.com/SciML/MathML.jl This is our parser for hooking into Julia's CAS system. FYI

[agarny]

Here is an OMEX file for a DAE model. Its corresponding MATLAB version can be found here (as generated by the CellML API). As can be seen, the CellML API somehow uses 0.1 as an initial guess for all four variables involved in the non-linear algebraic system.

[agarny]

Closing this issue as a result of incorporating some DAE models directly in PR #942.