Analysis of the model-based corrector approach for explicit cosimulation

In many industries, the focus of testing is currently shifting away from classical hardware tests to the virtual verification and validation of products. To this end, cosimulation has become a common tool for the simulation and analysis of complex systems that span multiple engineering domains and usually involve multiple, heterogeneous and application-specific simulation environments. In particular, the so-called explicit cosimulation allows a widespread application since it has minimal requirements regarding the capabilities of the tool interfaces. However, explicit cosimulation also poses a numerical challenge, especially when the system includes stiff coupling loops. The model-based corrector approach presented in Haid et al. (The 5th Joint International Conference on Multibody System Dynamics, 2018) provides a method for the efficient cosimulation of such systems. In this article, this model-based corrector approach is extended to additional extrapolation methods. By modeling the cosimulation process through a linear recurrence equation and applying it to the two-mass oscillator test model, the influence of model-based correction on the underlying extrapolation methods in terms of stability, accuracy, and error convergence is analyzed. It is shown that adding model-based correction can significantly improve the overall cosimulation, allowing >10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$>10$\end{document} times larger macrostep sizes or reducing the cosimulation error by a factor of 10 or more in some cases.


Introduction
Nowadays, many industries are facing increasing product complexity, where a large number of hardware and software components must work together across different engineering domains.While on the one hand, this increases the necessary development effort, on the other hand, development costs have to be reduced, for example by reducing the number of real prototypes and shifting the focus of testing away from classical hardware tests towards the virtual verification and validation of products.In addition, new products must be brought to market in less and less time.To tackle these challenges, a holistic simulation approach that covers all necessary aspects of the final product is necessary.Since in recent years the use of simulation technology has become a common development tool within the individual engineering domains or departments, there are already established, application-specific simulation environments that are usually tailored towards the needs of the individual domain or department.So-called solver coupling or cosimulation offers the possibility to integrate the subsystems of each domain, which are modeled in their respective application-specific simulation environments, into an overall system simulation.This allows a holistic analysis of the overall product behavior, while still exploiting the advantages of the application-specific simulation environments and avoids remodeling of the entire system in a more general, multiphysical simulation tool.
This article starts by giving a short introduction to cosimulation, as well as a rough classification of the general cosimulation approaches and an overview of common methods for explicit cosimulation in Sect. 1.The model-based corrector approach for explicit cosimulation is presented in Sect. 2 and derived for 0th-, 1st-, and 2nd-order polynomial extrapolation.Section 3 introduces the test model and the analysis methodology used in this article.An indepth analysis of the model-based corrector approach when used with 0th-, 1st-, and 2nd-order polynomial extrapolation with regard to stability, accuracy, and error convergence is conducted in Sect. 4. Section 5 takes a brief look at the application of the algorithm to a nonlinear example.Finally, Sect.6 gives conclusions and an outlook on future research regarding the model-based corrector approach.

