Handling ODE with constraints (drift elimination for index reduced DAE)

[masoud-najafi]

Applying the index-reduction algorithm (Pantelides,1988) to a high index DAE to convert it to an ODE, introduces hidden constraints. Since in FMI for ModelExchange, only ODE are supported, the index reduction algorithm should reduce the index to zero (ODE) and all hidden constraints may be ignored. Ignoring the hidden constraints, and integrating solely using the ODE causes the systems suffer from drift in hidden constraints on the states. In order to reduce the drift in states, we should take into account the constraints that have been ignored when integrating the ODE by ODE solvers. A very straightforward solution to bring back the constraints, is the state projection method. In a simple projection process, after each completed integration step, the FMU computes the changes required for each state variable so that the current values of the system lie on the constraint manifold. Ideally this should be computed as the minimum (or near minimum) change to accomplish this, as the constraint problem is typically under-determined, so many solutions are possible. The following API, computes the amount of projection of the current solution to bring the states back on the constraint manifold.

fmi3Status fmi3Projection(fmi3Component c, fmi3Float64 P[], size_t np, fmi3Float64 projectionTolerance)

The projection is applied to the current state until the states satisfy the constraints to within projectionTolerance. The length of array P is numberOfContinuousStates. This function returns fmi2Discard if it is unable to project onto the manifold, otherwise it returns fmi2OK.

On the XML side, a new attribute for FMU for ModelExchange is needed to indicate that the API fmi3Projection is supported.

If the importer does not support this API, or does not call it to readjust the state values after each completed integration step, the FMU works, but since the hidden constraints are not verified, drift may appear in solutions.

Reference: "Simulation of high-index DAEs and ODEs with constraints in FMI", Masoud Najafi, Modelica conference, 2018, Tokyo

[chrbertsch]

This change is "smaller" than the proposal #144 and the FCP support for DAEs (possibly by introducing a new FMU kind), right? And could be introduced in FMI3.1, right?

Is this similar to https://de.mathworks.com/help/simulink/sfg/mdlprojection.html for Simulink S-Functions?

[masoud-najafi]

Yes, I made it smaller and much easier to implement from importer's point of view. As far as I know, it should be possible to implement in 3.1.

[KarlWernersson]

The problem with introducing something like this is that it would enforce all importers to support it, as this is not optional for an FMU that uses this functionality.

I also don't see the need for ODE's as this can be solved perfectly well internally by the FMU with other methods like dummy derivatives.

https://www.researchgate.net/publication/235324214_Index_Reduction_in_Differential-Algebraic_Equations_Using_Dummy_Derivatives

[chrbertsch]

