Convenience feature: Include algebraic variables in result

[MartinOtter]

In Modia/Modia3D there are usually many auxiliary variables that are computed in the model function and that should be available in the result. For example, for mechanical systems, not only the states (e,g, relative coordinates of joints) are interesting, but also forces and torques and relative positions between different bodies. It would be convenient to include such variables already in the standard result structure of DifferentialEquations. Issues that might need to be considered:

1. Some algebraic variables are needed to compute the state derivatives and need therefore be computed whenever derivatives need to be computed. Mathematically: w = f1(x,t); dx/dt=f2(x,w,t).
2. Some algebraic variables are a pure function of the states and are only used for the output. Mathematically dx/dt=f1(x,t), w = f2(x,t). These algebraic variables need to be computed only at communication points and not for every model function evaluation. They could be also computed after the whole simulation run is finished (but see problem below). When computing only at communication points, it would be nice to have a helper function, such as, "isCommunicationPoint(integrator)" and compute these variables only in such a case (to enhance efficiency).
3. Some algebraic variables are clocked variables: Synchronous languages typically abstract sampled data systems with clocks: A variable is associated to a clock and has only a value, when the clock ticks (so at a special event). Note, a continuous-time variable has (at least) two values at an event instant: the left limit and the right limit to the event instant (so in the result, at the same time instant there are two values, as in DifferentialEquations, when I understand correctly). Contrary to this, a clocked variable has only one value at an event instant. In Modelica tools, clocked variables are represented in a plot as small circle at the clock instant and grey-dashed line between event instants (holding the last value). Clocked variables are partially implemented in Modia currently.

Problems:

• Algebraic variables can have any type, e.g., it is common to have Boolean variables. For example, in a simple bouncing ball model there could be a Boolean variable "flying" that defines whether the ball is flying or lying on the ground and it would be convenient to have this variable in the standard result.
• The algebraic variables of the second and third item above, could also be computed after the simulation run, as a postprocessing step. There is then the complication, that DifferentialEquations need to store not only the continuous states of the integrator, but also all discrete states of the model function (so variables that change their values only at event instants and are constant during continuous integration; if the continuous functions depend on such variables, the discrete states must be stored in the result as well in order that a post-processing function can correctly compute all desired variables). Mathematically: x_discrete = f1(previous(x_discrete), previous(x)); dx/dt = f2(x,t,x_discrete), where previous(..) is the value of a variable from the last event and x_discrete is a variable that is only computed at event instants.

In order to make the whole thing not too complicated, it might be sufficient in a first step, to just provide support for algebraic variables of 1. (and then the user has to somehow handle cases 2./3. with this basic support).

[ChrisRackauckas]

How are clocked variables stored in Modelica? Does that not disrupt SIMD?

Out of curiosity, what of this is not possible with a DEDataArray or some saving discrete call back?

[MartinOtter]

How are clocked variables stored in Modelica?

To my understanding, most Modelica simulators just store the clock variable at every (continuous-time) instant, even if the value does not change. It would be better to have different time axes for that (one time axis for continuous-time variable and then one time-axis for every clock, and then clocked variables are only stored with the time-axis of the respective clock).

Does that not disrupt SIMD?

This I do not know.

Out of curiosity, what of this is not possible with a DEDataArray or some saving discrete call back?

I think all of this can be implemented by the modeler with the current version of DifferentialEquations.jl (by implementing an own result structure). It would be just a convenience feature in order that the nice result/plot features of DifferentialEquations could be used in such a case as well.

[ChrisRackauckas]

To my understanding, most Modelica simulators just store the clock variable at every (continuous-time) instant, even if the value does not change. It would be better to have different time axes for that (one time axis for continuous-time variable and then one time-axis for every clock, and then clocked variables are only stored with the time-axis of the respective clock).

But the integrator still acts on vectors of numbers, not a heterogeneous vector?

I think all of this can be implemented by the modeler with the current version of DifferentialEquations.jl (by implementing an own result structure). It would be just a convenience feature in order that the nice result/plot features of DifferentialEquations could be used in such a case as well.

I think the main thing for DifferentialEquations.jl is that it should be efficient yet flexible and serve as the core to modeling tools, while not forcing the user down some DSL-like path so that other front-ends can be created.

The clock variables would require far too much upfront information and I don't it would not be type-stable. The type-stable choice is to either completely specialize the clock variables out of the array into a special array type, which would require user-input as to the discontinuous variables and break SIMD optimizations, or would require making the state silently a vector of tuples which would double all saves. The cheapest option seems to me to be just saving twice at the discontinuity (if you do want to save and interpolate properly, otherwise you just turn off the saving). Or just have the user use the integrator interface and like Sundials let them choose exactly how to save if it's model-specific.

For the algebraic variables, I think the right thing to do may be to create a nice plot recipe on DEDataArray or make a fully broadcasting array type (like Flatten.jl). Any features missing from that? What issues did you have using DEDataArray types for this?

[MartinOtter]

But the integrator still acts on vectors of numbers, not a heterogeneous vector?

