DAE

ModiaMath.DAE — Module.

module ModiaMath.DAE

Interface between the ModiaMath.SimulationEngine and the index 1 DAE model. A DAE model is a struct that has a required field simulationState::ModiaMath.SimulationState in which the main properties of the DAE model are reported to the simulation engine:

# DAE model ModelName
mutable struct ModelName <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # other definitions (e.g. parameters of model)
end

The following functions can be called in the DAE model to inquire information about the simulation state:

The following functions can be called in the DAE model to set properties in the simulation engine:

The following functions can be either called in the DAE model or they can be called on a simulation model (before or after ModiaMath.simulate!(simulationModel, ...) is called).

Main developer

Martin Otter, DLR - Institute of System Dynamics and Control

ModiaMath.DAE.EventRestart — Type.

@enum EventRestart NoRestart Restart FullRestart Terminate

Define how to continue or restart integration after an event. Usually, Restart should be used. Only in special cases, the other flags are useful.

NoRestart, continue integration without restarting the integrator
Restart, restart integrator
FullRestart, restart integrator and optionally exchange simulationState (so dimensions may change)
Terminate, terminate integration

ModiaMath.DAE.SimulationState — Type.

simulationState = SimulationState(name, getModelResidues!, x_start,
                                  getVariableName; kwargs...)

Return a simulationState object that is described by a DAE with one of the supported structures defined with enumeration StructureOfDAE.

A model that shall be simulated with function ModiaMath.simulate!(model, ...) is required to be defined as:

mutable struct ModelName <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # other definitions (e.g. parameters of model)
end

Keyword arguments	defaults
structureOfDAE	DAE_LinearDerivativesAndConstraints (see `StructureOfDAE`)
is_constraint	fill(false, length(x_start))
has_constraintDerivatives	false
nz	0
nw	0
zDir	fill(0, nz)
x_fixed	fill(false, length(x_start))
x_nominal	fill(NaN, length(x_start))
x_errorControl	fill(true, length(x_start))
jac	nothing (not yet supported)
maxSparsity	0.1
hev	1e-8 (only for DAE_NoSpecialStructure)
scaleConstraintsAtEvents	true (not used; only for backwards compatibility)
nc	0 (not used; only for backwards compatibility)
getResultNames	`ModiaMath.getResultNames`
storeResult!	`ModiaMath.storeRawResult!`
getResult	`ModiaMath.getStringDictResult`
defaultTolerance	1e-4
defaultStartTime	0.0
defaultStopTime	1.0
defaultInterval	NaN
defaultMaxStepSize	NaN
defaultMaxNumberOfSteps	missing

Required arguments

name::Union{AbstractString,Symbol}: Name of model
getModelResidues!::Function: Function with arguments (model,t,x,derx,r,w) to compute the residues r and auxiliary variables w from time t, vector x and its time derivative derx.
x_start::Vector{Float64}: Start values of x.
getVariableName::Function=ModiaMath.defaultVariableName: Function that returns the name of a variable, given its type and its index.

Optional (keyword) arguments:

structureOfDAE::StructureOfDAE: Structure of DAE.
is_constraint::Vector{Bool}: = true, residue[i] is a constraint equation fc or its derivative. = false, residue[i] is not a constraint equation. (is only used for structureOfDAE = DAE_LinearDerivativesAndConstraints)
has_constraintDerivatives::Vector{Bool}: if ModiaMath.compute_der_fc returns true and is_constraint[i] = true: = true , residue[i] is the derivative of an fc equation = false, residue[i] is an fc equation (is only used for structureOfDAE = DAE_LinearDerivativesAndConstraints)
nz::Int: Number of event indicators
nw::Int: Number of auxiliary variables (Float64 variables that are additionally computed and stored at communication points, and where start values can be provided for initialization)
zDir::Vector{Int}: Interpretation of event indictors: zDir[i] = 0: Root is reported for both crossing directions, = 1: Root is reported when crossing from negative to positive direction = -1: Root is reported when crossing from positive to negative direction
x_fixed::Vector{Bool}: = true, x_start[i] is fixed during initialization. = false, x_start[i] might be changed, e.g., due to an initial impulse.
x_nominal::Vector{Float64}: Nominal values of x. if x_nominal[i]=NaN no nominal value is defined for x[i] and a nominal value is computed (x_nominal[i] = abs(x_start[i]) > 1e-7 ? abs(x_start[i]) : 1.0).
x_errorControl::Vector{Bool}: = true, the absolute error tolerance is set to 0.1 * relativeTolerance * x_nominal[i]. = false, the absolute error tolerance is switched off (is set to a large value). This is recommended for variables that are basically not limited (for example the angle of a shaft that is permantently rotating in the same direction and therefore becomes larger and larger).
hev::Float64: Stepsize used during initialization and at event restart if structureOfDAE = ModiaMath.DAE_NoSpecialStructure. Otherwise hev is ignored.
scaleConstraintsAtEvents::Bool: Dummy argument. Only kept for backwards compatibility.
nc::Int: Dummy argument. Only kept for backwards compatibility.
jac: Sparse Jacobian datastructure (currently not supported).
maxSparsity::Float64: A sparse Jacobian is only used during simulation if sparseness of jac < maxSparsity (currently not supported)
getResultNames::Function: Function that returns the names of the variables to be stored in the result data structure.
storeResult!::Function: Function that stores the raw results.
getResult::Function: Function that resturns the result data structure after the simulation.
defaultTolerance::Float64: Model specific default relative tolerance, if not redefined in the call to ModiaMath.simulate!.
defaultStartTime::Float64: Model specific default start time in [s], if not redefined in the call to ModiaMath.simulate!.
defaultStopTime::Float64: Model specific default stop time in [s], if not redefined in the call to ModiaMath.simulate!.
defaultInterval::Float64: Model specific default interval in [s], if not redefined in the call to ModiaMath.simulate!. Result data is stored every defaultInterval seconds. If defaultInterval=NaN, the default interval is computed as interval = (stopTime - startTime)/500.0.
defaultMaxStepSize::Float64: Model specific default of the maximum absolute value of the step size. If defaultMaxStepSize=NaN, the default maximum step size of the integrator is used.
defaultMaxNumberOfSteps::Union{Int,Missing}: Model specific default of the maximum number of steps to be taken by the solver in its attempt to reach the next output time. If defaultMaxNumberOfSteps=missing, the default value of the integrator is used (= 500).

ModiaMath.DAE.StructureOfDAE — Type.

@enum StructureOfDAE
      DAE_LinearDerivativesAndConstraints
      DAE_ExplicitDerivatives
      DAE_NoSpecialStructure

Enumeration defining the structure of the DAE of the simulation model. The following DAE structures are supported (function getModelResidues!(model, t, x, derx, r, w; simulation=true) returns the residues r):

ModiaMath.DAE_LinearDerivativesAndConstraints (default)

\[\begin{align} z &= f_z(x,t) \\ 0 = r &= \left[ \begin{array}{l} M_d(x,t,z_i>0) \cdot \dot{x} + b_d(x,t,z_i>0) \;\; (= f_d) \\ f_c(x,t,z_i>0) \end{array} \right] \\ J &= \left[ M_d; \frac{\partial f_c}{\partial x} \right] \; \text{is regular (matrix is invertible)} \end{align}\]

where

\[\lim_{\epsilon \rightarrow 0} x(t_0 - \epsilon) = x_0^{-}\]

Equations $f_d$ are linear in the derivatives $\dot{x}$. It is required that the Jacobian $J$ is regular, that is the DAE has an index 1 (= by differentiating $f_c$ once, the system can be transformed to an ODE).

Equations $z=z(t)$ are zero-crossing functions. Whenever a $z_i(t)$ crosses zero, an event is triggered and simulation is halted. During an event, $z_i > 0$ can change its value. The equations above are solved with a fixed-point iteration scheme (= event iteration) until $z_i > 0$ does not change anymore. Afterwards, integration is restarted and the Boolean variable $z_{pos} = z_{i,ev} > 0$ keeps its value until the next event occurs.

At an event instant some $f_c$ equations might become $f_d$ equations and vice versa. The constraint equations $f_c$ can be at any position of the residue vector r and at an event instant $f_c$ equations might become $f_d$ equations and vice versa. When instantiating a SimulationState, the initial definition of the constraint equations is provided with vector argument is_constraint. Note, it is also possible to define time events, so triggering events at pre-defined time instants, for example to model sampled data systems (see ModiaMath.setNextEvent!).

