DAEs and specifying states

[B-LIE]

I'm currently using Modelica, which is a modeling language for DAEs -- Modia for Julia is similar in scope. I have used the DifferentialEquations.jl DAE solvers; they run much faster than the Modelica tools I've tested, but have some disadvantages:

1. In Modelica, it suffices to specify the state of the system, and then Modelica figures out the initial conditions for other variables. In Julia's DifferentialEquations package, I must specify the initial values for all differential variables, their derivatives, as well as the algebraic variables.

2. In Modelica, I can specify which initial values are fixed, and which are not. If I specify the initial values of states as fixed, Modelica adjusts all other initial values such that the initial states are as given. If I over-specify the fixed initial values, I get an error message in Modelica.

This Modelica property of knowing that the specified states are not changed by the solver initiation procedure is essential for doing, e.g., state estimation for DAE models. And for knowing that I get the intended solution.

I'm here using the standard definition of states as in systems theory: basically, the state is a minimal length vector which uniquely specifies the solution of the model (together with possible input functions).

Observe that state is not the same as differential variables. The following examples makes this clear: -- if I model two water pipes (incompressible water, inelastic walls) using Newton's second law, and add the constraint that the same water runs through both of them, I have 2 differential variables (e.g., the volumetric flow rate through the pipes), but with an algebraic constraint that the two volumetric flow rates are the same. Thus: two differential variables, but there is only one state -- the common volumetric flow rate.

-- higher index DAEs may have more states than differential variables. Using Pantelidis' algorithm, the DAE may be converted to higher order DAEs.

In summary, it is the states that I need to be able to specify in such a way that the initiation routine doesn't adjust their values, not the differential variables (for very simple systems, these are, of course, the same, but not in general).

Chris Rackauckas told me that the Sun dials solvers actually allows to fix the initial values of the differential values. Again, it is the initial values of the states I need to fix. As far as I know, Modelica tools (e.g., OpenModelica) symbolically reorganizes the DAEs into ODEs, and for that case, it is (in principle) "straight forward" to fix the initial values of the states.

QUESTION: would it be possible to add the possibility to fix the initial states of DAEs in Julia's DifferentialEquations.jl package (but allow for adjusting differential variables which are not states + the derivative of differential variables + algebraic variables to achieve consistency), so that I can be sure that I get the intended solution, and so that I can use the DAE solvers of the DifferentialEquations.jl package for state estimation, etc.?

CR suggested I ping @kanav99 on this.

[ChrisRackauckas]

In Modelica, it suffices to specify the state of the system, and then Modelica figures out the initial conditions for other variables. In Julia's DifferentialEquations package, I must specify the initial values for all differential variables, their derivatives, as well as the algebraic variables.

You always have to give an initial value to the rootfinder. IIRC, Modelica tries 0, so you could always do the same.

In Modelica, I can specify which initial values are fixed, and which are not. If I specify the initial values of states as fixed, Modelica adjusts all other initial values such that the initial states are as given. If I over-specify the fixed initial values, I get an error message in Modelica.

We could add that. Sundials doesn't natively do that.

http://docs.juliadiffeq.org/latest/solvers/dae_solve.html#Sundials.jl-1

init_all = false

If that's set to true then all are init'd, otherwise only the states of the differential variables are fixed while the states of the algebraic ones are changed. This is a limitation of Sundials and will be a limitation until we have a native Julia method. And that's precisely what we'll be working on after the nlsolver changes give us a native DAE method in OrdinaryDiffEq.

This is our route to native DAEs in OrdinaryDiffEq:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/305#issuecomment-504809270

• [x] Expand the signature of OrdinaryDiffEq to allow DAEProblem
• [x] Allow W-operator to work on DAEs somehow? In progress
• [x] Newton needs to get setup for a DAE. In progress
• [x] Initial dt method for DAEs
• [x] Consistent IC finding

Once we have that for an Euler method, all of the stuff like interpolations and event handling comes for free, then all we need to do is add the higher order methods.

[B-LIE]