Yes (continous real variables can be arrays, e.g., the absolute vector to the center of mass of a body. However, the generated residue-function of Modia maps all datastructures to vectors and returns a Vector{Float64} for the residues (to be used by the integrator) and a Vector{Float64} for continuous algebraic variables (to be stored in the result data structure)).

The clock variables would require far too much upfront information and I don't it would not be type-stable.

Yes, I completely agree that clock variables should not be handled specially (and is not needed from Modia point of view; especially since this is only an improvement to reduce the size of the result data).

For the algebraic variables, I think the right thing to do may be to create a nice plot recipe on DEDataArray or make a fully broadcasting array type (like Flatten.jl). Any features missing from that? What issues did you have using DEDataArray types for this?

I had a look at the DifferentialEquations.jl documentation and DEDDataArrays seems to be perfect for variables that change only at events (but not for continuous algebraic variables).

[ChrisRackauckas]

How are you handling continuous algebraic variables when using Sundials.jl? Do you only allow them in a DAE?

[MartinOtter]

How are you handling continuous algebraic variables when using Sundials.jl? Do you only allow them in a DAE?

They are also allowed in an ODE: This just means that an ODE integrator can be used and whenever results are stored, then additionally the algebraic variables are stored too.

With a DAE integrator, it is additionally possible to handle algebraic variables that cannot be explicitely solved for.

So to summarize: In ModiaMath the continuous variables are split into two classes:

• The explicit solvable algebraic variables are just stored whenever the states are stored (otherwise, the algebraic variables are just internal to the model).
• The implicit algebraic variables (variables that cannot explicitely solved for) are included in vector x, if the DAE is described as f(dx/dt, x, t) = 0.

[ChrisRackauckas]

So then the saving value isn't defined by the input value? Is this better done via the SavingCallback? I assume you have to add a hook over the interpolation as well to make it also compute the explicit algebraic variables after each interpolation?

[MartinOtter]

So then the saving value isn't defined by the input value? Is this better done via the SavingCallback? I assume you have to add a hook over the interpolation as well to make it also compute the explicit algebraic variables after each interpolation?

Sorry, I overlooked this (or is it new?): Yes, this is completely sufficient.

O.k., you are right, I have to look at the interpolation as well.

Thanks for the discussion. I understand now that I can implement everything regarding algebraic variables in Modia with the available functionality of DifferentialEquations.

One minor question: How to associate symbols/names to the variables that are defined with SavingCallback.

[asprionj]

1. Some algebraic variables are needed to compute the state derivatives and need therefore be computed whenever derivatives need to be computed. Mathematically: w = f1(x,t); dx/dt=f2(x,w,t).

What's the solution for those now? Just re-calculate them in a SavingCallback? In the model equations I encountered over the years, there was often the case that multiple w's (that are then used in the expression for the state derivatives) are i) of interest (i.e. need to be available for analysis after the simulation) and ii) make up almost all the calculation time of the model function.

Convenience is not the point here; I'd just have to define a method to calculate these intermediate quantities and then call it from the actual model function and in the SavingCallback. My concerns lie with the performance of the simulation. In my use cases, there will often my a vast amount of simulations (e.g. to test some change on the entire set of possible scenarios). So I don't want to calculate a large part of the model again at each saving point.

To assess this in more detail: What usually is the fraction of saving points compared to overall model evaluations for the different types of ODE solvers?

[ChrisRackauckas]

What usually is the fraction of saving points compared to overall model evaluations for the different types of ODE solvers?

That completely depends on the application. In some applications I have seen far more saving points that time steps (DynamicalSystems.jl has a lot of examples of this). But PDE applications might save 0 intermediate points and just want the final solution. Then there's everything in the middle.

Convenience is not the point here; I'd just have to define a method to calculate these intermediate quantities and then call it from the actual model function and in the SavingCallback. My concerns lie with the performance of the simulation. In my use cases, there will often my a vast amount of simulations (e.g. to test some change on the entire set of possible scenarios). So I don't want to calculate a large part of the model again at each saving point.

The point that you're saving at likely has never calculated this value before, so there's no way it could've been cached. Saving occurs via an interpolation, where the interpolation is done so that the steps are not constrained by where the user wants to save (this is standard in ODE solver software and not DifferentialEquations.jl specific). The reason this is done is because then you maximize the number of time steps, and thus decrease the total number of f calls and increase the performance. Each interpolation does not use f, only a linear combination of precomputed derivatives, so those would avoid ever having to build the algebraic variables. So they would only ever be instantiated at the solver steps (which are not the save points, specifically to reduce the number of instantiations) and whenever the user wanted to use them in post-solution processing, which is why it's just as efficient to hand it off to the user (or the SavingCallback, which is essentially the user just asking for more save values).

[asprionj]

Ok, thanks for the explanation. Everything getting clearer. So there will be no to negligibly small overhead (if accidentally a solver point and a save point coincide).

[ChrisRackauckas]

Exactly. If the save and the solver point collide, then you get one extra algebraic variable instantiation. A 5th order RK method uses 7 f calls, so that's 1/8 are unnecessary. But in most other cases, you will step as large as possible, and each time you cut a step by allowing interpolation to be the saving routine, you cut 7 instantiations. The worst case then is that you do step to every save point, in which case you add 1/8 to the cost, which still isn't terrible.