Energy Conservation and Power Bonds in Co-Simulations: Non-Iterative Adaptive
Step Size Control and Error Estimation
Severin Sadjina, Stian Skjong, Eilif Pedersen
Department of Marine Technology, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway

Lars T. Kyllingstad

arXiv:1602.06434v2 [] 6 Oct 2016

SINTEF Fisheries and Aquaculture, NO-7465 Trondheim, Norway
Here, we study the flow of energy between coupled simulators in a co-simulation environment
using the concept of power bonds. We introduce energy residuals which are a direct expression of
the coupling errors and hence the accuracy of co-simulation results. We propose a novel EnergyConservation-based Co-Simulation method (ECCO) for adaptive macro step size control to improve
accuracy and efficiency. In contrast to most other co-simulation algorithms, this method is noniterative and only requires knowledge of the current coupling data. Consequently, it allows for
significant speed ups and the protection of sensitive information contained within simulator models.
A quarter car model with linear and nonlinear damping serves as a co-simulation benchmark and
verifies the capabilities of the energy residual concept: Reductions in the errors of up to 93 % are
achieved at no additional computational cost.

I.

INTRODUCTION

Engineering systems continue to increase in complexity, comprising of all sorts of physical phenomena, and
often characterized on vastly different time scales. At
the same time, development and manufacturing are constrained by an ever increasing pressure to keep costs low
and time-to-market periods short. Because of this, there
is a strong virtualization trend in engineering. Simulation methods play a crucial role in early product and
system development where they allow to identify and
avoid design flaws and risk potentials and replace timeintense and costly real testing. The automotive and the
maritime industries1–4 have recently, amongst others, increased their efforts to incorporate virtual prototyping
into all stages of the development cycle.
Simulator coupling, or co-simulation, has gained significant interest due to its attractiveness for scientific and
industrial applications: It allows for independent and parallel modeling, the efficient use and reuse of domain-specific models, software tools, and expert knowledge, and
significant simulation speed ups through parallelization
and the use of suitable solvers within each simulator.5
Co-simulation further enables the protection of intellectual property, which is often a strong requirement in industrial applications. Its basic idea is the inclusion of a
solver with each model in a loose coupling environment.
There are, however, two major challenges for this fullsystem simulation approach: stability and accuracy. A
general-purpose, easy-to-use and robust framework to estimate and control errors and ensure sufficient accuracy
and stability is lacking. Constant co-simulation step sizes
(macro step sizes) and simple synchronization algorithms
can still be considered state-of-the-art. In practice, the
step size is often chosen manually by trial-and-error, sacrificing both accuracy and efficiency, and uncertainty estimates are simply not available. Clearly, there is a
need for more elaborate concepts that also provide easy-

enough feedback about the quality of the co-simulation
results. This is especially true for industrial applications:
The use of commercial software will in most cases prohibit iterative coupling schemes or at least make them
computationally unfeasible. Furthermore, the protection
of simulator-internal data is typically a major concern
and implementation details about the model or the internal solver method are unavailable.
Unfortunately, this excludes the use of almost all advanced co-simulation schemes that have been proposed so
far: There exist a number of adaptive step size algorithms
for co-simulation6–12 which, based on a specific error estimation method, aim to choose an optimal macro step
size to exchange coupling data efficiently. However, practically all of these require simulator-internal data (such
as the time integration order) or the ability to revert to
a previous time point and redo an entire macro time step
(rollback). Mostly, they require both. One exception is
the step size controller proposed in Refs. 8 and 9 which
adaptively controls the macro step size by use of a noniterative predictor/corrector error estimator.
A different approach is to correct for the coupling errors directly. In the Interface Jacobian-based Co-Simulation Algorithm 13 , coupling conditions are solved iteratively with Newton’s method. While allowing for accurate solutions, it generally requires Jacobians along with
the coupling variables from each simulator at each macro
time step. This, together with its iterative nature practically prohibits any industrial use. An interesting noniterative co-simulation method is the Nearly Energy Preserving Coupling Element 14 . It is based on the realization that the extrapolation of input signals itself can be
considered an artificially introduced subsystem with its
own error contributions. These effectively change system
behavior and violate the conservation of (generalized) energy. Corrections to the coupling variables are then introduced in order to improve accuracy.
In the present paper, we take this approach in a dif-

2
ferent direction and directly study the flow of energy
throughout a co-simulation in detail using the concept
of power bonds. We show that, if coupling variables are
given in quantities whose product is a physical power—
such as force and velocity, pressure and flow rate, or voltage and current—the energy exchanges between simulators are directly accessible. Moreover, we investigate the
concept of energy residuals which directly alter the total energy of the overall coupled system due to the fact
that the individual subsystems are solved independently
of each other between communication points. Successfully ensuring that all energy flows between simulators
are (reasonably well) balanced is a strong indication for
the validity of the co-simulation results and a prerequisite
to stability. These findings are exploited in a novel Energy-Conservation-based Co-Simulation method (ECCO)
which is generally applicable for simulator coupling due
to its non-iterative nature and the fact that it is solely
based on the current coupling data.
The paper is organized as follows: In Section II, we
introduce the reader to power bonds and how they can
be used to conveniently study the flow and conservation
of energy in co-simulations. We also present the concept
of power and energy residuals as a direct consequence of
the inherent violation of the conservation of energy. Next,
we demonstrate their use for global error estimation to assess the quality of co-simulation results in Section III. A
quarter car co-simulation benchmark model is used in Section IV to demonstrate the performance of a non-iterative
adaptive step size controller based on energy conservation. The interesting effects of nonlinear damping and
different system reticulations on accuracy and efficiency
are also investigated for this benchmark model. Section V
provides a comparison to the predictor/corrector step size
controller proposed by Busch et al.8,9 . Finally, a conclusion is given in Section VI.

eral equation of continuity,
∂(x, t)
+ ∇ · j (x, t) = σ (x, t),
∂t

where (x, t) is the local energy density at position x and
time t, j is the energy flux, and σ is the net rate of energy dissipation (or creation) per unit volume and unit
time. Eq. (1) states that energy is locally continuous: A
change in the local energy density is always accounted for
by appropriate local energy fluxes and energy dissipation
(or creation). With the modular structure of complex
engineering systems and co-simulation environments in
mind, we are interested in already reticulated systems
where the transport of energy between subsystems and
the dissipation, creation, and distribution of energy inside of subsystems are typically restricted to discrete regions and strongly heterogeneous. Under these assumptions, we can cast the energy continuity equation into a
form much more suited for our purposes,15
−

X
α

Pα (t) =

X dEβ (t)
dt

β

−

X

Σγ (t).

(2)

γ

Here, Pα is the total energy flux (that is, the power transmitted) through the power port kα of a subsystem, Eβ is
the energy contained by a discrete sub-volume Vβ inside
the subsystem, and, analogously, Σγ is the rate of energy
dissipation or creation inside a discrete sub-volume Vγ .
Eq. (2) expresses the fact that any net flux of energy
into a subsystem is either stored or dissipated. In other
words, if we correctly account for all net power transmitted between the subsystems as well as the energy stored
and dissipated within the subsystems, we inevitably guarantee energy conservation and continuity.

