Object-Oriented Modeling of Classic Physical Systems using Linear Implicit Equilibrium Dynamics

1. Introduction

Pairs of potential (effort) and flow variables have long been established for the equation-based, object-oriented modeling of physical systems. For instance, Table 1 lists the interface currently in use for the Modelica Standard Library [1] and Table 2 lists the interface typically used for Bond-graphs [2,3].

For bond-graphs, the product of these pairs always represents a flow of energy. Although Modelica is less dogmatic, the flow of energy is still contained in the information of the connection pair.

The sheer size of the Modelica Standard Library, its industry wide adoption and also the use similar interfaces in a wide array of commercial tools provide empiric evidence that these pairs form very useful interfaces. And indeed, these pairs provide an answer for what is necessary for the object-oriented modeling of classic physical systems:

The equations whose solution represent a physical system can be distributed among its components

Thanks to these pairs, we can indeed distribute the equations of physical systems among different components. Unfortunately, the very same pairs are also the root cause for many very persistent problems encountered in the object-oriented modeling of physical systems.

• Irregular systems can be composed and are hard to diagnose.
• Very difficult handling of variable structure systems [4] (change of equations at run-time) due to the high index of the Differential Algebraic Equation (DAE) system.
• Very difficult code generation for large-scale system simulation [5] since often the complete equation system is needed for structural analysis
• Highly complex generation of simulation code requiring compilers that have high development costs. (The decade spanning effort of OpenModelica [6] provides evidence)

All these problems will be diminished (as demonstrated in later chapters), when we fulfill what is sufficient for object-oriented modeling:

Any valid combination of components (under rules of limited complexity) shall have a solution representing a physical system.

This statement centers on the ability to make an a-priori statement on the solvability of the modeled system and sadly (as we will see later), we can cannot fulfill this statement with the presented pairs of effort and flow. Instead, this paper suggests to use triplets as in Table 3.

Before we can aim to find a sufficient form and understand these triplets, we first need to understand what actually led us to the necessary pairs of potential (effort) and flow variables. Hence, we need to revisit the very basic fundamentals of classic physics.

2. The Principle of Stationary Action

The central question is what minimum interface is needed so that we can distribute the equation of systems in classical physics. The natural starting point is hence the most fundamental law forming the basis of all classical physics: the principle of stationary action also known as Hamilton‘s Principle. It is based on the action $S$ for an element moving along a path $q\left(t\right)$:

$$S={\int }_{t_a}^{t_b}L(q\left(t\right),\dot{q}\left(t\right))\mathrm{}dt={\int }_{t_a}^{t_b}T\left(\dot{q}\left(t\right)\right)\mathrm{}-V\left(q\left(t\right)\right)\mathrm{}dt$$

The action is the integral of the Lagrangian $L$ over the timespan $[t_a,t_b]$. The Lagrangian itself can be expressed as the difference between kinetic energy $T$ and potential energy $V$. The principle of stationary action states that:

$$\frac{\partial S}{\partial q\left(t\right)}=0$$

This principle is the basis of classic physics because it is also fulfilled by the underlying quantum physics although in a statistical manner as described by Feynman [7]. Yet for a sufficient large number of quantum events, the expected value will be fulfilled with great precision and an assumption of continuity can be made for macroscopic systems.

The principle of stationary action is here written for a conservative system. It would require extensions for non-conservative systems but there is no need here to complicate matters. In any case, the following literature may be advised [8,9,10,11] (chapter 20). A few points are however often underrepresented and hence deserve discussion here.

The principle of stationary action is based on the Lagrangian $T-V$. We see that a different sign is applied for the kinetic energy $T$ than for all other forms of (potential) energy $V$. The kinetic energy has thus a special role. Not only is its role special but also its form. The gradient of the kinetic energy forms a continuous bijective function for the total domain (whether relativistic or not), meaning that it has a well-defined inverse over the complete domain. This property is important for the discussion of another question: is there always a solution for the principle of stationary action?

