Switching between ODE and DAE behavior is controlled by a structural parameter called enableDAEModeParameter.

A structural parameter (causality = structuralParameter) may only be set during Configuration Mode or Reconfiguration Mode. This constraint ensures that the importer can re-examine the model structure from the manifest and reinitialize its solver before any computation takes place.

The variable enableDAEModeParameter is declared in the modelDescription.xml as a Boolean under <ModelVariables>. The name and valueReference is assigned freely by the FMU exporter; the importer locates the parameter by its valueReference declared in the <EnableDAEParameter> element of fmi-ls-manifest.xml. The description attribute is optional.

An example of how the enableDAEModeParameter is declared:

Because enableDAEModeParameter is a structural parameter, the importer must change it inside a configuration or reconfiguration scope.

Using the semi-explicit formulation defined in Section 2.2, the differential equations are

with

and the algebraic equation is

with

The initial values in this notation are

The Jacobian of the residual equation with respect to the algebraic variable is

which is non-singular. Hence, the algebraic equation can be locally solved for \(\mathbf{a}_c\), and the Robertson problem is an index-1 DAE in the representation supported by this layered standard.

Assume the FMU has the following variables:

Here, \(y_1\) and \(y_2\) are the components of \(\mathbf{x}_c\), while \(y_3\) is the single component of \(\mathbf{a}_c\). The residual \(\mathbf{f}_{\mathit{res}}\) is represented by the variable \(r\).

The DAE-specific manifest elements are:

Because the residual equation is locally solvable for \(\mathbf{a}_c\), the algebraic variable can be eliminated to obtain an explicit ODE.

For fixed \(\mathbf{x}_c\), \(\mathbf{u}_c\), \(\mathbf{p}\), and \(t\), consider the mapping

Its inverse with respect to \(\mathbf{a}_c\), evaluated at residual zero, gives

For the Robertson problem,

and therefore

Substituting this inverse into the differential equations gives the explicit ODE

For the Robertson problem, this becomes

The quantity corresponding to the algebraic variable in the DAE formulation can still be computed in the ODE representation as

The model is then represented as an ODE Model Exchange FMU. Using the value references of the Robertson FMU, the relevant variables are:

In the ODE representation, \(y_3\) is computed from \(y_1\) and \(y_2\) and is not listed as an <AlgebraicVariable>.

The corresponding regular ODE ModelStructure available in the modelDescription.xml is:

The generic fully-implicit DAE, or implicit DAE, formulation:

where \(z(t)\) represents state and additional variables and \(t\) is time.

A similar and sometimes more useful representation:

with dynamic variables \(x(t)\), algebraic variables \(y(t)\) and time \(t\).

Note that \(x(t)\) is not necessarily the state variable, more on this in Section 4.1.7.2.

The semi-explicit index 1, also called Hessenberg index 1, DAE representation

with \(g\) solvable for \(y\), \(det \left( \frac{\partial}{\partial y}g\right) \neq 0\).

Here \(f\) is called the differential equations and \(g\) the algebraic equations.

A fully implicit DAE

can always be rewritten as a semi-explicit DAE by introducing a new algebraic variable \(x^{\prime}\)

The semi-explicit DAE representation

with \(g^1\) solvable for \(y\), \(det \left( \frac{\partial}{\partial y}g^1\right) \neq 0\) and \(\left\{ \frac{d}{dt} g_i^k \right\} \subseteq \{ g_i^{k-1}\}, \; k > 1\).

where mass matrix \(M\) is singular.

An equivalent explicit formulation:

See [cui2020mass].

where matrix \(E\) is singular.

There are multiple, sometimes closely related, definitions of indices for DAE systems. Usually these indices are a measure of how far from an ODE the DAE is.

A simple initial-value problem

with variables \(x,y,z \in \mathbb{R}\) and start condition

at time \(t_0\). For this specific example a start condition \(x_0\) is sufficient to calculate consistent start values. The independent variable (usually time) \(t\) is omitted for simplicity.

The Modelica model Modelica.Electrical.Analog.Examples.CauerLowPassAnalog from the Modelica Standard Library ^[3]

has structural singularities.

To produce a complete and consistent set of initial conditions for the DAE, exporting tools usually use a symbolic index reduction algorithm.

The model has nine state candidates

and four constraint equations

Modelica tools that want to export an ODE ModelExchange FMU need to reduce the perturbation index to index 1 to end up with an ODE representation of the DAE system.

One option to achieve this is the so called dummy derivative method. From the set of state candidates, five state variables are chosen (e.g. C1.v, C3.v, C5.v, L1.i, L2.i) and the constraint equations

are differentiated

which results in four additional dummy states ($DER.L1.v, $DER.L2.v, $DER.C2.v, $DER.C4.v). Here the Modelica keyword stateSelect influences how variables are picked.

For simulation models representing power grids on a national or larger scale, traditional ODE solvers could struggle with convergence and performance. Also, the time needed to perform symbolic transformation from DAE to ODE is growing with the number of equations.