S1
II.

(1)

Pk

S2

POWER BONDS AND ENERGY FLUXES

Considering energy balances directly, that is, working
with the Lagrangian or Hamiltonian equations is a poor
starting point for co-simulations because such approaches are based on the total energy of the complete system.
Instead, we shall follow the approach that is the theoretical foundation of bond graph theory15 and describe
the energy transactions between subsystems in terms of
local energy continuity equations.16 This way, we can balance the flow of energy for each subsystem separately and
thereafter connect them all in a modular fashion while,
in principle, satisfying the conservation of energy.

Figure 1. Two coupled subsystems, S1 and S2 , exchange energy through a power bond k at a rate Pk = ek fk which is
determined by the product of a flow fk and an effort ek

Naturally, we have little control over energy storage
and dissipation and whether they are correctly accounted
for inside simulators when they are coupled in a general
co-simulation environment. We can, however, safely conclude that no energy should get lost,17
− (Pk1 + Pk2 ) = 0,

A.

Energy Continuity

The transport of energy, that is, any energetic interaction between systems (or within systems) follows a gen-

(3)

while it is being transmitted between the power ports k1
and k2 of two coupled simulators S1 and S2 . The energetic coupling between two subsystems via power ports is
an example of a power bond, see Fig. 1. Such a bond is defined by a pair of power variables: a flow and an effort. If

3
we consider mechanical translation, for example, the flow
fk is a velocity and the effort ek is a force. The product
of two power variables is always a power, Pk = ek fk , and
gives the energy flow through the power bond (the rate
at which energy is being exchanged). This makes it practical to keep track of the flow and the conservation of
energy and is the basis for bond graph theory.

B.

The Power Bond in a Co-Simulation

In a co-simulation, the models include their own solvers
that evolve their internal states x between the discrete
communication time points ti and ti+1 ,

ẋ(t) = f x(t), ũ(t) , t ∈ (ti , ti+1 ],
(4a)
after having received inputs u from each other. The values of the input variables are generally unknown between
communication points, and an extrapolation ũ(t) ≈ u(t)
has to be used.18 Often, the inputs are simply held constant during the time integration, hence ũ(t) = u(ti ) for
t ∈ (ti , ti+1 ]. Outputs y are then computed from the
internal states,

y(ti+1 ) = g x(ti+1 ), ũ(ti+1 ) ,
(4b)
and used as inputs again.19 We can express the simulator
coupling as
u(ti+1 ) = Ly(ti+1 ),

uk 1 (ti )
S2
uk 2 (ti )

(5a)

where ũk1 and yk1 are the input and output, respectively.
Simulator S2 , on the other hand, will report that the
power it received from S1 through its power port k2 is
indeed
Pk2 (t) = ũk2 (t)yk2 (t),

(5b)

which leaves us with a predicament because generally
−(Pk1 + Pk2 ) , 0, in violation of Eq. (3).
Before we have a closer look at this discrepancy and
its consequences, let us first define an approximation for
the total power transmitted through the power bond k
from simulator S1 to simulator S2 using both simulator
outputs,

Pk12 (t) = σk12 yk1 (t)yk2 (t) .
(6)
Here, the sign σk12 ≡ (Lk12 − Lk21 )/2 is determined by
the corresponding elements of the connection graph matrix L. This is very useful in studying the flow of energy
throughout a co-simulation, as we shall see in Sec. IV.

C.

Residual Powers and Energies

The violation of energy conservation which we just uncovered in the coupling between two subsystems is, of
course, an inherent issue with co-simulation. It stems
from the fact that the simulators evolve their states independently of each other between communication time
instances. Luckily, it plays in our favor because it is also
a measure of how accurate the co-simulation solutions
are at any given time. In order to demonstrate this, let
us first define a residual power 20
− (Pk1 + Pk2 ) ≡ δPk

(a) Inputs are set at t = ti

yk 2 (ti+1 )
S1

Pk1 (t) = ũk1 (t)yk1 (t),

(4c)

where L is a connection graph matrix that relates outputs
and inputs at communication time points.

S1

power bond—this is easily done. Simulator S1 will conclude that the power it transmits to the other simulator
S2 through its power port k1 is given by

S2
yk 1 (ti+1 )

(b) Outputs are retrieved at t = ti+1 after time integrations

Figure 2. Two coupled simulators, S1 and S2 , exchange energy
through a power bond in a co-simulation

Now that we have set the stage, let us investigate the
energy flow through a power bond k between two simulators in a co-simulation setting, see Fig. 2. Assuming that
one variable (either the input or the output) represents
a flow and the other an effort—which is the case for a

(7)

for the power bond k, see Fig. 3 for an illustration. If we
construct input and output vectors such that
 
 
ũk1
yk1
ũk ≡
, yk ≡
,
(8a)
ũk2
yk2
the residual power can be conveniently calculated from
their scalar product,
δPk = −ũk · yk .

(8b)

If the bond between the simulators obeys energy conservation between communication points, the input and the
output vector are orthogonal to each other and the residual power vanishes, δPk = 0. As already mentioned, this
will not be the case in general, and energy either leaks
from (δPk < 0) or accumulates at (δPk > 0) the power

4
bond. Let us now calculate a corresponding residual energy
Z ti+1
δEk (ti+1 ) ≡
δPk (t) dt,
(9)

input extrapolation (i.e., the input values are held constant during the time integrations) for this reason.
Note that we can generalize Eqs. (8) and group inputs
and outputs together according to power bonds,

ti

which tells us how much energy was incorrectly created
due to the independent time integrations of the simulators during the macro time step ti → ti+1 = ti + ∆ti .
When introduced into the total coupled system, these
residual energy fluxes inevitably lead to altered energy
densities and dissipation inside the subsystems due to
the fact that energy is locally continuous, as expressed
by Eq. (2). They will, in other words, distort the dynamics of the system under consideration and in general
decrease the accuracy of the co-simulation and challenge
its stability. This is a profound statement as it tells us
how much energy is wrongfully added to the total energy
of the overall coupled system and provides an intuitive
explanation of why a co-simulation produces inaccurate
results or diverging solutions.
δPk
S1

Pk 1

Pk 2

S2

Figure 3. A residual power δPk = −(Pk1 + Pk2 ) emerges
and distorts the dynamics of the full system when energy is
exchanged between two simulators, S1 and S2 , in a co-simulation

Let us for now assume constant extrapolation of the
input values. The residual power for the power bond k
at the communication point ti+1 then becomes
δPk (ti+1 ) = −uk (ti ) · yk (ti+1 ).

(10a)

When calculating the residual energy, it is usually sufficient to make use of the rectangle quadrature rule,
δEk (ti+1 ) ≈ δPk (ti+1 )∆ti .

(10b)

