On meeting Energy Balance Errors in Cosimulations

TL;DR

This paper addresses the challenge of energy balance errors in coupled simulations, proposing a method to accurately balance energy across subsystems to improve simulation accuracy.

Contribution

It introduces a novel method to precisely balance energy in co-simulations, overcoming limitations of previous balance correction techniques.

Findings

The proposed method effectively reduces energy balance errors.

It improves the accuracy of coupled simulation results.

The method maintains system-wide energy consistency.

Abstract

In engineering, it is a common desire to couple existing simulation tools together into one big system by passing information from subsystems as parameters into the subsystems under influence. As executed at fixed time points, this data exchange gives the global method a strong explicit component, and as flows of conserved quantities are passed across subsystem boundaries, it is not ensured that systemwide balances are fulfilled: the system is not solved as one single equation system. These balance errors can accumulate and make simulation results inaccurate. Use of higher-order extrapolation in exchanged data can reduce this problem but cannot solve it. The remaining balance error has been handled in past work with balance correction methods which compensate these errors by adding corrections for the balances to the signal in next coupling time step. Further past work combined smooth extrapolation of exchanged data and balance correction. This gives rise to the problem that establishing balance of one quantity a posteriori due to the time delay in general cannot establish or even disturbs the balances of quantities that depend on the exchanged quantities, usually energy. In this work, a method is suggested which allows to choose the quantity that should be balanced to be that energy, and to accurately balance it.

Tables

Table 1. Table 1: Cosimulation scheme

General scheme
$S_1$	$S_2$
System States
${𝒙}_1$	${𝒙}_2$
Outputs
${𝒖}_{21}$	${𝒖}_{12}$
Inputs
${𝒖}_{12}$	${𝒖}_{21}$
Equations
$\dot{{𝒙}_1}={𝒇}_1({𝒙}_1,\mathrm{Ext}\left({𝒖}_{12}\right))$	$\dot{{𝒙}_2}={𝒇}_2({𝒙}_2,\mathrm{Ext}\left({𝒖}_{21}\right))$

Table 2. Table 2: Cosimulation scheme with balance correction

Cosimulation scheme
$S_1$	$S_2$
System States
${𝒙}_1$	${𝒙}_2$
Outputs
${𝒖}_{21}$	${𝒖}_{12}$
Inputs
${𝒖}_{12}$	${𝒖}_{21}$
Equations
$\dot{{𝒙}_1}={𝒇}_1({𝒙}_1,\mathrm{Ext}\left({𝒖}_{12}\right))+\Delta E_1^{j-1}{\varphi }_j\left(t\right)$	$\dot{{𝒙}_2}={𝒇}_2({𝒙}_2,\mathrm{Ext}\left({𝒖}_{21}\right))+\Delta E_2^{j-1}{\varphi }_j\left(t\right)$

Table 3. Table 3: Cosimulation Schemes for the spring-mass system, left constant, right linear extrapolation

Spring	Mass
System States
$x_1:=s=x$	$x_2:=v=\dot{x}$
Outputs
$u_{21}:=F=-cx$	$u_{12}:=v=\dot{x}$
Inputs
$u_{12}$	$u_{21}$
Equations
$\dot{x_1}=\mathrm{Ext}\left(u_{12}\right)=v$	$\dot{x_2}=-\frac{1}{m}\mathrm{Ext}\left(u_{21}\right)$
	$=-\frac{F}{m}$

Equations

$$S_1:$$

$${\dot{\mathbf{x}}}_1$$

$$0$$

$$S_2:$$

$${\dot{\mathbf{x}}}_2$$

$$0$$

$$S_1:$$

$${\dot{\mathbf{x}}}_1$$

$$0$$

$$S_2:$$

$${\dot{\mathbf{x}}}_2$$

$$0$$

$${\dot{\mathbf{x}}}_{{\mathcal{D}}_k}=f_{{\mathcal{D}}_k}(t,\left[{\mathbf{x}}_{{\mathcal{D}}_k},{\mathbf{x}}_{{\mathcal{I}}_k}\right]),$$

$$\Vert {\mathbf{ϵ}}_{\Delta ,S}(T_{j+1})\Vert \le \frac{1}{L}\left(\widetilde{\stackrel{~}{C}}{h}_{\Delta }^p+C{H}^{P+1}\right)\left(e^{L(T_{j+1}-T_0)}-1\right)$$

$$\Delta E_i^j:={\int }_{T_{j-1}}^{T_j}{\mathbf{u}}_{i}dt-{\int }_{T_{j-1}}^{T_j}{\overline{\mathbf{u}}}_{i}dt$$

$$\dot{\mathbf{x}}=\text{tensor}[]A\mathbf{x}$$

$$\{\begin{array}{l}{\dot{x}}_1=a_{1,1}{x}_1+a_{1,2}\mathrm{Ext}(x_2)\\ {\dot{x}}_2=a_{2,2}{x}_2+a_{2,1}\mathrm{Ext}(x_1)\end{array}.$$

$$\dot{\mathbf{x}}=\text{tensor}[]A\mathbf{x}=\left(\begin{array}{cc}0& 1\\ -\frac{c}{m}& \left(-\frac{d}{m}\right)\end{array}\right)\mathbf{x},\mathbf{x}=\left(\begin{array}{c}x\\ \dot{x}\end{array}\right)$$

$$\text{tensor}[]A=\left(\begin{array}{cc}0& 1\\ -\frac{c}{m}& 0\end{array}\right).$$

$$\mathbf{z}={\mathbf{z}}_0\exp \left(\text{tensor}[]A(t-t_0)\right)$$

$$z_1=x=\sqrt{\frac{m}{c}}{\dot{x}}_0\sin \left(\sqrt{\frac{c}{m}}(t-t_0)\right)+x_0\cos \left(\sqrt{\frac{c}{m}}(t-t_0)\right).$$

$$\begin{array}{r}x=x_0\cos \left(\sqrt{\frac{c}{m}}t\right),\dot{x}=-\sqrt{\frac{c}{m}}{x}_0\sin \left(\sqrt{\frac{c}{m}}t\right),\\ \text{and }\ddot{x}=-\frac{c}{m}{x}_0\cos \left(\sqrt{\frac{c}{m}}t\right).\end{array}$$

$${\hat{P}}_{12}(P_{21},P_{12})=-{\hat{P}}_{21},$$

$$P_{21}({\mathbf{x}}_1(t),{\mathbf{u}}_{12\setminus k},(u_{12}{)}_k(t))=\mathrm{Ext}({\hat{P}}_{12})$$

$$(u_{12}{)}_k(t)=P_{21}^{-{\mathbf{u}}_k}(\mathrm{Ext}({\hat{P}}_{12}),{\mathbf{x}}_1,{\mathbf{u}}_{12\setminus k})$$

$$P_{21}^{-v}$$

$$P_{12}^{-f}$$

$$\dot{\mathbf{x}}=-\text{tensor}[]M{\nabla }_{\mathbf{x}}{\mathcal{P}}^T.$$

$$\dot{\mathcal{P}}(x)=⟨{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}),\dot{\mathbf{x}}⟩=⟨{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}),-\text{tensor}[]M{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}{)}^T⟩,$$