OK -- let me clarify on the part of Modelica only needing initial values for states -- perhaps I should have been more precise. Modelica is a language (like C, etc.). Actual tools to solve the DAE models are, e.g., Dymola (commercial; currently owned by Dassault but initiated by Hilding Elmquist who is also behind Modia), OpenModelica (free; developed at Linköping University in Sweden, with co-operation with a German research group, etc.), etc. -- these are more like Visual C, etc.

In OpenModelica, the traditional solution strategy has been to start with a DAE: F(dy,y,u; p, t) = 0 and translate it into a form dy = f(y,u;p,t) which is an ODE, and AEs of type g(z;x,u,p,t) = 0 where the AE part is converted to a (block) triangular set of algebraic equations where the "state" (x), inputs (u), parameters (p), etc. are known.

Thus, if the state is specified and if the AE part has unique solutions, then only the initial state is needed. This is, of course, the case if the AEs are linear, or if the AEs are nonlinear and really only have a single realistic solution (e.g., suppose one knew that the solution to the nonlinear AEs have to be positive real, that might in some cases make the solution unique).

As you indicate, in the general case, nonlinear AEs may have multiple solutions, and in that case it is -- of course -- vital to give a good initial guess.

When I did a comparison of OpenModelica, as run from Python via the interface API OMPython ("OpenModelica Python") against Julia DAE solvers some 1 1/2 years ago, I observed the following:

• In my case study, for OpenModelica, I only needed to specify the initial state. [Yes, initial guesses for the algebraic variables defaults to zero.]
• For the Julia DAE solver, I needed to specify initial values for all variables, i.e., for differential variables, for algebraic variables, and for derivatives. And it didn't work if I specified the same initial values for the states as in OpenModelica, and zero for all the others. I actually had to find initial values for the algebraic variables in Julia that were relatively close to the true values (which OpenModelica finds automatically by solving g(z;x,u,p,t) = 0 -- given that the solution is unique).
• I did set the initial guess of derivatives to zero, though.

Anyway, my point is that it is easier to use OpenModelica than the Julia DAE solver because normally, OpenModelica finds a decent solution even if when I don't bother about the initial AE values -- Julia DAE solver appeared to be more sensitive to this. I am pretty sure the reason is the way OpenModelica decomposes the DAE into an ODE part and a (block) triangular DAE part.

I should say that with this approach, I think OpenModelica actually uses an ODE solver and adds some nonlinear solvers (including continuation methods). Newer versions of OpenModelica also offers solution by proper DAE solvers. I don't know whether this advantage is maintained in that case.

In Modelica, the initial values are specified as Real T(start=3,fixed = true) for a real variable T starting at value 3. For a variable where I want to suggest the initial value and leave some freedom to the initiation procedure, the syntax is Real T(start=3,fixed=false).

Anyway, I should also mention that the Julia DAE solver was some 10 times faster than the OpenModelica solver. Which I was somewhat disappointed by -- after all, OpenModelica translates the equation based model into assignments in C + compiles the code into an exe file. There is probably some overhead in the OMPython API (using ZMQ), but still...

[ChrisRackauckas]

Putting this on DiffEq7. From implementing the initialization algorithms, I realized that indeed all you really need are the states. The algebraic equation initial values and the du0 are just used in some (not even all...) initialization algorithms, even then they are only used as the initial guess to the nonlinear solve for the initial values. So they are essentially initial guess throw aways and thus we shouldn't require that the user gives us them: they should be extra niceties but not required.

[ramboldio]

Maybe to get a first prototype going of being able to simulate a modelica model with DiffEq, we could use the xml representation of the DAE system that openmodelica generates (Example). It even comes with initial values AFAIK.

Would this be interesting/useful to try out?

(Working on "TrussFab" and also looking for solutions to speedup OpenModelica simulations and doing more advanced stuff with them)

[kanav99]

Sounds like a nice project @ChrisRackauckas

[ChrisRackauckas]

Yes, it would be really nice to have tools that read modelica formats into ModelingToolkit so we can do all sorts of nice codegen on it.

[ChrisRackauckas]

This is completed via the the ModelingToolkit work and all, and OpenModelica.jl stuff reads Modelica into ModelingToolkit, etc.