In principle, higher-order corrections to the residual powers and energies can easily be calculated if the simulators
also output (and input) higher-order derivatives along
with the base values. This, however, depends on the individual simulator implementation and it may be very
impractical or even impossible to retrieve Jacobians. If,
on the other hand, inputs are approximated inside a simulator in terms of previous values and without the use
or knowledge of derivatives (e.g. by use of Lagrange or
Hermite polynomials9 , or Newton series), this will generally be unknown at the co-simulation level. Consequently,
Eqs. (10) are a fair choice to calculate residual energies
on a general basis, and the forms provided here should
therefore constitute a sufficient and practical approximation. We will continue under the assumption of constant

 
ũa1
ũa2 
 ũ 
b 
ũ = 
 ũ 1  ,
 b2 
..
.

 
ya1
ya2 
y 
b 
y=
y 1  ,
 b2 
..
.

(11a)

such that the total residual power of all power bonds
{a, b, . . . } of the co-simulation is obtained from
δP =

X

δPk = −ũ · y.

(11b)

k∈{a,b,... }

Consequently, energy conservation is satisfied throughout the entire co-simulation if the input and output vectors (11a) are orthogonal to each other. On the other
hand, Eqs. (11) allows to easily calculate the total residual energy which is incorrectly added to the full coupled
system due to co-simulation coupling errors.
The concept of inaccurate energy transactions discussed here is closely related to the theoretical footing of the Nearly Energy Preserving Coupling Element
(NEPCE) described in Ref. 14. There, the local error of a
simulator input and the resulting local error in the transmitted (generalized) power are discussed, albeit with sequential time integrations in the coupled simulators in
mind. In the present section, we aimed to cast a different light on the issue focusing in on parallel executions of
coupled simulators and energy conservation throughout
the entire coupled system.
We also intend to promote the use of power bonds
for simulator coupling: Power and energy transactions
and their accuracies then become directly accessible with
nothing more than the knowledge of the coupling variable
values. In this work, we only focus on the energetic coupling between two simulators. The entire formalism along
with all statements and results can easily be extended to
direct couplings between a general number of simulators.
This is most effectively done by use of bond-graph-like
junctions, and also keeps the modular character of cosimulations intact.

III.

ERROR ESTIMATION AND ADAPTIVE
STEP SIZE CONTROL

In Section II, we have seen how the fact that subsystems evolve their states independently of each other between discrete communication instances leads to the violation of energy conservation on the co-simulation level.
The flow of energy between coupled simulators was discussed and the idea of power and energy residuals as a
manifestation of this violation was introduced. It is now

5
time to reap the fruits of our efforts and apply these concepts to introduce a novel co-simulation global error estimator without the need to repeat macro time steps and
without any knowledge about the simulator implementations. With it, we shall further propose a non-iterative
adaptive macro step size controller based on energy conservation. But first, let us review how global errors are
usually estimated for co-simulations.
A.

Local Error Estimators

For simulator coupling without algebraic loops the global error is bounded in terms of local errors.11 This allows
us to find global error estimates in terms of local errors
as a direct extension of ODE and DAE integration methods. Error estimates are commonly based on comparing
˜ (ti+1 ), for the initial
two approximations, ỹ(ti+1 ) and ỹ
0
value problem, y = f (y, t), y(t0 ) = y0 , at time ti+1
in terms of an already obtained approximate solution at
times t ∈ {t0 , t1 , . . . , ti }. The difference between the two
˜ can then be used to estimate local
approximations ỹ − ỹ
errors and, consequently, derive a global error estimator.
Once a suitable estimator is found, it can also be used for
adaptive control of the step size to obtain a balance between accuracy on one hand and efficiency on the other.
In a co-simulation, there is naturally no access to the
subsystem equations and the state variables are unavailable for local error estimation. Instead, co-simulation
error estimators typically try to give an estimate of all
numerical errors in the simulator outputs when solving
Eqs. (4) with approximated inputs. Commonly, classic
error estimation approaches are adapted directly for cosimulation,21 and two approximations for the simulator
˜ , have to be obtained by carrying out
outputs, ỹ and ỹ
a macro step twice: An adaptive step size control by
Richardson extrapolation10,11 is realized by carrying out
a larger step ti → ti+2 with step size 2∆t and two smaller
steps ti → ti+1 → ti+2 with step size ∆t and comparing
the resulting output vectors. In an embedded methods
approach6 , the integrations are performed twice with input polynomials of varying degrees. Another approach7
involves performing a compound macro step for two simulators and repeating the same step for the stiffer of the
two subsystems using the updated coupling variables obtained from the first step. Again, the difference between
both results yields an error estimate. Finally, error estimation based on Milne’s device9 compares a predicted
with a corrected solution, both of the same polynomial
order. In practice, all of these approaches suffer from one
or several of the following shortcomings, and are not easily implemented for coupling with commercial simulation
tools:
1. They all require rollback, that is the ability to redo
a macro step. This is often either not supported at
all or complicated to realize in practice. Rollback
is especially prohibited for real-time applications,
such as Hardware-in-the-Loop.

2. In addition, re-stepping is a time-consuming procedure, especially when performed at every communication point. Consequently, it somewhat mitigates
the efficiency gain from adaptive step size control.
3. Lastly, information essential to the error estimator
method (for example, the internal time integration
order) is often not available. Moreover, this requirement may expose sensitive information whose concealment is one of co-simulation’s strong suits, at
least from an industrial perspective.
The error estimator introduced by Busch et al.8,9 is an
exception worth mentioning because it requires neither
re-stepping nor internal simulator information. We will
discuss it further in Sections III D and V.
B.

Error Estimation Based on Energy
Conservation

As we have seen, co-simulation coupling errors are directly expressed in terms of residual energies when power
bonds are used. These energy residuals give precisely the
amount of energy that is incorrectly added to the total energy of the overall coupled system, and should naturally
be contained in order to ensure co-simulation accuracy
and stability. Because of this, they are well suited as
versatile energy-based error estimators.
Consider then the local error with respect to the power
output of power port kα for the macro time step ti →
ti+1 = ti + ∆ti ,
∆Pkα (ti+1 ) = Pkα (ti+1 ) − Pk0α (ti+1 ),

(12)

where Pk0α and Pkα are the exact solution and the cosimulation result, respectively, where Pkα (ti ) = Pk0α (ti ) is
assumed. For input extrapolation of order m, the local error in the inputs is22,23 ∆uk = O(∆tm+1 ) for sufficiently
smooth problems and constant step sizes ∆ti = ∆t.
Then, from considering Eq. (4), the local error in the
states is ∆xk = O(∆tm+2 ) and, consequently, the local
error in the outputs is ∆yk = O(∆tm+2 ) if there is no
direct feed-through (i.e., no direct dependence on any inputs), and ∆yk = O(∆tm+1 ) otherwise. Using Eqs. (5)
and (12) then yields ∆Pkα = O(∆tm+1 ), irrespective of
the presence of direct feed-through.
Let us now take the sum of the local errors for both
power ports of a power bond k. We obtain δPk =
−(∆Pk1 + ∆Pk2 ), where we used Eqs. (3) and (7). Therefore, the residual power gives the average of the local
errors in the power,
1
∆Pk = − δPk .
(13a)
2
Similarly, the residual energy gives the average local error
in the energy transmitted during the macro time step
ti → ti+1 ,
1
∆Ek (ti+1 ) = − δEk (ti+1 ),
2