$$\dot{\mathcal{P}}(x)=⟨{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}),\dot{\mathbf{x}}⟩=⟨{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}),-\text{tensor}[]M{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}{)}^T⟩<0,$$

$$⟨\dot{\mathbf{x}},{\nabla }_{\mathbf{x}}\mathcal{P}⟩=⟨{\nabla }_{\mathbf{x}}\mathcal{P},-\text{tensor}[]M{\nabla }_{\mathbf{x}}\mathcal{P}(\mathbf{x}{)}^T⟩=0,$$

$$\text{tensor}[]M=\text{tensor}[]M{}_{\text{skew}}+\text{tensor}[]M{}_{\text{pos.semidef}},$$

$$\dot{\mathcal{P}}({\Phi }^t\mathbf{x})=⟨\nabla \mathcal{P},\dot{{\Phi }^t\mathbf{x}}⟩\le 0.$$

$$\dot{{\mathcal{P}}_i}(x):=P_i$$

Full text

The remaining balance error has been handled in past work with balance correction methods which compensate these errors by adding corrections for the balances to the signal in next coupling time step. Further past work combined smooth extrapolation of exchanged data and balance correction.This gives rise to the problem that establishing balance of one quantity a posteriori due to the time delay in general cannot establish or even disturbs the balances of quantities that depend on the exchanged quantities, usually energy. In this work, a method is suggested which allows to choose the quantity that should be balanced to be that energy, and to accurately balance it.

Keywords:

Cosimulation, coupled problems, simulator coupling, explicit coupling, stability, convergence, balance correction

1 Introduction

Engineers are increasingly relying on numerical simulation techniques. Models and simulation tools for various physical problems have come into existence in the past decades. The desire to simulate a system that consists of well described and treated subsystems by using appropriate solvers for each subsystem and letting them exchange the data that forms the mutual influence is immanent.

The situation usually is described by two coupled differential-algebraic systems $S_{1}$ and $S_{2}$ that together form a system $S$:

[TABLE]

The $(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{1},\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{2})$ are the differential states of $S$, their splitting into $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i}$ determines the subsystems $S_{i}$ together with the choices of the $\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i}$. In Co-Simulation the immediate mutual influence of subsystems is replaced by exchanging data at fixed time points and subsystems are solved separately and parallely but using the received parameter:

[TABLE]

where $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{i}$ are given by coupling conditions that have to be fulfilled at exchange times $T_{k}$

[TABLE]

and are not dependent on subsystem $i$’s states any more, so are mere parameters between exchange time steps.

Full row rank of $d_{\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{i}}\mathchoice{\displaystyle\boldsymbol{g}}{\textstyle\boldsymbol{g}}{\scriptstyle\boldsymbol{g}}{\scriptscriptstyle\boldsymbol{g}}_{i}$ can be assumed, such that the differential-algebraic systems are of index 1. This description of the setting is widespread ([1]). With the $\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{ij}$ being solved for $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}$ inside the $S_{j}$ (let solvability be given), for systems with more than two subsystems it is more convenient to write output variable ${\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}}_{j}$ and now redefine $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}$ as the input of $S_{i}$, consisting of some components of the outputs $\mathchoice{\displaystyle\boldsymbol{y}}{\textstyle\boldsymbol{y}}{\scriptstyle\boldsymbol{y}}{\scriptscriptstyle\boldsymbol{y}}_{j}$ [2]. This structure is defined as kind of a standard for connecting simulators for cosimulation by the Functional Mockup Interface Standard [3]. It defines clearly what information a subsystems implementation provides.

In Co-Simulation the variables establishing the mutual influence of subsystems are exchanged at fixed time points. This results in continuous variables being approximated by piecewise constant extrapolation, as shown in the following picture:

If one does not want to iterate on those inputs by restarting the simulations using the newly calculated inputs, one just proceeds to the next timestep.

This gives the calculations an explicite component, the mutual influence is now not immediate any more, inducing the typical stability problems, besides the approximation errors.

But for good reasons, explicit co-simulation is a widely used method: It allows to put separate submodels, for each of which a solver exists, together into one system and simulate that system by simulating each subsystem with its specialised solver - examination of mutual influence becomes possible without rewriting everything into one system, and simulation speed benefits from the parallel calculation of the submodels. Usually it is highly desirable that a simulation scheme does not require repeating of exchange time intervals or iteration, as for many comercial simulation tools this would already require too deep intrusion into the subsystems method and too much programming in the coupling algorithm.

The following fields of work on explicite co-simulation can be named to be the ones of most interest:

Improvement of the approximation of the exchanged data will most often improve simulation results [4]. This is usually done by higher-order extrapolation of exchanged data, as shown in this plot, where the function plotted with dots is linearly extrapolated:

When the mutual influence between subsystems consists of flow of conserved quantities like mass or energy, it turns out that the improvement of the approximation of this influence by extrapolation of past data is not sufficient to establish the conservation of those quantities with the necessary accuracy. The error that arises from the error in exchange adds up over time and becomes obvious (and lethal to simulation results many times). In a cooling cycle example ([5, Section 6.3]), a gain of 1.25% in coolant mass occurs when simulating a common situation.

It has been tried to meet this challenge by passing the amount of exchanged quantity for the past timestep along with the actual flow on to the receiving system, where then the error that has just been commited is calculated and added to the current flow to compensate the past error. For well damped example problems in fluid circles this method has fulfilled the expectations [6]. It has been labelled balance correction. 3. 3.

There is good reason to prevent jumps in exchanged data by smoothing. Higher order extrapolation polynomials cannot make extrapolated data at the end of the exchange timestep match the newly given value.

While under all circumstances it is desirable to guess the influence information at the subsystem boundaries from past data as well as possible, balance correction techniques bear the profound problem that for establishing conservation * a posteriori* over the whole time an instantanious error in the exchanged data has to be accepted. More precisely, balance correction means making an error in the exchanged signal $u$, for lowering the accumulated error in the amount of that quantity and thus lowering the error in states $x$ that would be caused by this now persistent error in amount. The purposefully commited error which the refeed of past errors actually is of course bears an error in the derivatives of the exchanged data in it, which might cause dynamics in the receiving system and its neighbours, especially if subsystems act on quicker time scales. Our past work [7] lowered the error in the derivatives by construction and use of suitable functions for smoothing during switching and adding of correction terms.

Although in [1] and [8] the convergence of explicite cosimulation methods for ODE and index one PDE was proven, and the result in [8] was easily extended to balance correction methods, which thus are proven justified – it can be easily shown that those methods are not stable, neither with balance correction [8]. This is the most severe restriction to explicite simulator coupling.

From seeing instability as a rise in systems energy as mentioned in [8, Sec.4.2.2] and [7, Sec.3.2] and considering the power that is given by variables at the subsystems interfaces, in this contribution a stable explicit coupling scheme is derived.

1.1 Some preliminaries and notation