Initial conditions $x_{ev}^{-}$ must be provided before simulation can start ($x_{ev}^{-} = x_0^{-}$) or at an event restart. They need not to fulfill the constraint equations, so $f_c (x_{ev}^{-},t_{ev} ) \neq 0$ is allowed. If this is the case, initialization/re-initialization will simulate for an infinitesimal small time instant so that $x_{ev}^{-}$ changes discontinuously to $x_{ev}^{+}$ with $f_c (x_{ev}^{+},t_{ev} )=0$. This is performed by analytically integrating over the initial time or the event time and might imply to integrate over Dirac impulses (in case $x$ is discontinuous at this time instant). Under certain conditions a numerical approximation of the mathematical (exact) solution is computed.

The derivative of the constraint equations $f_c$ can be provided at event restart. In this case $\frac{\partial f_c}{\partial x}$ is computed with the provided $\dot{f}_c$, instead of computing it with finite differences (which is numerically less reliable). In both cases $\dot{x}_{ev}^{+}$ is then computed by solving a linear equation system with Jacobian $J$.

ModiaMath.DAE_ExplicitDerivatives

\[\begin{align} z &= f_z(x, t) \\ 0 = r &= - \dot{x} + b_d(x,t,z_i>0) \;\; (= f_d) \end{align}\]

where

\[x(t_0) = x_0\]

This is a special case of the first form, where no constraint equations are present and the derivatives are explicit. This form is also called ODE (Ordinary Differential Equations in state space form). It has the advantage that many ODE integration methods can be used to solve this system.

Initialization and re-initialization is trivial, because $x_{ev}^{+}$ is provided as initial value or at event restart from the model and then:

\[der(x_{ev})^{+} := f_d(0, x_{ev}^{+}, t_{ev})\]

Note, if a Dirac impulse occurs in the model, then this property has to be handled inside the model and the result of analytically integrating over the event instant must be returned as $x_{ev}^{+}$ from the model.

ModiaMath.DAE_NoSpecialStructure

\[\begin{align} z &= f_z(x,t) \\ 0 = r &= f(\dot{x},x,t,z_i>0) \end{align}\]

This is a general DAE. Initialization and re-initialization is performed by using an implicit Euler step. When appropriately scaling r, using a step size that tends to zero and under further assumptions, then this solution can be interpreted as analytically integrating over the time instant. This might mean to integrate over Dirac impulses (in case x is discontinuous at this time instant). Since the selected step size is not close to zero, the implicit Euler step will give a very rough approximation of the analytical integral. A much better approximation is achieved with option DAE_LinearDerivativesAndConstraints above where a step size of zero is used.

This structure is only provided for backwards compatibility. It is numerically not reliable and should no longer be used.

ModiaMath.DAE.change! — Method.

ModiaMath.change!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Trigger an event, whenever crossing > 0 changes from false to true or from true to false. The function returns always false.

ModiaMath.DAE.compute_der_fc — Method.

ModiaMath.compute_der_fc(m::ModiaMath.[AbstractSimulationModel|SimulationState])

If true is returned, return the derivative of the constraint equation in residues r[i]. Otherwise return constraint equations.

ModiaMath.DAE.defaultVariableName — Method.

defaultVariableName(model, vcat, vindex)

Return default names for the variables (e.g. x[1])

ModiaMath.DAE.edge! — Method.

ModiaMath.edge!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Trigger an event, whenever crossing > 0 switches from false to true. The function returns always false.

ModiaMath.DAE.getResultNames — Method.

getResultNames(model)

Return a vector of result names.

ModiaMath.DAE.getStartTime — Method.