(13b)

6
where we carried out the time integration according to
Eq. (9). Note that Eqs. (13) are exact, irrespective of
macro or micro step sizes, and irrespective of the simulator implementations.
1000

C.

Adaptive Step Size Controller

Let us now define a scalar error indicator,
v
u
2
N 
u1 X
δEk (t)
 ,
(t) ≡ t
N
rk E0k + |Ek (t)|

(14)

k=1

DP

100

10

1

0.1
0.0001

0.001

0.01

0.1

Dt

Figure 4. Absolute error in the power (dashed) and power
residual estimate (solid) as a function of the macro step
size for the power bond in the linear benchmark model from
Sec. IV

The suitability of energy residuals as energy-conservation-based error estimators is exemplified in Fig. 4 for
the power bond in the co-simulation benchmark model
studied in Sec. IV. Shown is Eq. (13a) averaged over the
entire simulation time T ,
1 X
1
|δPk | =
|δEk (ti+1 )|,
2
2T i
against the average error in the power

which contains the residual energies δEk and the energies Ek (ti+1 ) ≈ Pk12 (ti+1 )∆ti transmitted per time step
for all N power bonds. Here, E0k is a typical energy
scale and rk a relative tolerance, respectively. Both can
be chosen individually per power bond and effectively
determine the bond’s energy resolution. If  ≤ 1, the
error is sufficiently small with respect to the defined tolerances and energy scales. Values of  > 1 indicate that
the tolerances are exceeded and the co-simulation result
is potentially unreliable. Applying, for example, a PIcontroller27,28 determines a new optimal step size based
on Eq. (14),
∆ti+1 = αs (ti )−kI −kP (ti−1 )kP ∆ti ,

(15)

where kI and kP are the integral and the proportional
gain, respectively, and αs ∈ [0.8, 0.9] is a safety factor.
We choose kI = 0.3/(m + 2) and kP = 0.4/(m + 2), where
m is the extrapolation order (m = 0 for constant extrapolation). It is beneficial to introduce user-defined bounds
for the step size, ∆tmin ≤ ∆ti+1 ≤ ∆tmax , as well as
limit its rate of change, Θmin ≤ ∆ti+1 /∆ti ≤ Θmax . Typical parameter choices are Θmin ∈ [0.2, 0.5] and Θmax ∈
[1.5, 5.0].
Table I. Configuration of the adaptive step size controller for
the benchmark model in Sec. IV
Value Unit

1X
|∆Pk | =
|∆Pk (ti+1 )|∆ti ,
T i
where ∆Pk (ti+1 ) = Pk12 (ti+1 ) − Pk012 (ti+1 ), and Pk012 is
the exact solution.
Using energies as error metrics as opposed to nonenergetic quantities (e.g. positions) has the advantage
that energy considerations are directly taken into account
as well and offers a more holistic approach.24 Outside the
realm of co-simulations, energy errors are used as a measure of quality in hybrid earthquake simulations25,26 or in
molecular dynamics simulations, for example. Additionally, the method proposed here naturally solves the issue
of the numerical values of the outputs lying on very different scales for force-displacement9,10,23 coupling, or the
force-velocity coupling we shall investigate in the next section: The outputs representing forces will typically have
much larger numeric values than the ones representing
displacements or velocities. As a result, local errors from
some simulators may be given weightings which are much
too large compared to others.

αs
0.8
∆tmin
0.1 ms
∆tmax 10.0 ms
Θmin
0.2
Θmax
1.5
E0
750.0 J

We have thus defined a novel co-simulation method for
error estimation and adaptive step size control based on
energy conservation considerations. It does not require
the repetition of macro time steps, nor does it demand
any knowledge about simulator implementations. In Section IV, we will demonstrate how it can be used to make
co-simulations more accurate and more efficient.

D.

Predictor/Corrector Error Estimator

As mentioned previously, another error estimation
and step size control co-simulation algorithm which is
non-iterative in nature is the one proposed by Busch

7
et al.8,9 . It makes use of polynomial extrapolation of
(r)
order r = m + 1 to construct a prediction ỹk (ti+1 ) of
the output vector (4b) using the previously obtained values {yk (ti ), . . . , yk (ti−r )}. The difference between this
predictor and the actually obtained output vector—the
corrector—gives a local error estimate,8,9
(r)

δ ỹk (ti+1 ) ≡ yk (ti+1 ) − ỹk (ti+1 ) ≈ ∆yk (ti+1 ),

(16)

and can be used for adaptive control of the macro step
size: A PI-controller (15) with kI = 0.3/(m + 1) and
kP = 0.4/(m + 1) proposes a new step size based on the
scalar error indicator8,9

equation for stability analysis of integration methods29
is a linear 1-DOF oscillator. Because it can be considered two coupled Dahlquist equations, the use of a linear
2-DOF oscillator as a co-simulation test case has been
proposed.8–11,23,24,30–32 This test system is also a simplification of numerous relevant real engineering systems.
One such example is the quarter car model33 depicted in
Fig. 5. We too shall adopt it as a co-simulation benchmark model and slightly modify it to include nonlinearity.

mc

S1

zc (t)

(r)

ykα (t) − ỹkα (t)
1
,
˜(t) ≡ max
α TOL 1 + ρ max |yk (t)|, |ỹk (t)|
α
α

where ρ sets the weight between the relative and absolute
errors and TOL is a user-defined tolerance.
This predictor/corrector approach and the energy-conservation-based method (ECCO) discussed in the previous sections are both easily implemented for non-iterative
co-simulation error estimation and adaptive step size control. But they are several noteworthy differences to be
found between them:
• ECCO gives the exact value of the average local
error in the power outputs of two connected power
ports, see Eqs. (13). This is true for any macro
or micro step size, and irrespective of the simulator implementations. By contrast, the predictor/corrector method may be able to provide good
estimates only for relatively small macro time step
sizes and is generally sensitive to the polynomial
order of the predictor.8,9
• ECCO is based on the use of power bonds, and
while we firmly believe that these are generally a
great choice for co-simulation—especially from an
engineering and industrial perspective—it is also
safe to assume that the vast majority of available
models and tools does not use power bonds as of
yet. In that sense, the predictor/corrector method
is less demanding and likely more applicable given
the status quo.
• As mentioned in Section III B, the use of powers
and energies also means that ECCO’s global error
estimator (14) is insensitive to the scaling of simulator output values, while Eq. (17) is not. This can
lead to practical difficulties in the implementation
of the predictor/corrector method, as we shall see
in Section V.
IV.

kc

(17)

CO-SIMULATION BENCHMARK MODEL

A generic test case is needed in order to examine and
compare the performance of various co-simulation methods. The mechanical interpretation of Dahlquist’s test

S2

dc

mw

zw (t)
kw

z(t)

Figure 5. The quarter car benchmark model is split into the
subsystems S1 and S2 for co-simulation

Table II. Parameters for the linear quarter car benchmark
model according to Ref. 11
Value

Unit