Always, a time interval or quantities belonging to one are indexed with the index of its right boundary: $\Delta t_{i}:=[t_{i-1},t_{i})$. Indexing of times begins with 0. Big letters refer to time exchange steps, e.g $\Delta T_{k}:=[T_{k-1},T_{k})$ is an interval between exchange times, whereas above interval might denote a subsystems step. Consistently, $k=1,...,N$, and above $i$ ranges from $1$ to $n$ too.

Let $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})_{ji}^{k}(t)$ denote the value of the input variable as it is assumed to be by $S_{j}$, some extrapolation of ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}$, and $\overline{\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})}_{ji}^{k}(t)$ be the flow as it is used for calculation constructed on $\left[T_{k-1},T_{k}\right)$ – it may be different from $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})_{ji}^{k}(t)$, for example a smooth combination of $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})_{ji}^{k-1}(t)$ and $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})_{ji}^{k}(t)$. The error $\Delta E_{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}},k}=\int_{T_{k-1}}^{T_{k}}\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ji}-\overline{\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}})}_{ji}^{k}(t)dt$ is the balance error if $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ji}$ is a flux of a conserved quantity, but it is defined for arbitrary quantities. The expression $(\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}})_{i}$ denotes the $i$-th component of the vector $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}$, similarly, $(\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}})_{\mathcal{I}}$ denotes all component of the vector $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}$ whose indices are in the set $\mathcal{I}$.

1.2 Explicit Simulator Coupling

In industrial environment, explicit simulator coupling cannot always be avoided. There have been efforts to standardize model interfaces in order to enable coupling into one monolitic system and solve them using one ODE/DAE solver. Such an effort is the FMI standard [3]. But those efforts so far have not led to replacement of simulator coupling. This is, of course, due to the numerous given legacy codes, due to the fact that parallelization is also advantageous, especially if subsystems softwares are equipped with solvers that are customized to the problem, among other reasons. For example, including the residual of an FEM equation system into the monolithic global system is not sensible due to the amount of data that has to be passed to the monolitic solver, and the build-in solver usually is highly adapted to the problems needs. The solver for the global system cannot be optimal for all subsystems.

So, simulator coupling remains and will remain a field of work.

1.2.1 Iterative Methods

One can repeat the calculation for timestep $k$ using the inputs $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ji}(t)$ determined from the subsystems states numerical solution $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{\Delta}(t)$ just calculated - doing so all at the same time would correspond to Jacobi iteration or waveform iteration, doing so one after the other, preferably in the order of the dependency of the systems, would correspond to the Gauss-Seidel-Scheme [4]. These schemes usually converge [1]. One then has to program an external iteration procedure including a convergence criteria, and information exchange now concerns data that is time-dependent everywhere in $T_{i-1},T_{i}$. Implementing this requires skills that are frequently not sufficiently available in the environments in question, and sometimes it is computationally cheaper to lower the exchange time step size instead of iterating.

1.2.2 Making Inputs consistent

The coupling equations (9) - (10) together with the algebraic equations (6) and (8) of the subsystems form a global system of algebraic equations. If the graph of the information flow through the subsystems has no loops and one applies an according Gauss-Seidel scheme, then the exchanged data is consistent in the sense that the equations (9), (10), (6) and (8) are fulfilled - in the general case they are not if one just solves (9), (10) with respect to the input.

To avoid those errors, it has been suggested to solve the system (9), (10), (6) and (8) at exchange times. While this is obviously desirable for iterative and explicite schemes, it requires that subsystems solvers provide an the residual of their algebraic equations to callers. As mentioned above repeatedly, this need not be the case. Further, this again requires considerable programming effort, including a nonlinear solver call and convergence criteria.

1.2.3 Trying to tackle some drawbacks in explicit coupling

Finally, with the reasons given above, frequently the setting allows only for explicit simulator coupling. The drawbacks of using piecewise explicit extrapolation are:

• •

discontinuity,

• •

disbalance in amounts of conserved inputs but also in quantities depending on them, for example energy. Balance correction as given by Table 2, first suggested in [6], provides relief to disbalances in conserved inputs to a degree that it makes simulations possible that would be useless without.

Balance correction methods were applied to nonconserved quantities in [7], as even such quantities have some conservation properties in space and time due to their continuity, and examination will go on here.

• •

high derivatives in signals which are induced by the method, especially the balance correction recontributions, which are added as product of the missed amount and a hat function with integral 1, can cause high derivatives.

In [7] the issues of balance correction and smoothing of values as well as smoothing of derivatives were considered together, showing ways to minimize unwanted behavior and side effects of balance correction and smoothing. This was mainly done by extending the recontribution of the correction over several time intervals. But any time delayed correction, as all balance corection methods are, cannot establish balance in quantities depending on inputs as energy. Translating energy gain into mathematical terminology, we state that nothing can make an explicit method stable.

1.3 Aim of this work

The state of the art offers relief for many problems that arise when one extrapolates data during cosimulation, also to errors in balances as consequence of extrapolation errors. But finally, cosimulation is an explicit method and is not stable ([7], [8]), even when balance correction is applied. We see the instability as a rise in energy, which usually is a norm or half-norm for the system, that stems from errors in power acting on subsystems interfaces due to extrapolation errors in the factors of that power. Even if balance correction is applied on the factor, the time delay with which the refeed is applied in general has side effects, as the exchanged quantity then arrives at a time when the systems states already changed. For example the force $f$ is a factor of the mechanical energy. If a balance correction method is applied to the impulse $p=\int f\,dt$, it is not established that the energy balance is fulfilled. This is described in Section 3 and an example is given in Section 3.3. In Section 4.1 an exchange scheme is presented that balances exchanged energy by first choosing the power as exchanged variable and then making subsystems agree on which amount of energy should be exchanged, making the algorithm stable.

2 Exchanging Factors of conserved Quantities

2.1 Convergence of cosimulation schemes

In our preceding work [8] convergence rates for the standard cosimulation scheme

as given in table 1 and before in [8, Section 3.1] were derived. Shifting the solving of the equations (9) ff., $0={\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}}_{ij}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{j},\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}_{j},\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij})$ to the receiving system establishes that the exchanged quantities $u_{ij}$ are given by some of the subsystems states, see [8, Section 3.1]. This restricts generality slightly, as one has to assume $\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{ij}=\mathchoice{\displaystyle\boldsymbol{h}}{\textstyle\boldsymbol{h}}{\scriptstyle\boldsymbol{h}}{\scriptscriptstyle\boldsymbol{h}}_{ij}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{j},\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij})$ independent of any algebraic variable for that shifting of solving, but it simplifies the analysis.

Furthermore, we have to assume that the system, if DAE, is of order 1, such that the algebraic part can be solved and $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}$ can be used when evaluating ${\mathchoice{\displaystyle\boldsymbol{f}}{\textstyle\boldsymbol{f}}{\scriptstyle\boldsymbol{f}}{\scriptscriptstyle\boldsymbol{f}}}_{i}$, equivalently, that the system can be solved by the state space method [8, Section 3.1]. Such, arguments from the theory of ODE can be applied. This given, for scheme 1 the result [8, theorem 3.6] was proven:

Theorem 2.1**.**

Let S be a set of ODE which is split into disjoint subsystems $S_{k}$ of the shape

[TABLE]