A short introduction to cosimulation
In the solver coupling or cosimulation approach, each subsystem is solved independently for a small time interval (the so-called macrostep size T [7]) by the solver of the respective simulation environment using the subsystem-specific solver step size (the so-called microstep size δT [7]).Therefore, cosimulation is a multirate, multimethod approach.To simulate the behavior of the overall system, the individual subsystems are connected via inputs u and outputs y, where data is exchanged at and only at the end of each macrostep.As the outputs of a subsystem are generally unknown before the subsystem completes its macrostep, but are usually necessary inputs for the calculation of other subsystems, at least some if not all of these inputs have to be estimated via a suitable extrapolation method (depending on the execution scheme, see Sect.1.2).As these estimated inputs usually differ from the actual inputs an error is introduced during each macrostep, leading to inaccurate simulation results.Effectively minimizing this error is the basic goal of all the different cosimulation methods available today.
Figure 1 shows a block diagram of a simple cosimulation scenario, where 3 nonlinear subsystems represented in state-space form are interconnected via inputs u and outputs y.It is easy to see that, for example, in order to calculate subsystem I I (denoted by superscript Roman numerals) the outputs y I of subsystem I ( y I Fig. 1 Block diagram of a simple cosimulation scenario with 3 nonlinear subsystems represented in statespace form subsystem I in turn also needs outputs from the other two subsystems (y I I 1 and y I I I 1 ) this is not possible without estimating at least some of the corresponding inputs (e.g., u I 1 and u I 2 ).Coordination and interfacing of the subsystems (scheduling, synchronization, and data exchange) is facilitated by a so-called master algorithm.This master algorithm may either be directly included in one of the simulation environments used or provided by a so-called middleware, which is dedicated software that acts as a neutral layer between the simulation environments.Using a dedicated middleware offers advantages in terms of compatibility, usability, and robustness, especially as systems grow more complex (more subsystems and/or a larger number of inputs and outputs per subsystem).
The coupled system shown in Fig. 1  The connections between the subsystems are represented by the algebraic equations formed by the coupling matrix L, which is a special selection matrix with exactly one nonzero entry per row (usually 1), since each input of a subsystem is connected to exactly one output (obviously an output can in turn be connected to multiple inputs).For the coupled system shown in Fig. 1 this coupling equation therefore results in: 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 Note that eqn.(1) represents the coupled system of Fig. 1 without any of the errors introduced by cosimulation, therefore the solution of this DAE would give the true results (the monolithic reference solution without the error introduced by cosimulation).However, since in practical applications the subsystems are simulation models solved by the simulation environments solver this is usually not possible, hence the cosimulation approach.As stated before, the inputs of a subsystem are usually outputs of the other subsystems and therefore unknown before the macrostep is completed and vice versa.In cosimulation, this causality problem is solved by replacing the subsystem inputs u for the duration of the macrostep with estimated inputs ũ calculated by a suitable extrapolation function e that only depends on quantities known before the macrostep, for example, previous outputs or output derivatives.This decouples the DAE eqn.( 1) into an ODE for the duration of the macrostep, for example for the nth macrostep T n → T n+1 : (3) By applying the different cosimulation methods to eqn.(3) their influence on the overall system behavior can be analyzed.
When the subsystems are linear (or linearized), the DAE of the coupled system (see eqn. ( 1)) can be written as (5)

Classification of cosimulation methods
The different cosimulation approaches available today can be roughly classified by a few key features.These key features are whether a cosimulation method is implicit or explicit, sequential or parallel, and fixed step or variable step.When each macrostep is calculated only once this leads to a noniterative coupling scheme called explicit cosimulation [20] or weak or inconsistent approach [23].On the other hand, iteratively repeating macrosteps to improve the solution is commonly referred to as implicit cosimulation [20], Waveform-Iteration [4,22] or a strong or a consistent approach [21].Since in explicit cosimulation each macrostep is calculated only once, the error induced due to the extrapolation of the inputs is carried over into the next macrostep.In contrast, the repetition of macrosteps in implicit cosimulation reduces this error with every iteration.As a result, it is quite easy to achieve high stability and accuracy with an implicit cosimulation approach, but it also requires advanced interface capabilities like the reversion of macrosteps.Explicit cosimulation, on the other hand, has minimal requirements regarding interface capabilities and is therefore applicable for most simulation environments, but,depending on the coupling problem, advanced coupling methods might be necessary to achieve good stability and accuracy.
Depending on whether all subsystem solvers are executed at the same time or one after another, cosimulation methods can also be divided into the parallel Jacobi type and the sequential Gauss-Seidel type [2].While in a parallel Jacobi-type scheme all inputs have to be extrapolated, in a sequential Gauss-Seidel-type scheme some inputs are already known, since the respective subsystems have already calculated the macrostep and the respective outputs are known.As a result, sequential cosimulation is usually more accurate since it requires less extrapolation, which is the major source of error but also requires more time to calculate since only one subsystem is performing its macrostep at any given time.
Finally, the macrostep size used can be either fixed throughout the simulation or dynamically adjusted by a macrostep-size controller (see, e.g., [5,6,16,29,30]).Using a fixed macrostep size is a simple and robust approach, but can also lead to unnecessarily small macrostep sizes in phases with low system dynamics.By using a variable macrostep size larger macrostep sizes can be used during such phases, which usually improves the simulation performance.On the downside, there is a risk that the step-size controller cannot adjust the macrostep size in time if the system dynamics suddenly change, causing large errors or instability in the simulation.
For the purpose of this article, an explicit parallel cosimulation with a constant macrostep size is assumed in order to reduce the complexity of the analysis.Nevertheless, it is easily possible to apply the coupling approaches and methodologies presented in this article to sequential and/or variable step coupling schemes.

Cosimulation methods for explicit cosimulation
Probably the most popular cosimulation approach for explicit cosimulation is to extrapolate the inputs via a polynomial of sufficient order p, where a range from zero to two is most common [10,11].While higher-order polynomials can achieve higher accuracy and better error convergence, they are also more susceptible to oscillations, also known as Runge's phenomenon, which leads to stability problems for relatively large macrostep sizes with respect to the signal frequency and decreases the region of both numerical and zero stability [11].In addition, the subsystem interfaces involved must be able to process these higher-order polynomials by adjusting the input values for each microstep (and each intermediate step in the case of a multistep method) of the subsystem solver according to this extrapolation polynomial.For polynomial extrapolation either macrosteps, microsteps, exact derivatives, or a combination of these can be used.While access to macrosteps is a basic capability of all subsystem interfaces used for cosimulation, access to microsteps or even exact derivatives of the outputs is rarer, but if available, the use of microsteps or exact derivatives can often increase the accuracy of the extrapolation [7].
While (higher-order) polynomial extrapolation adjusts the input values only as a function of time, more advanced approaches like Semi-Implicit Coupling (SIC) [32] or Linear Implicit Stabilization (LIS) [3] also adjust the inputs as a function of the subsystem outputs.To this end, SIC uses a first-order multivariable Taylor polynomial with the corresponding directional derivatives to obtain a first-order approximation of the subsystem inputs with respect to the outputs.On the other hand, LIS uses a local linear representation of the residual system in ODE form, which is solved in conjunction with the subsystem itself by the subsystem solver.Both, SIC and LIS, require model knowledge in the form of directional derivatives, and in the case of LIS, access to the subsystem states is also necessary.Since both methods must be embedded in the subsystems, the implementation of the subsystem interfaces can be challenging.
Another method that uses model knowledge is the so-called Model-Based Coupling approach (MBC) [33] originally introduced for real-time cosimulation with test benches or hardware-in-the-loop systems.Since real-time applications involve additional delay due to signal latency or signal loss, MBC uses a local linear model to predict the subsystem inputs of multiple future macrosteps.The linear model is identified during runtime using either Kalman Filters or Recursive Least Squares [34].
The coupling error caused by the difference between the extrapolated inputs and the actual inputs can also be interpreted as residual energy, which is introduced during each macrostep.Since after each macrostep the actual inputs are known, this residual energy can be estimated.Energy-conservation methods like the Nearly Energy Preserving Coupling Element (NEPCE) [7] compensate this residual energy by adding an appropriate correction signal to the extrapolated inputs during the next macrosteps.An extension to NEPCE that also takes the direct feedthrough of the subsystems into account is given in [28].Naturally, these methods can correct the residual energy only after the fact.The Model-Based Pre-Step Stabilization (MBStab) [13] on the other hand uses model knowledge to approximate the monolithic reference solution and then optimizes the extrapolated inputs such that the difference between the approximated reference solution and the cosimulation becomes minimal.Therefore, this approach already compensates the residual energy during the macrostep in which it is produced, which leads to very good performance in terms of stability and accuracy.
Since the master algorithm is responsible for exchanging the subsystem outputs after each macrostep, these outputs can be adjusted during the data exchange to correct for the error in the extrapolated inputs.Parallel Linear Correction (PLC) [1] estimates and corrects for the error that was transferred directly to the outputs due to direct feedthrough.The Model-Based Corrector (MBCorr) [17] takes an extended but similar approach and is explained in detail in the next section.

The model-based corrector approach for explicit cosimulation
As stated before, the extrapolation of the inputs required for explicit cosimulation induces an error to the overall system with every macrostep.The model-based corrector approach for explicit cosimulation addresses this problem by estimating this induced error and correcting the outputs of the subsystems accordingly after every macrostep.Since the states of the subsystems usually cannot be accessed by the master algorithm in practical applications, the error introduced to the states is corrected during the next macrostep by applying a constant correction signal δu c to the subsystem inputs.This two-step process is shown in Fig. 2.

Correcting the error in the outputs
For the correction of the outputs y c = y + y, again consider the nonlinear coupled system eqn.(1).For this section, all true inputs, outputs, and states of the coupled system are denoted by an overline.The exactly solved system that is therefore unaffected by any cosimulation approach can be written as: Fig. 2 The two-step process of the model-based corrector approach This nonlinear coupled system is then linearized for the duration of the current macrostep Note that in practice it is necessary that the subsystems supply their Jacobians to the master algorithm after every macrostep.Alternatively, system identification methods like Multivariable Output Error State Space identification (MOESP) [19,25] or Recursive Least Squares (RLS) [24] can be used in order to acquire the subsystem Jacobians.The master algorithm then combines the subsystem Jacobians into the block diagonal matrices A n , B n , C n and D n , respectively, (see eqn (5)).
With the extrapolated inputs ũ = Le(y n , . . ., y n−p ) the corresponding cosimulated system can be written in linearized form as well by using the same Jacobians: Subtracting eqn.(8) from eqn. (7) gives the error in the states x = x − x and the outputs y = y − y due to the cosimulation approach: The error in the outputs at the end of the macrostep T n → T n+1 can directly be calculated as: The error in the states at the end of the macrostep (x n+1 − x n+1 ) can be calculated by solving the state error differential equation eqn.(9) for the duration of the macrostep T n → T n+1 .The solution of this type of ODE is given, for example, in [26], for eqn.(9) this results in: Since the model-based corrector approach only corrects for the error induced during the current macrostep (the local error), the state error at the beginning of the macrostep x n − x n is assumed to be 0, which eliminates the first term of eqn.(11).To solve the convolution integral in the second term, both the true inputs u and the extrapolated inputs ũ are approximated by a Taylor-series expansion of sufficient order p.For example, for the extrapolated inputs this Taylor series is given as: As a result, eqn. ( 11) can be written as: By inserting eqn.( 13) back into eqn.( 10) and using the Taylor-series expansions for u and ũ at T n+1 , the estimate for the local error in the outputs can be calculated as where the extrapolated inputs and their derivatives depend on the actual cosimulation method used, while the true inputs and their derivatives must to be approximated via a suitable method.For now, let us assume that the extrapolated inputs are constant for the duration of the macrostep and equal to the corresponding true outputs of the previous macrostep plus the constant correction offset ( ũn = Ly n + δ δ δu c,n and un . . .ũ(p) n = 0), which results in the simplest cosimulation approach, the so-called Zero-Order Hold (ZOH).By assuming the true inputs as piecewise constant as well and equal to the corresponding true outputs of the current macrostep (u n = Ly n+1 and un . . .u (p)  n = 0), eqn.( 14) can be written as which can be solved for the true outputs as where G ZOH,n denotes the correction matrix for ZOH for the nth macrostep.By looking at the definitions of the output error y = y − y and the output correction y c = y + y it can be seen that the true outputs calculated by eqn.( 16) are actually the desired corrected outputs for the macrostep.
If instead of ZOH the extrapolation function is a 1st-order polynomial this leads to a so-called First-Order Hold (FOH) approach.By constructing this 1st-order extrapolation polynomial through the outputs of the two previous macrosteps y n and y n−1 (which is the most common approach) and by approximating the true outputs also with a 1st-order polynomial through y n+1 and y n , the extrapolated and true inputs plus their derivatives result in: Inserting both into eqn.( 14) and again solving for the true outputs (aka the corrected outputs) at the end of the current macrostep y n+1 yields where G F OH,n denotes the correction matrix for FOH for the nth macrostep.For a higherorder polynomial extrapolation approach the same process can be repeated in order to obtain the appropriate calculation rule for the corrected outputs.For example, for Second-Order Hold (SOH) this calculation rule can be given as: Note that the model-based corrector approach is only applicable if eqn. ( 14) can be solved for the corrected outputs y n+1 for the actual cosimulation method, which might not be possible in all cases.Furthermore, the calculation rule for the corrected outputs does not need to depend exclusively on the previous outputs, it can also depend on any known quantity.

Correcting the error in the states
As can be seen from eqn. (11), the difference between the true inputs u and the extrapolated inputs ũ accumulates an error in the states during the macrostep.As stated before, the accumulated error in the states at the end of the macrostep cannot be corrected directly, since in practical applications the states of the subsystems are usually not directly accessible to the master algorithm.Instead, the error in the states is corrected similar to an energy-correction scheme, by applying a constant correction signal δ δ δu c to the extrapolated inputs during the next macrostep (see [7] and [28]).To this end, the accumulated error E r induced by the extrapolated inputs is tracked over the course of the simulation, which can be calculated for the current macrostep T n → T n+1 as: A certain fraction α of this accumulated error is then compensated during the next macrostep by choosing a correction signal δ δ δu such that Since the model-based corrector approach uses a constant correction offset δ δ δu c for each macrostep (therefore T δ δ δu c,n+1 = αE r n+1 , T δ δ δu c,n = αE rn , etc.), inserting eqn.(20) into eqn.( 21) yields the constant correction offset for the next macrostep δ δ δu c,n+1 : Note that setting α = 1 is the most common choice that results in instant correction of the accumulated error during the next macrostep, but it might also be helpful to over-or undercompensate in some cases.The true and extrapolated inputs in eqn.(22) are again approximated via a Taylor-series expansion (see eqn. ( 12)) yielding where again the extrapolated inputs and their derivatives depend on the actual cosimulation method used, while the true inputs and their derivatives must be approximated via a suitable method.
For example, for ZOH the extrapolated inputs and their derivatives are again given as ũn = Ly n + δ δ δu c,n and un . . .ũ(p) n = 0, while in contrast to the output correction the true inputs are approximated with a 1st-order polynomial though y n+1 and y n (see eqn. ( 17)), which results in: Similarly, for FOH the extrapolated inputs are given by a 1st-order polynomial, while the true inputs are approximated by a 2nd-order polynomial, which results in: Again, deriving the correction offset for other extrapolation methods follows the same process.

Summary of the model-based corrector algorithm
Finally, the model-based corrector approach consists of the following steps, which are performed by the master algorithm after every macrostep: 1. Collect the current outputs of each subsystem in the output vector y n+1 .
Since the model-based corrector approach only requires access to the inputs and outputs of the subsystems, at least when used in conjunction with a system-identification method, it is applicable to a wide variety of practical applications, as this is a basic capability that virtually all cosimulation interfaces provide.If a subsystem provides advanced capabilities like access to the subsystem Jacobians, these may be used as well.Either way, exotic interface capabilities such as direct access to the subsystem states are not required.

Test model and analysis methodology
To analyze the stability, accuracy, and error convergence of the model-based corrector approach and compare it to other coupling methods, the well-known two-mass oscillator with force-displacement coupling is used as a test model.This test model can be interpreted as an extension of Dahlquists test equation [18] commonly used for the analysis of numerical integration methods to cosimulation [11].Finally, such a linear test model allows for highly efficient analysis, since in this case the cosimulation process for most cosimulation methods can be described by a linear recurrence relation of the form where the so-called propagation matrix Φ Φ Φ advances the solution of the cosimulated system from one macrostep to the next, resulting in a discrete formulation of the cosimulated system's behavior, where the discrete steps represent the macrosteps of the cosimulation.Note Fig. 3 Two-mass oscillator test model separated by a force-displacement coupling approach (see, e.g., [12]) that the subsystem ODEs are solved exactly and therefore the cosimulation approach is the only source of error.As expected, the propagation matrix changes depending on the cosimulation method, the subsystem parameters, and the macrostep size T .All necessary data is collected in the information vector z, which represents the overall state of the cosimulated system after each macrostep.Using this propagation-matrix approach makes the analysis of the cosimulation stability, accuracy, and error convergence highly efficient since the cosimulation process is fully described by a series of matrix multiplications in which the subsystems are always solved exactly without the need to choose a sufficiently small microstep size as in classical numerical cosimulation.

The two-mass oscillator test model
The linear two-mass oscillator with force-displacement coupling shown in Fig. 3 is probably the most commonly used test model in cosimulation related contributions in the literature (see [1,7,9,12,14,15,31] among others).As shown in Fig. 3, the two-mass oscillator is split into two subsystems at the connection of mass m 1 and the coupling spring damper element c k , d k .Therefore, subsystem I provides the displacement s 1 = y I 1 and velocity v 1 = y I 2 of mass m 1 to subsystem I I , which returns the coupling force F k = y I I 1 (see Sect. 1.1 for details on the notation).The DAE for the coupled two-mass oscillator then results in (see eqn. (4)): For the purpose of this article the parameters of the two single-mass oscillators are set to  7 and T ∈ 3e −3 , 1e −6 .

Cosimulation stability in the time domain
As mentioned before, cosimulation induces an error that can lead to the instability of the system.Therefore, a fundamental requirement for the applicability of a cosimulation method for a given problem is the stability of the cosimulated system.
Given a set of initial conditions (an initial information vector z 0 ), the linear recurrence relation eqn.(26) completely defines the solution of the cosimulated system, for example for n macrosteps (from T 0 to T n ) as z n = Φ Φ Φ n z 0 .If the resulting information vector is bound by some real value ε for an infinite number of macrosteps, the cosimulated system is stable and, since the coupled system itself is assumed to be stable, therefore the cosimulation approach is stable, which leads to the following condition: In order for the power sequence of the propagation matrix Φ Φ Φ to be bound, the spectral radius of Φ Φ Φ has to be either less than or equal to 1, meaning all eigenvalues of Φ Φ Φ have to be inside the unit circle: Note that multiple eigenvalues exactly on the unit circle (λ 1..n = 1) must be nondefective, meaning their algebraic and geometric multiplicity must be the same.By evaluating eqn.( 29) it can be directly seen whether a cosimulation method is stable for the given coupled system and macrostep size.A basic requirement for any cosimulation method to be feasible is to fulfill eqn.( 29) at least as the macrostep size approaches 0. This property is called zero stability [21].

Cosimulation accuracy in the time domain
While stability is a basic prerequisite for any cosimulation to be feasible, stability alone does not automatically imply good results.For example, consider a cosimulation where all subsystems themselves are stable.Then, a cosimulation method that just sets all coupling signals to 0 would provide a stable cosimulation for arbitrarily large macrostep sizes.Therefore, the accuracy of the cosimulated system must also be considered.A straightforward approach is to calculate a reference solution and compare it to the results of the cosimulation.In the case of a linear test model such as the two-mass oscillator this reference solution can also be written as a linear recurrence equation similar to eqn.(26) as where Φ Φ Φ * denotes the propagation matrix and z * denotes the information vector of the reference solution (denoted by a superscript * ).By assuming a common information vector z n at the beginning of the macrostep, the local error of the current macrostep ε ε ε l,n+1 can be directly calculated: By using separate information vectors for the cosimulated and reference solutions the global error for the current macrostep ε ε ε g,n+1 can be calculated as: Equations ( 31) and ( 32) produce an error vector with an error value for every element of the information vector (all states, outputs, previous outputs, etc.) for every macrostep.The resulting k × n error matrix produced by a cosimulation with n macrosteps has to be reduced by an appropriate error norm.In this article, a normalized root mean square error (NRMSE) is used to combine the macrosteps into a single normalized error vector, where the empirical standard deviation is used for normalization.The resulting normalized error vector is then combined into a scalar error indicator via a root mean square: Since in practical applications the states of the subsystems are not available, only the outputs of the current macrostep are included in the calculation of the error indicator in this article.
To obtain the propagation matrix for the reference solution Φ Φ Φ * , the coupling equation is inserted into both the state differential and the output equation in the linear coupled system DAE eqn.(4), resulting in: Inserting this modified output equation in the state differential equation yields the coupled system DAE in the form of an autonomous linear system ODE, since the coupling equation is now implicitly included: The solution of this ODE for the duration of a macrostep T n → T n+1 is once again given in [26] (see Sect. 2.1) as which can be written in propagation matrix form as: Note that in the reference solution the next macrostep depends solely on the previous states x n and therefore the columns corresponding to the previous outputs y n are zero but are kept to match the dimension of the information vector shared with the cosimulation approach.As expected, since ρ(Φ Φ Φ * ) = ρ(A * d ), the exact solution is stable for all macrostep sizes if the underlying system is stable.

Analysis of the model-based corrector approach
As shown in Sect.2, the model-based corrector approach is always used in conjunction with a specific extrapolation method.In this section, the influence of model-based correction on the stability, accuracy, and error convergence of polynomial extrapolation (ZOH, FOH, and SOH) is analyzed.To apply the methodology presented in Sect. 3 for the analysis of the respective polynomial extrapolation approaches with and without model-based correction, the corresponding propagation matrices have to be determined.
Again, starting with the linear coupled system DAE eqn.( 4), replacing the inputs u with extrapolated inputs ũ and solving the resulting state ODE for the duration of a macrostep T n → T n+1 in the same way as the state error differential equation eqn.(9) (see Sect. 2.1) yields: (38) Again, the extrapolated inputs ũn and their derivatives un . . .ũ(p) n depend on the actual cosimulation method used.
For example, for the model-based corrector approach when used with ZOH, the extrapolated inputs for ZOH ũn = Ly n + δ δ δu c,n and their derivatives un . . .ũ(p) n = 0 are inserted into eqn.(38): Inserting the resulting output equation into the output correction equation for ZOH (eqn.( 16)) yields which can then be inserted into the calculation rule for the correction offset δ δ δu c,n+1 for ZOH (eqn.( 24)) resulting in: The equations above can be combined in the propagation matrix for the model-based corrector approach used with ZOH as: By setting the correction matrix G ZOH to 0 and eliminating the correction offset δ δ δu by setting H x , H y and H δu to 0, the propagation matrix for ZOH without model-based correction can directly be seen from eqn. (42).
Note that in the case of linear subsystems with perfect model knowledge K y = G ZOH and therefore M y = M δu = 0, which effectively eliminates the influence of the previous outputs y n and the correction offset δ δ δu c,n in the calculation of the current outputs y n+1 .Therefore, any error in the previous outputs is not directly propagated into the next macrostep.The distinction between K y and G ZOH is kept due to the difference in meaning and to include imperfect model knowledge.The propagation matrices for FOH or SOH with and without model-based correction are derived in the same way by using the respective formulations for the extrapolated inputs, their derivatives plus the output correction and the correction offset.

Stability of the model-based corrector approach
As described in Sect.3, the coupling stiffness c k and damping d k of the two-mass oscillator (see Sect.T are uniformly distributed in logarithmic space such that there is a factor of √ 10 between adjacent macrostep sizes, with a corresponding gray scale starting from black (largest) to light gray (smallest).To cover the wide range of values all three axis (c k , d k and T ) are in logarithmic scale.
Figure 4 shows the resulting contour plots that allow an easy analysis of the influence of model-based correction on the stability regions of the different extrapolation methods.In addition, Table 1 lists the macrostep size T , where ρ(Φ Φ Φ) = 1 for several combinations of c k and d k .
As can be seen in the lower-right quarter of Fig. 4(b), adding model-based correction to ZOH largely improves stability for combinations of high coupling stiffness c k and low coupling damping d k .Furthermore, there is a large improvement in stability for relatively large macrostep sizes along both the c k and d k axis (see lower-left quarter).Finally, there is also some improvement in stability for large damping, low stiffness combinations.
Figure 4(d) shows the effect of model-based correction on FOH extrapolation, which is in principle very similar to ZOH (large improvement for high c k , low d k or large T ; some improvement for high d k ), but overall the improvement is not as large as for ZOH.
Finally, it can be seen from Fig. 4(f) that the stability improvement for high coupling stiffness, low coupling damping combinations due to model-based correction has almost disappeared for SOH, but there is still some improvement for high coupling damping scenarios.
As a result, adding model-based correction to polynomial extrapolation can greatly improve stability for high coupling stiffness, low coupling damping combinations, or when relatively large macrostep sizes are used, especially for low polynomial orders.Furthermore, model-based correction also improves stability for high coupling damping and higher polynomial orders.
It can also be seen in Fig. 4 that in some cases there are small stable areas for high stiffness, high damping scenarios for relatively large macrostep sizes (most prominent in Fig. 4(c)).Furthermore, in the case of ZOH with correction (Fig. 4(b)) the stability for high damping is better for the largest macrostep size used than for the two next smaller macrostep sizes.Since the local error introduced by the cosimulation does not necessarily work towards the destabilization of the system, but might also have a stabilizing effect in some cases, the spectral radius for a given combination of c k and d k does not monotonically decrease as the macrostep size decreases.For example, the spectral radius ρ(Φ) with respect to the macrostep size T for FOH without correction (see Fig.  , which is the same macrostep size that exhibits a small stable area in Fig. 4(b) for this combination of c k and d k .The practical consequence of this effect is that the assumption that given a stable macrostep size for a certain case all smaller macrostep sizes are also stable is false in these cases.ZOH without correction 3.5e −3 1.8e −3 4.9e −6 7.2e −5 1.7e −3 4.8e −6 N/A N/A N/A ZOH with correction 1.0e −2 6.7e −3 1.4e −5 1.5e −3 7.0e −3 1.4e −5 4.9e −6 3.7e −6 6.3e −7 ZOH macrostep gain 2.9 3 .7 2 .9Also note that since the model-based corrector approach corrects the local error induced during each macrostep, the amount of correction and therefore the improvement due to model-based correction depends on the effectiveness of the underlying extrapolation method.The better the underlying extrapolation method performs for the actual case, the less local error it induces during each macrostep, and the less correction has to be applied by the model-based corrector leading to a smaller improvement.For example, it can be seen in Fig. 4 that for low coupling damping d k the stability improves with higher polynomial order, therefore the improvement due to model-based correction gets smaller as the polynomial order increases in this case.

Accuracy of the model-based corrector approach
For the analysis of the influence of model-based correction on the accuracy of polynomial extrapolation methods, again the coupling stiffness c k and damping d k of the two-mass oscillator (see Sect. 3.1) are varied in the ranges c k ∈ 1, 10 10 , d k ∈ 1, 10 7 for several macrostep sizes T in the range T ∈ 10 −2.5 , 10 −6 in the same way as for the stability analysis (see Sect. 4.1).
Unlike the stability limit, the accuracy limit is not fixed, but a matter of definition.In this analysis, the cosimulation is considered accurate if the global normalized root mean square error indicator e g,NRMSE (see Sect. 3.3) is below a threshold of e g,NRMSE ≤ 0.01.The error threshold was chosen in such a way that the signal plots of the subsystem outputs of both the cosimulation and the reference solution are indistinguishable without zooming in, which represents a practical choice for the acceptable error in many industrial applications (assuming the limits on the y-axis are the signal min and max values and the x-axis is reasonably scaled, e.g., from 0 s to 0.3 s in this case).The accuracy limit where e g,NRMSE = 0.01 is marked by a corresponding contour line for each macrostep size T .Figure 6 shows the As can be seen from Fig. 6(b), similar to stability adding model-based correction to ZOH leads to a large accuracy improvement for high coupling stiffness, low coupling damping combinations (lower-right quarter), and some accuracy improvements for high coupling damping (upper-left quarter).Furthermore, model-based correction greatly improves accuracy for large macrostep sizes.It can be seen that without model-based correction there are no accurate solutions for macrostep sizes of T ≤ 3.16 • 10 −4 s, while there is already a substantial accuracy region for T = 3.16 • 10 −4 s when model-based correction is applied.To achieve the same accuracy region without model-based correction a macrostep size of about T = 1 • 10 −5 s is necessary.Figure 6(d) shows the effect of model-based correction on FOH extrapolation.Again, there is a large improvement for high c k , low d k combinations or large macrostep sizes and also some improvement for high damping.Similar to stability, the accuracy improvement for FOH is not as large as for ZOH Finally, it can be seen from Fig. 6(f) that in contrast to stability, there is still some accuracy improvement for high coupling stiffness, low coupling damping combinations when model-based correction is applied to SOH, although the effect is further reduced compared to FOH or ZOH extrapolation.Also, there is an accuracy improvement for large macrostep sizes and high damping scenarios, but again the effect is reduced compared to lower-order polynomial extrapolation in both cases.
It can be seen in Fig. 6 (e.g., (a), (b), (d), and (f)) that for certain macrostep sizes there are bulges, dents, or branches in the contour line.This is caused by the fact that for each coupling method analyzed in this section, the global error reaches a local minimum for a certain coupling stiffness c k .Changes in the coupling damping d k or the macrostep size T move this local minimum mostly up and down as they change the magnitude of the global error.Therefore, in most cases, this local minimum is either completely above or below the defined accuracy limit and thus invisible, but when it is close to the defined accuracy limit, it causes the observed effects, depending on the exact shape and position.Figure 7 depicts the global error e g,NRMSE with respect to the coupling stiffness c k for SOH with correction (see As a result, adding model-based correction to polynomial extrapolation can greatly improve accuracy for high coupling stiffness, low coupling damping combinations, or when relatively large macrostep sizes are used, especially for low polynomial orders.Furthermore, model-based correction also improves accuracy for high coupling damping and higher polynomial orders.This is the same observation as for stability, although in general, the improvement in accuracy due to model-based correction is larger than the improvement in stability.Also, note that the correlation between higher polynomial order and smaller improvement due to model-based correction is much clearer for accuracy than for stability, since the correlation between higher-order polynomial extrapolation and improved accuracy is much clearer in this case.

Error convergence of the model-based corrector approach
Finally, the influence of model-based correction on the convergence of the coupling error as the macrostep size is reduced shall be analyzed.To this end, the coupling stiffness c k and the coupling damping  3 shows the local error for several given macrostep sizes as well as the necessary macrostep size for several desired local errors for each method.
Figure 8 shows the convergence of the local error for polynomial extrapolation for 0th order to 4th order with and without model-based correction.It can be seen that adding modelbased correction increases the error convergence order of the cosimulation such that the convergence of the next higher polynomial is achieved.For example, ZOH with modelbased correction achieves the same error-convergence order as FOH without model-based correction, and so on.In addition, applying model-based correction to polynomial extrapolation also achieves a constant offset regarding the local error (in log scale) compared to the next higher polynomial extrapolation order.As a result, adding model-based correction to a  Finally, model-based correction also stabilizes the cosimulation for large macrostep sizes, especially for high polynomial extrapolation orders (e.g., 3rd-order hold or 4th-order hold), which otherwise suffer instability due to Runge's phenomenon.It can be seen that, for example, for 4th-order hold the local error increases rapidly for macrostep sizes of T >∼ 1.8 • 10 −4 s.

Application to a nonlinear two-mass oscillator
While the use of linear test models for the analysis of a cosimulation method is common practice since they represent a well-understood and well-defined subset of possible test models, the validity with respect to nonlinear systems is highly dependent on how well these systems can be linearized.Since the model-based corrector approach is derived by assuming a linearized coupled system for the duration of the current macrostep (see Sect. 2) the applicability of the algorithm to nonlinear systems depends on how well this linearization captures the behavior of the actual system in the current macrostep.This raises the question of the robustness of the model-based corrector regarding imperfect model knowledge, e.g., due to nonlinearity in the system, and is beyond the scope of this article.Therefore, the model-based corrector approach is applied to a nonlinear test model as an example in this section.
To this end, the two-mass oscillator introduced in Sect.3.1 is used as a basis, but the constant coupling damping d k is replaced by the following calculation rule resulting in a nonlinear spring damper characteristic: Figure 9 shows the displacement s 1 of mass m1 for polynomial extrapolation from 0th order to 2nd order with and without model-based correction for several values of the nonlinearity factor n d .For these cases the global error e g,NRMSE is also listed in Table 4.It can be seen that for ZOH the improvement due to model-based correction becomes smaller with increasing nonlinearity factor n d , up to the point where ZOH without correction becomes unstable and adding model-based correction stabilizes the cosimulation.Therefore, model-based correction significantly improves the stability for high nonlinearity factors in this case.The same effect can be observed for higher-order polynomial extrapolation, where model-based correction can stabilize the cosimulation for higher nonlinearity factors up to a certain point.Also, note that the global error grows with increasing nonlinearity factor n d for all cosimulation methods presented.In addition, it can be seen that increasing nonlinearity proves more problematic for higher polynomial orders in this example, since polynomial extrapolation also assumes some degree of linearity in the system.
As a result, model-based correction seems perfectly applicable also for nonlinear problems, but of course a detailed analysis in this area still needs to be carried out.

Conclusion and outlook
This article demonstrates that the model-based corrector approach presented in [17] for ZOH can be easily applied to other extrapolation methods, for example to higher-order polyno-  mial extrapolation like FOH or SOH.By modeling the cosimulation process via a linear recurrence equation of the form z n+1 = Φ Φ Φz n where the coupled system with the respective cosimulation method is represented by a propagation matrix Φ Φ Φ, the stability, accuracy and error convergence of this cosimulation method are efficiently analyzed without any side effects from numerical solvers.Through the application of this methodology to the wellknown two-mass oscillator test model, which is coupled with ZOH, FOH, and SOH extrapolation with and without model-based correction, respectively, it can be seen that modelbased correction can significantly improve the stability, accuracy, and error convergence of the underlying extrapolation method (see key takeaways below).While the model-based correction algorithm itself consists only of a series of matrix operations and is therefore highly efficient (provided the number of outputs is not too large), obtaining the necessary subsystem Jacobians can be computationally expensive since either the subsystem has to calculate its Jacobians itself or they have to be calculated through a system-identification approach like MOESP.Depending on how fast the subsystem dynamics change, it might be possible to keep the Jacobians for several macrosteps, which can significantly reduce the computation cost.One possible approach is to check how well a linear model using the "old" Jacobians predicts the current macrostep, for example by using a so-called extended state-space model (ESSM [27]), and refresh or keep the Jacobians based on this prediction error (see [17] for details).
Nevertheless, in some cases, there can be a significant overhead compared to simpler coupling algorithms such as plain polynomial extrapolation.Therefore, model-based correction is best used when either the macrosteps of the subsystems require a high computational effort themselves, in which case model-based correction can speed up the simulation if it allows larger macrostep sizes, or the subsystems have technical limitations regarding minimum macrostep size or maximum input extrapolation order, in which case model-based correction can still allow accurate results although the trivial solution (smaller macrostep sizes, higher extrapolation order) is not feasible.To conclude, the overhead due to the model-based corrector approach has to be justified, but if it is, the approach can significantly improve the quality of the cosimulation.
The current work focuses on an indepth robustness analysis of the method with respect to imperfect model knowledge due to system identification and/or underlying nonlinear subsystems.A case study using different examples from the literature and industry will also be conducted in the near future.In addition, the state error correction is based on a correction offset that is calculated only from the values of the past and current corrected outputs.Future research will therefore try to incorporate model knowledge, which is required for the output correction in any case, into the correction of the states as well.

Key takeaways:
-Model-based correction can be easily extended to other extrapolation methods.
-Modeling the cosimulation process as a linear recurrence equation allows the efficient analysis of cosimulation methods.-Model-based correction offers significant stability and accuracy improvement in some cases (mainly high stiffness, low damping combinations, or large macrostep sizes).-In other cases, model-based correction still offers some stability and accuracy improvement (mainly low stiffness, high damping combinations, or balanced stiffness and damping).-Lower-order polynomial extrapolation generally benefits more from model-based correction.-The improvement in accuracy due to model-based correction is generally more pronounced than the stability improvement.-Error convergence is increased by one order when model-based correction is used.-Model-based correction offers a significant improvement of the necessary macrostep size or produced error when high accuracy requirements or small macrostep sizes are given.-Future research focuses on robustness with respect to imperfect model knowledge due to system identification and/or underlying nonlinear subsystems and on improving the state error correction.
Funding Note Open access funding provided by Graz University of Technology.

1 y I I 2 T 3 T 3 T
can also be written as a set of differential algebraic equations (DAEs) I , u I f I I x I I , u I I f I I I x I I I , u I I I I , u I g I I x I I , u I I g I I I x I I I , u I I I f (x, u) y = g (x, u) u = Ly , (1) where for example y I I denotes the output vector y I I , u I I denotes the input vector , and x I I denotes the state vector x I I 1 of subsystem I I .
where the system matrices A, B, C and D are block diagonal matrices that consist of the respective subsystem matrices:

1
Ns m and d 2 = 2 Ns m , while the parameters of the coupling spring damper element c k , d k and the macrostep size T are varied in the range c k ∈ 1, 1e 10 , d k ∈ 1, 1e 3.1) are varied in ranges c k ∈ 1, 10 10 , d k ∈ 1, 10 7 for several macrostep sizes T in the range T ∈ 10 −2.5 , 10 −6 .The stability limit where ρ(Φ Φ Φ) = 1 (see Sect. 3.2) is marked by a corresponding contour line for each macrostep size T .The macrostep sizes 4(b)) and c k = 10 9 , d k =

Fig. 4
Fig. 4 Stability regions for ZOH, FOH and SOH with and without model-based correction

Fig. 6
Fig. 6 Accuracy regions for ZOH, FOH and SOH with and without model-based correction

Fig. 7 4 Fig. 6 (
Fig. 7 Global error e g,NRMSE with respect to the coupling stiffness c k for SOH with correction at d k = 40, T = 3.16 • 10 −4 d k are fixed to moderate values of c k = 2 • 10 5 and d k = 5 • 10 2 , respectively.The macrostep size T is varied again in the range T ∈ 3 • 10 −3 , 1 • 10 −6 and the resulting local normalized root mean square error indicator e l,NRMSE (see Sect. 3.3) is observed in log-log scale.In addition, Table

Fig. 8
Fig. 8 Error convergence for different polynomial extrapolation orders with and without model-based correction where n d represents a "nonlinearity factor" (n d = 1 results in a linear system).The values for c k and d k are fixed to moderate values of c k = 2 • 10 5 and d k = 5 • 10 2 , respectively, and the macrostep size is set to T = 1e − 3, while the nonlinearity factor n d is varied in the range n d ∈ [1 2].The limits of the coupling damping are set to d k,min = 1 • 10 1 and d k,max = 1 • 10 4 .

Fig. 9 4
Fig. 9 Displacement s 1 of mass m 1 for different nonlinearity factors

Table 1
Largest stable macrostep size for ZOH, FOH, and SOH with and without model-based correction for several combinations of c k and d k

Table 3
Local error for a given macrostep size and necessary macrostep size for a desired local error for 0th-order to 4th-order polynomial extrapolation with and without model-based correction