See for example Modelica models ScalableTestSuite.Electrical.TransmissionLine.ScaledExperiments.TransmissionLineModelica_N_10 to ScalableTestSuite.Electrical.TransmissionLine.ScaledExperiments.TransmissionLineModelica_N_1280 from ScalableTestSuite ^[4].

The simulation was performed using OpenModelica v1.27-dev-92 with the following settings:

Since OpenModelica does not support FMI 3.0 or FMI-LS-DAE at the time of writing, the export and simulation times were measured using OpenModelica’s default C export and simulation.

When a mathematical model is derived for a physical system, the fundamental conservation laws (such as energy or momentum) or geometric properties are directly embedded into the differential equations. These conservation laws define the system’s invariants (quantities that remain constant over time). Examples include total energy in a closed system, total mass in a chemical reaction, or angular momentum in a gravitational field.

The analytical solution to the differential equation respects these invariants perfectly. However, when numerical solvers (like Euler’s method, Runge-Kutta, or BDF methods) are used to simulate the differential equations, numerical drift is introduced.

Numerical solvers solve the ODEs/DAEs by approximating the continuous-time variable with a series of discrete steps. At each step, they calculate an update based on the current state, inevitably accumulating a small, local truncation error. While high-order methods can make this error very small for the state variables themselves, they are generally not designed to preserve arbitrary invariants of the system.

This means that even if the initial condition perfectly satisfies an invariant (e.g., total energy = \(E_0\)), after hundreds or thousands of time steps, the computed solution may slowly but surely "drift" away from the manifold defined by that invariant. The solver is effectively solving a nearby differential equation, one whose true solution does not possess the same conserved quantities.

As a consequence of numerical drift, the simulation can produce results that violate fundamental physical laws. A satellite’s orbit might slowly decay or escape when it should be stable. In molecular dynamics, the temperature (related to kinetic energy) might artificially rise or fall, rendering the simulation useless. Also, small drifts can accumulate, leading to a complete breakdown of the solution over long time horizons. The numerical solution may blow up or behave chaotically even when the true solution is well-behaved.

Consider a physical system consisting of a mass \(m\) attached to a spring with a spring constant \(k\). When the mass is displaced by a distance \(x\) from its equilibrium position, the spring exerts a restoring force given by Hooke’s law:

Applying Newton’s second law (\(m\ddot{x} = F\)) gives:

To solve this numerically, we convert it into a system of two first-order ODEs by introducing velocity (\(v = \dot{x}\)):

The total mechanical energy \(I(x, v)\) (the Hamiltonian) of this closed system is the sum of kinetic and potential energy, which must remain strictly constant:

In the phase space of \((x, v)\), the analytical solution traces out a perfect, closed ellipse, forever staying on this fixed physical "energy manifold".

Now, let’s simulate this system using the explicit Euler method with a time step \(h\):

If we substitute this numerical update step back into our energy invariant function \(I(x, v)\), we find:

Expanding this expression reveals:

The energy is not conserved. Because the error term on the right is strictly positive (\(h^2 > 0\), \(k > 0\), \(m > 0\)), the total calculated energy artificially increases at every single time step. Instead of tracing a stable, closed ellipse, the numerical solution spirals outward indefinitely. The solver has no built-in knowledge of the system’s physics and lacks a mechanism to constrain its steps to the true energy manifold.

The constraints are not always related to conservation laws. If the index of a DAE is reduced to zero to obtain an ODE, some equations may be omitted. In order to ensure that these equations are taken into account and used to avoid drift, these equations should be exploited.

The connection between an ODE solved with invariants and a DAE is fundamental, particularly when dealing with index reduction. Many physical systems, such as multibody dynamics, are naturally modeled as high-index DAEs. Reducing the DAE index to zero (transforming it into an ODE) simplifies the solution process but removes the original algebraic constraints from the solved equations. These constraints, which define the system’s geometric invariants (e.g., the set of physically allowed positions), are then subject to numerical drift.

After sufficient differentiation, we often obtain a pure ODE \(\dot{y}=F(y)\) for which the original algebraic constraint \(g(y)=0\) is no longer part of the differential equations. Instead, it becomes an invariant:

Since the consistent initial condition satisfies \(g(y(0))=0\), the reduced ODE should preserve \(g(y(t))=0\). Thus, the algebraic constraints of the original DAE become invariants of the reduced ODE.

Consider the planar pendulum, an index-3 model constrained by the length of the rod. The full index-3 DAE in Cartesian coordinates is:

Index reduction eliminates \(\lambda\) and produces an ODE in \((x, y, u, v)\), but removes the original constraints. The resulting ODE no longer enforces:

If the invariants are ignored, numerical drift in the integration of the ODE alone causes, for example, the pendulum to stretch or shrink over time. To prevent numerical drift, both constraints must be enforced via projection after each integration step (or every few steps).

This is where invariant and manifold checking (and correction) becomes essential. To preserve physical fidelity, simulations must actively monitor and correct for drift using one of two strategies.