Specify which flavour(s) of DAE residuals are available in the FMU-DAE

[casella]

The collection of use cases that we have gathered corresponds to DAEs with different structural properties. The FMI-DAE format should be agnostic from this point of view and support as many as it is reasonably feasible. However, from a user's point of view, an FMU is only usable if it provides residuals with some structural properties that match the requirements of the solver one wants to use. So, a specific solver should be able to enquire the FMU to know if residuals corresponding to a formulation with certain properties are available, and if they are, actually get them.

This applies to FMUs that were generated with only one specific formulation in mind, as well as to FMUs that may contain multiple different formulations and can make a specific one available on demand.

For starters, here is a preliminary, incomplete and not particularly orderly list of such properties

• DAEs for collocation-based optimal control must be index-1, contain no discrete variables, no events, and only continuous and twice-differentiable constraint equations, so that Jacobians and Hessians are available.
• DAEs for simulation using certain solvers must be index-1, meaning that the tool generating them should take care of performing symbolic index reduction before generating the code for the residuals; however, they could contain discrete variables, discontinuous expressions and state or time events.
  • @MartinOtter mentioned that Brenan and Petzold in their book describe some cases of DAE systems that are index-1 but yet problematic for BDF solvers like DASSL. We may want to categorize them explicitly.
  • We may also want to identify semi-explicit index-1 DAEs, w.r.t. generic fully implicit DAEs
• In some application contexts, having index-2 or index-3 DAEs is instead fine, if you have solvers that can handle this case numerically.
• In some cases the FMU is required to provide code for consistent initialization and return initial conditions that can safely assumed to be consistent, so they can be used to initialize solvers that need this property to work correctly,
• In other cases, the FMU just provides the raw residual equations and it's up to the solver using the FMU to compute consistent initial conditions, using suitable algorithms.
• In some cases, the DAE residuals correspond to each and every differential and algebraic equation of the original system model, possibly after symbolic index reduction and trivial simplifications like alias elimination.
• In some other cases, the FMU could just expose to an external DAE solver the residuals of differential equations and implicit algebraic equations, whereas variables and equations that can be solved explicitly after causalisation are computed internally in the FMU code and not be exposed to the external solver. In this case, additional functions should be made available to compute all such variables, or a selection of them which are relevant as outputs, once a step has been accepted.

We should complete this list based on all the collected use cases, and then try to come up with a systematic classification of properties that can be declared in the FMU manifest file, so a solver can enquire the FMU, make sure that it can provide residuals with certain structural properties, and eventually get them by calling the appropriate functions. Specific FMI-DAE functions may need to be defined to support some of these structural options.

This classification will also be useful to formalize what kind of requirements could be asked to an FMU-generating tool when generting the FMU-DAE code for a specific system.

[jaeandersson]

@ALL: Feel free to add your use case

[casella]

There's been some discussion about FMU-DAEs that provide representations of the same system with different indeces. I'd like to post a trivial index-2 example to make sure we are all on the same page, at least with the terminology. Never mind that this system can be easily turned into index-1 by trivial alias elimination, that's irrelevant for the purposes of this post.

Original formulation

$C_1\frac{dx_1}{dt} = -y$
$C_2\frac{dx_2}{dt} = y - 1$
$x_1 = x_2$

Case 1. An FMU-DAE could just provide the system as-is:

$r_1 = C_1\frac{dx_1}{dt} + y$
$r_2 = C_2\frac{dx_2}{dt} - y + 1$
$r_3 = x_1 - x_2$

without any structural analysis. It would be entirely up to the solver to figure out how to solve it.

Case 2. An FMU-DAE could provide the system as-is:

$r_1 = C_1\frac{dx_1}{dt} + y$
$r_2 = C_2\frac{dx_2}{dt} - y + 1$
$r_3 = x_1 - x_2$

with the additional indication that the system is index-2. This tells the importing tool that it needs a solver that can handle index-2 systems.

Case 3. An FMU-DAE could provide the system as-is:

$r_1 = C_1\frac{dx_1}{dt} + y$
$r_2 = C_2\frac{dx_2}{dt} - y + 1$
$r_3 = x_1 - x_2$

with the additional indication that $r_3$ is a state constraint equation that should be differentiated in order to get an index-1 system. This gives additional information to the solver regarding the "problematic" algebraic equations.

Case 4. An FMU-DAE could provide the system as index-1 implicit DAE

$r_1 = C_1\frac{dx_1}{dt} + y$
$r_2 = C_2 x^{der}_2 - y + 1$
$r_3 = x_1 - x_2$
$r_3 = \frac{dx_1}{dt} - x^{der}_2$

This is obtained by applying Pantelides' algorithm to figure out which equations must be differentiated, adding them to the system, and then applying the dummy derivatives algorithm to turn some states into dummy states.

As I understand, all these cases (and possibly more, e.g. involving semi-explicit formulations, or models with index > 2) can be useful in some context. Hence, the new standard should be able to provide all of them, each with clearly specified properties.

Of course this doesn't mean that an FMU-DAE must provide all this formulations, that it is up to the exporting tool and to how it is configured by the user generating the FMU. But it should support all of them and describe them in the manifest file, so that a tool can pick the best formulation for a given task, or complain that the required formulation is missing.

Do I miss something?

@elmir-nahodovic, @HansOlsson I'm not sure I understand the notation of the old slide shown in #22, nor the subsequent comments. Would you mind posting how the simple system of this post would be mapped into that notation? Thanks!

