Stability is a widely used term and stability in systems and solutions are often separated because they usually do not correlate that much. In general, dynamics stability in systems is a measure of convergence to a stable value, or set, and is related to the dynamics in the system. Dynamic stability is well documented in the field of modelling [1] [2] as well as in the field of Control [3] [4] . On the other hand, stability in simulations is a measure of numerical convergence and is more or less only dependent on the eigenvalues of the system to be solved and the solver dynamics.

In distributed systems and simulations, these stability terms seem to be more connected in comparison to “continuous” systems, and the theory of distributed systems can be associated with the theory of sampled systems [5] and discrete dynamical systems [3], although ordinary stability theory of sampled systems does not state anything about numerical convergence regarding solver dynamics.

When distributing systems in a simulation, the distributed systems exchange data at a given rate, while all signals are held constant in between. This is also a source of error that can cause instabilities in the total distributed system, because the amount of energy exchanged between the subsystems is not correct since signals are held constant in between each global communication step.

Both dynamic stability, simulation stability and energy balance in a distributed system are discussed in the following.

Dynamical System stability can be analysed by various tools available in the literature. However, one of the most complete stability theories for both linear and non-linear dynamical system is the classical Lyapunov stability analysis [4], even though it focuses on continuous systems. However, extensions to this theory, that opens for analysing discrete- and sampled dynamical systems, is available such as hybrid dynamical systems theory [6].

Even though many analysis tools and theories are available for analysing dynamical systems, both distributed or not, the amount of work required for analysing the dynamical system stability increases drastically with the distributed system. Neither does it help that such analysis often includes iterations, e.g. for Lyapunov stability analysis a Lyapunov function candidate is chosen, tested, and if nothing can be concluded a new Lyapunov function candidate should be considered. Moreover, if these Lyapunov function candidates are chosen for local subsystems in the total distributed system, the sum of all Lyapunov function candidates must also hold for the Lyapunov stability criteria. Thus, a more practical and less time consuming way of estimating the stability condition in a distributed system would be beneficial.

One such practical approach for estimating stability of the dynamics in a distributed system is to map the energy flow in the total system in a sink-source study. In other words, if all energy sources and all energy dissipations can be established, it is possible to conclude something about the rate of energy, the power, in the total system. Hence, if more energy is being dissipated than added, dynamic stability in the total system can in most cases be concluded. This approach is related to passivity theory in the field of control, and has its roots from classical Lyapunov stability analysis.

In distributed systems, both dynamic stability and numerical stability must be assured at the same time because both depend on each other through the local and global time steps. This is elaborated in the following.

The stability of the local solvers in each subsystem is closely related to the eigenvalues in the subsystem itself. By knowing the eigenvalues it is possible to choose a solver and a time step, if the solver is a fixed step solver, that, in general, stabilizes the local solution. For example, a system given as \begin{equation} \dot{x}=f(x) \end{equation} having the eigenvalue $\lambda$, can be solved with the Euler integration method such that \begin{equation} x(t+\Delta t)=x(t)+f(x(t))\Delta t \end{equation} Then, $\Delta t$ must be chosen such that \begin{equation} |1+\lambda \Delta t|\leq 1,\; \Delta t>0 \end{equation} However, even though each subsystem in a distributed simulation are locally numerically stable, it is not guaranteed that the total distributed system is stable. Hence, a combined stability requirement is needed which takes into account both the dynamical stability and the numerical stability and relates them through the local solver time steps and the global communication time step. This is elaborated in the following.

An analysis tool for combined stability in linear distributed dynamical systems has been proposed by combining dynamic stability and solver stability [7] [8]. The method propose to use a numerical solver with a low order, such as the Euler integration method, Runge Kutta 2 or Runge Kutta 4, which often are basis for more complex numerical solver methods, and step through local time steps until the global simulation time step is reached. Then, it is possible to express the algebraic solution of the local system, including the solver dynamics, in between global time steps. By combining such local algebraic solutions into a global algebraic solution of the total distributed system, it is possible to check whether the amplitude of the eigenvalues of the global algebraic solution are less than 1 or not. If the former is a fact, the global numerical solution will converge. For detail matters the reader is referred to [7][8]. The drawbacks with this method are that it can be quite time consuming to calculate the local algebraic solutions in between each global time step, and that exact results can only be obtained for linear dynamical systems. Extensions to non-linear dynamical systems have been studied in [8], but these results are more conservative.

