# F3ORNITS : a flexible variable step size non-iterative co-simulation method
handling subsytems with hybrid advanced capabilities

Yohan ÉGUILLON1 , Bruno LACABANNE1 and Damien TROMEUR-DERVOUT2
1Siemens Industry Software, Roanne, France
2Institut Camille Jordan, Université de Lyon ,UMR5208 CNRS-U.Lyon1,
Villeurbanne, France
{yohan.eguillon<EMAIL_ADDRESS>damien.tromeur-dervout@univ-
lyon1.fr
https://orcid.org/0000-0002-9386-4646 https://orcid.org/0000-0003-1790-3663
https://orcid.org/0000-0002-0118-8100

###### Abstract

This paper introduces the F3ORNITS non-iterative co-simulation algorithm in
which F3 stands for the $3$ flexible aspects of the method: flexible
polynomial order representation of coupling variables, flexible time-stepper
applying variable co-simulation step size rules on subsystems allowing it and
flexible scheduler orchestrating the meeting times among the subsystems and
capable of asynchronousness when subsystems’ constraints requires it. The
motivation of the F3ORNITS method is to accept any kind of co-simulation
model, including any kind of subsystem, regardless on their available
capabilities. Indeed, one the major problems in industry is that the
subsystems usually have constraints or lack of advanced capabilities making it
impossible to implement most of the advanced co-simulation algorithms on them.
The method makes it possible to preserve the dynamics of the coupling
constraints when necessary as well as to avoid breaking $C^{1}$ smoothness at
communication times, and also to adapt the co-simulation step size in a way
that is robust both to zero-crossing variables (contrary to classical relative
error-based criteria) and to jumps. Two test cases are presented to illustrate
the robustness of the F3ORNITS method as well as its higher accuracy than the
non-iterative Jacobi coupling algorithm (the most commonly used method in
industry) for a smaller number of co-simulation steps.

##  1 INTRODUCTION

Co-simulation consists in processing a simulation of a modular model, that is
to say a model composed of several dynamical subsystems connected together.
This is usually motivated by the need to simulate a model with multiphysical
parts. Designing a subsystem representing the physics of a given field
(electricity, mechanics, fluids, thermodynamic, etc) allows the use of a
specific and adapted solver for this field, or even the modelling with a
specific third party software. Regarding industrial applications, a modular
model is prefered because subsystem providers can focus on a part of the
global system without taking the rest into account. Nonetheless, gathering
different subsystems is not straightforward: a simulation of the equations of
the global system cannot be retrieved. Co-simulation is the field
investigating the ways to solve such systems, based on regular data
communications between the subsystems solved separately.

The co-simulation method (or co-simulation algorithm) is the rule used to
process the simulation on such modular systems. It namely deals with: the
determination of the times of the data communications, the way the inputs each
subsystem should use at each step are computed, and the way the outputs are
used and so on. Many co-simulation algorithms have been established until now
[Kübler and Schiehlen, 2000] [Arnold and Unther, 2001] [Gu and Asada, 2004]
[Bartel et al., 2013] [Sicklinger et al., 2014] [Busch, 2016] and studied [Li
et al., 2014] [Schweizer et al., 2016], with various complexity of
implementation, subsystems capabilities requirements, or physical field-
specific principles (see also the recent state of art on co-simulation of
[Gomes et al., 2018]). It may appear in some co-simulation algorithms that a
time interval should be integrated more than once on one or more subsystems.
This is called iterative co-simulation. The systems that need to do so must
have the capability to replay a simulation over this time-interval with
different data (time to reach, input values, …). Depending on the model
provider, this capability, so-called ”rollback”, is not always available on
every subsystem. When it is available, an iterative algorithm such as
[Éguillon et al., 2019] can be used. Nevertheless, if at least one subsystem
cannot rollback, a non-iterative method has to be chosen. The rollback
capability is not available in most industrial models. As the F3ORNITS method
presented in this paper is non-iterative, it can be applied to configurations
with rollback-less subsystems.

Besides the rollback, other plateform-dependent capabilities may lead to an
impossibility of use of a given co-simulation method on certain subsystems.
Amongst them, the ability to provide inputs with $n$th order time-derivatives
and the ability to obtain the time-derivatives of the outputs at the end of a
macro-step can be mentioned. The motivation of the F3ORNITS method is to
accept any kind of modular model, including any kind of subsystem, regardless
on their available capabilities. The first consequence of this specification
is that F3ORNITS is non-iterative (in order to accept rollback-less
subsystems) and asynchronous (in order to accept modular models with
subsystems with an imposed step, even when several subsystems with imposed
steps non multiple of one another). Such a method is presented on figure 1.
F3ORNITS also represents the inputs as time-dependent polynomials (when
supported by the subsystems concerned), and adapts the communication step size
of the subsystems which can handle variable communication step sizes. This
both allows accuracy win when frequent exchanges are needed, and saves time
when coupling variables can be represented with a high enough accuracy.

Figure 1: Visualization of the behavior of a non-iterative asynchronous co-
simulation method on $2$ subsystems

This algorithm is based on a variable order polynomial representation of every
input (the polynomial degree might be different for each input variable at a
same time) determined by an a posteriori criterion and redefined at each
communication step. A smoother version can be triggered by interpolating on
the extrapolated values, as done in [Busch, 2019] with the so-called
”EXTRIPOL” technique. This technique avoids the non-physical jumps on coupling
variables at each communication time, and may help the solvers of the
subsystems to restart faster after a discontinuity since we can guarantee the
$C^{1}$ smoothness of the input variables. The F3ORNITS method adapts this
smoothening with flexible order polynomials, variable step size subsystems,
and asynchronous cases.

The paper is structured as follows. Section 2 presents the mathematical
formalism we used to develop the co-simulation algorithm. This formalism is
made of general notations, of subsystems topologies detailing and of common
polynomial techniques used in F3ORNITS method. Section 3 presents the F3ORNITS
algorithm. On the one side, it shows the way time-dependent inputs are
determined. On the other side the way the step size determination is handled.
Section 4 gives the results of F3ORNITS method both on a controlled speed
model and on the classical linear two-mass oscillator test case [Schweizer et
al., 2016] [Éguillon et al., 2019]. Comparison of F3ORNITS algorithm with
different options and basic non-iterative Jacobi method (most basic non-
iterative co-simulation technique with fixed step size) is also achieved. The
conclusion is given in section 5.

##  2 MATHEMATICAL FORMALISM AND MOTIVATIONS

A subsystem that communicates with other subsystems to form a co-simulation
configuration will be represented by its equation. As the context is about
time integration, these equations will be time differential equations. First
of all, we present the general form of a monolithic system, in other words a
closed system (which neither has inputs nor outputs) represented by a $1$st
order differential equation which covers, among others things, ODE, DAE and
IDE as follows:

$\left\\{\begin{array}[]{lll}F(\frac{d}{dt}x,x,t)&=&0\\\
x(t^{[\text{init}]})&=&x^{[\text{init}]}\\\ \end{array}\right.$ (1)

where

$\begin{array}[]{l}n_{st}\in\mathbb{N}^{*},\
x^{[\text{init}]}\in\mathbb{R}^{n_{st}},\
t\in[t^{[\text{init}]},t^{[\text{end}]}],\\\
{[}t^{[\text{init}]},t^{[\text{end}]}]\subset\mathbb{R}\text{ so that
}t^{[\text{end}]}-t^{[\text{init}]}\in\mathbb{R}_{+}^{*},\\\
F:\mathbb{R}^{n_{st}}\times\mathbb{R}^{n_{st}}\times[t^{[\text{init}]},t^{[\text{end}]}]\rightarrow\mathbb{R}^{n_{st}},\\\
\end{array}$ (2)

are given, and

$x:[t^{[\text{init}]},t^{[\text{end}]}]\rightarrow\mathbb{R}^{n_{st}}$ (3)

is the solution state vector whose $n_{st}$ components are called the state
variables.

This paper will only cover the ODEs case, which are of the form (4). The terms
”equation”, ”equations”, and ”differential equations” will refer to the ODE of
a system throughout this document.

$\frac{d}{dt}x=f(t,x)\\\ $ (4)

where

$f:[t^{[\text{init}]},t^{[\text{end}]}]\times\mathbb{R}^{n_{st}}\rightarrow\mathbb{R}^{n_{st}}$
(5)

###  2.1 Framework and notations

In a cosimulation context, the principle of linking the derivatives of the
states to themselves (and potentially with time) is the same as in the
monolithic case, yet the inputs and outputs also have to be considered.

Let $n_{sys}\in\mathbb{N}^{*}$ be the number of subsystems. Please note that
the case $n_{sys}=1$ corresponds to a monolithic system. The cases that will
be considered here are connected subsystems, that is to say $n_{sys}\geqslant
2$ subsystems that need to exchange data to one another: each input of each
subsystem has to be fed by an output of another subsystem. A subsystem will be
referenced by its subscript index $k\in[\\![1,n_{sys}]\\!]$ so that subsystem-
dependent functions or variables will have an index indicating the subsystem
they are attached to.

Considering the inputs and the outputs, the co-simulation version of (4) - (5)
for subsystem $k\in[\\![1,n_{sys}]\\!]$ is:

$\left\\{\begin{array}[]{lll}\displaystyle{\frac{d}{dt}}x_{k}&=&f_{k}(t,x_{k},u_{k})\\\
y_{k}&=&g_{k}(t,x_{k},u_{k})\\\ \end{array}\right.$ (6)

where

$\begin{array}[]{l}x_{k}\in\mathbb{R}^{n_{st,k}},\
u_{k}\in\mathbb{R}^{n_{in,k}},\ y_{k}\in\mathbb{R}^{n_{out,k}}\\\
f_{k}:[t^{[\text{init}]},t^{[\text{end}]}]\times\mathbb{R}^{n_{st,k}}\times\mathbb{R}^{n_{in,k}}\rightarrow\mathbb{R}^{n_{st,k}}\\\
g_{k}:[t^{[\text{init}]},t^{[\text{end}]}]\times\mathbb{R}^{n_{st,k}}\times\mathbb{R}^{n_{in,k}}\rightarrow\mathbb{R}^{n_{out,k}}\\\
\end{array}$ (7)

Equations (6) - (7) are the equations representing a given subsystem. Ther are
the minimal data required to entirely characterize any subsystem, yet they do
not define the whole co-simulation configuration (connections are missing and
can be represented by extra data, see link function [Éguillon et al., 2019]).

To stick to precise concepts, we should write $x$ as a function
$x:[t^{[\text{init}]},t^{[\text{end}]}]\rightarrow\mathbb{R}^{n_{st}}$
(respectively for $y$, and $u$) as it is a time-dependent vector. That being
said, we will keep on using $x$, $y$ and $u$ notations, as if they were simple
vectors.

This representation enables us to evaluate derivatives at a precise point that
has been reached. It is compliant to methods that need only the subsystems’
equations, such as Decoupled Implicit Euler Method and Decoupled Backward
Differentiation Formulas in [Skelboe, 1992].

When $n_{in,k}$ and/or $n_{out,k}$ are $0$ for a given subsystem
$k\in[\\![1,n_{sys}]\\!]$, we will tag the topology of the subsystem $(S_{k})$
with a special name: NI, NO, NINO, or IO depending on the case (see figure 2).
This will be usefull to treat different behaviors when scheduling.

Figure 2: The four different topologies of subsystems, depending on the number
of their coupling variables

As a co-simulation implies communications at discrete times, let’s introduce
discrete notations. Let $n\in\mathbb{N}$ be the time index. Subsystem
$(S_{l})$ (with $l\in[\\![1,n_{sys}]\\!]$) communicates at times
$t_{l}^{[0]},t_{l}^{[1]},t_{l}^{[2]},...$. At these times, the $n_{out,l}$
outputs of $(S_{l})$ are known: the value of the $j$th output (with
$j\in[\\![1,n_{out,l}]\\!]$) of subsystem $(S_{l})$ at time $t_{l}^{[n]}$ will
be written $y_{l,j}^{[n]}$.

Let $(S_{k})$ be a subsystem with $n_{in,k}>0$. The $i$th input is given by
$u_{k,i}^{[n]}$ on macro-step $[t_{k}^{[n]},t_{k}^{[n+1]}[$.

Please note that $\forall n\in\mathbb{N}$:

$\forall l\in[\\![1,n_{sys}]\\!],\forall
j\in[\\![1,n_{out,l}]\\!],y_{l,j}^{[n]}\in\mathbb{R}$ whereas

$\forall k\in[\\![1,n_{sys}]\\!],\forall i\in[\\![1,n_{in,k}]\\!],$
$u_{k,i}^{[n]}:[t_{k}^{[n]},t_{k}^{[n+1]}[\rightarrow\mathbb{R}$.

###  2.2 Polynomial calibration

For $q\in\mathbb{N}^{*}$, let $\mathbb{A}(q)\subsetneq\mathbb{R}^{q}$ be the
set of elements that do not have any common values at different coordinates.
In other words:

$\begin{array}[]{l}\mathbb{A}(q)=\big{\\{}(t^{[r]})_{r\in[\\![1,q]\\!]}\in\mathbb{R}^{q}\
\big{|}\\\ \hskip 14.22636pt\forall(r_{1},r_{2})\in[\\![1,q]\\!]^{2},\
r_{1}\neq r_{2}\Rightarrow t^{[r_{1}]}\neq t^{[r_{2}]}\big{\\}}\end{array}$
(8)

For a given set of $q$ points $(t^{[r]},z^{[r]})_{r\in[\\![1,q]\\!]}$ whose
abscissas satisfy $(t^{[r]})_{r\in[\\![1,q]\\!]}\in\mathbb{A}(q)$, we define
the two following polynomials:

$\Omega_{q-1}^{Ex}:\left\\{\begin{array}[]{lcl}\mathbb{A}(q)\times\mathbb{R}^{q}\times\mathbb{R}&\rightarrow&\mathbb{R}\\\
\left(\begin{array}[]{c}(t^{[r]})_{r\in[\\![1,q]\\!]}\\\
(z^{[r]})_{r\in[\\![1,q]\\!]}\\\
t\end{array}\right)&\mapsto&\Omega_{q-1}^{Ex}\left(t\right)\\\
\end{array}\right.$ (9)

$\Omega_{q-2}^{CLS}:\left\\{\begin{array}[]{lcl}\mathbb{A}(q)\times\mathbb{R}^{q}\times\mathbb{R}&\rightarrow&\mathbb{R}\\\
\left(\begin{array}[]{c}(t^{[r]})_{r\in[\\![1,q]\\!]}\\\
(z^{[r]})_{r\in[\\![1,q]\\!]}\\\
t\end{array}\right)&\mapsto&\Omega_{q-2}^{CLS}\left(t\right)\\\
\end{array}\right.$ (10)

respectively called the extrapolation and the constrained least squares
polynomials111For the sake of readability, we will sometimes write only the
last variable of $\Omega_{q-1}^{Ex}$ and $\Omega_{q-2}^{CLS}$..

These polynomials are defined in order to have
$\Omega_{q-1}^{Ex}\in\mathbb{R}_{q-1}[t]$ and
$\Omega_{q-2}^{CLS}\in\mathbb{R}_{q-2}[t]$, where $\forall p\in\mathbb{N},\
\mathbb{R}_{p}[t]$ is the set of polynomials of the variable $t$ with
coefficients in $\mathbb{R}$ and with a degree lower or equal to $p$. We have:

$\forall r\in[\\![1,q]\\!],\Omega_{q-1}^{Ex}(t^{[r]})=z^{[r]}$ (11)

and

$\begin{array}[]{c}\Omega_{q-2}^{CLS}:t\mapsto\sum_{i=0}^{q-2}a_{i}t^{i}\
\text{where}\\\ (a_{i})_{\begin{subarray}{c}\\\
i\in[\\![0,q-2]\\!]\end{subarray}}=\\!\\!\\!\\!\\!\\!\\!\underset{\begin{subarray}{c}(\bar{a}_{i})_{i\in[\\![0,q-2]\\!]}\\\
\Omega_{q-2}^{CLS}(t^{[1]})=z^{[1]}\end{subarray}}{\arg\min\vspace{-0.1cm}}\\!\\!\\!\\!\\!\\!\Big{\\{}\underset{r=1}{\overset{q}{\sum}}\\!\\!\Big{(}\\!\\!z^{[r]}\\!-\\!\\!\underset{i=0}{\overset{q-2}{\sum}}\\!\bar{a}_{i}(t^{[r]})^{i}\\!\Big{)}^{2}\Big{\\}}\end{array}$
(12)

In practice, the $z$ variables will be either inputs $u$ or outputs $y$, and
index will correpond to a time index where $z^{[1]}$ is the latest one, so
that the constraint on the least squares (12) is the equality on the most
recent point.

When we will consider one of these polynomials of degree $p$ (either
extrapolation on $p+1$ points or constrained least squares on $p+2$ points)
without specifying which one is used, we will use the generic notation
$\Omega_{p}$.

The coefficients of $\Omega_{p}^{Ex}$ can be computed by several methods
(Lagrange polynomials, Newton’s formula, barycentric approach [Berrut and
Trefethen, 2004]), and the coefficients of $\Omega_{p}^{CLS}$ can be obtained
with a cosntrained linear model (formula $(1.4.11)$ page $22$ of [Amemiya,
1985]).

Finally, we define the Hermite interpolation polynomial $\mathcal{H}$ in the
specific case with two points and first order derivatives:

$\mathcal{H}:\left\\{\begin{array}[]{lcl}\mathbb{R}^{2}\times\mathbb{R}^{2}\times\mathbb{R}^{2}\times\mathbb{R}&\rightarrow&\mathbb{R}\\\
(\left(\begin{subarray}{c}t^{[1]}\\\
t^{[2]}\end{subarray}\right),\left(\begin{subarray}{c}z^{[1]}\\\
z^{[2]}\end{subarray}\right),\left(\begin{subarray}{c}\dot{z}^{[1]}\\\
\dot{z}^{[2]}\end{subarray}\right),t)&\mapsto&\mathcal{H}(t)\end{array}\right.$
(13)

where

$\begin{array}[]{l}\mathcal{H}\in\mathbb{R}_{3}[t]\ \text{and}\ \forall
r\in[\\![1,2]\\!],\\\ \mathcal{H}(t^{[r]})=z^{[r]}\ \text{and}\
\textstyle{\frac{d\mathcal{H}}{dt}}(t^{[r]})=\dot{z}^{[r]}\\\ \end{array}$
(14)

The coefficients of $\mathcal{H}$ can be computed with the square of the
Lagrange polynomial basis or by using divided differences [Hildebrand, 1956].

Hermite interpolation will be used for smoothness enhancement in 3.1.4.

##  3 F3ORNITS ALGORITHM

We introduce here the F3ORNITS method, standing for Flexible Order
Representation of New Inputs, including flexible Time-stepper (with variable
step size, when applicable) and flexible Scheduler (asynchronous-capable, when
applicable).

The method stems from the desire to keep the dynamical behavior of the
coupling variables, what zero order hold (ZOH) does not do. At a given
communication time, the outputs of the past data exchange will be reused in
order to fit polynomial estimations for the future (the upcoming step). This
is done in several co-simulation methods [Kübler and Schiehlen, 2000] [Busch,
2016], yet usually the polynomial order is decided in advance. F3ORNITS has a
flexible order and will decide of an order for every coupling variable at each
communication time. Moreover, we will focus on the way to use these
estimations in an asynchronous context (as it is not always possible to pass
information at the same simulation time from one system to another). The
error, depending on the polynomial order, will also be used to decide the
evolution of the macro-step size (see figure 1). Regarding subsystems with
limited capabilities for holding time-dependent inputs, strategies are
proposed to take advantage of the variable order thanks to an adaptation of
the data to the capabilities of such subsystems. The latter strategies fit the
specification to handle any modular model regardless the missing capabilities
in every subsystem.

Finally, the time-stepping strategy (including the time-stepper itself and the
scheduler) will be presented. The time-stepper is based on the error made by
the estimation described in 3.1. The normalization of this error has a strong
impact on the time-stepping criterion, so several normalization methods will
be described and a new one will be introduced: the damped amplitude
normalization method. The scheduler is also a part of the F3ORNITS method, it
occurs once the time-stepper produced an estimation of the upcoming step
sizes, yet it shall not be fully detailed for the sake of space.

The smoothness enhancement [Busch, 2019] [Dronka and Rauh, 2006] will also be
presented as it is compliant with the F3ORNITS method. Nevertheless, we
adapted it to the context of a flexible order and variable step size method.
The motivation is similar to the one in [Éguillon et al., 2019]. In the case
where subsystems do not have sufficient capabilities (up to $3$rd order
polynomial inputs), the F3ORNITS method can still run without smoothness
enhancement.

###  3.1 Flexible polynomial inputs

We arbitrary set the maximum degree for polynomial inputs to:

$M=2$ (15)

Let $(m_{k})_{k\in[\\![1,n_{sys}]\\!]}$ be the maximum degrees for polynomial
inputs supported for each subsystem. As we want to support every kind of
subsystem, we cannot assume anything on $m_{k}$ (we only know that $\forall
k\in[\\![1,n_{sys}]\\!],m_{k}\geqslant 0$).

Let’s define the effective maximum degree for each subsystem, by adding the
constraint (15):

$\forall k\in[\\![1,n_{sys}]\\!],\ M_{k}:=\min(M,m_{k})$ (16)

The time-dependent inputs that will be generated will always satisfy (17), and
for each degree in these ranges, the maximum supported degree for each
subsystem will never be exceeded thanks to (16).

$\forall k\in[\\![1,n_{sys}]\\!],\ \forall n\in\mathbb{N},\
u_{k}^{[n]}\in\left(\mathbb{R}_{M_{k}}[t]\right)^{n_{in,k}}$ (17)

The determination of the order to use for polynomial inputs is quite
straightforward. However, as this order may be different for each variable and
as it is determined in the subsystem holding the corresponding variable as an
output, some complications due to the asynchronousness may appear on
subsystems having an input connected to this variable.

In order to clarify the process and to deal with properly defined mathematical
concepts, we will split this explanation into three parts:

To begin with, we will define the order function for every output variable and
see how the values of this function are found. Then, we will define the
estimated output variables based on order function. Finally, we will see how
these estimated output variables are used from the connected input’s
perspective.

####  3.1.1 Order function

Let’s consider a subsystem $(S_{l})$ with $l\in[\\![1,n_{sys}]\\!]$. Let’s
consider we have already done $n$ macro-steps with $n\in\mathbb{N}^{*}$, id
est at least one.

For $j\in[\\![1,n_{out,l}]\\!]$ and $q\in[\\![0,\min\left(M,n-1\right)]\\!]$,
we define ${}^{(q)}err_{l,j}^{[n+1]}$ the following way:

$^{(q)}err_{l,j}^{[n+1]}=\left|y_{l,j}^{[n+1]}-\Omega_{q}^{Ex}\left(t_{l}^{[n+1]}\right)\right|$
(18)

where $\Omega_{q}^{Ex}$ is calibrated on $(t_{l}^{[n-r]})_{r\in[\\![0,q]\\!]}$
and $(y_{l,j}^{[n-r]})_{r\in[\\![0,q]\\!]}$ as explained in (9) and (11). In
other words the value obtained at $t_{l}^{[n+1]}$ is compared to the
extrapolation estimation based on the $(q+1)$ exchanged values before
$t_{l}^{[n+1]}$ (excluded), that is to say at
$t_{l}^{[n-q]},t_{l}^{[n-q+1]},t_{l}^{[n-q+2]},...,t_{l}^{[n]}$ as shown on
figure 3.

We now define the order function $p_{l}$ for subsystem $(S_{l})$, which is a
vectorial function for which each coordinate corresponds to one output of the
subsystem $(S_{l})$. It is a step function defining the ”best order” of
extrapolation for each macro-step, based on the previous macro-step. This
includes a delay that makes it possible to integrate the subsystems
simultaneously.

$p_{l}:\left\\{\begin{array}[]{lcl}[t^{[\text{init}]},t^{[\text{end}]}[&\rightarrow&[\\![0,M]\\!]^{n_{out,l}}\\\
t&\mapsto&\left(p_{l,j}(t)\right)_{j\in[\\![1,n_{out,l}]\\!]}\\\
\end{array}\right.$ (19)

For $j\in[\\![1,n_{out,l}]\\!]$, the $p_{l,j}$ functions are defined this way:

$\begin{array}[]{l}p_{l,j}:\left\\{\begin{array}[]{lcl}[t^{[\text{init}]},t^{[\text{end}]}[&\rightarrow&[\\![0,M]\\!]\\\
t&\mapsto&p_{l,j}(t)\ \text{with:}\\\ \end{array}\right.\\\
p_{l,j}(t)=\displaystyle{\sum_{n=1}^{n_{\max}}}\left(\mathds{1}_{[t_{l}^{[n]},t_{l}^{[n+1]}[}(t)\\!\\!\\!\\!\\!\\!\\!\underset{q\in[\\![0,\min(M,n-1)]\\!]}{\arg\min}\\!\\!\left\\{{}^{q}err_{l,j}^{[n]}\right\\}\right)\end{array}$
(20)

where $n_{\max}$ satisfies $t_{l}^{[n_{\max}]}=t^{[\text{end}]}$.

An illustration of the ”$\arg\min$” choice in (20) is presented on figure 3,
and the order function itself can be visualized on its plot in figure 4.

Figure 3: Order determination

####  3.1.2 Estimated outputs

We still consider the subsystem $(S_{l})$ and its $n_{out,l}$ outputs. Let’s
admit the step $[t_{l}^{[n]},t_{l}^{[n+1]}[$ is computed. Therefore, for all
$j\in[\\![1,n_{out,l}]\\!]$, the value of $p_{l,j}(t)$ for $t$ in
$[t_{l}^{[n+1]},t_{l}^{[n+2]}[$ is known.

We can now determine the estimated outputs
$\hat{y}_{l}^{[n+1]}\in\left(\mathbb{R}_{q}[t]\right)^{n_{out,l}},\
t\in[t_{l}^{[n+1]},t_{l}^{[n+2]}[$. We have two choices for the way this
estimation is made: Extrapolation mode and Constrained Least Squares (CLS)
mode.

For the sake of genericity, we will define the abstraction of this choice
$\Omega_{q}$ introduced in 2.2.

Estimated outputs are defined the following way on the step
$[t_{l}^{[n+1]},t_{l}^{[n+2]}[$:

$\hat{y}_{l}^{[n+1]}:\left\\{\begin{array}[]{lcl}[t_{l}^{[n+1]},t_{l}^{[n+2]}[&\rightarrow&\mathbb{R}^{n_{out,k}}\\\
t&\mapsto&\left(\hat{y}_{l,j}^{[n+1]}(t)\right)_{j\in[\\![1,n_{out,l}]\\!]}\\\
\end{array}\right.$ (21)

where

$\hat{y}_{l,j}^{[n+1]}:\left\\{\begin{array}[]{lcl}[t_{l}^{[n+1]},t_{l}^{[n+2]}[&\rightarrow&\mathbb{R}\\\
t&\mapsto&{\Omega_{q}}_{\big{|}_{\big{[}t_{l}^{[n+1]},t_{l}^{[n+2]}\big{[}}}\left(t\right)\\\
\end{array}\right.$ (22)

where $\Omega_{q}$ is calibrated on
$(t_{l}^{[n+1-r]},y_{l,j}^{[n+1-r]})_{r\in[\\![0,p]\\!]}$ in extrapolation
mode (see (11)), on
$(t_{l}^{[n+1-r]},y_{l,j}^{[n+1-r]})_{r\in[\\![0,p+1]\\!]}$ in CLS mode (see
(12)), and where $q:=p_{l,j}(t_{l}^{[n+1]})$ so that
$\deg(\hat{y}_{l,j}^{[n+1]})=p_{l,j}(t_{l}^{[n+1]})$ (order determined
previously is respected).

Let $\breve{y}_{l,j}^{[n+1]}$ be the extension of the polynomial
$\hat{y}_{l,j}^{[n+1]}$ on the whole $\mathbb{R}$ domain.

The degree of $\breve{y}_{l,j}^{[n+1]}$ is given by $p_{l,j}(t_{l}^{[n+1]})$
as shown on figure 4.

Figure 4: Order function and corresponding output estimation with
corresponding orders

As in (22), for the rest of 3.1.2, we will use the following notation (for the
sake of readability):

$q=p_{l,j}(t_{l}^{[n+1]})$ (23)

The definition of $\Omega_{q}$ varies depending on the mode. In
”Extrapolation” mode, the oldest point is forgotten (id est not taken into
account in the extrapolation): $(t_{l}^{[n-q]},\ y_{l,j}^{[n-q]})$. This
point, taken into account to choose the order $q$ but not used in
extrapolation calibration, is represented as striped in figure 5.

On the other hand, the idea of ”CLS” mode is to take into account the point
described above in the estimation of the output. An idea can be to forget the
most recent point, that is to say: $(t_{l}^{[n+1]},y_{l,j}^{[n+1]})$, but this
would mean that the value given by the subsystem’s integrator would be unused,
so we will introduce a delay in the coupling process. Thus, the strategy is to
take into account all the $(q\\!+\\!1)$ points that have been used in the
determination of the chosen order, and the most recent point as well. We thus
have $(q\\!+\\!2)$ points to adjust a polynomial of degree at most $q$: an
extrapolation process cannot be made, but the ”best fitting” polynomial can be
found for the CLS criterion (12). Please note that removing the constraint
$\Omega_{q-2}^{CLS}(t^{[1]})=z^{[1]}$ in (12) corresponds to the relaxation
technique on the past refered to as ”method 1” in [Li et al., 2020] in the
particular case of $q=0$.

Figure 5: Extrapolation mode vs. CLS mode, on step
$[t_{l}^{[n+1]},t_{l}^{[n+2]}[$ once step $[t_{l}^{[n]},t_{l}^{[n+1]}[$ is
done (so that $p_{l,j}(t^{[n+1]})$ could be computed)

####  3.1.3 Estimated inputs

From the inputs perspective, the order function of the connected output is
used. Let’s consider a subsystem $(S_{k})$ with $k\in[\\![1,n_{sys}]\\!]$,
which have $n_{in,k}>0$ inputs. As the case $n_{in,k}=0$ should not be
excluded, let’s admit that nothing is done from the inputs perspective in this
case (because there is no input). From here and for the whole 3.1.3 section,
we will consider $n_{in,k}\in\mathbb{N}^{*}$.

We will consider the input $i$ with $i\in[\\![1,n_{in,k}]\\!]$ and, to
properly consider the connected output, we stand $l\in[\\![1,n_{sys}]\\!]$ and
$j\in[\\![1,n_{out,l}]\\!]$ so that input $i$ of subsystem $k$ is fed by
output $j$ of subsystem $l$ (with the link function notation of [Éguillon et
al., 2019], we would write $(l,j)=L(k,i)$).

As asynchronousness should be supported, a special care should be made when
the step $[t_{k}^{[n]},t_{k}^{[n+1]}[$ does not fit into one single definition
of $\hat{y}_{l,j}^{[m]}$ for an $m\in\mathbb{N}$. In this case, we will use:

$u_{k,i}^{[n]}:\left\\{\begin{array}[]{lcl}[t_{k}^{[n]},t_{k}^{[n+1]}[&\rightarrow&\mathbb{R}\\\
t&\mapsto&\breve{y}_{l,j}^{[m]}(t)\end{array}\right.$ (24)

where $\breve{y}$ is defined in 3.1.2 and where

$m=\max\left\\{m\in\mathbb{N}\ \big{|}\ t_{l}^{[m]}\leqslant
t_{k}^{[n]}\right\\}$ (25)

(24) (25) can be visualized figure 6.

Figure 6: Estimated input using corresponding extended estimated output

Some subsystems cannot hold polynomial inputs, and some other can but not
necessary at any order.

As F3ORNITS method only requires to hold up to order $2$ (except when
smoothness enhancement is triggered (see 3.1.4), which is a particular case),
we will consider $9$ cases represented in the table below.

Figure 7: Alternatives to decrease polynomial input degree for subsystems with
limited capabilities

####  3.1.4 Smoothness enhancement

Smoothness enhancement can be triggered to enable $C^{1}$ inputs. This mode
will only be applicable on subsystems supporting at least $3$rd order
polynomial inputs.

Let’s consider a system $(S_{k})$ with $k\in[\\![1,n_{sys}]\\!]$. Let’s
consider the step $[t_{k}^{[n]},t_{k}^{[n+1]}[$ with $n\in\mathbb{N}^{*}$. We
will also consider only the input $i$ of this system, with
$i\in[\\![1,n_{in,k}]\\!]$. As the smoothness enhancement process should be
applied on every input separately, it will be detailed only on $u_{k,i}$ here.

The idea is to guarantee that $C^{1}$ smoothness is not broken at time
$t^{[n]}$. In other word, we will remove the jump at the communication time.
Moreover, as several consecutive steps are concerned, the $C^{1}$ smoothness
won’t be broken on the whole time interval (union of the steps).

Regardless of the degree of the polynomial input, we will extend it to a third
order polynomial using Hermite interpolation, as shown on figure 8.

Figure 8: Redefinition of the time-dependent input on a co-simulation step in
the case of smoothness enhancement

As $(S_{k})$ reached the time $t_{k}^{[n]}$, we know the expression of the
input $u_{k,i}^{[n-1]}$ that has been used for $[t^{[n-1]},t^{[n]}[$, so we
can compute the left constraints.

$\begin{array}[]{c}u^{[n],left}:=\underset{\begin{subarray}{c}t\rightarrow
t^{[n]}\\\ t<t_{[n]}\end{subarray}}{\lim\ \ }u^{[n-1]}(t)\\\
\dot{u}^{[n],left}:=\underset{\begin{subarray}{c}t\rightarrow t^{[n]}\\\
t<t_{[n]}\end{subarray}}{\lim\ \ }\dot{u}^{[n-1]}(t)\end{array}$ (26)

Moreover, before performing the simulation on the step $[t^{[n]},t^{[n+1]}[$,
we compute the polynomial input (of degree $0$, $1$, or $2$, computed by
extrapolation or by CLS) $u_{k,i}^{[n]}$. We thus compute the right
constraints.

$\begin{array}[]{c}u^{[n],right}:=u^{[n+1],left}:=\underset{\begin{subarray}{c}t\rightarrow
t^{[n+1]}\\\ t<t_{[n+1]}\end{subarray}}{\lim\ \ }u^{[n]}(t)\\\
\dot{u}^{[n],right}:=\dot{u}^{[n+1],left}:=\underset{\begin{subarray}{c}t\rightarrow
t^{[n+1]}\\\ t<t_{[n+1]}\end{subarray}}{\lim\ \ }\dot{u}^{[n]}(t)\end{array}$
(27)

Finally, instead of using $u_{k,i}^{[n]}$ on $[t^{[n]},t^{[n+1]}[$, we will
use the ”smooth” version of it:

$^{smooth}u_{k,i}^{[n]}:t\mapsto\mathcal{H}\Big{(}\left(\begin{subarray}{c}t^{[n]}\\\
t^{[n+1]}\end{subarray}\right),\left(\begin{subarray}{c}u^{[n],left}\\\
u^{[n],right}\end{subarray}\right),\left(\begin{subarray}{c}\dot{u}^{[n],left}\\\
\dot{u}^{[n],right}\end{subarray}\right),t\Big{)}$ (28)

where $\mathcal{H}$ denotes the Hermite interpolation polynomial described in
section 2.2.

Using Hermite interpolation on values known by extrapolation is sometimes
refered to as ”extrapolated interpolation” [Busch, 2019] [Dronka and Rauh,
2006].

###  3.2 Flexible time management

The other major aspect of F3ORNITS algorithm (the first one being the flexible
polynomial inputs) is the time management, which includes a time stepper and a
scheduler.

The time-stepper defines the next communication time after a macro-step is
finished. The scheduler ensures the coherence of the rendez-vous times of all
subsystems based on their connections, topologies and constraints.

####  3.2.1 Time-stepper

Let’s consider a system $(S_{l}),l\in[\\![1,n_{sys}]\\!]$ that is either IO or
NI (see figure 2). We have: $n_{out,l}>0$. Let’s consider an output
$y_{l,j},j\in[\\![1,n_{out,l}]\\!]$ of $(S_{l})$.

The aim of the time-stepper is to determine $t_{l}^{[n+2]}$ once step
$[t_{l}^{[n]},t_{l}^{[n+1]}[$ have been computed.

For the sake of readability, we will use the notation
$p:=p_{l,j}(t_{l}^{[n]})$ (different from $q$ in (23) as here we focus on the
lastly computed step, and not the upcoming one). Let’s introduce the macro-
step size and the dilatation coefficient, respectively (29) and (30).

$\delta t_{l}^{[n]}=t_{l}^{[n+1]}-t_{l}^{[n]}$ (29)

$\rho_{l}^{[n+1]}=\nicefrac{{\delta t_{l}^{[n+1]}}}{{\delta t_{l}^{[n]}}}$
(30)

They will be detailed later in this subsection, but for now we only need to
know that the dilatation coefficients are bounded:
$\exists(\rho_{\min},\rho_{\max})\in\mathbb{R}_{+}^{*},\forall
l\in[\\![1,n_{sys}]\\!],\forall
n\in\mathbb{N},\rho_{l}^{[n]}\in[\rho_{\min},\rho_{\max}]$.

In extrapolation mode, it exists
$\zeta_{t}\in\Big{[}t_{l}^{[n-p]},t_{l}^{[n+1]}\Big{]}$ so that

$\begin{array}[]{l}y_{l,j}^{[n+1]}-\breve{y}_{l,j}^{[n]}(t_{l}^{[n+1]})\\\
=\underbrace{\displaystyle{\frac{1}{(p+1)!}}\displaystyle{\frac{d^{(p+1)}y_{l,j}}{dt^{(p+1)}}}(\zeta_{t})}_{\text{independent
of }\delta
t_{l}^{[n]}}\displaystyle{\prod}_{r=0}^{p}(t_{l}^{[n+1]}-t_{l}^{[n-r]})\\\
=c_{1}\cdot\displaystyle{\prod}_{r=0}^{p}\Big{(}\sum_{s=0}^{r}\delta
t_{l}^{[n-s]}\Big{)}\\\ =c_{1}\cdot\displaystyle{\prod}_{r=0}^{p}\Big{(}\delta
t_{l}^{[n]}\Big{(}1+\sum_{s=1}^{r}\Big{(}\prod_{c=1}^{s}\textstyle{\frac{1}{\underbrace{\rho_{l}^{[n+1-c]}}_{\in[\rho_{\min},\rho_{\max}]}}}\Big{)}\Big{)}\Big{)}\\\
\end{array}$

$\begin{array}[]{l}\leqslant
c_{1}\cdot\displaystyle{\prod}_{r=0}^{p}\Big{(}\delta
t_{l}^{[n]}\Big{(}1+\sum_{s=1}^{r}\Big{(}\prod_{c=1}^{s}\textstyle{\frac{1}{\rho_{\min}}}\Big{)}\Big{)}\Big{)}\\\
\leqslant c_{1}\cdot\Big{(}\delta
t_{l}^{[n]}\Big{)}^{p+1}\underbrace{\displaystyle{\prod}_{r=0}^{p}\Big{(}\sum_{s=0}^{r}\Big{(}\textstyle{\frac{1}{\rho_{\min}}}\Big{)}^{s}\Big{)}}_{\text{independent
of }\delta t_{l}^{[n]}}\\\ \leqslant c_{1}\cdot\Big{(}\delta
t_{l}^{[n]}\Big{)}^{p+1}\cdot c_{2}\\\ \end{array}$ (31)

where $c_{1}$ and $c_{2}$ are random constants.

Final expression of error in (31) shows that the error is of order $(p+1)$ on
the mean macro-step sizes. Analogously, we will consider that the error is of
the same order error in CLS mode: $(p+1)$ for a polynomial of degree $p$
(generated with a constrained least-square fitting on $(p+2)$ points).

The time-stepper uses this known error and error order to adapt the step size
accordingly. We use a formula similar to the one in [Schierz et al., 2012]
(also mentioned in [Gomes et al., 2018]) to define a dilatation coefficient
candidate per output $j\in[\\![1,n_{out,l}]\\!]$.

$\rho_{l,j}^{[n+1]}:=\sqrt[p+1]{\frac{1}{error_{l,j}^{[n+1]}}}$ (32)

In (32), the $error$ term is expected to be a relative error either relative
to the values of the concerned variable and to a relative tolerance given
which will determine the error threshold over which the step size is expected
to decrease.

Given such a tolerance $tol_{rel}$ and an absolute tolerance $tol_{abs}$, a
first approach refered to as Magnitude relies on the order of magnitude of the
variable at the moment of the communication. In that case, the error is
defined as follow.

$^{Magn.}error_{l,j}^{[n+1]}:=\frac{\left|y_{l,j}^{[n+1]}-\breve{y}_{l,j}^{[n]}(t_{l}^{[n+1]})\right|}{tol_{abs}+tol_{rel}\cdot
y_{l,j}^{[n+1]}}$ (33)

The problem in this approach is that when values of $y_{l,j}$ are close to
zero, the trend is to give a big error (and the step size will then be
reduced). It is particularly problematic with variables with a great order of
magnitude and periodically crossing the value zero, such as sinusoids.

An approach that might reduce this effect is to normalize the error according
to the amplitude (observed since the beginning of the simulation) instead of
the order of magnitude. This approach will be refered to as Amplitude and
defines the error as follow.

$\begin{array}[]{l}{}^{Ampl.}error_{l,j}^{[n+1]}:=\\\ \\\
\lx@intercol\hfil\displaystyle{\frac{\left|y_{l,j}^{[n+1]}-\breve{y}_{l,j}^{[n]}(t_{l}^{[n+1]})\right|}{tol_{abs}+tol_{rel}\\!\left(\\!\underset{m\in[\\![0,n+1]\\!]}{\max}\\!\\!\left(y_{l,j}^{[m]}\right)-\\!\\!\\!\\!\underset{m\in[\\![0,n+1]\\!]}{\min}\\!\\!\left(y_{l,j}^{[m]}\right)\\!\\!\right)}}\lx@intercol\end{array}$
(34)

The problem in this approach is that when the values of $y_{l,j}$ undergoes a
great jump at one single moment (e.g. at initialization), the trend will be to
produce artificially small errors. Therefore, big step sizes will be produced
by the time-stepper, and the accuracy may dwindle.

Solving this problem can be done by damping the amplitude, in order to
progressively erase the effects of jumps while keeping the local amplitude of
the variable $y_{l,j}$. The error produced with this principle will refer to
the Damped Amplitude strategy, and will be defined as follow.

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!\\!{}^{\begin{subarray}{c}Damped\\\
Ampl.\end{subarray}}error_{l,j}^{[n+1]}:=\\\ \\\ \lx@intercol\hfil\hskip
28.45274pt\displaystyle{\frac{\left|y_{l,j}^{[n+1]}-\breve{y}_{l,j}^{[n]}(t_{l}^{[n+1]})\right|}{tol_{abs}+tol_{rel}\\!\left(\\!{}^{\begin{subarray}{c}damp\\\
max\end{subarray}}\\!\\!y_{l,j}^{[n+1]}-^{\begin{subarray}{c}damp\\\
min\end{subarray}}\\!\\!y_{l,j}^{[n+1]}\\!\\!\right)}}\lx@intercol\end{array}$
(35)

Expression (35) refers to the damped minimal and damped maximal sequences
which are recursively defined as follow.

$\left\\{\begin{array}[]{>{\hspace{-0.2cm}}l<{\hspace{-0.3cm}}c<{\hspace{-0.3cm}}l}\alpha^{[0]&=&0\\\
^{\begin{subarray}{c}damp\\\
max\end{subarray}}\\!\\!y_{l,j}^{[0]&=&y_{l,j}^{[0]}\\\
^{\begin{subarray}{c}damp\\\
min\end{subarray}}\\!\\!y_{l,j}^{[0]&=&y_{l,j}^{[0]}\\\
\alpha^{[m]&=&{}^{\begin{subarray}{c}damp\\\
max\end{subarray}}\\!\\!y_{l,j}^{[m]}-^{\begin{subarray}{c}damp\\\
min\end{subarray}}\\!\\!y_{l,j}^{[m]}\\\ ^{\begin{subarray}{c}damp\\\
max\end{subarray}}\\!\\!y_{l,j}^{[m]&=&\max\\!\Big{\\{}y_{l,j}^{[m]},^{\begin{subarray}{c}damp\\\
max\end{subarray}}\\!\\!y_{l,j}^{m-1}\\!-\textstyle{\frac{\nu\cdot\delta
t_{l}^{[m-1]}}{2}}\cdot\alpha^{[m-1]}\Big{\\}}\\\ ^{\begin{subarray}{c}damp\\\
min\end{subarray}}\\!\\!y_{l,j}^{[m]&=&\min\\!\Big{\\{}y_{l,j}^{[m]},^{\begin{subarray}{c}damp\\\
min\end{subarray}}\\!\\!y_{l,j}^{m-1}\\!+\textstyle{\frac{\nu\cdot\delta
t_{l}^{[m-1]}}{2}}\cdot\alpha^{[m-1]}\Big{\\}}\\\ \end{array}\right.}}}}}}$
(36)

where $\forall n\in\mathbb{N},\ \delta t_{l}^{[n]}=t_{l}^{[n+1]}-t_{l}^{[n]}$
denotes the step size.

The damping coefficient $\nu\geqslant 0$ has to be defined. The greater it
will be, the faster an event such as a jump will be ”forgotten”; the smaller
it will be, the closer ${}^{\begin{subarray}{c}Damped\\\
Ampl.\end{subarray}}error_{l,j}^{[n+1]}$ will be to
${}^{Ampl.}error_{l,j}^{[n+1]}$.

We can now provide a proper definition of the subsystem’s dilatation ratio
mentioned in (30) as the safest $\rho$ among the candidates (32).

$\rho_{l}^{[n+1]}:=\underset{j\in[\\![1,n_{out,l}]\\!]}{\min}\left(\rho_{l,j}^{[n+1]}\right)$
(37)

With this dilatation ratio, we obtain an estimation of the next communication
time for subsystem $(S_{l})$:

$\begin{array}[]{lcl}\delta t_{l}^{[n]}&=&t_{l}^{[n+1]}-t_{l}^{[n]}\\\
\widetilde{\delta t}_{l}^{[n+1]}&=&\rho_{l}^{[n+1]}\cdot\delta t_{l}^{[n]}\\\
\tilde{t}_{l}^{[n+2]}&=&t_{l}^{[n+1]}+\widetilde{\delta
t}_{l}^{[n+1]}\end{array}$ (38)

This next estimation temporary solves the problematic introduced at the
beginning of 3.2.1, the determination of $t_{l}^{[n+2]}$ once
$[t_{l}^{[n]},t_{l}^{[n+1]}[$ is computed. Nonetheless, this is only an
estimation as the scheduler presented in 3.2.2 may modify it and determine the
$t_{l}^{[n+2]}$ to use.

Remark 1 (ratio bounds): In practice, a minimal and a maximal value for the
dilatation ratio (37) can be define in order to avoid unsafe extreme step size
reduction/increase. In section 4, the dilatation coefficient is projected in
the interval $[10\%,105\%]$.

Remark 2 (subsystems without outputs): The above procedure only works for
subsystems having at least one output variable. For the other subsystems (NO
and NINO), we will use the following rule.

$\begin{array}[]{l}\forall l\in\left\\{l\in[\\![1,n_{sys}]\\!]\ |\
n_{out,l}=0\right\\},\\\ \hskip 56.9055pt\forall
n\in\mathbb{N}\cap[2,+\infty[,\
\tilde{t}_{l}^{[n]}=t^{[\text{end}]}\end{array}$ (39)

Remark 3 (co-simulation start): The rules presented above do not enable to set
the two first times. We will obviously set the first one with $\forall
l\in[\\![1,n_{sys}]\\!],\tilde{t}_{l}^{[0]}=t^{[\text{init}]}$ and the next
one by using an initial step $\delta t_{l}^{[0]}$ given as co-simulation
parameter for every subsystem. This initial step size will have a limited
influence on the whole co-simulation. We then set
$\tilde{t}_{l}^{[1]}=t^{[\text{init}]}+\delta t_{l}^{[0]}$ for each subsystem.

####  3.2.2 Scheduler

The scheduler is intended to adjust the communication times of the subsystems
according to the times they reached individually. It is based on topological
data (an adjacency matrix representing the connections between the subsystems
and the set of constraints regarding step sizes on each subsystem) and can
handle asynchronousness in a way that makes it robust to co-simulation
involving subsystems with imposed step sizes. Several stages happen
sequentially once every subsystem has produced an estimated next communication
time (see 3.2.1). Among them, some are dedicated to avoid to use a variable
further than the time it is supposed to be used according to the time-stepper
(trying to avoid phenomena like the one happening on figure 6), and others act
as optimisations, as they avoid to communicate data when this has no effect
(NI systems feeding a NO subsystem with a large imposed step size, for
instance).

The stages of the scheduler will not be detailed here for the sake of space.
It can be seen as a black box producing the effective next communication steps
once the time-stepper produced the estimated ones.

In other words, once $\forall l\in[\\![1,n_{sys}]\\!],\
[t_{l}^{[n_{l}]},t_{l}^{[n_{l}+1]}[$ is computed ($n_{l}$ might be different
per subsystems as some subsystems may be idle when other keep on going, so
they might do a different amount of steps), the time-stepper has produces
$(\tilde{l}^{[n_{l}+1]})_{l\in[\\![1,n_{sys}]\\!]}$, and the scheduler can be
seen as the black box $Sch$ below.

$(t_{l}^{[n_{l}+2]})_{l\in[\\![1,n_{sys}]\\!]}:=Sch\left(\begin{array}[]{c}(t_{l}^{[n_{l}+1]})_{l\in[\\![1,n_{sys}]\\!]}\\\
(\tilde{t}_{l}^{[n_{l}+2]})_{l\in[\\![1,n_{sys}]\\!]}\end{array}\right)$ (40)

##  4 RESULTS AND BEHAVIOR ON TWO TEST CASES

Two test cases will be presented in this section. Both of them have been made
with Simcenter Amesim.

The first model will present the importance that keeping the dynamics on the
coupling variables may have. This is actually done with the polynomial inputs
in F3ORNITS method. The second model is a variation of the 2-masses-springs-
dampers [Éguillon et al., 2019] [Busch, 2016] and performances will be
measured on several co-simulations.

The most widely used co-simulation method is the non-iterative Jacobi
(parallel) fixed-step algorithm, with zero-order-hold inputs. Therefore, the
comparisons will be made between the latter (refered to as ”NI Jacobi”) and
F3ORNITS algorithm.

###  4.1 Car with controlled speed

The model consists in two subsystems, one corresponds to a $1000$kg car
(simplified as a mass) moving on a $1$D axis (straight road) modelled with
Newton’s second law, which takes a force on entry and gives its position as
output, and the other corresponds to a controller, producing a force and using
the car position on input.

Figure 9: Test case 1: car with controlled speed - Subsystems and monolithic
reference in Simcenter Amesim

As it can be seen on figure 9, the car subsystem adds a random perturbation to
the input force (this may be seen as a 1D wind, for instance). Moreover, on
the first $10$s of the (co-)simulation, the controller will output a force
that is predetermined and that can be seen on figure 10.

Figure 10: Preset output force from the controller on $[0,10[$

At $t=10$s, the output of the controller subsystem becomes the output of a ”P”
controller based on the vehicle speed and which is designed to make the
vehicle reach (or maintain) a target speed of $16$m s$-1$. The velocity needed
by the ”P” controller is computed from the input position, using an explicit
equation $\dot{x}=v$ and a constraint equation $x-u=0$ where $v$ denotes the
computed speed, $x$ a hidden variable representing the position, and $u$ the
position input in the controller subsystem.

Due to the zero order hold, the NI Jacobi method does not allow $(S_{2})$ to
properly retrieve the vehicle speed. Indeed, on every co-simulation step, the
input position will be constant and the DAE will produce a null speed (see
figure 11, as it can be seen on the zoom on $[13,13.2]$.

Figure 11: Vehicle speed computed in controller $(S_{2})$ using its input
(vehicle position) with a ZOH input

The consequence of this is that the controller will put force to push the
vehicle in order to increase the speed (attempting to make it reach $16$m
s$-1$) without even realizing that the vehicle has a non-null velocity.
Therefore, the velocity will keep on increasing indefinitely (curve with
triangles on figure 12). Inversely, F3ORNITS represents the dynamics of the
coupling variables (Force and vehicle position) when needed (this is the
”flexible order” part presented in 3.1). The velocity computed from such a
position input is more relyable, and the controller $(S_{2})$ will send proper
force to make the vehicle maintain its target velocity.

The effective vehicle velocity inside of $(S_{1})$ across the time is
presented in figure 12 in the three cases of (co-)simulation.

Figure 12: Vehicle speed in $(S_{1})$ in the monolithic reference and in co-
simulations with respectively the NI Jacobian (ZOH) method and the F3ORNITS
algorithm

###  4.2 Two masses, springs and dampers

The model consists of two masses coupled with force, displacement and velocity
as shown on figure 13. It is a variant of the test cases presented in
[Éguillon et al., 2019] or in [Busch, 2016] in the sense that, after a
transition time (set to the middle of the simulation time interval), the
behavior of one of the models changes. This generates a discontinuity on
coupling variables and different behaviors before and after.

$\begin{array}[]{|lrr|lrr|lrr|}\hline\cr c_{1}&=&1\ \text{kN/m}&\ d_{1}&=&1\
\text{kN/(m/s)}&\ x_{1}(t^{[\text{init}]})&=&-1\ \text{m}\\\ c_{2}&=&1\
\text{kN/m}&\ d_{2}&=&0\ \text{kN/(m/s)}&\ x_{2}(t^{[\text{init}]})&=&0\
\text{m}\\\ c_{3}&=&1\ \text{kN/m}&\ d_{3}&=&1\
\text{kN/(m/s)}&v_{1}(t^{[\text{init}]})&=&0\ \text{m/s}\\\ m_{1}&=&1000\
\text{kg}&\ m_{2}&=&1000\ \text{kg}&v_{2}(t^{[\text{init}]})&=&0\
\text{m/s}\\\ \vrule\lx@intercol\hfil[t^{[\text{init}]},t^{[\text{end}]}]=[0\
\text{s},\ 200\ \text{s}]\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\end{array}$

Figure 13: Two masses model with coupling on force, displacement and velocity.
Behavior of the right mass changes after $100$s

First of all, for the sake of performance and accuracy comparisons, the model
have been co-simulated using the NI Jacobi method with five fixed step size
values: $0.01$s, $0.05$s, $0.1$s, $0.2$s and $0.4$s. Then, F3ORNITS method was
used with $\delta t_{1}^{[0]}=\delta t_{2}^{[0]}=0.01$s (also used as ”minimal
step size value”), for different strategies regarding polynomial calibration
(extrapolation / CLS, see 2.2), smoothness enhancement (disabled / make
interfaces $C^{1}$, see 3.1.4) and error normalization (w.r.t. amplitude /
magnitude / damped amplitude with $\nu=5\%$, see 3.2.1).

Figure 14: Evolution of $y_{1,1}$ (coupling variable corresponding to the
position of the left mass $x_{1}$) and its damped bounds sequences a defined
in (36) Figure 15: RMSE/#steps trade-off (on $x_{1}$ and $v_{1}$) for various
sets of parameters of F3ORNITS algorithm on 2 masses model

All co-simulations run instantly (measured elapsed time: $0.0$s), but a good
performance indicator can be the number of steps proceeded. Indeed, the main
cost of a co-simulation run is the frequent restart of subsystems’ solvers
after each discontinuity introduced by a communication time. On the other
side, bigger co-simulation steps are expected to produce inaccurate results:
the larger a co-simulation step is, the older an input value might be used.
The precision criterion will be the RMSE (root mean square error) compared to
the monolithic reference.

The error normalization criterion in the time-stepper has a strong impact on
the results: as expected, ${}^{Magn.}error$ generates larger errors and
smaller co-simulation steps. The total amount of steps is bigger than the
other normalization methods, but the results are also more accurate. The
trade-off is: results are either more accurate than with NI Jacobi for the
same total number of steps, or obtained in less steps than NI Jacobi to
achieve the same accuracy (see right-hand dashed circled sets of points on
figure 15).

In order to reduce the effect of artificial step size reduction when coupling
variables cross zero, the error normalization ${}^{Ampl.}error$ method may be
used in the time-stepper. The steps are larger, and the total numbers of steps
are the smallest ones on every case (compare the $3$rd subtable of table 1 to
the $2$nd and the $4$th ones). However, as the model is damped
($d_{1}=d_{3}>0$), the initial highest and lowest values of the variables are
not forgotten throughout the whole co-simulation (this can be observed on the
blue curve of figure 14). The trade-off is still better than NI Jacobi on
every case.

In order to improve the accuracy of the results obtained using the
${}^{Ampl.}error$, the damped amplitude alternative for error normalization
has been proposed in subsubsection 3.2.1. Indeed, damping the amplitude with a
factor $\nu=5\%$ allows the coupling to progressively forget the initial high
values, as the minimum and maximum values will ”follow” the order of magnitude
of the variables amplitude. The ${}^{\begin{subarray}{c}damp\\\
max\end{subarray}}\\!\\!y_{1,1}^{[n]}$ and ${}^{\begin{subarray}{c}damp\\\
min\end{subarray}}\\!\\!y_{1,1}^{[n]}$ sequences can be seen on figure 14,
$y_{1,1}$ being the coupling variable associated to $x_{1}$.

Table 1 compiles the results of figure 15 and shows precisely the rmse values.
Among others, it shows that the polynomial inputs calibration method (either
with an extrapolation $\Omega_{q}^{Ex}$ or a CLS fitting $\Omega_{q}^{CLS}$)
has minor incidence on the final results on the two masses test case.
Nonetheless, it can be observed on this table that the smoothness enhancement
(presented in 3.1.4) slightly damages accuracy, except when the error is
normalized using ${}^{Magn.}error$. The reason of this is the fact that the
value of the inputs is not the one given as output by its connected subsystem
at the beginning of the step. In other word, this phenomena (that can be
observed at $t_{k}^{[n]}$ on figure 8) generates a small extra error at the
beginning of every co-simulation step, but the value of the input is quickly
adapted to match the one it would have had without the smoothness enhancement.
The smaller the steps are, the less this phenomena can be observed: that is
the reason why the subtable regarding ${}^{Magn.}error$ does not show it. The
benefit of smoothness enhancement is mainly the possibility to restart faster
the solvers inside of the subsystems due to the $C^{1}$ smoothness of the
input variables. The consequence of it would be a faster co-simulation run for
the same number of co-simulation steps.

Table 1: Results on 2 masses test model - comparing number of co-simulation
steps and rmse on $2$ state variables NI Jacobi (ZOH)
---
$\delta t$ | #steps | $\begin{subarray}{c}\text{rmse}\\\ \text{on}\ x_{1}\end{subarray}$ | $\begin{subarray}{c}\text{rmse}\\\ \text{on}\ v_{1}\end{subarray}$
$0.01$ | $20000$ | $0.030\%$ | $0.074\%$
$0.05$ | $4000$ | $0.197\%$ | $0.462\%$
$0.1$ | $2000$ | $0.412\%$ | $0.981\%$
$0.2$ | $1000$ | $0.772\%$ | $2.044\%$
$0.4$ | $500$ | $1.928\%$ | $5.194\%$
F3ORNITS
$\\!\Omega_{q}\\!$ | Smoothness | #steps | $\begin{subarray}{c}\text{rmse}\\\ \text{on}\ x_{1}\end{subarray}$ | $\begin{subarray}{c}\text{rmse}\\\ \text{on}\ v_{1}\end{subarray}$
time-stepper using ${}^{Magn.}error$
$\\!\Omega^{Ex}_{q}\\!$ | disabled | $1935$ | $0.017\%$ | $0.024\%$
enhanced | $1848$ | $0.014\%$ | $0.019\%$
$\\!\Omega^{CLS}_{q}\\!$ | disabled | $1768$ | $0.017\%$ | $0.024\%$
enhanced | $1786$ | $0.017\%$ | $0.024\%$
time-stepper using ${}^{Ampl.}error$
$\\!\Omega^{Ex}_{q}\\!$ | disabled | $669$ | $0.083\%$ | $0.183\%$
enhanced | $674$ | $0.122\%$ | $0.269\%$
$\\!\Omega^{CLS}_{q}\\!$ | disabled | $636$ | $0.083\%$ | $0.188\%$
enhanced | $672$ | $0.104\%$ | $0.227\%$
time-stepper using damped amplitude error
$\\!\Omega^{Ex}_{q}\\!$ | disabled | $989$ | $0.019\%$ | $0.030\%$
enhanced | $1007$ | $0.024\%$ | $0.048\%$
$\\!\Omega^{CLS}_{q}\\!$ | disabled | $1071$ | $0.020\%$ | $0.037\%$
enhanced | $911$ | $0.025\%$ | $0.048\%$

##  5 CONCLUSIONS

The F3ORNITS coupling algorithm used with a time-stepping criterion based on
the local error estimation normalized with regard to magnitude gives a better
error/#steps trade-off that the non-iterative Jacobi method. The number of co-
simulation steps can even be reduced by tuning the error normalization method,
and the involved coefficient (refered to as $\nu$ damping coefficient). In
other words, a safe approach could be to start with a co-simulation with
F3ORNITS using the error normalized with regard to the order of magnitude and,
if the total number of steps is not satisfactory, the strategy can be to use
the damped amplitude with a decent $\nu$ coefficient and to progressively
decrease it ($\nu=0\%$ corresponds to the classical amplitude approach).

On the second test case, with error normalized regarding amplitude damped with
a coefficient $\nu=5\%$ (shapes in the middle dashed circles on figure 15), we
can achieve, with approximately $20$ times less co-simulation steps (and as
many avoided discontinuities), an accuracy similar to the one obtained with
the non-iterative Jacobi method (or ever a bit more accurate results, as non-
iterative Jacobi method does not reach an rmse lower than $0.03\%$ for $x_{1}$
even with $20\ 000$ steps). From another point of view, a similar amount of
communication times produces an error $38$ times smaller with F3ORNITS than
with the non-iterative Jacobi method. Moreover, a co-simulation step size does
not need to be chosen in advance as it is automatically adapted with F3ORNITS
on the subsystems allowing it.

The robustness of the method comes from the motivation at the starting point:
F3ORNITS can be applied regardless of the subsystems and their capabilities.
It is therefore possible to integrate the method on an industrial product that
must handle subsystems coming from a wide range of platforms.

The CLS approach, designed in order to avoid forgetting data used for
calibration, does not show any significant difference compared to the
classical extrapolation approach on the second test case, yet the behavior on
a wider set of models has to be studied.

Regarding the enhancements, we can point out a faster restart of embedded
solvers inside the subsystems when $C^{1}$ smoothness of inputs is garanteed
(presented in 3.1.4) at the communication times: rough discontinuities do not
necessarily occur at every communication time. This can only be considered for
subsystems supporting $3$rd order polynomial inputs.

## REFERENCES

  * Amemiya, 1985 Amemiya, T. (1985).  Advanced econometrics.  Cambridge, Mass.: Harvard University Press.
  * Arnold and Unther, 2001 Arnold, M. and Unther, M. G. (2001).  Preconditioned dynamic iteration for coupled differential-algebraic systems.  Bit, 41(1):1–25.
  * Bartel et al., 2013 Bartel, A., Brunk, M., Günther, M., and Schöps, S. (2013).  Dynamic iteration for coupled problems of electronic circuits and distributed devices.  SIAM J. Sci. Comp., 35(2):315–335.
  * Berrut and Trefethen, 2004 Berrut, J.-P. and Trefethen, L. N. (2004).  Barycentric lagrange interpolation.  SIAM Review, 46(3):501–517.
  * Busch, 2016 Busch, M. (2016).  Continuous approximation techniques for co-simulation methods: Analysis of numerical stability and local error.  ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 96(9):1061–1081.
  * Busch, 2019 Busch, M. (2019).  Performance Improvement of Explicit Co-simulation Methods Through Continuous Extrapolation.  In IUTAM Symposium on solver-coupling and co-simulation, volume 35 of IUTAM Bookseries, pages 57–80. IUTAM.  IUTAM Symposium on Solver-Coupling and Co-Simulation, Darmstadt, Germany, September 18-20, 2017.
  * Dronka and Rauh, 2006 Dronka, S. and Rauh, J. (2006).  Co-simulation-interface for user-force-elements.  In Proceedings of SIMPACK User Meeting, Baden-Baden.
  * Éguillon et al., 2019 Éguillon, Y., Lacabanne, B., and Tromeur-Dervout, D. (2019).  IFOSMONDI: a generic co-simulation approach combining iterative methods for coupling constraints and polynomial interpolation for interfaces smoothness.  In Science, S. and Publications, T., editors, Proceedings of the 9th International Conference on Simulation and Modeling Methodologies, Technologies and Applications, pages 176–186. INSTICC.
  * Gomes et al., 2018 Gomes, C., Thule, C., Broman, D., Larsen, P. G., and Vangheluwe, H. (2018).  Co-simulation: a survey.  ACM Computing Surveys (CSUR), 51(3):1–33.
  * Gu and Asada, 2004 Gu, B. and Asada, H. H. (2004).  Co-Simulation of Algebraically Coupled Dynamic Subsystems Without Disclosure of Proprietary Subsystem Models.  Journal of Dynamic Systems, Measurement, and Control, 126(1):1.
  * Hildebrand, 1956 Hildebrand, F.-B. (1956).  Introduction to Numerical Analysis.  Dover Publications, 2nd edition.
  * Kübler and Schiehlen, 2000 Kübler, R. and Schiehlen, W. (2000).  Two methods of simulator coupling.  Mathematical and Computer Modelling of Dynamical Systems, 6(2):93–113.
  * Li et al., 2014 Li, P., Meyer, T., Lu, D., and Schweizer, B. (2014).  Numerical stability of explicit and implicit co-simulation methods.  J, 10(5):051007.
  * Li et al., 2020 Li, P., Yuan, Q., Lu, D., Meyer, T., and Schweizer, B. (2020).  Improved explicit co-simulation methods incorporating relaxation techniques.  Archive of applied mechanics, 90:17–46.
  * Schierz et al., 2012 Schierz, T., Arnold, M., and Clauß, C. (2012).  Co-simulation with communication step size control in an fmi compatible master algorithm.  pages 205–214.
  * Schweizer et al., 2016 Schweizer, B., Li, P., and Lu, D. (2016).  Implicit co-simulation methods: Stability and convergence analysis for solver coupling approaches with algebraic constraints.  ZAMM Zeitschrift fur Angewandte Mathematik und Mechanik, 96(8):986–1012.
  * Sicklinger et al., 2014 Sicklinger, S., Belsky, V., Engelman, B., Elmqvist, H., Olsson, H., Wüchner, R., and Bletzinger, K.-U. (2014).  Interface Jacobian-based Co-Simulation.  International Journal for Numerical Methods in Engineering, 98:418–444.
  * Skelboe, 1992 Skelboe, S. (1992).  Methods for Parallel Integration of Stiff Systems of ODEs.  BIT Numerical Mathematics, 32(4):689–701.