ModiaMath.getStartTime(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return startTime of the actual simulation run.

ModiaMath.DAE.getStopTime — Method.

ModiaMath.getStopTime(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return stopTime of the actual simulation run (when the simulation will be terminated).

ModiaMath.DAE.getTime — Method.

ModiaMath.getTime(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return actual simulation time.

ModiaMath.DAE.getTolerance — Method.

ModiaMath.getTolerance(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return (relative) tolerance of the actual simulation run.

ModiaMath.DAE.getVariableName — Function.

getVariableName(model, vcat, vindex, nx=0;
                xNames   =nameVector("x", nx),
                derxNames=fcNameVector("der", xNames),
                wNames   =String[])

Given category vcat and index vindexof category, return the full path name of the respective variable.

ModiaMath.DAE.get_is_constraint — Method.

ModiaMath.get_is_constraint(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return reference to is_constraint vector (this vector can be modified in getModelResidues! if ModiaMath.isEvent returns true.

ModiaMath.DAE.isAfterSimulationStart — Method.

ModiaMath.isAfterSimulationStart(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if after start of simulation and false if during initialization.

ModiaMath.DAE.isEvent — Method.

ModiaMath.isEvent(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if event phase of simulation (including initialization)

ModiaMath.DAE.isInitial — Method.

ModiaMath.isInitial(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if initialization phase of simulation.

ModiaMath.DAE.isStoreResult — Method.

ModiaMath.isStoreResult(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if communication point and variables shall be stored in the result data structure.

ModiaMath.DAE.isTerminal — Method.

ModiaMath.isTerminal(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if terminal phase of simulation.

ModiaMath.DAE.isZeroCrossing — Method.

ModiaMath.isZeroCrossing(m::ModiaMath.[AbstractSimulationModel|SimulationState])

Return true, if simulation model shall compute zero-crossing functions (as required by the integrator).

ModiaMath.DAE.negative! — Method.

ModiaMath.negative!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Return crossing < 0 such that a state event is triggered whenever crossing < 0 changes its value.

ModiaMath.DAE.positive! — Method.

ModiaMath.positive!(sim, nr, crossing, crossingAsString; restart = ModiaMath.Restart)

Return crossing > 0 such that a state event is triggered whenever crossing > 0 changes its value.

Note, a small hysteresis is added to the crossing function during continuous integration, in order to reduce issues during event iterations due to small numerical errors. However, at an event instant, crossing > 0 is returned without any hysteresis.

Arguments

sim::ModiaMath.SimulationState: (Internal) simulation state provided by simulation engine.
nr::Int: (> 0 required) Every call of positive!(..) must be identified by a unique value of nr This value is used as index in a vector that holds the internal memory for crossing > 0.
crossing::Float64: Zero crossing function.
crossingAsString::String: crossing as string representation. This string is used for log messages.
restart::ModiaMath.EventRestart: Restart action after the crossing > 0 event occurred.

Example

function computeVariables!(m::Model, sim::SimulationState)
   ...
   # f = s > 0 ? fMax : 0.0
   m.sPos.value = ModiaMath.positive!(sim, 1, m.s.value, "s")
   m.f.value    = m.sPos.value ? m.fMax.value : 0.0
   ...
end

ModiaMath.DAE.setNextEvent! — Method.

ModiaMath.setNextEvent!(sim, nextEventTime; integratoToEvent=true, restart = ModiaMath.Restart)

At an event instant (isEvent(sim) == true) trigger the next time event.

Arguments

sim::ModiaMath.SimulationState: (Internal) simulation state provided by simulation engine.
nextEventTime::Float64: Time instant of the next time event.
integrateToEvent::Bool: If true, the integrator integrates exactly to this event. If false, the integrator might integrate beyond nextEventTime (so the step size of the integrator is not influenced by this event). This option is only useful, if information is inquired at the event and restart=ModiaMath.NoRestart is used. The benefit is that the integrator is practically not influenced by this event.
restart::ModiaMath.EventRestart: Restart action after the event at which setNextEvent! was called.

ModiaMath.DAE.setNominal! — Method.

ModiaMath.setNominal!(sim::ModiaMath.SimulationState, x_nominal::Vector{Float64})

At initialization provide x_nominal, the nominal values of the x-vector, and store them insim`, the simulation state. These values are used to compute the absolute tolerances of the x-vector for the integrator.

Example

function getModelResidues!(m::Model, t, x, derx, r, w)
   sim = m.simulationState
   if ModiaMath.isInitial(sim)
      ModiaMath.setNominal!(sim, [1.0, 1e-5, 1e8])   # if length(x)=3
   end
   ...
end