the $\mathcal{D}_{k}$ and $\mathcal{I}_{k}$ denoting index sets. Disjoint subsystems means that the $\mathcal{D}_{k}$ are mutually disjoint and ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\cup\mathcal{D}_{k}}={\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$, and $\mathcal{D}_{k}$ and $\mathcal{I}_{k}$ are disjoint. Let ${T_{i}}$ be a time grid with width $H$, and let the inputs ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}_{k}}$ of all subsystems be extrapolated at $T_{i}$ with polynomial order $P$ and then be solved with an one-step method of order $p$ and maximal stepwidth $h\leq H$. Then for the error of the numerical solution ${\mathchoice{\displaystyle\boldsymbol{\epsilon}}{\textstyle\boldsymbol{\epsilon}}{\scriptstyle\boldsymbol{\epsilon}}{\scriptscriptstyle\boldsymbol{\epsilon}}}_{\Delta,S}(T_{j+1})$ the estimate

[TABLE]

holds.

The cosimulation scheme with balance correction reads as in Table 2,

using upper indices for the time interval and

[TABLE]

are added in the time step $j+1$. The correction that is applied in the $j$-th interval is then $\Delta E_{i}^{j-1}\phi_{j}(t)$.

For this method, [8, Theorem 4.1] tells that the convergence result theorem 2.1 still holds:

Theorem 2.2**.**

Let $S$, $S_{k}$, ${{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}_{\mathcal{D}_{k}}$, ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\mathcal{I}_{k}}$ and ${T_{i}}$ be as in theorem 2.1, but balance correction contributions be added to the extrapolated variables. Then the estimate from theorem 2.1 still holds.

Numerical tests and considerations indicate that balance correction methods improve the convergence order of explicit cosimulation schemes by one, but no proof has been written down for this yet.

2.2 Linear Problem for convergence examination

These results were confirmed numerically using the twodimensional linear problem

[TABLE]

which written as a cosimulation problem is

[TABLE]

The matrix entries are chosen such that

• •

an unidirectional dependency on input is given: $a_{11},a_{21},a_{22}\neq 0$, $a_{12}=0$,

• •

a mutual dependency is given: $a_{12},a_{21}\neq 0$, $a_{11}=a_{22}=0$.

Explicit and implicit extrapolation was used.

As ODE solver on subsystems, any solver that does not dominate the convergence and stability behavior of the cosimulation scheme could be used. Choice was vode and zvode from the numpy Python numerics library, which both implement implicit Adams method if problem is nonstiff and BDF if it is.

Figure 5 shows the convergence result for the four situations. To show the predictions made by (12) in terms of $H$, it is necessary that the subsystems methods contribution $\tilde{\tilde{C}}h_{\Delta}^{p}$ is of higher order than the extrapolation and that the method used is a one-step method, as many multistep methods use one-step methods of unknown order at start, and that would become visible here as at each exchange step requires one restart [8]. Thus dopri5, an explicite Runge-Kutta method, was chosen, using the built-in stepsize control with default absolute tolerance $10^{-12}$.

The figure shows that convergence is of order 1 for constant extrapolation and of order 2 for linear extrapolation, as predicted by (12). As discussed, there is no order loss for linear extrapolation and circular dependency of inputs, but not even an higher error, in spite of the negative effects that should occur.

The error of the first component of the lower triangular, thus unidirectionally coupled system is very low as it has no extrapolation contribution, thus indicating that the error made by $dopri$ is low enough to allow for judgement of the effect of extrapolation error.

2.3 Stability

But concerning stability, in [7] and [8] it is shown that cosimulation schemes are not stable for linear problems. The stability for linear problems replaces the notion of A-stability, as the one-component equation used there cannot be split. Linear stability is shown numerically as well as by arguing that the ODE (15) which is induced by extrapolating inputs is not stable [8, Section 3.5].

Balance correction cannot settle this issue: the linear spring-mass oscillator

[TABLE]

with mass $m$, spring constant $c$, and in which the damping constant $d$ shall vanish, was solved to numerically examine the stability of the method. This problem is an implementation of the linear problem (14) and the most simple problem possible that is linear and can be splitted, as it has the least number of components required for coupling, and coupling is its only contribution to derivatives, and this moreover linear. The eigenvalues of the uncoupled problem are purely imaginary, so the problem is stable.

It is solved with the cosimulation scheme given as above in this section, which yields 3. Output of the spring is the force $F=-cx$, that of the mass is the velocity $v=\dot{x}$.

The numerical convergence examination backs up the results Theorem 2.1 and 2.2. But the method is unstable for its explicit contributions, as proven in [8, Section 3.5]. This means $\left\|{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}\right\|\longrightarrow\infty$ for $t\longrightarrow\infty$. The energy of our system is $E=\frac{1}{2}mv^{2}+\frac{1}{2}cs^{2}=\left<{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}},{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}\right>_{\frac{1}{2}\operatorname{diag}(m,c)}$, which is an equivalent norm, so lack of stability is equivalent to energy augmentation.

So here, this lack of stability can be interpreted in physics as a consequence of extrapolation errors in factors of power acting on subsystems boundaries: In [7] it was pointed out that errors in the force $y_{1}$ made during data exchange lead to errors in the power that acts on the mass. The system picks up energy and behaves unstable (See figure 6).

In the following sections, it will be shown that balance correction techniques applied to the impulse as the integral of the force do not prevent this effect, as the a posteriori refeed of force then acts at another system state than it should, as the states have changed meanwhile - here the mass has changed its velocity.

3 Reclassification of Balance Errors

3.1 Errors in balances while exchanging conserved Quantities

These errors are defined as those that can be calculated as $\Delta E_{u_{ij}}^{k}=\int_{t_{k-1}}^{t_{k}}\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}dt-\int_{t_{k-1}}^{t_{k}}\overline{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}}^{k}dt$, thus those resulting directly from the extrapolation error in flow of the exchanged conserved quantity. Classical balance errors of extrapolations of $k$-th order arise on intervals where $dt^{k+1}f_{m}$ does not change signs. The balance error is often partially compensated during passing through intervals with positive and negative sign. This does not imply it is small during all intervals of the simulation.

As a typical such simulation situation think of an automotive driving cycle: Using piecewise constant extrapolation of exchanged data, there usually is a loss in conserved quantity after the phase of rising system velocity has passed - this loss remains uncompensated during the ( typically significant) phase of elevated speed, see e.g. [5]. This originally motivated the balance correction method [6].

3.2 Errors in balances while exchanging factors of conserved quantities

If an exchanged quantity influences a conserved quantity, the balance of that quantity is disturbed by the approximation error of an input as described above, even if the input quantity is nonconserved. Furthermore, as the conserved quantity depends on other factors, its imbalance may persist even if the balance of the exchanged quantity is reestablished, e.g. by errors compensating each other.

The prominent example of a globally conserved quantity that is exchanged between subsystems via its factors is energy. Think of an exchange of mechanical energy which in a cosimulation context is exchanged by passing displacement $s$ from $S_{1}$ to $S_{2}$ and the force $f$ vice versa. There is a power $p=vf$ acting at the interface between the two systems. Approximation errors on the receiving side cause an error in this power: Into the one subsystem flows an amount of energy different from the amount of energy leaving the other.