Given that we can assume no strict restrictions on the complexity of the potential $V$ (other than having a derivative), it may seem unsolvable in the general case. Fortunately, the plot of Figure 1 indicates how to approach this problem.

This graph is a topological map of the action over time for paths of different constant velocity. We recognize that in this (completely arbitrary) example, the action has become a very complex function after a time-step of one second. However, for a small step of time, the action remains simple and well-natured. The explanation is straight forward: no matter how complex $V$ is shaped, it takes a motion through time and space to gather this complexity. Hence at $t=t_a$, the well-natured kinetic energy dominates the action.

This means we can reliably solve the principle of stationary action if the step in time $t_b-t_a$ is small enough. What is small enough? Under the assumption of continuity, there is only one answer: $dt$. This leads us ultimately to the famous Euler-Lagrange equation and is also the explanation why time derivatives are so useful in physics.

$$\frac{\partial L}{\partial q}-\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}=0$$

Now we see also why it is so important that the partial derivative of the kinetic energy has a special form and a well-defined inverse. This ensures that the second term can always compensate the first in an unambiguous way. Please note that such a general statement on solvability is extremely helpful.

For illustration, let us now apply the Euler-Lagrange to the following example system for the pressure wave of a fluid in a pipe. The derivative of our path: $\dot{q}\left(t\right)$ represents the volume flow $\dot{Q}$. Correspondingly the path $q\left(t\right)$then represents the shift in volume $Q$.

• The kinetic energy is then defined as: $\mathrm{T}=\frac{\mathrm{I}\mathsf{\rho }}{2}{\dot{\mathrm{Q}}}^2$with $\rho$ being the density and the inertance defined by the geometry length $s$ and cross section $A$: $\mathrm{I}=\int \frac{\mathrm{d}\mathrm{s}}{\mathrm{A}}$
• The potential energy can be formulated by $V=\frac{K}{{2Q}_{ref}}{Q}^2$ where $K$ is the bulk modulus and $Q_{ref}$ a reference volume.

The Lagrangian is now defined by:

$$L=\frac{I\rho }{2}{\dot{Q}}^2-\frac{K}{2Q_{ref}}{Q}^2$$

and we can apply the Euler-Lagrange equation which results in:

$$\frac{\partial \left(\frac{I\rho }{2}{\dot{Q}}^2-\frac{K}{2Q_{ref}}{Q}^2\right)}{\partial Q}-\frac{d}{dt}\frac{\partial \left(\frac{I\rho }{2}{\dot{Q}}^2-\frac{\mathrm{K}}{2Q_{ref}}{Q}^2\right)}{\partial \dot{Q}}=0$$

Fortunately, many terms vanish when taking the partial derivative and the resulting equation is a differential equation of second order that expresses the wave in the fluid based on its volumetric shift:

$$Q\frac{K}{Q_{ref}}+I\rho \ddot{Q}=0$$

The elegance of the Euler-Lagrange equation is that it provides a solution of $n$ dimensions with $n$ equations for an underlying variational problem that is actually a $2n$ dimensional problem (please note that one shall treat $\partial q$ and $\partial \dot{q}$as independent variables [11]).

This elegant “compression” is however working against the intent of object-oriented modeling: for examples of higher dimensions, it will be very difficult to redistribute the equations to individual components. Hence the result of Euler-Lagrange is mostly unsuited for object-oriented modeling.

Fortunately, there is an alternative form that can be achieved by the Legendre transformation that transforms the Lagrangian $T-V$ into the Hamiltonian $T+V$. The stationary action principle can now be reformulated for our example. For the Legendre transformation we require the generalized momentums that are defined as

$$p_i=\frac{\partial L}{\partial {\dot{q}}^i}$$

For illustration, let us reiterate our previous example. Based on the former Lagrangian, we can determine the generalized momentum $p$:

$$\mathit{p}=\frac{\partial \mathit{L}}{\partial \dot{\mathit{q}}}=\frac{\partial (\frac{\mathit{I}\mathit{\rho }}{2}{\dot{\mathit{Q}}}^2-\frac{\mathit{K}}{2{\mathit{Q}}_{\mathit{r}\mathit{e}\mathit{f}}}{\mathit{Q}}^2)}{\partial \dot{\mathit{Q}}}=\mathit{I}\mathit{\rho }\dot{\mathit{Q}}$$

Using $p$ to express the kinetic energy $T$ is then applied for composing the Hamiltonian $H=T+V$:

$$H=\frac{1}{2}\frac{p^2}{I\rho }+\frac{K}{2Q_{ref}}{Q}^2$$

Given the Hamiltonian, we can now apply the famous pair of Hamilton equations:

$$\frac{dq}{dt}=\frac{\partial H}{\partial p}$$

$$\frac{dp}{dt}=-\frac{\partial H}{\partial q}$$

In our example, these translate to:

$$\frac{dq}{dt}=\dot{Q}=\frac{\partial \left(\frac{1}{2}\frac{p^2}{I\rho }+\frac{K}{2Q_{ref}}{Q}^2\right)}{\partial p}=\frac{p}{I\rho }$$

$$\frac{dp}{dt}=-\frac{\partial \left(\frac{1}{2}\frac{p^2}{I\rho }+\frac{K}{2Q_{ref}}{Q}^2\right)}{\partial q}=-\frac{Q}{\kappa Q_{ref}}$$

We now get $2n$ first order differential equations, based on generalized position and generalized momentum. These equations express the same system as the previous solution derived by Euler-Lagrange. (simply differentiate the first equation and eradicate $dp/dt$ by substitution). Different from the Euler-Lagrange solution, these equations can be nicely distributed among components, especially if we reformulate them based on a corresponding pair. In our setting it is natural to choose the volume flow $\dot{Q}=dq/dt$ and the pressure $P=dp/dt$ as pair. The corresponding components are presented in Table 4, row 1 and 2.

We can then choose to model a different storage of potential energy such as gravitational pressure. Even though all our derivations were done on conservative systems, the same pair can be used to model non-conservative components such as the friction of a valve. These two additional components are presented in Table 4, row 3 and 4.

Given these 4 components of Table 4 and the pair of effort and flow, we can now assemble a more complex system and for instance model the example of 3 communicating vessels as in Figure 2. Classic object-oriented modeling becomes functional.

In the corresponding simulation results (Figure 3) we can see the oscillation that result from the natural parameters of water. As typical for many systems the exchange between potential and kinetic energy is of high frequency.

3. From Necessary to Sufficient

We have seen that the well-known pairs of potential and flow variables can be derived from the Hamiltonian form of the principle of stationary action (although not always in the exact same way as for the hydraulic example above).

Different form the Lagrangian $T-V$, the kinetic energy loses its special role in the Hamiltonian $T+V$. Treating all energy forms equal may deceivingly suggest that energy flows are all what the modelers needs (Bond graph modelling is a prime example of this fallacy). Yet, whereas each valid Hamiltonian represents a sum of energy terms, not all sums of energy terms represent a valid Hamiltonian. This leads to a serious problem for object-oriented modeling. Whereas the pairs are necessary to redistribute the equations resulting out of a valid Hamiltonian, only very few recombination possible by these pairs form a valid Hamiltonian. This implies that the modeler can formulate non-physical systems too and modelers even do this on purpose (albeit often without being fully aware of the consequences). For instance, let us consider the following recombination of components as depicted in Figure 4.

Here the modeler has now completely removed the inertia and compressibility of water from the model. This means that there is no representation of the kinetic energy anymore. We have removed exactly the component that is required to formulate the principle of stationary action and guarantee its solvability. Instead we now use the pairs to formulate a constraint problem in the domain of energy. Whether this can be solved is in general undefined. And indeed, we have to solve a non-linear equations system in this example to perform the simulation. In this example, the solution could reliably be found during simulation time but in the general case no statement on solvability can be given.