Example 1 illustrate the method applied to a linear system.

Example 1 Combined Stability Analysis

A linear dynamical system given as \begin{equation} \begin{split} \dot{\boldsymbol{x}}&=\boldsymbol{Ax}+\boldsymbol{Bu}\\ \boldsymbol{y}&=\boldsymbol{Hx} \end{split} \end{equation} where $\boldsymbol{x}$ is the state vector, $\boldsymbol{u}$ is the system input, $\boldsymbol{y}$ is the system output and $\boldsymbol{A}$, $\boldsymbol{B}$ and $\boldsymbol{H}$ are matrices containing the system dynamics, input mapping and output mapping, respectively. Assuming that this system is solved numerically by using the Euler integration method with local solver time step $\Delta t$, and that the global time step is given as $T_d$, the number of local solver steps between each global time step is defined as \begin{equation} n=\frac{T_d}{\Delta t} \end{equation} Note that $\Delta t$ and $T_d$ should be chosen such that $n\in \mathcal{N}_{>0}$ where $\mathcal{N}_{>0}$ is the set of all integers larger than zero. Then, it can be verified that the local solution between each time step can be expressed as \begin{equation} \boldsymbol{x}(t_0+T_d)=\boldsymbol{A}_n\boldsymbol{x}(t_0)+\boldsymbol{B}_n\boldsymbol{u}(t_0) \end{equation} where \begin{equation} \begin{split} \boldsymbol{A}_n&=(\boldsymbol{I}+\Delta t\boldsymbol{A})^{n}\\ \boldsymbol{B}_n&=\boldsymbol{A}^{-1}(\boldsymbol{A}_n-\boldsymbol{I})\boldsymbol{B} \end{split} \end{equation} Thus, by having $N$ number of subsystems in the distributed system and by utilizing the I/O mapping (e.g. $\boldsymbol{y}_1=\boldsymbol{H}_1\boldsymbol{x}_1$, $\boldsymbol{u}_1=\boldsymbol{y}_2$), the solution of the distributed system in between global time steps can be expressed as \begin{equation} \boldsymbol{x}_d(t_0+T_d)=\boldsymbol{S}_d\boldsymbol{x}_d(t_0) \end{equation} where $\boldsymbol{S}_d$ is called the global distributed simulation solution matrix. Hence, if $|\text{eig}(\boldsymbol{S}_d)|<1$, the total distributed numerical simulation is stable.

In general, energy is incorrectly transferred between two coupled simulators due to the fact that their states are evolved independently of each other between discrete communication time points. In effect, energy is either created or destroyed through the co-simulation coupling during each macro time step. This residual energy directly alters the total energy of the overall coupled system. It thus changes its dynamics, deteriorates simulation accuracy, and may pose a threat to numerical stability. For detail matter, the reader is referred to [9] and ECCO.

1. Dean C. Karnopp, Donald L. Margolis, and Ronald C. Rosenberg., 2006. System Dynamics: Modeling and Simulation of Mechatronic Systems.. Wiley & Sons, Inc., New York, NY, USA.
2. Francois E. Cellier and Ernesto Kofman., 2006. Continuous System Simulation. Springer-Verlag New York, Inc., Secaucus, NJ, USA.
3. Chi-Tsong Chen., 1998. Linear System Theory and Design.. Oxford University Press, Inc., New York, NY, USA, 3rd edition.
4. H. K. Khalil., 2002. Nonlinear Systems.. Prentice Hall, PTR.
5. T. Chen and B. A. Francis, 1991. Input-Output Stability Of Sampled-Data Systems.. IEEE Transactions on Automatic Control, 36(1), p.50-58.
6. R. Goebel and R. G. Sanfelice., 2012. Hybrid Dynamical Systems: Modeling, Stability, and Robustness.. Princeton University Press.
7. S. Skjong and E. Pedersen, 2016. The Theory of Bond Graphs in Distributed Systems and Simulations. International Conference on Bond Graph Modeling and Simulation (ICBGM'16), Montreal, Canada. Spciety for Modeling & Simulation International.
8. S. Skjong and E. Pedersen, 2016. On Numerical Stability in Dynamical Distributed Simulations.. To be submitted.
9. S. Sadjina, L. T. Kyllingstad, S. Skjong, and E. Pedersen., 2016. Energy Conservation and power bonds in co-simulations: non-iterative adaptive step size control and error estimation.. Engineering with Computers.