Of course, the received data influences the other factors of the energy, thus one cannot directly conclude from the input error to the energy error. It is correct, anyhow, that this influence is not a healing one in general.

3.3 Example for augmentation of energy

Oscillating systems are most suggestive to illustrate the nature of those energy-based balance errors as they accumulate energy rapidly and furthermore the effect of extrapolation errors is comparableto that of phase angles in electric engineering. Consider the simple linear mass-spring system (16) given by $\dot{{\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}}}=\tensor[]{A}{}{\mathchoice{\displaystyle\boldsymbol{z}}{\textstyle\boldsymbol{z}}{\scriptstyle\boldsymbol{z}}{\scriptscriptstyle\boldsymbol{z}}}$,

[TABLE]

The solution of an initial value problem with this system is

[TABLE]

or, equivalently,

[TABLE]

Thus the derivatives $\dot{x}$ and $\ddot{x}$ are sines and cosines. Choosing for simplicity $v_{0}=0$ and $t_{0}=0$ we have

[TABLE]

The mechanical force is $f=m\ddot{x}$.

This system was treated with a cosimulation scheme given in Table 2. Output of the spring is the force $f=-cx$, that of the mass is the velocity $v=\dot{x}$.

The power acting on the subsystem mass is $P=Fv=m\ddot{x}\dot{x}$. The energy received by the mass thus is $W=\int_{t_{0}}^{t_{1}}m\ddot{x}\dot{x}\,dt$, and as $P=\text{const}\cdot\operatorname{sin}\left(2\sqrt{c/m}t\right)$ energy is of the form $\text{const}\cdot\operatorname{cos}\left(2\sqrt{c/m}t\right)$. It vanishes at $t_{1}$ with $t_{1}-t_{0}=k\pi$, $k\in\mathbb{N}$, elsewhere it is bounded for all $t$.

If the split system is calculated with a cosimulation method and piecewise constant extrapolation of inputs, the force as it is seen by the mass is effectively shifted to the right. The analogy with the reactive power and the real power of an electrical network is apparent. As $\sin(\omega t)\cos(\omega t+\phi)=(\sin(2\omega t+\phi)+sin(-\phi))/2$, work from oscillating systems with phase shift contains the integral over a constant and thus grows unbounded. The plots in figure 6 clearly demonstrate this.

This kind of error occurs whenever a quantity inducing an energy on both sides of a subsystem boundary is exchanged across that boundary, making an approximation error, and thus energy appears on the one side without vanishing on the other side. Thus, it is also observed using degree 1 and higher extrapolation polynomials. For example, assuming linear Hermitean extrapolation, also derivatives of these quantities are exchanged. Let the displacement of the mass increase and the neighbouring system give a force against this movement as a response, i.e. let the mass move towards a spring and compress it. The force will decelerate the mass and reverse its movement. Now it is clear that due to the piecewise linear extrapolation of the data in this zone where it is convex (concave) a bigger displacement is assumed by the force delivering subsystem, making it respond with a higher force. The energy of the system is augmented - see the plots in figure 7.

Balance correction techniques applied to the impulse as the integral of the force do not prevent the system from picking up energy and behaving unstable, as the a posteriori refeed of force then acts at another system state than it should, as the states have changed meanwhile - here the mass has changed its velocity.

Mathematically, these balance errors are errors in arguments of a functional which would be conserved in the exact solution but is not in the cosimulation solution. This motivates our method that is functional conserving, which is presented in section 4.1 and which we finally examine for stability.

4 Enforcing Balance by sharing the view on Potential flow

It now is clear that balance correction methods can hardly stabilize systems that suffer from the effects described in 3 and 3.2 as it considers only an error made in the amount of a quantity, but the correction of amount is done at another, so wrong, time.

4.1 The method

The key feature to establish energy balance is exchanging the value of power and calculating the variable of interest from that power. Consider a cosimulation problem with subsystems $S_{1}$ and $S_{2}$ as given by equations (5) - (10) with states ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{1}$ and ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2}$ respectively and inputs ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21}$ and ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12}$. We suggest the following procedure to enforce energy balance between subsystems $S_{1}$ and $S_{2}$.

At data exchange timepoint $T_{n}$ the powers $P_{ij}$ as the flux of energy are calculated in both subsystems, using up-to-date input $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}^{n}_{ij}$. In general $P_{ij}\neq P_{ji}$, although in some situations equality holds. In the input vectors ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{ij}$ one component is replaced by $P_{ij}$ and that new vector is exchanged between the subsystems. Applied to $S_{1}-S2$ setting, the values $P_{21}$ and $P_{12}$ are exchanged, this means $S_{1}$’s point of view about the power is passed on to $S_{2}$ and vice versa. The ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12}$ and ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21}$ are also exchanged, but one component $(u_{ij})_{k}$ of each ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{ij}$ is omitted as it is now calculated from $P_{ij}$. 2. 2.

Now both subsystems have the same information and thus the opportunity to draw the same conclusion on what energy exchange should be assumed. We denote this assumed energy exchange as

[TABLE]

a straightforward choice is $\hat{P}_{21}=(P_{12}-P_{21})/2=-\hat{P}_{12}$, where now it is necessary to define flow directions: $P_{ij}$ shall be negative if it leaves $S_{j}$, so it is counted with opposite sign in $S_{i}$.

The former input $(u_{12})_{k}(t)$ is calculated subject to

[TABLE]

Indexing is like follows: $P_{21}$ is the power calculated in $S_{1}$ for passing to $S_{2}$, calculated using ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12}$, the input into $S_{1}$. Analogously $(u_{21})_{k}(t)$ s.t. $P_{12}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21\setminus l},(u_{21})_{l})=\operatorname{Ext}(\hat{P}_{21})$ is calculated. This requires that the maps $P_{ij}(.,.,(u_{ji})_{k})$ are monotone. The expression ${12\setminus k}$ in index is to say that the $k$-th component of the vector is left out.

Now it is established that the inputs of $S_{1}$ and $S_{2}$ are consistent in terms of energy conservation for all $t$.

Inversion of $P$ and its Notation

The solution of finding $({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{k}(t)$ subject to $P_{12}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{1}(t),{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12\setminus k},(u_{12})_{k}(t))=\operatorname{Ext}(\hat{P}_{12})$ will in the following be denoted as

[TABLE]

in analogy to the usual denotion of inverse functions by exponent $\cdot^{-1}$. This inversion in fact can be ill-conditioned in practice.

Precisely, the cosimulation scheme

[TABLE]

is replaced by

[TABLE]

Energy balance still holds when the method is extended to more than two subsystems as balance holds at each inter-subsystem boundary. Again, the expression ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{ij\setminus m}$ in index is to say that the $m$-th component of the vector ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{ij}$ is left out.

4.1.1 Idle negotiation in two subsystems situation

Given two subsystems with one interface between them, from physics of course $P_{21}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{1},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12}(x_{2}))=-P_{12}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{2},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21}(x_{1}))$ holds. In fact, $P_{12}$ and $P_{21}$ depend on the same variables $x_{1}$ and $x_{2}$ during cosimulation, and although the two powers may use different formula, they give equal result. In this case, the negotiating step begins with equal views and so is idle.