It is however important to recognize that modelers do such things (mostly) not out of malevolence or ignorance but to solve relevant problems. The slight additional cost of solving the non-linear equation system was rewarded in this example by the much higher gain of removing the high frequency from the system enabling a very efficient simulation of the system (albeit significantly altering the trajectory).

Figure 5. Simulation result of the model of Figure 4, showing the 3 volume flows going through the tank’s orifices. In this simplified model the high frequency has been eliminated and the flow rate is now dominated by the resistance.

Figure 5. Simulation result of the model of Figure 4, showing the 3 volume flows going through the tank’s orifices. In this simplified model the high frequency has been eliminated and the flow rate is now dominated by the resistance.

Hence, what would we need to do to ensure solvability and also enable powerful modeling approach? Essentially, we have to fulfill two criteria:

(a)

Ensure that the kinetic energy is always sufficiently represented so that we can guarantee solvability according to the steady-action principle.

(b)

Enable the modeler to reduce the interaction between kinetic and potential energies to key points in the system so that irrelevant high frequency behavior can be effectively suppressed.

For some domains, these two criteria can be met in a methodological sound manner by going from a pair to a triplet.

Part of this triplet is a pair that expresses the kinetic energy. This pair is accompanied by a signal that describes the non-linear configuration under which the kinetic energy can be exchanged. This rather abstract picture is best described by concrete examples. One for thermo-fluid systems and one for mechanical systems.

In the following we discuss the modeling of two domains using triplets.

3.1. Triplet Interface for Thermofluid Systems

Here is the triplet interface for thermo-fluid streams:

• $r$: inertial pressure (potential)
• $\dot{m}$: mass-flow rate (flow)
• $\widehat{\mathsf{\Theta }}$: Vector representing the thermodynamic state of the medium (signal)

-

$\widehat{p}$: steady-mass flow pressure

-

$\widehat{h}$: steady-mass flow enthalpy

-

$X$: mass fractions

Key to understand this interface is the dynamic of the inertial pressure $r$. It represents the pressure needed to accelerate a mass-flow which forms linear law based on the inertance [12,13] of the fluid $L$:

$$r=L\frac{d\dot{m}}{dt}$$

The inertance $L$ is thereby independent of the thermodynamic state of the fluid and just defined by the geometry of the flow with the length$ds$and cross section area $A$:

$$L=\int \frac{ds}{A}$$

If we now enforce that the inertance law is part of every single component, we have fulfilled criterion (a).

To fulfill criterion (b), we decompose the pressure into two components:

$$p=\widehat{p}+r$$

The complement to $r$ is hence the newly defined steady-mass flow pressure $\widehat{p}$ For the steady-state with a constant mass flow rate, $r$ will go to zero and $p=\widehat{p},$ which explains its name.

A very elegant way to reduce the interaction points between potential energies and kinetic energy is to use $\widehat{p}$ to form an approximation for the thermodynamic state $\widehat{\mathsf{\Theta }}$ of the fluid. Furthermore, where it is not reasonable to mix pressure, it is feasible to provide mixing laws for $\widehat{p}$. This enables an explicit downstream computation of $\widehat{\mathsf{\Theta }}.$The kinetic energy will then only interact at the boundaries of streams which are either volume elements or sources and sinks. At junctions, the kinetic energy of different streams may interact with each other. Ultimately this means that a different spatial resolution is applied for the effects of $r$ than for $\widehat{p}$. This assumption is in line with many modeling tasks and deviations are limited to rapid transients where model uncertainty is typically very high. At steady mass flows there is no difference which is important for most applications.

The resulting computational scheme is that all the non-linear computations regarding the thermodynamic state including properties like viscosity, heat conductivity or density are performed in explicit form downstream. This means that the signal $\widehat{\mathsf{\Theta }}$ is always computed from the source to the sink. To ensure this, all circular flows need to contain at least one volume. A linear system then computes the time derivative of all mass-flows. This system is formed by the equation for the pairs of $r$ and $\dot{m}$. It is guaranteed to be linear because of the characteristics of $L$ and it will be regular since all inertances are strictly greater than zero.