[casella]

Then, we have an additional interesting case, adding invariants to models that are already solvable, to improve the accuracy of the simulation results. For example:

$\frac{dx}{dt} = v$
$m \frac{dv}{dt} = - kx$
$m\frac{v^2}{2} +kx = m\frac{v_0^2}{2} +kx_0$

where the last equation is not strictly necessary to solve the problem, but provides additional constraints on the state variables that are known to hold at all time, e.g. energy conservation. Solvers can use this information to implement projection methods that keep the solution on the manifold defined by these additional constraints.

As I understand, these constraint equations will have to be identified, as they need to be handled in a specific way. Of course this issue is completely orthogonal to the index issue.

[MartinOtter]

General DAE solvers such as DASSL or Sundials.IDA cannot solve every index-1 DAE. The limitations are, for example, described in pages 51-54 of the book:

• K.E. Brenan, S.L. Campbell, and L.R. Petzold (1996): Numerical Solution of Initial Value Problems in Differential-Algebraic Equations. SIAM.

In particular, several DAE structures are shown in this book where BDF methods of order k with fixed or variable step-size h converge with O(h^k). The following special index 1 DAE does not fall into this proof:

(1)  
a11*der(x1) + a12*x2 = b1
a21*der(x1) + a22*x2 = b2
 
where matrix [a11, a12; a21, a22] is regular

When transforming this DAE to an ODE, everything is uncritical, because just the linear, regular algebraic equation system

(2)
[a11 a12;
 a21 a22]*[der(x1),x2] = [b1;b2]

has to be solved for [der(x1), x2] and der(x1) is returned to the ODE solver that solves der(x1) = f(x1).

Assume implicit Euler is used to discretize the residual form:

(3)  
residual = [a11*der(x1) + a12*x2 - b1;
                  a21*der(x1) + a22*x2 - b2]

This results in the following Jacobian ("h" is the step size):

(4)
J = [a11/h   a12;
      a21/h   a22]

Its simpler to analyze h*J:

(5)   
lim_h->0:  h*J = [a11  a12*h;
                            a21  a22*h]
                        = [a11 0;
                             a21 0]

So, for h->0 the Jacobian J becomes singular. Variable step integrators reduce the step size h if the error estimate is too large. This works under the assumption that reducing the step-size makes the results more accurate because in the limit for h->0, the mathematical exact solution is achieved. For the discretized system above, reducing the step-size will make the solution worse, if h is too small. So a step-size selection algorithm will have difficulties to decide what to do.

One way to fix the problem above is to introduce a transformation proposed by Gear:

(6)
residual = [a11*der(x1) + a12*der(x2_int) - b1;
                  a21*der(x1) + a22*der(x2_int) - b2]
with x2_int(0) = 0 (and x2 = der(x2_int))

Discretizing der(x1) and der(x2_int) with implicit Euler gives a Jacobian of (6) in x1 and x2_int that is regular if h->0. In order to be able to do this, one has to provide the additional information that x2 needs to be specially handled.

Another solution is provided by Matlab where ode15s can solve systems of the form M(y,t)*der(y)=f(y,t) with singular M (see also Shampine et al. (1999): Solving Index-1 DAEs in MATLAB and Simulink, provided this is an index-1 system. This means, not the residual needs to be provided, but the matrix M and the right hand-side vector separately - so additional information.

In the paper Otter/Elmqvist (2017): Transformation of Differential Algebraic Array Equations to Index One Form a quite general approach is proposed, by transforming a DAE to the following special index 1 DAE:

(7)
0 = fd(der(x), x, t) 
0 = fc(x,t)

 where J = [der(fd;der(x)); 
                   der(fc,x)       ] is regular

This index 1 DAE falls into the proof given by Brenan et al. (1996), so BDF methods converge on this DAE.

The standard structural algorithms of Modelica tools (such as the Pantelides algorithm) are used to differentiate equations until the highest derivative equations are structurally regular with respect to the highest derivative unknowns. Under the assumption that these highest derivative equations are also numerically regular with respect to the highest derivative unknowns, the differentiated equations can be transformed to (7) without solving algebraic equations and use generalizations of some proposals from Gear et al. that have been developed for multibody systems. Since no algebraic equations need to be solved, this approach is suited for large systems (with potentially large algebraic systems of equations). Note, multibody systems are contained as a special case, see equations (14,15) in Otter/Elmqvist (2017).

One idea could be to use (7) as a DAE form supported by fmi-ls-dae. The benefit would be that a large class of systems are described in one DAE format. Basically this means, that the specification of fmi-ls-dae defines that DAEs in residual form (7) are supported. If dummy derivatives, such as x2=der(xs_int) need to be introduced, these equations can be provided as algebraic output variables (or they are specially defined in the xml file, to handle them specially, e.g. always have initial values of zero). Therefore, equations need not be marked as having an index or be the derivative of another equation (could be still additional useful information that is provided). Examples could be provided to show how to transform special DAEs to this form, e.g., multibody equations. Initialization is not yet clear to me. Probably, for initialization the information about which equation is the derivative of another equation would be useful.

The alternative is to provide several different DAE structures that are supported by fmi-ls-dae, such as (a) semi-explicit DAEs, (b) systems with mass-matrix as for ode1d, (c) multibody systems, ...