4.2 Example

To apply the scheme given in 4.1 above to the model of a spring-mass system as (16), replacing the standard cosimulation scheme from table 3, one first calculates the energies of the systems parts, powers acting on subsystems boundaries, and their derivatives. As $P_{i}=\dot{W}_{i}$, $P_{i}<0$ indicates that energy leaves $S_{i}$.

[TABLE]

$\dot{f}$ is available as output of spring, as usually serves as derivative of input. Now the following systems are treated:

Spring

Mass

System States

$x_{1}:=s=x$

$x_{2}:=v=\dot{x}$

Outputs

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21,\text{Std}})_{1}:=f=-cx$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12,\text{Std}})_{1}:=v=\dot{x}$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21,\text{Std}})_{2}:=\dot{f}=-cv$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12,\text{Std}})_{2}:=\dot{v}=f/m$

(intermediately exchanging ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{ij,\text{Std}}$)

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{1}=P(x_{1},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})=cxv=cx_{1}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{1}$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{1}=P(x_{2},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})=fv=x_{2}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{1}$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{2}=\dot{P}(x_{1},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})=c(v^{2}+xa)=c(({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{1}^{2}+x_{1}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{2})$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{2}=\dot{P}(x_{2},{\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})=m(a^{2}+v\frac{\dot{f}}{m})=m\left(\frac{({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{1}}{m}^{2}+x_{2}\frac{({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{2}}{m}\right)$

Inputs

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{1}:=\hat{P}$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{1}:=-\hat{P}$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{2}:=\hat{\dot{P}}$

$({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{2}:=-\hat{\dot{P}}$

Variables (inputs of standard method ${\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{\text{std}}$ ) depending on Inputs (Power)

$v=\frac{\operatorname{Ext}(\hat{P})}{cs}=\frac{\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{12})_{1}}{cx_{1}}$

$f=-\frac{\operatorname{Ext}(\hat{P})}{v}=\frac{\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{21})_{1}}{x_{2}}$

Equations

$\dot{x}_{1}=v$

$\dot{x}_{2}=\frac{f}{m}$

Remark that here

[TABLE]

have to be calculated. For exact values, those are well defined and bounded as $P\longrightarrow 0$ if $s\longrightarrow 0$ and if $v\longrightarrow 0$. On a computer, they are neither defined nor bounded. One has to switch to d’Hopitals rule for calculation near denominators zeros.

5 Stability of power balanced schemes

As discussed in Section 2.3 and shown in [8] , stability for linear systems of a partly explicite scheme is not given. This section shall relate energy conservation of the method to stability. We give an outline of this section:

• •

Switch to gradient flow view

• •

Introduce split system

• –

identify coupling contributions

• –

Characerize potential conservation/dissipation properties

• •

See method as decoupling ODE – Insert calculation of inputs from power into orig. equations

• •

relate decoupled ODEs stability properties to stability of original systems

• –

show that negotiated exchange conserves $\dot{\mathcal{P}}\leq 0$.

• –

use Lyapunovs direct method

• –

additionaly, one can argue that maximum stable stepwidth for dissipative systems is augmented (method is closer to B-stable than Extrapolation of Inputs).

• •

if such stable subsystems ODEs are solved with methods preserving that stability, overall solution will be stable.

5.1 Structure of coupled system

Many system’s behavior is governed by conservation of some energy or by entropy related potential minimization. This behavior is described by relating states time derivative to the functional gradient of the potential, $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\sim\nabla_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}^{T}$, usually linearly:

[TABLE]

The properties of the mobility matrix$\tensor[]{M}{}$ are, as

[TABLE]

connected to the properties of the system as follows:

for near-equilibrium potential driven processes, $\tensor[]{M}{}$ is symmetric due to the Onsager reciprocal relations, and has positive spectrum due to vicinity to a stable equilibrium point, so is positive definite. Such a system moves towards the potentials minimum, and

[TABLE]

accordingly $\leq 0$ if $\tensor[]{M}{}$ is pos. semi-def.. See e.g. [9] for a detailed description and example. 2. 2.

for systems that preserve a total energy. For conservation, the systems evolution has to be perpendicular to the functional gradient, which is

[TABLE]

so skew or implementing a dirac structure for conserving/Hamiltonian processes [10], [11]. Examples for those are mechanical systems, for example the micromechanical force balance, see [9, App 3], or all Hamiltonian systems and those that can be seen as such, e.g. spring-mass systems as equation (62) . 3. 3.

Usually

[TABLE]

as systems have conserving as well as dissipative properties.

Stability of methods is the conservation of the stability of the numerical solution of some IVP by the method. The class of gradient flow problems that are in some sense stable is relevant:

Definition 5.1** (stabilities).**

Let $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}$ be an equilibrium point of the ODE $\dot{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}=f(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})$ and $\phi^{t}\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$ the solution for the IVP with $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}(t_{0})=\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$. Then $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}$ is

• •

stable if $\forall\epsilon$ $\exists\delta>0$ $\left\|\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}-\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}\right\|<\delta$ $\Rightarrow\quad\left\|\phi^{t}\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}-\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}\right\|<\epsilon$ $\forall t$

• •

asymptotically stable if $\exists r:\left\|\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}-\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}\right\|<r$ $\Rightarrow\quad\lim_{t\longrightarrow\infty}\phi^{t}\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}=\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}$.

Stability of an IVP can be proven using

Theorem 5.2** (Lyapunovs direct method).**

*Let $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}$ be an eq. point of ODE.

Let $\mathcal{P}:V\longrightarrow\mathbb{R}^{+}:$ $\mathcal{P}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*})=0$, $\mathcal{P}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})>0$ $\forall\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}\neq\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}$ such that*

[TABLE]

Then $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}^{*}$ is a stable equilibrium point as defined in 5.1. If $\dot{\mathcal{P}}<0$, $x^{*}$ is an asymptotically stable eq. point.

which is proven in numerous higher analysis textbooks. It immediately follows by $\dot{\mathcal{P}}(x)\leq 0$ that gradient flow problems of potentials that are convex around a minimum and $\tensor[]{M}{}$ positive definite are stable in the sense of Lyapunov. Vice versa, it can be stated that near a stable equilibrium point, any ODE can be approximated by a gradient flow problem. Of course, all linear ODE can be seen as gradient flow problem with respect to a quadratic functional.

So gradient flow problems are not only widespread, but Lyapunov stable systems can be approximated by gradient flow problems. For all this, this problem class is useful for examining stabilities of the power balancing method.

5.1.1 Properties of Interaction in gradient flow problem seen as subsystems

Into the potential and its gradient flow description of the system $S$ the notion of subsystems is introduced: With $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}(t)\in\mathbb{R}^{N}\times\mathbb{R}$ being the states of $S$, let $I_{i}\subset\left\{1,2,...N\right\}$ for $i=1...k<N$ such that $I_{i}\cap I_{j}=\emptyset$ for $i\neq j$ and $\bigcup_{i}I_{i}=\left\{1,2,...N\right\}$.