The reader is advised to literature [14,15] to follow the above statements in detail. Revisiting our example of section 2, there is fortunately a less complex replacement model. Since our models assume constant density and viscosity, there is no need for the decomposition of $p$ and for the signal $\widehat{\mathsf{\Theta }}$. The solution would (in this particular case) be equivalent to represent the inertance for each mass-flow but we are free to neglect the compressibility. The corresponding diagram is shown in Figure 6 and the simulation result in Figure 7. Also because each mass flow is represented by the sum of all of its inertances, the frequency of the response is low.

3.2. Triplet Interface for Mechanical Systems

It is possible to formulate a similar triplet to model multi-body mechanics:

• $v$: velocity (potential)
• $f$: force (flow)
• $s$: position (signal)

Where $v$ and $s$ represent positions and velocities in a generalized way for both translatory and rotational motion.

The product of velocity and force thereby represents a flow of energy which will contain also the kinetic energy. To ensure criterion (a), we simply must enforce that there are no massless components or (by the means of connection rules) that there is at least a mass representing all degrees of freedom of a kinetic chain.

To fulfill criterion (b), we have to address the elements of a multi-body system that define the degrees of freedom for the system: the joint elements. These are typically revolute joints or prismatic joints but exist in various combinations that define a different set of positional states. Naturally, the velocity is defined as direct derivative of the position and the acceleration as second derivative:

$$\frac{dv}{dt}=\frac{d^{2}s}{d{t}^2}$$

In section 3.1, we have reduced the interaction of kinetic and potential energy by lowering the spatial resolution of the kinetic part. For multi-body mechanics we have to reduce the temporal resolution of the kinetic part. We can do so by applying a first-order filter:

$$\frac{dv}{dt}{T}_D=\frac{ds}{dt}-v$$

Using the time constant $T_D$, we can now effectively limit the occurring frequencies of the mechanical system (at the cost of becoming more energy dissipative) without modifying the steady-state solution.

Just for the sake of quick illustration: Figure 8 from [16] shows the dynamics of a lightweight object moved in clamp modeled by a very stiff spring. The figure simply illustrates how the oscillatory dynamics is approximated with a replacement dynamics leading to the same (quasi) steady-state solution.

Using the triplet, we can formulate a difference between the motion in the elastic regime $ds/dt$ and the velocity in the kinetic regime $v$. Because of this separation of two regimes, the author has denoted it as dialectic mechanics. The interface enables again that all non-linear equations regarding the geometrical configuration can be formulated in explicit form from the root of the kinematic chain to its end by the signal $s$. To ensure this tree structure, all kinematic loops must contain a flexible element. The pair of $v$ and $f$ then sets up a linear system of equations for the kinetic energy.

First implementations and analysis of this modeling approach are presented in [16,17]. Models using this interface are especially suitable for the simulation of contacts and limited joints also under hard-real time constraints. We expect many future applications in robotics and related tasks.

4. Triplets and Object-Oriented Modeling

The two triplets presented in section 3 can be directly transformed into code of the equation-based, object-oriented language Modelica. Table 5 contains the corresponding code. For both domains, corresponding Modelica libraries have been developed.

The DLR ThermoFluid Stream Library [15] is a fully functional open-source that is already in use both in industry and academia. Figure 9 shows an exemplary model diagram of the thermal architecture of an electric car. We can clearly see that a complex architecture is nicely composed out of its individual components. The end-user is not directly confronted with the triplets. The modeler simply needs to take care about the connection rules that all loops contain at least one volumetric element.

The library for Dialectic Mechanics is still under internal development and testing. Figure 10 shows the model diagram of a kinematic loop for a hypothetical landing gear kinematic. This demonstrates that also here, complex systems can be composed in an object-oriented fashion and the end-user simply has to uphold the connection rules. The kinematic loops a closed by flexible elements and each chain contains at least one mass.

Using Modelica for these examples, a complete DAE system is built based on the composed models. The complete system of equations is then flattened (meaning that all equations are collected in a global set) and then state variables are selected and symbolic index reduction (as for example in [18]) is performed to transform the DAE into a set of Ordinary Differental Equations (ODEs). This typically results in a larger piece of simulation code where all variables are globally collected.