mc
400.0 kg
mw
40.0 kg
kc
15 000.0 N/m
kw 150 000.0 N/m
dc
1000.0 Ns/m
nd
0.5

The equations of motion describe the displacements zc
and zw of the masses mc and mw of the chassis and the
wheel, respectively, according to
mc z̈c (t) = −Fc (t),
mw z̈w (t) = −Fw (t) + Fc (t).
The spring-damper forces are given by11,34

Fc (t) = kc zc (t) − zw (t)

+ dc sign żc (t) − żw (t) żc (t) − żw (t)
Fw (t) = kw (zw (t) − z(t)),

(18a)
(18b)

2/(1+2nd )

,
(18c)

(18d)

where kc and kw are the spring constants, dc is the damping constant, and nd tunes the linearity of the damping
force. The parameter values are chosen as listed in Ta-

8
ble II. The system is excited with the function
(
0,
t < 0,
z(t) =
0.1, t ≥ 0,

A.

(19)

and the reference solution for the displacements of chassis
and wheel with the initial conditions zc = zw = żc =
żw = 0 at t = 0 is shown in Fig. 6.
0.15

zHtL

0.10

Adaptive Step Size Control

Let us now demonstrate the performance of the adaptive step size controller as described in Section III to improve co-simulation accuracy and efficiency. The residual
energy concept is used to propose an optimal macro step
size for the next time step in order to minimize residual energies and hence increase accuracy. As an example, consider the previous case of the linear quarter car
model. We use the scalar error indicator (14) and configure the PI-controller (15) according to the parameters
listed in Table I, with a starting step size ∆t0 = ∆tmin .
The energy scale is set to the initial energy excitation
E0 = 750 J of the system, and the balance between accuracy and computational efficiency can conveniently be
tuned with the relative tolerance r.

0.05

Table III. Linear quarter car benchmark results with constant
step size and residual-energy-based adaptive step size control
0.00
0

1

2

3

Algorithm

4

tolerance

t

constant
Figure 6. Displacements of wheel (solid) and chassis (dashed)
for the linear quarter car benchmark

S1 models the chassis only, described by Eq. (18a) with
the coupling variables
u1 (ti ) = −Fc (ti ),

y1 (ti+1 ) = żc (ti+1 ).

(20a)

S2 contains the remainder of the system and is described by Eq. (18b) with the coupling variables
u2 (ti ) = żc (ti ),

y2 (ti+1 ) = Fc (ti+1 ).

(20b)

The connections (20) define a power bond and we can
apply the residual energy concept. We use the forward
Euler method to carry out the time integration in S2 with
micro step sizes of ∆tS2 = ∆t/10, where ∆t is the macro
step size.35
The energy transactions between the simulators are
plotted in Fig. 7 for a constant macro step size of ∆t =
1 ms: Energy is first transmitted from S2 to S1 as the
potential energy stored in kw is transformed into kinetic
energy of mc . However, these energy transfers are not
accurate, and residual energy is accumulated over the
power bond and added to the total
R energy of the overall
system. This way, a total of δP (t) dt ≈ 6.4 J is incorrectly added to the coupled system during the entire
simulation. By comparison, the total amount of energy
initially stored in the spring kw and finally entirely dissipated in the damper dc is E0 = 750 J.

Error

P12
W

|∆P |
W

∫ δP
J

1.0

0.4

1.3

6.4

−6

1.0

0.0

0.4

1.6

−5

2.9

0.1

1.3

5.0

adaptive 2.8 × 10
adaptive 3.1 × 10

For the co-simulation, the system is split into two simulators, see Fig. 5:

Power

∆t
ms

Let us first demonstrate the improvements in accuracy
and set the tolerance to r = 2.8 × 10−6 such that the
average macro step size is ∆t ≈ 1 ms. While keeping the
calculation time approximately the same, the mean absolute error in the power |∆P | is reduced by 70 % relative to
the result with constant step size, see Table III. As Fig. 8
exemplifies, the residual energy is kept significantly better
R contained with adaptive step size control: A total of
δP (t) dt ≈ 1.6 J is added to the coupled system during
the entire simulation, only 25 % of the previous value.
Next, we increase the relative tolerance such that the
mean absolute error in the power is approximately the
same as with a constant step size of ∆t = 1 ms. With
r = 3.1 × 10−5 , a value of |∆P | ≈ 1.3 W is now reached
with only about 1/3 of the total time steps required without the algorithm. Energy flows are still described more
accurately than with constant step sizes, and the total
residual
energy accumulated during the simulation run
R
is δP (t) dt ≈ 5.0 J.

B.

Nonlinear Damping

Let us now investigate the effects of nonlinear damping
on the system dynamics and co-simulation accuracy. To
this end, we modify the parameters in Eq. (18c) according
to Table IV. The reference solution for the displacements
of chassis and wheel is shown in Fig. 9. Note that both
masses return to their equilibrium positions much faster

9
6

500

5
0
P HtL

Ù ∆P HtLât

4

-500

3
2

-1000

1
-1500
0.0

0.1

0.2

0.3

0.4

0.5

0
0.0

0.1

0.2

0.3

0.4

0.5

t
t

Figure 7. Energy transactions for the linear quarter car benchmark with a constant macro step size ∆t = 1 ms: power
transmitted from S1 to S2 (left) and the residual energy accumulated over time (right)

6

Table IV. Parameter changes to include nonlinear damping
forces in the benchmark model according to Ref. 34

5

Ù ∆P HtLât

Value
4

Unit

dc 900.0 N(s/m)1/2
nd
1.5

3
2
1
0
0.0

0.1

0.2

0.3

0.4

0.5

t

Figure 8. Accumulation of residual energy for the linear quarter car benchmark: Adaptive step size control with ∆t ≈ 1 ms
(solid) and constant step size ∆t = 1 ms (dashed)

than in the previous case with linear damping. Because
of this, we only simulate the system up to t = 2 s.
0.15

As was the case with linear damping, using the energybased adaptive step size control significantly improves
accuracy. Table V lists all results discussed here. With
the tolerance set to r = 7.5 × 10−6 , the computational
cost is kept about the same as with constant step sizes
(∆t ≈ 1 ms) while the errors are minimized and energy flows described more accurately, as demonstrated
by Fig. 10. The mean absolute error in the power and
the total residual energy wrongfully added to the coupled system during simulation time are both reduced by
70 % with step size control. Increasing the relative tolerance to r = 1.0 × 10−4 keeps the mean absolute error in
the power approximately the same but improves computational efficiency: It only takes about 1/3 of the time
steps to obtain a value of |∆P | ≈ 4 W.

zHtL

0.10

Table V. Nonlinear benchmark results with constant step size
and residual-energy-based adaptive step size control

0.05

Algorithm
tolerance

0.00
0.0

0.5

1.0

1.5

2.0

t

Figure 9. Displacements of wheel (solid) and chassis (dashed)
for the nonlinear quarter car benchmark

constant

adaptive 1.0 × 10

Error
|∆P |
W

∫ δP
J

P12
W

1.0

1

4

5

0.0

1.1

1.6

0

4

6

adaptive 7.5 × 10−6 1.0
−4

Power

∆t
ms

3.1

10
500
4
0
P HtL

Ù ∆P HtLât

3

2

-500

1
-1000
0.0

0.1

0.2

0.3

0.4

0.5

0
0.0

0.1

0.2

0.3

0.4

0.5

t
t

Figure 10. Energy transactions for the nonlinear benchmark: Power transmitted from S1 to S2 with a constant macro step size
∆t = 1 ms (left), and the residual energy accumulated over time for adaptive control of the macro step size with ∆t ≈ 1 ms
(right, solid) and with a constant step size ∆t = 1 ms (right, dashed)

mc

The connections (21) define a power bond and the residual energy concept can be applied. This time, time integration needs to be carried out in both simulators, and we
initially choose micro step sizes of ∆tS1 = ∆tS2 = ∆t/10.
Again, the forward Euler method is used in both simulators.

zc (t)

S1
kc

dc

mw

zw (t)
Table VI. Linear benchmark results for alternative system
reticulation with constant step size and residual-energy-based
adaptive step size control

S2
kw

z(t)

Algorithm
tolerance

C.

Alternative System Reticulation

The way a system is reticulated can have a substantial
influence on accuracy, stability, and efficiency of a cosimulation. To exemplify this and further demonstrate
the applicability and performance of the adaptive step
size controller described in Sec. III, we consider the quarter car model with a different subsystem splitting, see
Fig. 11:
S1 is described by Eq. (18a) but now also includes the
spring-damper element (18c) with the coupling variables
u1 (ti ) = żw (ti ),

y1 (ti+1 ) = Fc (ti+1 ).

(21a)

S2 is described by Eq. (18b) with the coupling variables
u2 (ti ) = −Fc (ti ),

y2 (ti+1 ) = żw (ti+1 ).

(21b)

Power
P12
102 W

1.0 −1.92

constant
Figure 11. Alternative reticulation of the quarter car benchmark model with subsystems S1 and S2

∆t
ms

−7

adaptive 9.1 × 10

Error
|∆P |
102 W

0.12

∫ δP
102 J

0.23

1.0 −1.879 0.013 0.016

As can be seen from the simulation results listed in
Tables VI and VII, this system reticulation is much less
favorable than the one discussed previously, giving larger
errors for linear and nonlinear damping alike. The co-simulation is also less stable: With linear damping, for example, the previous reticulation became unstable around
a constant step size of ∆t ≈ 58.5 ms, while instability
already sets in at around ∆t ≈ 11.3 ms with the reticulation discussed here. The effects of such an ill-chosen
reticulation can be mitigated and the accuracy enhanced
substantially by employing the adaptive step size controller. The errors are then reduced by between 80 %
and 93 % at no additional computational cost.
Another interesting observation can be made if we
greatly lower the accuracy of the time integration in S2
by choosing ∆tS2 = ∆t. As expected, the errors with
constant step sizes are considerably larger yet. Again,
however, using the adaptive step size control based on
the residual energy improves the situation significantly
and helps keep the accuracy of the results within reasonable bounds. As an example, consider the simulation

11
Table VII. Nonlinear benchmark results for alternative system
reticulation with constant step size and residual-energy-based
adaptive step size control
Algorithm
tolerance

∆t
ms

Power
P12
102 W

1.0 −3.9

constant

Table IX. Linear quarter car benchmark results with predictor/corrector (pred./corr.)
and residual-energy-based
(ECCO) step size control

Error
|∆P |
102 W

∫ δP
102 J

0.3

0.5

adaptive 2.4 × 10−5 1.0 −3.77 0.05 0.05

Algorithm

Algorithm
tolerance

∆t
ms

Power
P12
102 W

1.0 −2.2

constant
−6

adaptive 1.0 × 10

Error
|∆P |
102 W

∫ δP
102 J

0.4

0.3

|∆P |
W

∫ δP
J

1.0

0.4

1.3

6.4

pred./corr. 6.7 × 10−1 1.0

0.3

0.7

2.9

0.0

0.4

1.6

constant

2.8 × 10

−6

Algorithm

1.0 −1.90 0.04 0.02
constant
pred./corr. 2.1
ECCO

7.5 × 10

−6

Now that we have investigated the performance and
capabilities of a non-iterative energy-conservation-based
approach to error estimation and step size control, let us
in the present section deploy the quarter car benchmark
model one last time to draw a comparison with the predictor/corrector method proposed by Busch et al.8,9 . To
this end, we use the PI-controller (15) with the error indicator (17) as outlined in Sec. III D. The weight between
the relative and absolute errors is set to ρ = 1.0 × 10−4 ,
and the tolerance TOL is chosen according to the desired accuracy. All remaining parameters are configured
according to Table I. Because our main focus lies on cosimulation with constant input extrapolation, we set the
polynomial extrapolation order to r = 1.
The thus defined predictor/corrector method generally
yields significant error reductions throughout the various
benchmarks: Errors are reduced by 40 % to 55 % compared to the cases with constant step sizes, see Tables IX
and X. For the alternative system reticulation discussed
in Sec. IV C, a reduction of between 40 % and 93 % is
achieved, see Tables XI and XII. In comparison, the residual-energy-based adaptive step size controller generally
leads to a more substantial reduction in the energy errors,
though both approaches achieve relatively very high accuracies for the linear benchmark model with alternative

Power

Error
|∆P |
W

∫ δP
J

∆t
ms

P12
W

1.0

1

4

5

1.0

0.4

1.9

3.1

1.0

0.0

1.1

1.6

Table XI. Linear benchmark results for alternative system
reticulation with predictor/corrector (pred./corr.) and residual-energy-based (ECCO) step size control
Algorithm
tolerance

V. COMPARISON TO
PREDICTOR/CORRECTOR METHOD

1.0

Table X. Nonlinear benchmark results with predictor/corrector (pred./corr.)
and residual-energy-based
(ECCO) step size control
tolerance

results for the linear quarter car benchmark model as
listed in Table VIII, where the residual energy is reduced
by 90 %. The situation is similar with nonlinear damping
forces. This is rather beneficial for industrial applications
where one does not necessarily have access to internal simulator settings in order to satisfy accuracy and stability
demands on co-simulation results.

Error

P12
W

ECCO
Table VIII. Linear benchmark results for alternative system
reticulation with a low accuracy time integration in S2

Power

∆t
ms

tolerance

Power
P12
102 W

1.0 −1.92

constant
pred./corr. 6.0 × 10
ECCO

∆t
ms

−1

Error
|∆P |
102 W

0.12

∫ δP
102 J

0.23

1.0 −1.877 0.013 0.017

9.1 × 10−7 1.0 −1.879 0.013 0.016

Table XII. Nonlinear benchmark results for alternative system reticulation with predictor/corrector (pred./corr.) and
residual-energy-based (ECCO) step size control
Algorithm
tolerance
constant
pred./corr. 6.5
ECCO

2.4 × 10

−5

∆t
ms

Power
P12
102 W

Error
|∆P |
102 W

∫ δP
102 J

1.0 −3.9

0.3

0.5

1.0 −3.92

0.18

0.21

1.0 −3.77

0.05

0.05

system reticulation, see Table XI.
It should be noted that the error indicator (17) is sensitive to the scale of the output values: Recalling the
corresponding discussion in Sec. III B, the output values representing the spring force Fc (18c) are typically
several orders of magnitude larger than the ones representing the velocities żc and żw . This effectively completely disregards the error contributions of the simulator outputting the velocity in Eq. (17), and leads to an
inaccurate representation of the global co-simulation cou-

12
pling error. By using individual simulator-specific values
TOLα and ρα , this situation can be mitigated and higher
accuracies achieved. Doing so increases the number of
tunable parameters to four, however.
The particular benchmark cases at hand also reveal
that the actual step size selection suggested by the predictor/corrector method as implemented according to Ref. 8
and 9 is highly oscillatory. This is likely due to improper
tuning of the PI-controller, and has a detrimental impact
on the performance of the method for the problem set at
hand.36 It is also not unlikely that these issues are tied
to the fact that the error contribution of one simulator is
effectively disregarded.
In any case, we expect the predictor/corrector approach to yield a performance comparable to ECCO if
i.) the sensitivity issue with different scales of the outputs is resolved, and ii.) the PI-controller is correctly
tuned for the problem at hand.

