1. Introduction
In practice, a network of co-simulation slaves is often used to model complex systems by coupling subsystems on the behavioral description level [
1]. The usefulness of this approach comes from the fact that a co-simulation slave can be exported from a simulation tool. Multi-disciplinary teams can use co-simulation to combine information from already developed models from multiple domains. A recent overview of existing co-simulation research can be found in [
2]. A Functional Mock-up Interface (FMI) [
3,
4] has been introduced to standardize the co-simulation interface. This effort allows the coupling of a growing number of commercial simulators [
5].
A co-simulation master is an algorithm responsible for the simulation of a co-simulation network. The master calculates the approximation of input signals and controls the execution of slaves. Each slave has a solver of the internal model to calculate its own state and output updates. The responsibility for the quality of co-simulation results is shared between the master and the solvers. The main objective of this article is to develop a practical co-simulation quality assessment and to illustrate its use in a variable-step co-simulation.
An implicit co-simulation master repeats the simulation steps of the slaves until the coupled inputs and outputs match. An explicit co-simulation master executes a single step and continues the execution regardless of the connection error. A comparison of implicit and explicit masters using the example of a two-mass oscillator can be found in [
6]. The comparison shows that implicit approaches have larger regions of stability than explicit approaches. Furthermore, Ref. [
1] states implicit co-simulation is zero-stable if zero-stable [
7] solvers are used. However, an implicit co-simulation requires the option to roll back a step of a co-simulation slave. That option is defined in the FMI standard but is rarely supported in practice. For this reason, only the explicit co-simulation is analyzed in this article.
Another very important reason to focus on explicit co-simulation is that hardware-in-the-loop simulation [
8,
9] is explicit co-simulation. Hardware co-simulation slaves included in the simulation loop cannot repeat a simulation step. Furthermore, there are very few techniques that can be applied to assess the quality of such a co-simulation. An energy-based quality assessment of an engine-in-the-loop simulation is presented in [
10]. In comparison, this article tries to provide an analysis of the whole co-simulation including integration and output equations, and not just the connections. The quality assessment in [
10] shows how the quality of power bonds [
11] can be analyzed. This article provides a quality assessment applicable to a wider range of co-simulation systems. Connections do not have to be pairs of effort and flow signals. One example of systems that can have connections that are not power bonds is kinematic models in robotics [
12].
Error estimation techniques used in ordinary differential equations provide suggestions for evaluating the simulation quality. Most algorithms for solving ordinary differential equations attempt to control either the local error or the defect of a numerical solution [
13]. Local error estimation techniques based on Richardson extrapolation for the co-simulation have been presented in [
14]. Assuming perfect subsystem integration, that article shows that the global error is bounded in terms of extrapolation error. That technique, however, requires the option to roll back a step of a co-simulation slave.
This article presents a co-simulation quality assessment based on the defect calculation [
15,
16,
17]. The research presented is the continuation of the work presented in [
15]. There, the numerical defect of the co-simulation is analyzed. That analysis showed that for co-simulation, when numerical defects are limited, numerical errors are limited. The main parts of this analysis are referenced in the next section.
The analysis in this article and [
15] is based on coupled ordinary differential equations. Coupled ordinary differential equations are a special case of differential and algebraic equations. An example of such a system is shown in
Figure 1. Ordinary differential equations are used to represent the state equations of systems modeled by co-simulation slaves. Algebraic equations are divided into output and connection equations. This is performed to reflect co-simulation practice, where co-simulation slaves are typically black boxes. State and output equations are not available to the co-simulation master and are therefore colored gray in
Figure 1.
This article proposes an explicit variable step co-simulation method based on numerical defect control. A co-simulation slave is responsible for generating its output signals, while the co-simulation master is responsible to solve the connection equation (
Figure 1). The state and output equations are grayed out to indicate that they are usually not available to the co-simulation master. However, the output defect depends on the co-simulation step size controlled by the co-simulation master. This is why this article assumes that the output defect should be estimated by the co-simulation master and why the output equations in
Figure 1 are lighter colored.
Ref. [
15] provides the basis for the defect analysis of the co-simulation. This article continues that work and proposes how to adapt a step size controller [
18] to co-simulation defect control. Non-iterative co-simulation in [
15] is a fixed-step, multi-rate co-simulation that uses zero-order hold. This article presents a non-iterative single-rate variable-step co-simulation using higher-order interpolation.
The next section shows the basis for the analysis in this work, carried over from [
15]. The defect control section introduces an explicit co-simulation variable-step method. That section shows how to calculate the connection defect and estimate the output defect. The proposed method controls the error by varying the communication step size using the PI controller. A simple application of the method is presented in the example section. This application serves to highlight some of the effects of using different orders of interpolation polynomials in co-simulation. The final section of the article contains conclusions and ideas for future work.
2. Error Bounds
This article extends the work performed in [
15] with variable step co-simulation presented in the next section and experimental analysis in the following section. There, coupled ordinary differential Equation (1) and their numerical solution (2) are used to analyze numerical errors of the co-simulation. An important result of [
15] is Theorem 1. It states that the global error of the co-simulation is limited when the co-simulation defect is limited. It is worth noting that a similar statement is proved in [
19] for a system of differential and algebraic equations. The main difference is that algebraic equations in this article are divided into output and connection equations (
Figure 1). This section repeats the expressions for coupled ordinary differential equations and the error bounds theorem from [
15] (Theorem 1). Coupled ordinary differential equations are the basis for the step size analysis and Theorem 1 justifies the error control presented in the next section.
A co-simulation models a system partitioned into
subsystems connected by the connection function L
where
i is the subsystem index,
is the state signal,
is the output signal,
is the input signal, and
is the initial state of the subsystem. The numerical solution of the system satisfies the following equations
where the numerical solution of the state, output, and input signals is denoted as
,
, and
, respectively. The signals found by the numerical solution are assumed to be piecewise continuous. The defect introduced to the numerical solution is partitioned into integration
, output
, and connection defect
. The system (1) represents the equations solved by the co-simulation and the system (2) represents the behavior of the solution obtained by the co-simulation.
The numerical solution (1) can be rewritten to
where signals of all subsystems in (2) are grouped into large column vector signals
The error of the numerical solution (3) is defined as the difference between the numerical and the analytic solution of the system (1)
where
is the integration error,
the output error and
the input error of the numerical solution.
Definition 1 (Lipschitz Continuity).
A function is said to be Lipschitz continuous if there exist constant such that for all :The constant is called the Lipschitz constant of the function . Definition 1 introduces Lipschitz continuity, which is used to describe the conditions for the uniqueness of the solution for (1) [
15]. It is used to formulate Assumption 1 and Theorem 1.
Assumption 1. Assume that there exists a Lipschitz continuous function that explicitly calculates the input signalsAssume that the aggregated state transition function (
3a)
is Lipschitz continuous. Assume that the aggregated output function (
3b)
is Lipschitz continuous. Assume that the numerical solution (2)
is continuous in every subinterval.
Theorem 1 (Error Bounds = Theorem 2.12 in [
15]).
Suppose Assumption (1)
holds. Then, the integration error is limitedthe input error satisfies is limitedand the output error satisfies is limited In [
15], it is shown that Theorem 1 requires the same conditions as the uniqueness of the solution for (1). In addition, the numerical solution obtained by co-simulation (2) must be piecewise continuous. This theorem provides the justification for the defect control presented next.
3. Defect Control
Theorem 1 suggests that the step size control of the co-simulation defect can be used to limit the global error of the co-simulation. According to [
3], a communication step size is the “distance between two subsequent communication points (also known as sampling rate or macro step size)”, where communication points are “time grid for data exchange between master and slaves in a co-simulation environment (also known as sampling points or synchronization points)”. This section describes how to control the communication step size of a slave using a non-iterative variant of the Jacobi co-simulation method (Algorithm 1). The proposed co-simulation method is shown in
Figure 2. The variable-step variant of the method generates a sequence of communication step sizes
with the explicit version of the PI control procedure introduced in [
18]. The defect controlled by the proposed method is based on the connection defect (
2d) and/or the output defect (
2b).
The numerical defect must be limited to ensure that the co-simulation error is limited. Theorems 2 and 3 show that the defect can be limited by reducing the communication step size. This section shows how to calculate the connection defect and estimate the output defect. Theorem 4 shows that the output defect estimate used is asymptotically correct.
The sequence of communication points
is determined by the communication step sizes
The sequence of points
at which a numeric solution discontinuity can occur is marked differently than the sequence of communication points
. The reason is that discontinuities in the numeric solution may occur during slave integration. Each slave can perform the integration with different individual integration steps (sometimes referred to as a micro step size ). It is assumed that the internal solver of the slave takes over the responsibility for the control of the integration step size. Control of the integration step size may be used to reduce the integration error [
16,
17]. Theorem 1 suggests that integration defect reduction is an important factor for the quality of co-simulation. This article focuses on controlling the size of the communication step with respect to output and connection defects. This choice respects the black box character of co-simulation slaves (
Figure 1) and leaves the integration defects in the responsibility of internal slave solvers.
Co-simulation is the simulation of a continuous time model with a computer, a discrete system. The co-simulation generates samples of the signal and its derivatives at communication points. This description is consistent with [
3]. The reconstruction of output signals from such samples can be obtained using the Taylor polynomial
where the samples of the output signal and its derivatives are determined by the co-simulation slave. The output signal of the
co-simulation slave is a piecewise polynomial signal
The inputs of the
co-simulation slave are extrapolated with the following polynomial
where the samples of the input signal and its derivatives are determined by the co-simulation master. A non-iterative Jacobi co-simulation master (
Figure 2) determines the input signals in the
step with the connected output signals from the
step
The input signal of the
co-simulation slave is a piecewise polynomial signal
The connection defect (
2d) can be calculated by comparing the extrapolation polynomials for the output (12) and input (14) signals, i.e.,
where individual scalar signals are selected according to the connection function.
This article assumes that the output signal and all its derivatives are perfectly sampled, i.e., the output defect can only deviate from zero between the communication points
Lemma 1. Assume that the numerical solution for the input signals is bounded , the numerical solution for the state signal is bounded , the integration defects are , and the state transition function is Lipschitz continuous (Definition 1). Then Proof. Since
there exists
such that for all
By the integration of (
2a)
the following inequality is obtained
The statement of the lemma (19) follows from the previous inequality
□
Theorem 2 (Connection defect).
Assume that the numerical solution for the input signals is bounded , the numerical solution for the state signal is bounded , the integration defects are , the state transition function is Lipschitz continuous (Definition 1), the output function of the subsystem isand the step sizes of connected simulators . The connection defect (
2d)
of a numerical solution converges in terms of the communication step size Proof. From (14)–(18) and (24) it follows that
where
. Since
is Lipschitz continuous it follows that
The statement of the theorem follows from Lemma 1, (26) and (27). □
Theorem 2 specifies the conditions under which the connection defect converges with respect to the communication step size. A number of simplifications are adopted to prove Theorem 2. In (26), the higher order terms are ignored to simplify the proof. The goal was to prove that the connection defect was converging in terms of communication step size, rather than finding the smallest limit. The next section shows the order of convergence of the connection defect in a simple example. Bounded inputs and state signals are reasonable assumptions for a stable system and a stable co-simulation. The assumptions (24), however, restrict the models used to those that have no direct output dependency on input signals. If this simplification is not applied, it is possible to construct a system with an algebraic loop that makes the co-simulation unstable. The subject of future work will be to analyze how some or all of these simplifications can be discarded or at least relaxed.
Theorem 3 (Output defect).
Suppose the function is continuously differentiable and the calculated state signal is continuously differentiable. Then, the output defect (
2b)
is Proof. From (
14) and the fact that
and
are continuously differentiable, it follows that
for
The expression (28) follows from (
2b), (12) and (18). □
Theorem 3 shows that the output defect can be controlled by reducing the communication step size. Theorem 2 and Theorem 3 justify the communication step size control. Numerical solvers are expected to minimize the remaining component of the numerical defect, the integration defect (
2a). The rest of the section analyzes how to estimate the output defect.
For the purpose of estimating the output defect between the communication points, a Hermite interpolation polynomial [
20] is łintroduced
A Hermite interpolation polynomial is consistent with multiple samples of the signals and their derivatives. The polynomial used in this paper is consistent with signal values at two communication points and signal derivatives at the later point
The Hermite interpolation polynomial is used to obtain an asymptotically correct estimate of the output defect.
Theorem 4 (Estimate of the Output Defect).
The estimation of the output defect is defined as the difference between interpolation polynomials (30)
and (12)
Suppose the function is times continuously differentiable on the interval andThen, the estimate of the output defect (32)
is asymptotically correct, i.e., for each Proof. Since
is continuously differentiable, it is also Lipschitz continuous. Lipschitz continuity and (
33) imply
The output defect (
3b) is the difference between the numerical output signal (12) and the output signal without defect (35). The output defect on the interval
is equal to
The estimate of the output defect on the interval
is equal to
The Equation (
34) follows directly from (36) and (37). □
The output defect estimate (32) is obtained by adding communication points during co-simulation. The additional points are added in the middle of the interval
shown in
Figure 2. These points are used to obtain a higher order interpolation polynomial (30) for use in the output defect estimate. Theorem 4 gives the conditions (33) under which the output defect estimate (32) is asymptotically correct. The integration should have a higher order of local error than the order of the output interpolation polynomial (12).
The connection defect calculation (17) and the output defect estimate (32) are used to define the controlled co-simulation defect
The calculation of the co-simulation defect is used in a step size control approach similar to the one introduced in [
18]. The approach uses a PI controller for the logarithm of an error measurement
where
is the controlled error approximation. Such a method has already been used for co-simulation [
21]. The difference to the method presented in this article is the controlled error estimate. In [
21], the authors have used an explicit step size control procedure to control the local co-simulation error of the output signal. In this article, the method controls the maximum of the connection defect and the output defect estimate (38).
The co-simulation method used for demonstrating the communication step size control is defined in Algorithm 1. It is also shown in
Figure 2 to simplify the introduction. Like any other co-simulation master, the method solves the connection Equation (
1d) and controls the execution of the co-simulation slaves. The connection equation is solved by the assignment of input signal derivatives (15). The presented co-simulation method allows a parallel execution with a distribution of the calculation performed in for-loops. The communication step sizes are controlled using the PI controller (39) to keep the co-simulation defect (38) constant.
4. Numerical Example
This section presents an example of using the proposed variable-step co-simulation method (Algorithm 1). The example system is a two mass oscillator that is commonly used to benchmark co-simulation master algorithms [
6,
21,
22,
23,
24,
25]. The slaves are implemented in the C programming language using Functional Mock-up Interface (FMI) [
3,
4]. The proposed connection defect calculation (17) and the output defect estimation (32) are illustrated in this example.
Algorithm 1 The pseudocode describes a non-iterative Jacobi co-simulation method. In the variable-step variant, the step size is computed with the PI controller (17), (32), (38) and (39). In the fixed-step variant, the step size is constant . |
- Require:
a partitioned system (1) without algebraic loops, an initial step size , defect tolerance - 1:
, , - 2:
calculate the initial output signals by solving the Equations ( 1b) and ( 1d) - 3:
repeat - 4:
- 5:
for to N do - 6:
assign the input signals (15) - 7:
for to N do - 8:
integrate the Equation ( 1a) on the interval - 9:
for to N do - 10:
obtain the output signal samples at - 11:
for to N do - 12:
integrate the Equation ( 1a) on the interval - 13:
for to N do - 14:
obtain the output signal samples at - 15:
compute - 16:
- 17:
until
|
The algorithms presented in this article and the code used to generate the results in this section are published at [
26]. The models are implemented in C, the algorithms are implemented in C++, and the figures are created in Python. In the repository [
26] an interested reader can find
a C++ implementation of Algorithm 1 in the template function
fmi::jacobi_co_simulation (src/fmi.h),
an implementation of the step size controller (39) in the method VariableStepSize::next (src/fmi.cpp),
an implementation of the Hermite polynomial calculation (30) in the method FMU::get_hermite (src/fmi.cpp),
an implementation of the output defect calculation in the function fmi::calculate _output_defects (src/fmi.cpp),
an implementation of connection defect calculation in the function fmi::calculate _connection_defects (src/fmi.cpp),
an implementation of co-simulation slaves according to the FMI 2.0 standard [
3] in the directories
src/OscillatorOmega2Tau and
src/OscillatorTau2Omega,
an implementation of the reference solution in the directory
src/TwoMassRotationalOscillator,
and the code to create the images presented in this section in the Python script scripts/results_analysis.py.
The implementation of the co-simulation slaves (
Figure 3b) follows the FMI standard [
3]. Since slaves are compiled together with the master, the shared library and the interface description are not packaged together for ease of implementation. Models solved by the slaves are described in the following equations.
The example system consists of two slaves
that solve the following system of equations
connected by
where
and
Model parameters are set to
The analytic solution is approximated by a monolithic solution of the system (
Figure 3a) solved using the CVODE solver [
27] with a tight tolerance bound. The absolute tolerance limit for the solver used is set to
. This solution is used to give a reference solution for the co-simulation and to approximate the numerical error.
The internal equations of co-simulation slaves are solved using the forward Euler method. The forward Euler method is a first-order numerical solver for the same equations
The state derivative values of the Euler solver are equal to
Output polynomials are Taylor polynomials (12) with output derivative values calculated using (18) and
Input polynomials are Taylor polynomials (14) with input derivative values calculated using (15).
Figure 4a,b shows the piecewise response for both fixed and variable-step co-simulation using Algorithm 1.
Figure 4a shows the responses obtained by fixed-step co-simulation with the step size
Figure 4b shows the responses obtained by variable-step co-simulation with the initial step size and tolerance
The step size is controlled by the step controller (
39) with the control parameters set to
The comparison of figures shows how the step size controller reduces the step size during the co-simulation. In
Figure 4a,b
Theorems 2 and 3 show that the connection and the output defect can be limited by reducing the communication step size. In order to verify this statement, the root mean square value of the connection
and output defect
is calculated for fixed-step co-simulations with different step sizes
The results are shown in
Figure 5a,b.
The output defect (
2b) is estimated using (32). The order of convergence of the output defect is given with Theorem 3.
Figure 5a confirms this theorem. It is interesting to observe the output defect of the Euler solver (45) in
Figure 5a. From (12), (
40b), (45), (46) and (47) it follows that
and
for the Euler solver. This result agrees with the asymptotic upper limit from Theorem 3. This is interesting because it shows that the output defect for such a solver setup is zero for any output interpolation order greater than 0. However, the numerical error is larger than that of an analytic solver (
Figure 6).
Theorem 2 shows that the connection defect converges (by assuming that the output equations are independent of the input signal), but does not show the order of convergence. The input assignment (15) suggests that the input signal depends on the order of the input polynomial as well as the connected output polynomial. However,
Figure 5b shows a relationship between connection defects and an internal solver. The connection defect when using the analytic solver seems to correlate with the interpolation order and appears to be
. In the case of the Euler solver, the connection defect seems to be limited to
. The latter seems to be a consequence of (46).
The focus of future work will be on the investigation of the order of convergence of connection defects. Theorem 2 uses a simplification to avoid the analysis of direct influences of input signals on output signals. This simplification prevents the analysis of algebraic loops, which could make the system unstable. Direct influences of input signals on output signals have an effect on the local co-simulation error shown in [
24]. This indicates that a direct output input dependency could affect the order of the connection defect. This will be a topic for future work, along with the analysis of the influence of the solver on connection defects.
Theorem 1 suggests that by limiting the overall co-simulation defect, the global co-simulation error should be limited. Theorems 2 and 3 show that by limiting the step size, connection and output defects can be limited.
Figure 5 confirms this. The three theorems suggest that by limiting the step size, the global co-simulation error can be limited. This is confirmed by
Figure 6. It is worth noting that the conclusions given apply to a coupled ordinary differential system 1.
Reducing the step size limit is not the only way to reduce connection and output defects. The order of the polynomials used to transmit input and output signals also has an effect. This effect can be observed in
Figure 6. Theorem 1 shows that the co-simulation error is bounded by connection, output, and integration defects. If defects are
then the errors are
All co-simulation errors are limited by the worst co-simulation defect.
Figure 6 shows how the integration defect of the Euler solver limits the co-simulation error to
.
Co-simulation errors can be limited by limiting the co-simulation defects. The rate of convergence of the output defect is given by Theorem 3. It is important to note that the order of convergence for the output defect estimate is influenced by the Euler solver. The Euler solver brakes the assumption (33). This is an example where an integration defect can affect the output error estimate.
It is interesting to observe the effect of the internal solver on the connection defect. In
Figure 5b, it can be seen that for the Euler solver and larger orders of extrapolation for input and output signals the connection defect is
. The connection defect (17) in the explicit co-simulation is influenced by the difference of the state signal at different co-simulation steps (26). This observation suggests that a connection defect can be used to detect numerical errors introduced by solving state, output, and connection equations. The authors plan future work to rigorously analyze whether there are conditions under which this hypothesis is true.
The previous experiments use fixed-step co-simulation to confirm Theorems 1, 2 and 3. Theorem 1 shows that defect control can be used to limit the co-simulation error. Theorems 2 and 3 show that output and connection defects can be limited by limiting the step size. The previous experiments and theorems justify the use of variable step size co-simulation using defect control. The next experiments show the results of applying the variable-step Jacobi co-simulation method to the system (40) and (41). The method is presented in Algorithm 1 and the step size is calculated with (39). The numerical experiment was performed with the reference
, controller parameters (50) and the initial step size
. The comparison of the obtained output signals with the monolithic solution (
Figure 3a) is shown in
Figure 4b.
Next, the tolerance was varied
to demonstrate that such a controller can limit output and connection defects.
Figure 7 shows that the output defect is limited by the tolerance.
Figure 8 shows that the connection defect is limited by the tolerance. This may not always be the case for explicit co-simulation. In the presented experiments, the initial step size was set to small
to ensure that the numerical defect produced in the initial step stays within tolerance. The experiment shows that by controlling (38) both connection and output defects can be controlled.
Furthermore,
Figure 9 shows that by reducing the tolerance, the co-simulation error can be reduced. It is interesting to observe that plots of the output defect (
Figure 7) and connection defect (
Figure 8) are similar shapes to the error plot (
Figure 9). This comparison suggests that there are cases where such variable step co-simulation can be used successfully.
It should be noted that the integration defect (
2a) is not directly controlled in this example. In the case of the analytic solver, the defect is completely eliminated. In the case of the Euler solver (45), the integration defect is
. In practice, co-simulation slaves are black boxes without the ability to monitor the integration. This is why the integration defect is not included in this analysis.
In
Figure 9, it can be seen that co-simulation errors are similar in order of magnitude for different extrapolation orders.
Figure 10 shows the benefit of increasing the extrapolation order. It shows how an average step size during co-simulation depends on the requested tolerance. The step size controller takes smaller steps to achieve the same tolerance if higher extrapolation order is used. This conclusion may not be generalized to more complex subsystems. It shows the idea that an extrapolation order could be used to decrease the CPU and communication network load during co-simulation.
5. Conclusions and Future Work
This article presents an analysis of the co-simulation defect for a system of coupled ordinary differential equations. The analysis is motivated to deepen the understanding of the co-simulation configuration. In practice, co-simulation slaves are black boxes coupled with connection equations (
Figure 1). A quality measure that does not require knowledge of the slave’s internal equations can facilitate the co-simulation configuration. The defect analysis was only applied to the co-simulation in [
15]. This article continues the application of defect analysis and applies it to variable-step co-simulation with different orders of interpolation polynomials.
The main contribution of this article is a non-iterative co-simulation method with variable communication step size (
Figure 2, Algorithm 1). Theorem 1 states that the co-simulation error is bounded if the co-simulation defect is bounded. Theorem 2 and Theorem 3 show that the connection and the output defect can be limited by reducing the communication step size. These theorems justify the use of variable step co-simulation based on defect control.
Section 4 shows an application of the proposed method to an example of a two-mass oscillator and gives a verification of the above statements.
Such a method is valuable in practice because it requires little configuration. The parameters for the procedure are the initial communication step size and the required tolerance . The method does not require a co-simulation slave to repeat a communication step. This relaxes the implementation requirements for co-simulation slaves. In addition, like any variable step method, it can save computation time by calculating the step size for the results of the desired quality.
One goal of future work would be to see if there is a way to eliminate the need to perform additional sampling of the communication points to estimate output defects. It is worth considering under what conditions the calculations of the connection defect are sufficient to assess the quality of the co-simulation.
Another goal of future work would be to focus on the properties of a model and try to estimate the correct initial step size for co-simulation. This would reduce the configuration effort even further and achieve an almost ideal configuration. In this case, only the required quality of the co-simulation is requested by a user .