5. Triplets and Linear Implicit Equilibrium Dynamics

We have derived the triplet interface out of the motivation to have a robust modeling interface that ensures solvability and yet enables effective modeling by the suppression of high frequency behavior (that occurs in real physical systems). It is remarkable that for both cases the resulting computational structure of the complete system is of a particular form that we define here as Linear Implicit Equilibrium Dynamics (LIED).

A DAE system with potential state derivatives $\dot{\mathit{x}}$, time $t$ and algebraic variables $\mathbf{w}$

$$0=\mathit{F}(\dot{\mathit{x}},\mathit{x},\mathbf{w},t)$$

is defined as LIED system when it can be transformed into the following form:

$$\left[\begin{array}{c}{\mathbf{w}}_E\\ {\dot{\mathit{x}}}_E\end{array}\right]=\mathrm{g}({\mathit{x}}_I,{\mathit{x}}_E,t)$$

$$\mathbf{A}({\mathit{x}}_I,{\mathit{x}}_E,{\mathbf{w}}_E)\left[\begin{array}{c}{\mathbf{w}}_I\\ {\dot{\mathit{x}}}_I\end{array}\right]=\mathrm{f}({\mathit{x}}_I,{\mathit{x}}_E,{\mathbf{w}}_E,t)$$

We see that both the algebraic variables as well as the state derivatives can be split into a fully explicit part (${\dot{\mathit{x}}}_E;{\mathbf{w}}_E$) and a part (${\dot{\mathit{x}}}_I;{\mathbf{w}}_I$) with a linear system in implicit form expressed by the matrix $\mathbf{A}$. Furthermore, the following conditions shall hold true:

• ${\dot{\mathit{x}}}_E\cap {\dot{\mathit{x}}}_I\subseteq \dot{\mathit{x}}$
• ${\mathbf{w}}_E\cap {\mathbf{w}}_I\supseteq \mathbf{w}$
• ${\dot{\mathit{x}}}_E\cap {\dot{\mathit{x}}}_I\cap {\mathbf{w}}_I\supseteq \dot{\mathit{x}}$
• ${\dot{\mathit{x}}}_E,{\dot{\mathit{x}}}_I,{\mathbf{w}}_E,{\mathbf{w}}_I\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{j}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}$

These conditions essentially mean that it is allowed to perform certain symbolic mechanism of index reduction such as the dummy derivative method [18] originating from Pantelides [19]. Using this method, states variables of $\mathit{x}$ can be transformed to algebraic variables in ${\mathbf{w}}_I$ and further derivatives may be added to ${\mathbf{w}}_I$ or ${\mathbf{w}}_E$. In practice this is important because it means that the linear implicit dynamics can be expressed by far fewer states than suggested in the original DAE formulation.

To get a practical feeling for this LIED form, let us discuss the decomposition of the linear and non-linear part, for the presented triplets.

• For thermofluid streams, the signal for the thermodynamic state $\widehat{\mathsf{\Theta }}$ will be part of${\mathbf{w}}_E$. States associated with this state such as the temperature of a volume will be part of ${\mathit{x}}_E$. A subset of the flow variables $\dot{m}$ will form ${\mathit{x}}_I$. The subset is chosen in such a way that all streams are described by ${\mathit{x}}_I$ in a non-redundant way. The redundant mass flows as well as the inertial pressures $r$ become part of ${\mathbf{w}}_I.$
• For dialectic mechanics, the signal for the position will be part of ${\mathbf{w}}_E.$ The potential and flow variables of the interface $v$ and $f$ will be part of ${\mathbf{w}}_I$. The state vector ${\mathit{x}}_E$ is typically empty. The members of ${\mathit{x}}_I$ are chosen from the internal variables in such a way that the motion of the kinematic chain is described in a non-redundant way. For each degree of freedom its kinetic velocity and its position form state variables in ${\mathit{x}}_I$.