VI.

CONCLUSION

When simulators are coupled via power bonds their inputs and outputs are given as so-called power variables
whose product is a physical power. Then, the flow of
energy throughout a co-simulation can be conveniently
studied. Moreover, because the subsystems are solved
independently of each other between discrete communication points, energy residuals emerge. These directly
alter the total energy of the overall coupled system and
distort its dynamic behavior.
In the present paper, we demonstrated that such residuals are easily computable with the simulator input and
output values only. Because they are a direct expression
of the coupling errors, they constitute a novel and versatile energy-based error estimation method. We showed
how energy residuals can be used for adaptive control of
the co-simulation step size. The performance and applicability of this non-iterative Energy-Conservation-based
Co-Simulation algorithm (ECCO) was investigated using
a quarter car benchmark model, both with linear and nonlinear damping characteristics. The proposed method
ensures that approximately the correct amount of energy is exchanged between subsystems. Consequently,
significant improvements in accuracy and efficiency were
demonstrated in comparison with constant co-simulation
step sizes. Additionally, using the proposed adaptive step
size control makes the accuracy of the global result less
dependent on the system reticulation or the accuracies of
the time integration methods inside the subsystems.
Unlike almost all other proposed co-simulation algorithms, ECCO makes do without two major restrictions:
i.) It does not require the repetition of entire co-simulation steps. Such rollback is often prohibited by use
of commercial software tools and computationally expensive. ii.) It does not require knowledge of any simulatorinternal information and, consequently, supports the protection of sensitive information and intellectual property