Could such a new capability with a new API interface function be realized in a layered standard? (in the sense of https://github.com/modelica/fmi-standard/issues/1643) Then perhaps not with the name "fmi3Projection", but with different name? (how could this look like?)

[masoud-najafi]

If You handle constraints and do not throw them away and ensure that you do not have drift, you need not use this API. But In some cases there are some extra constraints that purely algebraic and additional to the ODE, such as energy conservation equations. In current FMI you cannot exploit these extra constraints and drift will be present. This is a capability flag and can be ignored by any importer not supporting it. The naming may change.

[Maplesoft-fmigroup]

The projection method also has advantages from an FMU implementation standpoint. Unlike the dummy derivative method, a partitioning of the constraint set into any number of well conditioned sub-problems is not required. Instead the constraint set is viewed as an underdetermined subproblem on the states alone, and can be tackled purely numerically by the fmi3Projection implementation. This makes it more straightforward to construct an FMU for models with constraints.

Tools which utilize projection as their primary approach for their solvers may have difficulty adapting this to a dummy derivative method.

Use of fmi3Projection by an importer should not be too difficult, as it is primarily a matter of calling the projection upon completion of each integrator step.

If this additional function does pose difficulties, then an alternative approach could be used. The projection could be done automatically in the call to fmi3CompletedIntegratorStep, with a capability flag that indicates that this call can change state values, and that they need to be re-queried after this call. The fmi3Projection approach has advantages over this. One being that the calculation (projection) is directly exposed and can be measured by the importer if so desired.

[KarlWernersson]

But is this at all a backwards compatible change? An Importer cannot ignore these function calls this as the FMU would otherwise be faulty. Isn't it better to work on a proper DAE solution?

[Maplesoft-fmigroup]

A rough outline of how this could be made backward compatible as follows:

1. Add an attribute to the model description: applyProjection with default value "completedStep" which tells when projection should be applied. By default projection is then applied in the call to fmi3CompletedIntegratorStep

2. One other choice of applyProjection is "manual" in which case the fmi3Projection routine must be explicitly called by the importer, or not at all if they wish to ignore the constraints.

Moving forward, newer 3.1 FMUs could use applyProjection="none" to indicate that the FMU does not use projection.

This allows FMUs that require projection to run properly in 3.0 with no changes for the importer. If they require additional control of when projection is applied (for 3.1) then if applyProjection is specified as "completedStep" in the modelDescription, they can change it to "manual" (and if specified as "none" then no projection is in use).

The applyProjection option could be made into a pair of true/false flags to more closely match the existing capability flags instead.

[chrbertsch]

For the first case: why does the projection have to be exposed an cannot just be done silently within the FMU at the call to fmiCompletedIntegrorStep, with a change to the existing standards, that continuous states may be changed by a call of this function?

[Maplesoft-fmigroup]

I have only proposed exposure of the fmi3Projection interface for 3.1 to allow importers finer control over the fmu behavior.

For example, for a well behaved system it may not be necessary to project at every step.

In addition, if the import tool also uses projection, this step could be coupled.

[TorstenBlochwitz]

For the first case: why does the projection have to be exposed an cannot just be done silently within the FMU at the call to fmiCompletedIntegrorStep, with a change to the existing standards, that continuous states may be changed by a call of this function?

According to my memory we discussed this problem some time ago and decided not to introduce an extension since the projection can be done by the the FMU in Event Mode after requesting a step-event in fmi*ComletedIntegratorStep.

[masoud-najafi]

@TorstenBlochwitz The change in Event-mode can already be done in FMI-2. But that slows down the simulation a lot. Furthermore, for the constraint correction purpose, the change in continuous-time states, can only be done the numerical solvers. Just consider a variable step solver that needs to keep the Nordsieck polynomial of the solutions and provide the output dense (interpolation). Also it is possible that the numerical solver requires these corrections between steps, not when the step has already finished.

[Maplesoft-fmigroup]

As for the existing mechanism, in section 3.2.1 of the FMI3 specification the fmi3CompletedIntegratorStep function is described in detail. A couple examples of its usage are provided.

I would suggest that we add the following between the 1st and 2nd bullet points in the examples.

• Handle DAE projection by projecting the states onto the constraint manifold. Note that this occurrence is signalled to the importer by flagging that a state event has occurred.

This will formalize the current agreement on how projection should be handled, and provide a path for other FMU creators who use projection in their tools.

[andreas-junghanns]

I would suggest that we add the following between the 1st and 2nd bullet points in the examples.

* Handle DAE projection by projecting the states onto the constraint
  manifold. Note that this occurrence is signalled to the importer
  by flagging that a state event has occurred.

Could you provide a PR for this change?

[masoud-najafi]

I have no objection with adding this bullet to the FMI-3.0, but I would like to highlight that this is not the intention of this issue. As I explained above, we do not want to generate event each time the projection is applied.

[TorstenBlochwitz]

I have no objection with adding this bullet to the FMI-3.0, but I would like to highlight that this is not the intention of this issue. As I explained above, we do not want to generate event each time the projection is applied.

But your proposal generates an event for the ODE-Solver since the continuous states change discontinuously! Hence, I do not really see the benefit of the new mechanism.

[masoud-najafi]

I have no objection with adding this bullet to the FMI-3.0, but I would like to highlight that this is not the intention of this issue. As I explained above, we do not want to generate event each time the projection is applied.

But your proposal generates an event for the ODE-Solver since the continuous states change discontinuously! Hence, I do not really see the benefit of the new mechanism.

There are several ODE solvers on the market that can handle this situations. For example the latest sundials-CVODE (old CPODE) can handle this situation even for multi-step case. For single step solvers it is easy and can be implemented easily without need for generating any event.