6. How to Create Simulation Code for LIED Systems?

The original intention of the LIED approach was simply to ensure that no non-linear system is created in implicit form that spans across the components and hence a robust solution of the model evaluation could be taken for granted, given robust component models. When we started with it, we expected it to be the only notable change from other DAE models in Modelica and that all other features of a of code generation (as briefly outlined at the end of section 4) would basically remain untouched.

However, over time, it was gradually recognized that LIED systems are much simpler to transfer to simulation code than general DAEs. Let us go through these observed simplifications one-by one:

• Because we avoid the creation of non-linear equation systems, we do not need a non-linear equation system solver anymore.
• For the same reason, constraint equations among potential states cannot be non-linear and hence no dynamic state selection is needed [20]
• Even stronger: we can select the states on component level. For example, the stars in Figure 9 and 10 mark those components that define states in ${\mathit{x}}_I$.
• Because we can select the states on component level, this means that the dummy-derivative method can be applied also on component level before system composition.
• Since the goal of the linear equation system is to have a synchronized replacement dynamics towards the equilibrium, we know suitable tearing variables for this system. These will be the linear state derivatives: ${\dot{\mathit{x}}}_L$ or at least a subset of it.
• The residual for a tearing variable can be attributed to the same component as the tearing variable.

The terminology of the last two points may require further explanation. A tearing variable is a variable that can be used to probe the reaction of an algebraic system which is assessed by a corresponding residual. If the system is linear: $\mathbf{M}\mathbf{w}=\mathbf{b}$, the probing can be used to construct M and w and find the solution for w.

All of the points above represent observations resulting from modeling many components and system examples using the LIED approach. However, these observations have profound implications. For each component we know:

• the set of pairs of state-variables and their derivatives it adds to the system.
• the set of pairs of tearing and residual variables it adds to the system.

If this is the case, we can basically causalize everything already on the component level. In concrete terms, this means for each component:

• we stipulate the states
• we stipulate the tearing variables of the linear system and the corresponding residuals
• we perform the dummy derivative method on those equations where necessary.
• we define the causality of the interface variables
• we causalize all equations into assignments in a particular order
• we group the list of assignments depending on their dependence of the inputs.

Practical experience so far indicates that performing index reduction to transform to the LIED form can be performed in a very methodical and deterministic manner. It is thus far easier to generate simulation code for the LIED modeling approach than it would be for general higher-index DAEs. Neither there is a need for global flattening anymore nor are elaborate heuristics needed for the selection of state or tearing variables. Indeed, the generation of simulation code is so easy that a direct implementation in C++ becomes feasible. To demonstrate this, the following code excerpts illustrate the implementation for a thermofluid stream library (using idealized water) in C++.

First, we have to define the interface. This is naturally more tedious than in a language like Modelica because there is no direct support in the C++ language. Yet, it is feasible and after all, interfaces only need to be defined once:

To best understand the interface, let us look at the classes ThermalPlugOut for a nominal outlet flow and at ThermalPlugIn for a nominal inlet flow first. These contain the same 3 components as the corresponding Modelica connector of the DLR ThermoFluid Stream library.

There are two notable differences however. In Modelica, inertial pressure and mass flow were not causalized signals as in the C++ implementation. Also the mass-flow signal in the C++ library consists of the mass-flow rate and its derivative. In Modelica, this is not necessary since symbolic differentiation can be applied by the Modelica compiler. Using this interface, we can now implement a component such as the pressure drop:

First, we declare our interface for outlet and inlet, then we have to implement three methods. The first is evalState and computes the thermodynamic state downstream:

The second method is evalFlow to ensure what flows in is what flows out. However, this constraint is restated for the derivative. This is because the dummy derivative method is applied on the component level.

The third one is evalInertia that implements the law for the inertance as in the ThermoFluid Stream Library.