rights. With the traction that the Functional Mockup Interface (FMI)37,38 has been gaining over the past
years4,10–12,32 , the fulfillment of such properties is a major advantage. FMI is a tool independent standard for
the exchange and the co-simulation of dynamic models
that are, for most practical purposes, closed for inspection and modification. In addition, while FMI does support the repetition of macro steps and exposing derivatives along with the outputs, such traits are still rarely
found in models for co-simulation that are in use in the
industry. This makes ECCO especially attractive for industrial and engineering applications.
The same holds true for the predictor/corrector error estimation and step size control method proposed by
Busch et al.8,9 . It too showed significant improvements
in accuracy compared to using constant step sizes, albeit
less so than ECCO. This is likely due to its sensitivity to
the scaling of the output values and its improperly tuned
PI-controller, however. One major drawback of a predictor/corrector-based approach to error estimation is the
fact that it requires sufficiently small macro time steps
to begin with. Residual energies, on the other hand, are
an exact representation of the local errors in the power
outputs of two coupled simulators, irrespective of macro
or micro step size.
It should be noted that using power variables to realize
simulator coupling brings about three major advantages:
i.) For one, power and energy are the universal currencies
of physical systems, and power bonds make them directly
accessible. ii.) Power bonds provide a complete and universal, energy-flow-centered connectivity between mathematical models of different engineering and physical domains. iii.) Power variables typically represent higher
derivatives than the coupling variables of other coupling
schemes. In general, this is a favorable trait if we keep in
mind that numerical integration is much more desirable
than numerical differentiation. On the other hand, very
few models and tools use power bonds to date, which
may make other methods more advantageous in practice,
at least for the time being.
An issue we only touched upon briefly is the stability
of co-simulations. Due to the complexity of the topic
(different solver methods of various orders, different simulator coupling schemes of various orders, the presence of
direct feed-through and algebraic loops, etc.), a general
approach is far from trivial and theoretical treatments
lack behind practical implementations. For further details see, for example, Refs. 9, 22, 31, 39, and 40, and
references therein.

ACKNOWLEDGMENTS

This work was funded by the Research Council of Norway (grant number 225322), and the industrial partners
in the ViProMa project consortium (VARD, Rolls-Royce
Marine, and DNV GL). We are grateful for their financial
support.

13

1

G. Henry, I. Cox, P. Crossland, and J. Duncan, “Virtual
ships: Nato standards development and implementation,”
report, Systems Engineering and Assessment LTD Frome
(United Kingdom), 2009.
2
S. Harries, F. Tillig, M. Wilken, and G. Zaraphonitis, “An
integrated approach for simulation in the early ship design of a tanker,” in 10th International Conference on
Computer and IT Applications in the Maritime Industries,
pp. 411–425, Technische Universität Hamburg-Harburg,
2011.
3
H. Palm, J. Holzmann, R. Klein, S.-A. Schneider, and
D. Gerling, “A novel approach on virtual systems prototyping based on a validated, hierarchical, modular library.”
Lecture, 2013. Embedded World, Nuremberg, Germany.
4
C. Bertsch, E. Ahle, and U. Schulmeister, “The functional
mockup interface - seen from an industrial perspective,” in
10th International Modelica Conference, Linköping Electronic Conference Proceedings, pp. 27–33, Linköping University Electronic Press, 2014.
5
Note that we use the word ‘simulator’ throughout in the
sense of a subsimulator: a mathematical model of a subsystem coupled to other such models of subsystems to form a
full model of a total system.
6
M. Günther, A. Kværnø, and P. Rentrop, “Multirate partitioned runge-kutta methods,” BIT, vol. 41, no. 3, pp. 504–
514, 2001.
7
A. Verhoeven, B. Tasic, T. G. J. Beelen, E. J. W. ter Maten,
and R. M. M. Mattheij, “Bdf compound-fast multirate
transient analysis with adaptive stepsize control,” Journal
of Numerical Analysis, Industrial and Applied Mathematics, vol. 3, no. 3-4, pp. 275–297, 2008.
8
M. Busch and B. Schweizer, “An explicit approach for controlling the macro-step size of co-simulation methods,” in
Proceedings of the 7th European Nonlinear Dynamics Conference, 2011.
9
M. Busch, Zur effizienten Kopplung von Simulationsprogrammen. Dissertation, Universität Kassel, 2012.
10
T. Schierz, M. Arnold, and C. Clauß, “Co-simulation with
communication step size control in an fmi compatible master algorithm,” in 9th International Modelica Conference,
Linköping Electronic Conference Proceedings, pp. 205–214,
Linköping University Electronic Press, 2012.
11
M. Arnold, C. Clauss, and T. Schierz, “Error analysis and
error estimates for co-simulation in fmi for model exchange
and co-simulation v2.0,” Archive of Mechanical Engineering, vol. 60, no. 1, pp. 75–94, 2013.
12
A. Viel, “Implementing stabilized co-simulation of strongly
coupled systems using the functional mock-up interface 2.0,” in 10th International Modelica Conference
(H. Tummescheit and K.-E. Årzén, eds.), Linköping Electronic Conference Proceedings, pp. 213–223, Linköping
University Electronic Press, 2014.
13
S. Sicklinger, V. Belsky, B. Engelmann, H. Elmqvist,
H. Olsson, R. Wüchner, and K. U. Bletzinger, “Interface
jacobian-based co-simulation,” INT J NUMER METH
ENG, vol. 98, no. 6, pp. 418–444, 2014.
14
M. Benedikt, D. Watzenig, J. Zehetner, and A. Hofer,
“Nepce—a nearly energy preserving coupling element for
weak-coupled problems and co-simulation,” in V International Conference on Computational Methods for Coupled
Problems in Science and Engineering, Coupled Problems,