We write

• •

$(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})_{I_{i}}$ for vector of components (indexing operator),

• •

$\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{i}}$ for vector that contains all components of $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$ that are in $I_{i}$ and 0 elsewhere.

The subsystems $S_{i}$ of $S$ now are given by the $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{i}}$ being their states. Without loss of generality, all $I_{i}$ consist of subsequent numbers. Further the gradient flow context implies that what was modelled as input before e.g. in (6) depends on states such that it can be expressed by the states $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{j}}$ and applying the replacement $\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}_{ij}=\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{j}}$ is no further loss of generality.

With this, we can identify subsystems energy gain as

[TABLE]

Using this notation, to identify the interaction between subsystems it is split into $I_{i}$-block lines

[TABLE]

and then splitting the block lines into blocks, mentioning only $k$ and $l$-containing expressions,

[TABLE]

we identify

[TABLE]

as the potential production in $S_{k}$ by $S_{l}$s variables, or power acting from subsystem $l$ onto subsystem $k$ and

[TABLE]

as $S_{k}$s internal potential change.

5.1.2 Properties of Potential Production contributions

Let, as discussed above, $\tensor[]{M}{}=\tensor[]{M}{}_{\text{symm}}+\tensor[]{M}{}_{\text{skew}}$. Due to (40), $P_{kk}$ is

• •

energy absorbing iff $\rho(\tensor[]{M}{}_{\text{symm},kk})\subset\mathbb{R}^{+}$ and $\exists\lambda_{i}>0$

• •

energy conserving iff $\tensor[]{M}{}_{kk}=\tensor[]{M}{}_{\text{skew},kk}$: It is then a Hamiltonian system with inputs.

If $\tensor[]{M}{}$ is symmetric, then

[TABLE]

which can be expressed in words as: The states adjacent to the $S_{l}$–$S_{k}$ boundary influence Potential on both subsystems equally. For example, when coupling two spatial domains with a PDE-governed field on them, this means that DOFs that neighbour the other domain evolve such that the potential production (e.g. entropy production) is equal for both domains.

If $\tensor[]{M}{}$ skew, then analogously

[TABLE]

This implements a real potential flow across the boundary: $S_{k}$ takes what $S_{l}$ loses. Those flows cancel out in the overall potential production (35).

If the mobility matrix has both nonzero $\tensor[]{M}{}_{\text{symm}}$ and $\tensor[]{M}{}_{\text{skew}}$ contribution, the potential flow has real flow contributions that are no productions and balanced production contributions contributions .

Remark 5.3* (Relation to Port Hamiltonian systems).*

The setting in which $S$ is given above makes it a Port-Hamiltonian System according to [10], which is a system of the shape

[TABLE]

with skew $\tensor[]{J}{}$ that corresponds to $\tensor[]{M}{}_{\text{skew}}$, and $\tensor[]{R}{}$ is symmetric positive definite and corresponds to $-\tensor[]{M}{}_{\text{symm}}$. If a $P_{kk}$ is

• •

energy absorbing, i. e. $\tensor[]{M}{}_{kk}$ is symmetric postive definite, then it is a dissipative port

• •

energy conserving, i.e. $\tensor[]{M}{}$ skew, it is a Port-Hamiltonian subsystem.

5.1.3 Structure of decoupled and power balanced system

We now derive the ODE induced by cosimulation scheme from the original ODE by replacing $S_{i}$s inputs $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{j}}$, $j\neq i$, by the explicit expression replacing them: From directly inserting $\operatorname{Ext}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{j}})$ into $\dot{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}=-\tensor[]{M}{}\nabla_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}$

[TABLE]

results, which is no gradient flow any more and in general will be instable ( see section 2.3).

Analogously, power balanced method induces an ODE:

[TABLE]

The index $j$ is, without loss of generality, used to denote the component of ${\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{I_{j}}$ that was replaced by $P_{ij}$. The argument $(\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{I_{j}\setminus j}),P_{ij}^{-({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}})_{j}}(\operatorname{Ext}\hat{P}_{ij}))_{j\neq i}$ is a vector, as in fact as many $j\neq i$ as couplings appear in the argument. For readability, $\operatorname{Ext}({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{I_{j}\setminus j})$ will be omitted from now on.

5.1.4 Stability inheritance

Before it is sketched how from properties of the method induced ODE the stability properties of our scheme is derived, remember the stability concepts:

First, from the dissipativity of an equations right hand side

[TABLE]

follows the non-expansivity of its evolution ${\mathchoice{\displaystyle\boldsymbol{\Phi}}{\textstyle\boldsymbol{\Phi}}{\scriptstyle\boldsymbol{\Phi}}{\scriptscriptstyle\boldsymbol{\Phi}}}^{\tau}{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$

[TABLE]

(see [12, th. 6.49 ]), which if inherited to numerical solution is the B-stability of the method ([12, Section 6.3]).

Second, remember A-Stability and its vector valued generalization, the linear stability as discussed in Section 2.3. Remember further Definition of Lyapunov stability 5.1.

We now use the following obvious arguing: Let original ODE

[TABLE]

Then, if ODE induced by method inherits this property, which means it still has a stable point at $x^{*}$/ is still linear and stable/dissipative, all subproblems $\dot{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{I_{i}}$ have that property.

Then all subproblems numerical solution $\Psi{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{I_{i}}$ is stable if solved with a

[TABLE]

Then the overall solution is stable and so the split method is

[TABLE]

if the methods applied on the subsystems are. Stability of power balanced scheme now is shown using Lyapunovs direct method, Theorem 5.2.

By applying this method one gains the stability result for the power balanced method:

Theorem 5.4**.**

For a Lyapunov stable (asymptotically stable) initial value problem (IVP), the IVP resulting from the energy balancing method is also stable (asymptotically stable).

Proof.

Inserting (45) for $\dot{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$,

[TABLE]

results, the last equality by identifying subsystem internal production and exchange: Potential flow, as always, is:

[TABLE]

and by construction of the power negotiating method, without loss of generality $P_{\text{bal},ij}=\hat{P}_{ij}$ and $P_{\text{bal},ji}=-\hat{P}_{ij}$ - the $x_{j}=P^{-({\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}})_{j}}(\operatorname{Ext}\hat{P}_{i,j})$ are determined with respect to that. So the second sum cancels out,

[TABLE]

as $P_{ii}\leq 0$ for all $i$. So theorem 5.2 can be applied. ∎

Remark

The cancellation $\sum_{i}\sum_{j\neq i}\hat{P}_{ij}=0$ due to $\hat{P}_{ji}=-\hat{P}_{ij}$ obviously is consistent to the cancellation of real flows due to $\tensor[]{M}{}_{skew}$ in the original model - its flows also balance according to (42): $P_{ji}=-P_{ij}$.

Is it consistent with original model for the flows due to the symmetric part of $\tensor[]{M}{}$ like in (41)?

Symmetric part is the potential production on subsystem boundary DOFs, flowing equally into both subsystems. In the original ODE, $P_{ji}=P_{ij}$, and from negative semidefiniteness of $\dot{\mathcal{P}}$, the restriction

[TABLE]

holds. In power balanced scheme, it is replaced by

[TABLE]

This means: Using power balanced decoupling scheme here is a change in the model – a bit of damping is lost – leading to a small but $O(1)-$Modeling error in this case!

5.2 Towards B-stability - inheritance of dissipativity to the method induced ODE?

For a system that is governed by (26) for convex $\mathcal{P}$ one can evaluate the dissipativity relation (46) as follows

[TABLE]

This is lower than 0 for all $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$ and $\overline{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$ if and only if $\tensor[]{M}{}D^{2}_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}$ is positive definite for all $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$, lower or equal 0 for positive semidefiniteness. The Hesse matrix $D^{2}_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}$ is positive definite as of a convex functional. The result is stated formally:

Lemma 5.5**.**

Given a system of ODEs on whose states $\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}$ an energy functional $\mathcal{P}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}})$ is defined, let the ODEs be defined by its gradient flow $\nabla_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}$ by