In similar manner the other components of our introductory example can be implemented. Each of these components declares its interfaces, defines and implements methods representing the computational blocks. For all these components, meta-information must be collected by a dedicated virtual method to register state and tearing variables as well as to track the signal dependence of the computing blocks. Finally, we can compose the introductory example of section 2 and 3 directly in C++:

Regarding that C++ is a statically compiled imperative general-purpose language, the final result is astonishingly close to what a modeler is used to from a domain specific language like Modelica.

When an instance of the class is coupled to a simulator, the meta information of all its computational methods is collected as well as the structural information regarding the signals, the states and the tearing variables for the linear equation system. Then the methods are put into right order for complete or partial model evaluation. The full model evaluation is thus simply a list of method calls that are called in their respective order. The linear implicit system defined by $\mathbf{A}({\mathit{x}}_I,{\mathit{x}}_E,{\mathbf{w}}_E)$ is reconstructed by using the tearing variables and their residuals for probing. This requires a number of partial model evaluations.

This also means that all component code is statically compiled but the total system composition is performed at run-time. Advanced tasks such as variable structure systems would thus be comparably easy to achieve.

This is also the purpose of this code demonstration. It is not suggesting that we should use C++ as modeling language like Modelica but to highlight that for a certain class of models the object-oriented modeling code could (and maybe should) be translated to object-oriented imperative code that can be statically compiled even before total system composition. This would avoid the flattening of all equations before code generation and help to overcome many limitations of current Modelica compilers with respect to scalability.

7. Conclusions

This paper started with a very fundamental consideration about the solvability of general formulations for classic physical systems, then presented a new interface for the object-oriented modeling and ended up showing a simplified generation of simulation code potentially suitable to handle challenges such as large-scale system simulation or variable structure systems. It is quite remarkable that the two derived criteria for handling the kinetic energy have such practical consequences for object-oriented modeling.

Indeed, that a more restrictive class of modeling enables a simpler compilation scheme is neither new nor surprising. The same can be said about the many conventional signal-based modeling schemes or simple modeling schemes as Forrester’s System Dynamics [21]. Typically, the disadvantage is that the easier generation of simulation code has to be paid by an inferior modeling approach and indeed modeling complex mechanics or thermofluid streams is painful when using purely signal-based approaches (nevertheless this pain has been taken in industrial practice all too often).

The remarkable thing about the LIED approach is that a simulation engineer has a simple scheme for code generation but can also conveniently model both mechanics and thermo-fluid streams in a very robust manner. The corresponding Modelica Libraries prove this [15,16]. Both application domains are known to be rather difficult but LIED can even be applied to the challenging parts of these fields such as handling stiff contact mechanics or elaborate by-passes in complex thermal architectures. It is yet unclear for what other domains LIED is an attractive choice.

Figure 11 attempts to qualitatively depict the trade-off between computational complexity and algorithmic complexity. LIED forms a very exposed point on a hypothetical Pareto front. This means that for a large number of applications it is a very attractive choice.

It might be not unimportant to note that the actual discovery of the LIED approach happened not the way it was presented in this paper. Rather the opposite reflects reality: practical applications came first; generalization and theory came second. The approach was simply born out of the necessity (or desperation) to find robust solutions for aircraft environmental control systems and to handle gripping of stiff objects with robots.

For this paper, this means that most of its statements are to be regarded as empirical. For some specific LIED systems, formal statements on structure and solvability could be made [14] but in general, the statements on suitable system composition and code generation are merely observations from simulation practice and test implementations. Further solidification of the theoretical understanding is needed. This is a call to ourselves but also an invitation to the research community.

As the question marks in Table 3 indicate, there remains an open question whether the LIED approach can also be successfully applied for electrical systems so that it is of practical value. Clearly, kinetic energy plays a far less important role for electrical systems, also the formulation of potential energy requires a generalized form. Hence, the reasoning of this paper cannot be directly converted. Maybe extensions are needed, maybe it is infeasible? Again, this open question is an invitation to the research community.

The primary motivation of this paper remains to raise awareness on this class of models and the possibilities it enables for the generation of simulation code. Although it states mostly observations, these observations seem valuable to share.