pp. 1021–1032, International Center for Numerical Methods in Engineering, 2013.
15
H. M. Paynter, Analysis and Design of Engineering Systems: Class Notes for M.I.T. Course 2.751. Boston: M.I.T.
Press, 1961.
16
P. C. Breedveld, Physical systems theory in terms of bond
graphs. Dissertation, Twente University, 1984.
17
This is true if we assume an ideal bond that does not include physical energy sources or sinks. This is well justified
in a co-simulation setting where (relevant) energy creation
and dissipation ought to be accounted for inside subsystems. Our discussion is not limited by this choice, however,
and it is fully possible to work with non-ideal bonds.
18
Note that some coupling schemes interpolate input values.
In the Gauss-Seidel iteration pattern, for example, one simulator is stepped first using input value extrapolation, and
thereafter the second simulator is stepped using interpolation of the other system’s outputs.
19
The functions f and g may also depend on time explicitly.
20
Note that the sign convention for the residual power is
different than the one used for energy fluxes encountered
so far. A positive residual power means that the energy of
the total system increases.
21
For a more comprehensive overview see, for example, Ref. 9
and references therein.
22
M. Arnold, Multi-Rate Time Integration for Large Scale
Multibody System Models, vol. 1 of IUTAM BOOKSER,
book section 1, pp. 1–10. Springer Netherlands, 2007.
23
M. Arnold, S. Hante, and M. A. Köbis, “Error analysis for
co-simulation with force-displacement coupling,” Proceedings in Applied Mathematics and Mechanics, vol. 14, no. 1,
pp. 43–44, 2014.
24
F. González, M. A. Naya, A. Luaces, and M. González,
“On the effect of multirate co-simulation techniques in the
efficiency and accuracy of multibody system dynamics,”
MULTIBODY SYST DYN, vol. 25, no. 4, pp. 461–483,
2011.
25
C. Thewalt and M. Roman, “Performance parameters for
pseudodynamic tests,” J STRUCT ENG-ASCE, vol. 120,
no. 9, pp. 2768–2781, 1994.
26
G. Mosqueda, B. Stojadinovic, and S. Mahin, “Real-time
error monitoring for hybrid simulation. part i: Methodology and experimental verification,” J STRUCT ENGASCE, vol. 133, no. 8, pp. 1100–1108, 2007.
27
K. Gustafsson, M. Lundh, and G. Söderlind, “A pi stepsize
control for the numerical solution of ordinary differential
equations,” BIT, vol. 28, no. 2, pp. 270–287, 1988.
28
P. M. Burrage, R. Herdiana, and K. Burrage, “Adaptive
stepsize based on control theory for stochastic differential
equations,” J COMPUT APPL MATH, vol. 170, no. 2,
pp. 317–336, 2004.
29
G. Dahlquist, “Convergence and stability in the numerical integration of ordinary differential equations,” MATH
SCAND, vol. 4, pp. 33–53, 1956.
30
B. Schweizer and D. Lu, “Predictor/corrector cosimulation approaches for solver coupling with algebraic
constraints,” ZAMM-Z ANGEW MATH ME, vol. 95, no. 9,
pp. 911–938, 2014.
31
B. Schweizer, P. Li, and D. Lu, “Explicit and implicit cosimulation methods: Stability and convergence analysis for
different solver coupling approaches,” J COMPUT NON-

14

32

33

34

35

36

37

LIN DYN, vol. 10, no. 5, p. 051007, 2014.
C. Clauß, M. Arnold, T. Schierz, and J. Bastian, “Master zur simulatorkopplung via fmi,” in Tagungsband
der ASIM/GI-Fachgruppen STS und GMMS, pp. 57–70,
ARGESIM-Verlag, 2012.
K. Popp and W. Schiehlen, Ground Vehicle Dynamics.
Springer Berlin Heidelberg, 2010.
J. Bußhardt and R. Isermann, “Realization of adaptive
shock absorbers by estimating physical process coefficients
of a vehicle suspension system,” in American Control Conference, 1992, pp. 531–535, 1992.
S1 can be solved exactly with constant extrapolation of the
input Fc (t) = Fc (ti ) for t ∈ (ti , ti+1 ].
Setting kP = 0 and substantially lowering kI , for example, seems to be a much better choice to configure the
PI-controller in this case.
T. Blochwitz, M. Otter, M. Arnold, C. Bausch,
H. Elmqvist, A. Junghanns, J. Mauss, M. Monteiro, T. Neidhold, D. Neumerkel, H. Olsson, J.-V. Peetz, S. Wolf, and
C. Clauß, “The functional mockup interface for tool in-

38

39

40

dependent exchange of simulation models,” in 8th International Modelica Conference (C. Clauß, ed.), Linköping Electronic Conference Proceedings, pp. 105–114, Linköping
University Electronic Press, 2011.
T. Blochwitz, M. Otter, J. Åkesson, M. Arnold, C. Clauss,
H. Elmqvist, M. Friedrich, A. Junghanns, J. Mauss, and
D. Neumerkel, “Functional mockup interface 2.0: The standard for tool independent exchange of simulation models,” in 9th International Modelica Conference (M. Otter
and D. Zimmer, eds.), Linköping Electronic Conference
Proceedings, pp. 173–184, Linköping University Electronic
Press, 2012.
R. Kübler and W. Schiehlen, “Two methods of simulator
coupling,” MATH COMP MODEL DYN, vol. 6, no. 2,
pp. 93–113, 2000.
S. Skjong and E. Pedersen, “The theory of bond graphs
in distributed systems and simulations,” in 2016 International Conference on Bond Graph Modeling and Simulation, SIMUL SERIES, pp. 147–159, Society for Modeling
and Simulation International, 2016.