[TABLE]

If the right hand side $-\tensor[]{M}{}\nabla_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}$ is dissipative, i.e.

[TABLE]

*then $\tensor[]{M}{}\tensor[]{D}{}^{2}\mathcal{P}$ is positive semi-definite. Vice versa, if $\tensor[]{M}{}\tensor[]{D}{}^{2}\mathcal{P}$ is positive definite, the right hand side $-\tensor[]{M}{}\nabla_{{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}}\mathcal{P}$ is dissipative.

Inequality implies positive definiteness and vice versa.*

A big subclass of split gradient flow systems is given by systems $S=\cup{S_{i}}$ in which all couplings are due to $\tensor[]{M}{}_{\text{skew}}$, and all subsystems energies are due to subsystem variables only, which is $\mathcal{P}_{i}=\mathcal{P}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i})$. Example (63) gives a hint that the class of problems for which the assumptions are valid is relevant. For this class it holds that

[TABLE]

The first equality follows by $P_{kl}=-P_{lk}$ (42) from

[TABLE]

the second is equation (52) - saying, by $\mathcal{P}_{i}=\mathcal{P}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{i})$ and chancellation of flows, the potential production of the power balanced method is independent of extrapolations. Moreover, if written as scalar product using $P_{ii}=\left<\nabla_{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{i}}}\mathcal{P}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}),-\tensor[]{M}{}\nabla_{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}\mathcal{P}_{i}(\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}_{I_{i}})^{T}\right>$,

[TABLE]

which says that the projection of $\dot{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}$ and $\dot{\mathchoice{\displaystyle\boldsymbol{x}}{\textstyle\boldsymbol{x}}{\scriptstyle\boldsymbol{x}}{\scriptscriptstyle\boldsymbol{x}}}_{\text{bal}}$ onto the gradient of the potential is the same for the ODE and the balanced method induced ODE.

This projection property is independent of extrapolation, the property is given for all exchange step lengths if given for one.

As before, it here becomes obvious that for $\tensor[]{M}{}=\tensor[]{M}{}_{symm}+\tensor[]{M}{}_{skew}$ only $\tensor[]{M}{}_{symm}$ induces potential production, the skew part has no influence on it.

6 Comparison of numerical results and Discussion

6.1 Convergence

Simulations of the spring-mass problem (16) in its splitted form as described in Table 3, moreover with balance correction as in Table 2, [6], [7] and by negotiated power scheme (section 4.1) have been executed. Expressed as gradient flow problem, more precisely as a Hamiltonian system, as in (43) , it is written as

[TABLE]

where the energy of this system is

[TABLE]

and so

[TABLE]

As Hamiltonian system, this system has skew mobility matrix, is energy conserving and has a stable solution.

for the series of exchange step sizes $H={0.2,0.1,0.05...}$. As a subsystem solver dopri5 was used, a one-step method is necessary due to reasons that will follow. In this setting, the extrapolation error dominates $\epsilon_{\Delta}$.

Plot 9 shows that convergence rates predicted by (12) are met, also by the negotiated power method, which has extrapolation order 2. The error of the balance correction method behaves better than predicted by Theorem (2.2): it reduces approximately by 1/8 in each refinement step for linear extrapolation, and 1/4 for constant extrapolation. In the derivation of the estimate, the balance correction was treated as an arbitrary perturbation - in fact, it reduces the extrapolation error, and the numerical experiment suggests that this reduction leads to a gain of one in the order of convergence.

6.2 Stability

The numerical examination reveals that stability in practice is, albeit given for bigger $H$ than for all other coupled methods, not global. It turns out that this is due to the difficulties when calculating the inverse of Power, which is treated in section 6.3.1. There, it is also shown that without that phenomenon, the energy of the system would be conserved.

Figure 10 shows the solution of the spring-mass system simulated using $H=0.2$ and linear extrapolation without and with balance correction and (right) the power balanced scheme. The latter one is the only one that conserves the energy. This still holds for $H=0.15$, see Figure 11. For smaller $H$ all methods become stable.

6.3 Pitfalls

6.3.1 Inversion of $\hat{P}$

As mentioned in 4.2, the inverses $P_{i}^{-({\mathchoice{\displaystyle\boldsymbol{u}}{\textstyle\boldsymbol{u}}{\scriptstyle\boldsymbol{u}}{\scriptscriptstyle\boldsymbol{u}}}_{ij})_{j}}$ have to be calculated, which e.g. means finding $(u_{{\mathchoice{\displaystyle\boldsymbol{1}}{\textstyle\boldsymbol{1}}{\scriptstyle\boldsymbol{1}}{\scriptscriptstyle\boldsymbol{1}}}2,\text{Std}})_{j}$ such that

[TABLE]

• there, e.g. $P_{1}^{-v}=\frac{\operatorname{Ext}(\hat{P})}{cs}$ and $P_{2}^{-f}=\frac{\operatorname{Ext}(\hat{P})}{v}$. For exact values, those are well defined and bounded as $P\longrightarrow 0$ if $s\longrightarrow 0$ and if $v\longrightarrow 0$. On a computer, they might be undefined or arbitrary large. One has to switch to d’Hopitals rule for calculation near denominators zeros.

It turns out that the effect of the errors committed during solving (66) near zeros of the denominator is acting as an energy source for bigger step sizes and such is threatening stability, even if solved with d’Hopitals rule. Away from those zeros, energy of the system is conserved, and for systems where such divisions do not appear there will be no such unphysical sources. Furthermore, even for our case the power balanced method is stable for much bigger stepsizes than the formerly discussed Cosimulation methods.

7 Discussion, Conclusion and future work

The suggested method establishes balance of energy possibly at the cost of other balances. If we think of a system that is highly damped, this does not make sense. But if a system is undamped, the method enables applying cosimulation methods and implements an explicit but stable method. Moreover, the method has a clear interpretation in physics and can be implemented by anyone with understanding of the systems he wants to couple. For simulations in industrial research and development, the new method is a big step forward.

A future task could be to design explicit methods whose extrapolation is such that a variable connected to stability – as energy – is conserved.

8 *Acknowledgments

The author thanks Dirk Scharff for putting up the questions that led to this work.