Pre-step co-simulation method and device

Abstract

A pre-step co-simulation method and device for co-simulation of several subsystems is disclosed, utilizing partial derivatives of replacement subsystems. On this basis, the exact solutions of output variables of the individual subsystems are approximated by evaluation of an Error Differential Equation in order to estimate the future progress of the approximated exact output variables; corresponding input variables are determined by optimization prior to the current co-simulation time increment. The pre-step co-simulation method is realized within a computer-executable program and implemented in a network of computing nodes.

Claims

1. A method for simulating a real physical model on the basis of a pre-step co-simulation system, the method comprising: providing a first subsystem simulating at least one first parameter v of the real physical model on the basis of a first equation system, wherein the first equation system is configured for determining the first parameter v on a basis of at least one second parameter F, providing a second subsystem simulating the at least one second parameter F of the real physical model on the basis of a second equation system, wherein the second equation system is configured for determining the second parameter F on a basis of the at least one first parameter v, wherein the first input value v_in of the first parameter v in the second subsystem and the second input value F_in of the second parameter F in the first subsystem are estimated by extrapolation at a first coupling time point T, k−1, determining a first output value v_out of the first parameter v by the first subsystem and a second output value F_out of the second parameter F by the second subsystem, determining a first error between the estimated first input value v_in of the first parameter v and the determined first output value v out of the first parameter v and a second error between the estimated second input value F_in of the second parameter and the determined second output value of the second parameter F at a subsequent second coupling time point, defining a first local replacement model of the real physical model which is based on partial derivatives of the first subsystem, and defining a second local replacement model of the real physical model which is based on partial derivatives of the second subsystem, wherein on the basis of the first local replacement model and the second local replacement model an Error Differential Equation is defined, wherein the Error Differential Equation is configured for receiving the first error and the second error for determining a first characteristic of the first parameter v over time and/or a second characteristic of the second parameter F over time between the first coupling time point T, k−1 and a subsequent second coupling time point Tk, wherein on the basis of the first characteristic of the first parameter v and/or the second characteristic of the second parameter F the first input value v_in of the first parameter v and/or the second input value F_in of the second parameter F at the subsequent second coupling time point Tk is optimized for simulating the first parameter v in the first subsystem and/or the second parameter F in the second subsystem.

2. The method according to claim 1, defining a global replacement model of the real physical model which is based on partial derivatives of the first subsystem and on partial derivatives of the second subsystem, wherein the global replacement model is configured for estimating a future characteristic of the first parameter v and/or the second parameter F between the second coupling time point Tk and a third coupling time point T, k+1 such that a next first input value v_in of the first parameter v and a next second input value F_in of the second parameter F at the second coupling time point Tk can be determined.

3. The method according to claim 2, wherein the global replacement model is based on a transformation of the first subsystem and the second subsystem by time derivation of the first output value v_out of the first parameter v and the second output value F_out of the second parameter F.

4. The method according to claim 1, wherein the solution of the Error Differential Equation is approximated by an explicit or implicit numerical solver.

5. The method according to claim 1, wherein the estimation by extrapolation is a model-based extrapolation for prediction of the approximated exact output variables is utilized.

6. The method according to claim 1, wherein the estimation by extrapolation is a signal-based extrapolation for prediction of the approximated exact output variables is utilized.

7. The method according to claim 1, wherein the optimization of the input variables is based on a variational approach.

8. The method according to claim 1, wherein the partial derivatives of the first subsystem and the partial derivatives of the second subsystem are calculated instantaneously, in particular by using a subspace system identification algorithm.

9. The method according to claim 1, wherein a time increment between the first coupling time point T, k−1 and the second coupling time point Tk is adapted on the basis of the first error and the second error.

10. A simulating device for simulating a real physical model on the basis of a pre-step co-simulation system, the device comprising: a processor; a memory unit in communication with the processor, the memory unit having stored therein a first subsystem unit having a first subsystem simulating at least one first parameter v of the real physical model on the basis of a first equation system, wherein the first equation system when executed by the processor is configured for determining the first parameter v on a basis of at least one second parameter F, a second subsystem unit having a second subsystem simulating the at least one second parameter F of the real physical model on the basis of a second equation system, wherein the second equation system when executed by the processor is configured for determining the second parameter F on a basis of the at least one first parameter v, an estimating unit when executed by the processor configured for estimating the first input value v_in of the first parameter v in the second subsystem and the second input value F_in of the second parameter F in the first subsystem by extrapolation at a first coupling time point T, k−1, wherein a first output value v_out of the first parameter v is determinable by the first subsystem and a second output value F_out of the second parameter F is deteminable by the second subsystem, an error determination unit when executed by the processor configured for determining a first error between the estimated first input value v_in of the first parameter v and the determined first output value vout of the first parameter v and a second error between the estimated second input value F_in of the second parameter F and the determined second output value F_out of the second parameter F at a subsequent coupling time point Tk, an Error Differential Equation unit for defining a first local replacement model of the real physical model which is based on partial derivatives of the first subsystem and for defining a second local replacement model of the real physical model which is based on partial derivatives of the second subsystem, wherein on the basis of the first local replacement model and the second local replacement model an Error Differential Equation is definable by the Error Differential Equation unit, wherein the Error Differential Equation when executed by the processor is configured for receiving the first error and the second error for determining a first characteristic of the first parameter v over time and/or second characteristic of the second parameter F over time between the first coupling time point and a subsequent second coupling time point Tk, and an optimization unit when executed by the processor configured for optimizing the first input value v_in of the first parameter v and/or the second input value F_in of the second parameter F at the subsequent coupling time point at Tk on the basis of the first characteristic of the first parameter v and/or the second characteristic of the second parameter F for simulating the first parameter v in the first subsystem and/or the second parameter F in the second subsystem.

11. The simulating device according to claim 10, further comprising: a single or multiple computing nodes connected via a communication network to a network of computing nodes; a memory area to store the input values of the first and second parameter, the partial derivatives, output values of the first and second parameter and in particular output derivative variables, the computed partial derivatives of the subsystems and the first and second equations; the processor having access to the memory area and respective memory units of the computing nodes and configured to execute a method including, providing a first subsystem simulating at least one first parameter v of the real physical model on the basis of a first equation system, wherein the first equation system when executed by the processor is configured for determining the first parameter v on a basis of at least one second parameter F, providing a second subsystem simulating the at least one second parameter F of the real physical model on the basis of a second equation system, wherein the second equation system when executed by the processor is configured for determining the second parameter F on a basis of the at least one first parameter v, wherein the first input value v_in of the first parameter v in the second subsystem and the second input value F_in of the second parameter F in the first subsystem are estimated by extrapolation at a first coupling time point T, k−1, determining a first output value v_out of the first parameter v by the first subsystem and a second output value F_out of the second parameter F by the second subsystem, determining a first error between the estimated first input value v_in of the first parameter v and the determined first output value v out of the first parameter v and a second error between the estimated second input value F_in of the second parameter and the determined second output value of the second parameter F at a subsequent second coupling time point, defining a first local replacement model of the real physical model which is based on partial derivatives of the first subsystem, and defining a second local replacement model of the real physical model which is based on partial derivatives of the second subsystem, wherein on the basis of the first local replacement model and the second local replacement model an Error Differential Equation when executed by the processor is defined, wherein the Error Differential Equation when executed by the processor is configured for receiving the first error and the second error for determining a first characteristic of the first parameter v over time and/or a second characteristic of the second parameter F over time between the first coupling time point T, k−1 and a subsequent second coupling time point Tk, wherein on the basis of the first characteristic of the first parameter v and/or the second characteristic of the second parameter F the first input value v_in of the first parameter v and/or the second input value F_in of the second parameter F at the subsequent second coupling time point Tk is optimized for simulating the first parameter v in the first subsystem and/or the second parameter F in the second subsystem.

12. The simulating device according to claim 11, wherein the computing nodes represent real-time systems and the communication of data is preformed asynchronously.

13. A non-transitory computer-readable medium on which there is stored a computer program for simulating a real physical model, the computer program, when being executed by a data processor, is adapted for performing and/or controlling a method including, providing a first subsystem simulating at least one first parameter v of the real physical model on the basis of a first equation system, wherein the first equation system is configured for determining the first parameter v on a basis of at least one second parameter F, providing a second subsystem simulating the at least one second parameter F of the real physical model on the basis of a second equation system, wherein the second equation system is configured for determining the second parameter F on a basis of the at least one first parameter v, wherein the first input value v_in of the first parameter v in the second subsystem and the second input value F_in of the second parameter F in the first subsystem are estimated by extrapolation at a first coupling time point T, k−1, determining a first output value v_out of the first parameter v by the first subsystem and a second output value F_out of the second parameter F by the second subsystem, determining a first error between the estimated first input value v_in of the first parameter v and the determined first output value v out of the first parameter v and a second error between the estimated second input value F_in of the second parameter and the determined second output value of the second parameter F at a subsequent second coupling time point, defining a first local replacement model of the real physical model which is based on partial derivatives of the first subsystem, and defining a second local replacement model of the real physical model which is based on partial derivatives of the second subsystem, wherein on the basis of the first local replacement model and the second local replacement model an Error Differential Equation is defined, wherein the Error Differential Equation is configured for receiving the first error and the second error for determining a first characteristic of the first parameter v over time and/or a second characteristic of the second parameter F over time between the first coupling time point T, k−1 and a subsequent second coupling time point Tk, wherein on the basis of the first characteristic of the first parameter v and/or the second characteristic of the second parameter F the first input value v_in of the first parameter v and/or the second input value F_in of the second parameter F at the subsequent second coupling time point Tk is optimized for simulating the first parameter v in the first subsystem and/or the second parameter F in the second subsystem.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows a flowchart illustrating a complete cycle for co-simulation with respect to the herein described computer-implemented pre-step co-simulation method.

(2) FIG. 2 shows a diagram indicating the different steps to be performed for realizing the pre-step co-simulation method.

(3) FIG. 3 shows a schematic block diagram pointing out the herein described architecture of the co-simulation device as well as central and distributed implementation possibilities of the pre-step co-simulation method.

(4) FIG. 4 shows a simulation result indicating the performance improvement by application of the herein described computer-implemented co-simulation method.

DETAILED DESCRIPTION OF ILLUSTRATED EMBODIMENTS

(5) The aspects defined above and further aspects of the present invention are apparent from the examples of embodiment to be described hereinafter and are explained with reference to the examples of embodiment. The invention will be described in more detail hereinafter with reference to examples of embodiment but to which the invention is not limited.

(6) The illustrations in the drawings are schematically presented. It is noted that in different figures similar or identical elements are provided with the same reference signs.

(7) A pre-step co-simulation method of 100, 200 (see FIG. 1 and FIG. 2) and device 300 (see FIG. 3) is shown, that can be used for simulating a complex physical system and real physical models, which comprises an arbitrary number of subsystems 308. The subsystems 308 can represent models of real objects simulated in dedicated simulation tools or even real hardware components mounted on tailored testbeds in case of real-time applications. Due to the modular character of the co-simulation for synchronization purposes data is communicated via a network 302 (see FIG. 3) and is exclusively exchanged between subsystems at dedicated communication points in time (coupling time points). Subsystems 308 are executed independently between these communication point instances.

(8) FIG. 1 shows a flowchart illustrating a complete cycle for co-simulation with respect to the herein described computer-implemented pre-step co-simulation method. The pre-step co-simulation method further comprises the computation in step 104 (compute partial derivatives of (replacement) subsystems for all individual subsystems) of replacement system (Eq. 10-11) matrices and the assembly in step 105 (assemble required (global) matrices based on linear replacement subsystems) of derived matrices (Eq. 15-16) in step 101 (Compute partial derivatives of (replacement) subsystems for all individual subsystems), approximating the exact output variables/values of the respective parameters of the individual subsystems by evaluation of an Error Differential Equation in 106 (globally valve Error Differential Equation and approximate exact output), optionally adapting the (local) co-simulation time increment in step 107 (locally determined adequate step size), predicting the future progress of the exact output variables in step 108 ((globally) predict future progress of exact output) and optimizing corresponding input values of the first and second parameters in step 109 ((locally) optimise input variables based on predicted future progress) prior to the current co-simulation time increment in step 111 (local simulation of subsystems).

(9) FIG. 2 shows a diagram indicating the different steps to be performed for realizing the pre-step co-simulation method for simulating a real physical model according to the present invention. FIG. 2 shows graphs of simulations at certain coupling time points Tk of one subsystem 308. A respective first subsystem 308 is provided simulating at least one first parameter (e.g. the speed v of a vehicle) of the real physical model on the basis of a first equation system (e.g. a behavioral model) (local), wherein the first equation system is configured for determining the first parameter (v) on a basis of at least one second parameter (e.g. the torque of a shaft of the vehicle F).

(10) Furthermore, a second subsystem 308 simulates the at least one second parameter (F) of the real physical model on the basis of a second equation system (local), wherein the second equation system is configured for determining the second parameter (F) on a basis of the at least one first parameter (v). The first input value (v_in) of the first parameter (v) in the second subsystem and the second input value (F_in) of the second parameter (F) in the first subsystem are estimated by extrapolation (curve 201) at a first (e.g. an initial) coupling time point (T, k−1). The method further comprises determining a first output value (v_out) of the first parameter (v) by the first subsystem and a second output value (F_out) of the second parameter (F) by the second subsystem (curve 202).

(11) The method further comprises determining a first error between the estimated first input parameter (v_in) (curve 201) and the determined first output parameter (v_out) (curve 202) and a second error between the estimated second input parameter (F_in) (curve 201) and the determined second output parameter (F out) (curve 202) at a subsequent coupling time point (e.g. the Errors at T,k).

(12) A first local replacement model of the real physical model which is based on partial derivatives of the first subsystem is defined and a second local replacement model of the real physical model which is based on partial derivatives of the second subsystem is defined.

(13) On the basis of the first local replacement model and the second local replacement model an Error Differential Equation is defined, wherein the Error Differential Equation is configured for receiving the first error and the second error for determining a first characteristic of the first parameter over time and/or second characteristic of the second parameter over time between the first coupling time point and a subsequent second coupling time point (at coupling time point Tk, i.e. the curve 203 between Tk−1 and Tk). On the basis of the first characteristic of the first parameter and/or the second characteristic of the second parameter the first input value (v_in) of the first parameter (v) and/or the second input value (F_in) of the second parameter (F) at the subsequent coupling time point is optimized (at Tk) for simulating the first parameter (v) in the first subsystem and/or the second parameter (F) in the second subsystem.

(14) Thereby, optimization is based on the extrapolation of the first and second characteristic 204 at the second coupling time point till the third coupling time point by utilizing the global replacement model.

(15) FIG. 3 shows a schematic block diagram pointing out the herein described architecture of the co-simulation device as well as central and distributed implementation possibilities of the pre-step co-simulation method.

(16) The pre-step co-simulation method is implemented within a central or distributed co-simulation computing architecture 300, consisting of computing nodes 301 and a communication network 302. Each computing node 301 consist of a memory 305 for storing data, a network interface device 306 for data communication and a processor 303 for processing instructions of programs, where a single node can also host and execute several programs.

(17) The pre-step co-simulation method is implemented within a co-simulation program 304, 307, which coordinates all subsystem simulation processes 308 by starting, pausing, continuing, stopping and communicating data at communication point instances for synchronization purposes locally or to other relevant computing nodes 301. In particular, the co-simulation program 304, 307 is able to read from memory and write data to memory 305, such as relevant input and output variables and their time derivatives and optionally the corresponding partial derivatives of all subsystems. When executed by a processor 303 the processor conduct to compute (see step 104 in FIG. 1) the replacement system (Eq. 10-11) matrices and to assemble (see step 105 in FIG. 1) derived matrices (Eq. 15-16) in (step 101 in FIG. 1), to evaluate the Error Differential Equation (see step 106 in FIG. 1) for adjusting the co-simulation time increment (see step 107 in FIG. 1) and for approximation of the exact solution of output variables of the subsystems (see step 108 in FIG. 1) and to determine input variables by optimization (see step 109 in FIG. 1). Furthermore, when executed by a processor 303 the processor conduct to communicate the resulting inputs to corresponding subsystem simulation computing nodes prior to the current co-simulation time increment and to coordinate, i.e. start, continue, pause, stop, the individual subsystem simulation processes remotely.

(18) During non-iterative co-simulation the coupled subsystems 308 are independently solved by the use of tailored fixed or variable step-size solvers. In the context of co-simulation δT denotes the numerical steps performed by the underlying solver within the individual subsystems 308. Further, for synchronization purposes an additional step-size is introduced. Also this so-called co-simulation time increment ΔT may be time varying and determines the points in time t.sub.k+1.sup.i=t.sub.k.sup.i+ΔT.sub.k.sup.i at which coupling data is exchanged between subsystem simulations, where k denotes the k-th co-simulation time increment and i denotes the i-th subsystem. In the most general case, instead of updating only the value at the coupling time instant, all available values in the time interval of the last co-simulation time increment are exchanged. This functionality is supported by high performance co-simulation platforms and leads to improved synchronisation capabilities.

(19) Because of bidirectional dependencies between subsystems 308 and the modular (distributed) character of a co-simulation, several input signals (input values of first and second parameter) are unknown for the time interval of the current co-simulation time increment and thus, piecewise extrapolation of some coupling signals is necessary to solve the overall coupled system. Extrapolation is directly associated with an estimation error and thus, an inherent co-simulation discretization error is introduced. The invention described herein copes with this problem by disclosing a pre-step co-simulation method and the corresponding device.

(20) FIG. 1 shows a flowchart 100 which depicts steps of the herein described pre-step co-simulation method for a single co-simulation time increment. In the following the individual steps of the flowchart 100 are detailed.

(21) It is common practice to describe the behavior of physical real world objects via ordinary differential equations (ODEs). In general, in some cases additional algebraic constraints are stated, which directly leads to a subsystem 308 description via differential algebraic equations (DAEs) as shown in Equations (1-3); the explicit dependency on time is omitted:
{dot over (x)}.sup.i=f.sup.i(x.sup.i,u.sup.i),  (1)
y.sup.i=g.sup.i(x.sup.i,u.sup.i),  (2)
0=h.sup.i(x.sup.i,u.sup.i),i=1 . . . r.  (3)

(22) Where x.sup.i represents the state variables and parameters of the i-th subsystem 308; u.sup.i and y.sup.i represent the corresponding input 201 and output 202 (see FIG. 2) variables of the subsystem, respectively; and r denotes the number of subsystems. These subsystems are modelled within special simulation tools and are typically solved by tailored numerical solvers 308. For disclosing the invention a linear version of the DAE system is outlined below:
{dot over (x)}.sup.i=f.sup.i(x.sup.i,u.sup.i)≈A.sup.i.Math.x.sup.i+B.sup.i.Math.u.sup.i,  (4)
y.sup.i=g.sup.i(x.sup.i,u.sup.i)≈C.sup.i.Math.x.sup.i+D.sup.i.Math.u.sup.i,  (5)
0=h.sup.i(x.sup.i,u.sup.i)≈E.sup.i.Math.x.sup.i+F.sup.i.Math.u.sup.i.  (6)

(23) Depending on the overall system simulation at hand dedicated inputs and outputs of the individual subsystems are interconnected. The interconnections are described via a coupling matrix, where the coupling matrix is identified based on the adjacency matrix of the directed graph (nodes represent subsystems and edges represent interconnections), which is generated from the overall co-simulation system. For a system comprising two fully coupled subsystems, i.e. u.sup.1=y.sup.2 and u.sup.2=y.sup.1, the subsystems interconnections are described by the coupling matrix L:

(24) $\begin{array}{cc}\left(\begin{array}{c}{y}^1\\ y^2\end{array}\right)=\underset{\underset{L}{︸}}{\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)}.Math.\left(\begin{array}{c}{u}^1\\ u^2\end{array}\right).& \left(7\right)\end{array}$

(25) In contrast to known approaches, the herein disclosed invention is based on replacement subsystems, whose can be derived from the classical DAE descriptions (Eq. 1-3) by time derivation of the output variables of all subsystems till the input variables appears on the right hand side and derivation of the algebraic constraint equations as shown in Equations (8) and (9):
{dot over (y)}.sup.i=S.sup.i(y.sup.i,u.sup.i,{dot over (u)}.sup.i),  (8)
0=H.sup.i(y.sup.i,u.sup.i,{dot over (u)}.sup.i).  (9)

(26) For sake of simplicity the approach w.r.t. the disclosed invention is described herein based on a single time derivation of the subsystems outputs; for relative degree smaller than one; generalization to multiple time derivations and higher relative degrees is straight forward. This way the original system (Eq. 1-3) is transformed into a new state space form (Eq. 8-9). Below, two options are described to determine the parameters of the replacement subsystems, both are based on approximated linear system descriptions:

(27) $\begin{array}{cc}{\stackrel{.}{y}}^i=S^i\left(y^i,u^i,{\stackrel{.}{u}}^i\right)\approx \frac{\partial S^i}{\partial y}.Math.y^i+\frac{\partial S^i}{\partial u}.Math.u^i+\frac{\partial S^i}{\partial \stackrel{.}{u}}.Math.{\stackrel{.}{u}}^i,& \left(10\right)\\ 0=H^i\left(y^i,u^i,{\stackrel{.}{u}}^i\right)\approx \frac{\partial H^i}{\partial y}.Math.y^i+\frac{\partial H^i}{\partial u}.Math.u^i+\frac{\partial H^i}{\partial \stackrel{.}{u}}.Math.{\stackrel{.}{u}}^i.& \left(11\right)\end{array}$

(28) Option 1: By utilization of an system identification algorithm (see below), where the inputs to the algorithm represent the input u.sup.i and output y.sup.i variables of the subsystems and their time derivatives {dot over (u)}.sup.i and {dot over (y)}.sup.i, respectively. In this case no explicit transformation of the matrices of the original subsystems has to be performed.

(29) Option 2: By realizing the system transformation (Eq. 12-14). Therefore, the partial derivatives of the state derivatives w.r.t. the state and input variables and the partial derivatives of the outputs w.r.t. the state and input variables and of all original subsystems (Eq. 1-6) has to be available. These matrices can be delivered by the subsystems itself via an enhanced co-simulation interface (e.g. the Functional Mockup Interface by Blochwitz) or identified online by appropriate system identification approaches. The replacement subsystem dynamic and input matrices are determined as follows:

(30) $\begin{array}{cc}\frac{\partial S^i}{\partial y}=C^i.Math.A^i.Math.C^{i^{-1}},& \left(12\right)\\ \frac{\partial S^i}{\partial u}=C^i.Math.B^i-C^i.Math.A^i.Math.C^{i^{-1}}.Math.D^i,& \left(13\right)\\ \frac{\partial S^i}{\partial \stackrel{.}{u}}=D^i,\mathrm{with}{A}^i=\frac{\partial f^i}{\partial x};B^i=\frac{\partial f^i}{\partial u};C^i=\frac{\partial g^i}{\partial x}\mathrm{and}{D}^i=\frac{\partial g^i}{\partial u}.& \left(14\right)\end{array}$

(31) Notably, matrices C.sup.i has to be regular to ensure an inversion. However, in other cases an inversion by computation of the well-known pseudo-inverse serves as work-around.

(32) Above, two times system identification were mentioned. For offline (non-real-time) and online (real-time) identification different approaches exist and can be applied by one skilled in the art. However, herein the well-known subspace method is mentioned for completeness. Based on the input and output variables of a subsystem the subspace method is able to identify the state-space representation, as stated in Equ. (10-11). This way, there is no need to specify the order of the replacement subsystems as the method automatically choose the best model order based on the singular-value decomposition, which depends on the input and output variables.

(33) In order to ensure a stable co-simulation of especially stiff systems the cross-coupling effects of the interconnected subsystems have to be taken into account; therefore the individual replacement system (Eq. 10-11) and derived matrices (Eq. 15-16) of the individual subsystems 308 are assembled to global system matrices (see step 101 in FIG. 1) at coupling time instances.

(34) Depending on the system at hand a reduced number of updates is possible, e.g. every third co-simulation increment or with each significant change of the system dynamics. For sake of simplicity, the following details are explained for a co-simulation consisting of two (r=2) subsystems 308; the generalization for a more complex system consisting of more subsystems is straight forward.

(35) The following (global) matrices will be used by the Error Differential Equation (see step 106 in FIG. 1) and the result thereof is utilized to adjust the co-simulation time increment (see step 107 in FIG. 1) and to approximate the exact output variables 108 of the subsystems:

(36) $\begin{array}{cc}\tilde{A}=\left(\begin{array}{cc}\frac{\partial S^1}{\partial y}& \frac{\partial S^1}{\partial u}\\ \frac{\partial S^2}{\partial u}& \frac{\partial S^2}{\partial y}\end{array}\right),\tilde{B}=\left(\begin{array}{cc}0& \frac{\partial S^1}{\partial \stackrel{.}{u}}\\ \frac{\partial S^2}{\partial \stackrel{.}{u}}& 0\end{array}\right),\tilde{C}=\left(\begin{array}{cc}\frac{\partial S^1}{\partial u}& 0\\ 0& \frac{\partial S^2}{\partial u}\end{array}\right),\tilde{D}=\left(\begin{array}{cc}\frac{\partial S^1}{\partial \stackrel{.}{u}}& 0\\ 0& \frac{\partial S^2}{\partial \stackrel{.}{u}}\end{array}\right).& \left(15\right)\end{array}$

(37) To generalize (global) matrices the couple matrix L essentially has to be taken into account, which is omitted here and can be derived by one skilled in the art.

(38) On the other hand, for prediction of the future progress (see step 108 in FIG. 1) of the output variables (global) model-based extrapolation is performed based on the following (global) matrices:

(39) $\begin{array}{cc}\hat{A}=\left(\begin{array}{cc}\frac{\partial S^1}{\partial y}& 0\\ 0& \frac{\partial S^2}{\partial y}\end{array}\right),\hat{B}=\left(\begin{array}{cc}\frac{\partial S^1}{\partial u}& 0\\ 0& \frac{\partial S^2}{\partial u}\end{array}\right),\hat{C}=\left(\begin{array}{cc}\frac{\partial S^1}{\partial \stackrel{.}{u}}& 0\\ 0& \frac{\partial S^2}{\partial \stackrel{.}{u}}\end{array}\right).& \left(16\right)\end{array}$

(40) The generalization works without knowledge of the coupling matrix, due to the fact that the matrices remain blockdiagonal matrices in any case.

(41) In FIG. 1 the flowchart 100 depicts the main steps of the invented co-simulation algorithm. FIG. 2 shows a diagram 200 which indicates exemplarily the steps for one single co-simulation time increment. The computed co-simulation result 202 until the coupling time instant t.sub.k are communicated to the master algorithm. At this point in time (t.sub.k) the determined and applied input u.sup.i 201 at t.sub.k−1 over the interval [t.sub.k−1, t.sub.k), the co-simulation result 202, i.e. the output y.sup.i at t.sub.k over the interval (t.sub.k−1, t.sub.k], as well as the individual partial derivatives for assembling the (global) matrices (Eq. 15-16) are available.

(42) In a first step the deviation between the applied past input u.sup.i 201 and the resulting output y.sup.i 202 are used to determine the co-simulation discretization error ∈.sup.i over the past co-simulation time increment:
∈.sup.i=y.sup.i−u.sup.i.  (17)

(43) On this basis the exact output {tilde over (y)}.sup.i 203 over the past step is approximated by evaluation of the (global) Error Differential Equation 106:
{dot over (δ)}=(I−{tilde over (B)}).sup.−1.Math.(Ã.Math.δ+{tilde over (C)}.Math.ε+{tilde over (D)}.Math.{dot over (∈)}),  (18)

(44) with δ=(δ.sup.1,δ.sup.2).sup.T and ∈=(∈.sup.1,∈.sup.2).sup.T, leading to an additional error measure δ.sup.i:
δ.sup.i={tilde over (y)}.sup.i−y.sup.i,  (19)

(45) where Ã, {tilde over (B)}, {tilde over (C)} and {tilde over (D)} are assembled matrices (Eq. 15) for the overall system. The approximation of the exact output {tilde over (y)}.sup.i 203 in step 106 is then determined by:
{tilde over (y)}.sup.i=y.sup.i+δ.sup.i.  (20)

(46) Since Equation (18) represents an ordinary differential equation there is a variety of ways to solve it; focus is especially taken in linearly-implict methods. Due to the fact that y.sup.i typically represents a physical signal it has to be continuous and therefore it follows that δ.sup.L is continuous as well. Combining this property with the global initial value δ(t.sub.start)=0 with Equation (18) leads to an initial value problem. Utilization of the linearly-implicit Euler scheme, i.e. a numerical solver, leads to:

(47) $\begin{array}{cc}{\delta }_k={\delta }_{k-1}+\Delta T.Math.{\left(I-\frac{\Delta T}{2}{\left(I-\tilde{B}\right)}^{-1}.Math.\tilde{A}\right)}^{-1}.Math.{\left(I-\tilde{B}\right)}^{-1}.Math.\hspace{1em}\left[\tilde{A}.Math.{\delta }_{k-1}+0.5.Math.\tilde{C}.Math.\left({ϵ}_k+{ϵ}_{k-1}\right)+0.5.Math.\tilde{D}.Math.\left({\stackrel{.}{ϵ}}_k+{\stackrel{.}{ϵ}}_{k-1}\right)\right]& \left(21\right)\end{array}$

(48) Note: The Error Differential Equation approach can be utilized to determine the approximation of the exact output of the subsystems in a standalone manner without co-simulation discretization error.

(49) Due to errors δ.sup.i and ∈.sup.i a locally determination of an appropriate co-simulation time increment ΔT.sub.k.sup.i is possible 107. Due to fact that δ.sup.i and ∈.sup.i are independent from the other subsystems it is feasible to choose different constant co-simulation time increments for every subsystem. The co-simulation time increment can be updated at each communication point t.sub.k.sup.i for every subsystem, which leads to a local and adaptive co-simulation time increment control 107. Three possible options are outlined below:

(50) Option 1: Like in most well-known step size controls an error indicator:

(51) $\begin{array}{cc}{\mathrm{err}}_{\mathrm{indicator}1}^i=\frac{{\int }_{t_{k-1}^i}^{t_k^i}.Math.\delta \left(t\right).Math.\mathrm{dt}}{{\int }_{t_{k-2}^i}^{t_{k-1}^i}.Math.\delta \left(t\right).Math.\mathrm{dt}},& \left(22\right)\end{array}$
is utilized with the aim of decreasing the co-simulation time increment if err.sub.indicator1.sup.i is greater than a given threshold and increasing the co-simulation time increment if err.sub.indicator1.sup.i is less than a given threshold.

(52) Option 2: If numerical steps, based on δT, are available the strategy is to evaluate

(53) $\tilde{y}\left(t_{k-1}^i+\frac{\Delta T_k^i}{2}\right)={\tilde{y}}_{k-1,1/2}^i$
and then compute:

(54) $\begin{array}{cc}{\mathrm{err}}_{\mathrm{indicator}2}^i=\frac{\begin{array}{c}{\int }_{t_{k-1}^i}^{t_{k-1}^i+\frac{\Delta T_k^i}{2}}.Math.{\tilde{y}}_{k-1}^i+\frac{2.Math.\left({\tilde{y}}_{k-1,1/2}^i-{\tilde{y}}_{k-1}^i\right)}{\Delta T_k^i}.Math.t-{\tilde{y}}^i\left(t\right).Math.\mathrm{dt}+\\ {\int }_{t_{k-1}^i+\frac{\Delta T_k^i}{2}}^{t_k^i}.Math.{\tilde{y}}_{k-1,1/2}^i+\frac{2.Math.\left({\tilde{y}}_k^i-{\tilde{y}}_{k-1,1/2}^i\right)}{\Delta T_k^i}.Math.t-{\tilde{y}}^i\left(t\right).Math.\mathrm{dt}\end{array}}{{\int }_{t_{k-1}^i}^{t_k^i}.Math.\delta \left(t\right).Math.\mathrm{dt}},& \left(23\right)\end{array}$
as the required error indicator. The strategy to adapt the co-simulation time increment is the same as in Option 1 above, exclusively the error indicator changes to err.sub.indicator2.sup.i.

(55) Option 3: a model-based step-size control utilizes the slope of the extrapolated output variables and computes the error indicator

(56) 0$\begin{array}{cc}{\mathrm{err}}_{\mathrm{indicator}3}^i=\frac{.Math.s_{k-1}^i-s_k^i.Math.}{.Math.s_{k-1}^i.Math.+.Math.s_k^i.Math.},& \left(24\right)\end{array}$
where s.sub.k.sup.i denotes the slope of the extrapolation of the local output variable in t.sub.k.sup.i for the i-th subsystem. The strategy to determine the co-simulation time increment is the same as in Option 1 above, exclusively the error indicator changes to err.sub.indicator3.sup.i.

(57) Note: All three options above are based on local variables and therefore the co-simulation time increment can be determined locally, independently for every subsystem. For practical issues an additive combination of all options may be implemented.

(58) To ensure the implict nature of this pre-step co-simulation algorithm the prediction of the future progress of the approximated exact output {tilde over (y)}.sup.i (graph 203 in FIG. 2) over the next co-simulation time increment (see step 108 in FIG. 1), denoted as ŷ.sup.i (graph 204 in FIG. 2), is essential. Depending on the available information two options for prediction of the future progress are possible:

(59) Option 1: A (global) model-based extrapolation is performed, which is based on the (global) assembled dynamic matrix  and input matrices {circumflex over (B)} and Ĉ (Eq. 16)
{circumflex over ({dot over (y)})}=(I−Ĉ.Math.L).sup.−1.Math.(Â+{circumflex over (B)}.Math.L).Math.ŷ  (25)

(60) with ŷ=(ŷ.sup.1, ŷ.sup.2).sup.T. The initial condition ŷ.sub.k.sup.i for the solution of the initial value problem (Eq. 25) originates from the solution of the Error Differential Equation, it holds:
ŷ.sub.k.sup.i={tilde over (y)}.sub.k.sup.i.  (26)

(61) Option 2: A (evtl. local) signal-based extrapolation (graph 204 in FIG. 2) is performed based on the approximation of the exact output (graph 203 in FIG. 2) and the data available until the coupling time instant. Local means, that every subsystem can perform it's own extrapolation independently. There is wide variety of well-known extrapolation techniques like polynomial extrapolation using time derivatives of the signal. There are a lot of well-known methods for computing derivatives; herein the proposal is to use the recently proposed DELC-Method, see Genser-2017.

(62) Notably, Option 1 is preferred as the cross-coupling effects of the subsystems are considered, which are essential for handling of stiff system integration.

(63) After prediction of the future progress 204 based on the approximated output 203, the input 201 of the next co-simulation time increment has to be determined by global or local optimization 109. Local means, that this part of the pre-step co-simulation method can be computed locally 307 at the different computing nodes 301 and there is no need for globally assembled matrices in this part of the algorithm. The optimization is done by minimizing the following objective function over the next co-simulation time increment:
min.sub.u.sub.i.sub.,t∈[t.sub.k.sub.,t.sub.k+1.sub.]|ŷ.sup.i(t)−y.sup.i(t)|.  (27)

(64) As this optimization 109 and the resulting input variable 201 is also computed at the coupling time instant prior to the next co-simulation increment, this approach is called pre-step co-simulation method. Note, in the particular case of handling DAE subsystems the algebraic subsystem equations like (Equ. 11) have to be considered as an additional constraint within the input optimization. The five following options are possible.

(65) Option 1: The input 201 for the next co-simulation time increment can be determined by energy preservation means, where the integral of the prediction of the future progress 204 and the integral of the input 201 over the next co-simulation time increment equals:

(66) $\begin{array}{cc}{\int }_{t_k^i}^{t_{k+1}^i}{\hat{y}}^i\left(t\right)\mathrm{dt}={\int }_{t_k^i}^{t_{k+1}^i}{y}^i\left(t\right)\mathrm{dt}.& \left(28\right)\end{array}$

(67) Option 2: The input 201 to the next co-simulation time increment can be determined by explicitly solving the minimization problem (Equ. 27). The input is determined by discretization of the minimization problem (Equ. 27), e.g. by the co-simulation time increment ΔT, leading to the explicit algebraic equation:

(68) $\begin{array}{cc}{u}_k^i={\varphi }_{\frac{\partial S^i}{\partial u}}^{-1}.Math.\left({\hat{y}}_{k+1}^i-{\varphi }_{\frac{\partial S^i}{\partial y}}.Math.y_k^i\right).& \left(29\right)\end{array}$
where

(69) ${\varphi }_{\frac{\partial S^i}{\partial y}}\mathrm{and}{\varphi }_{\frac{\partial S^i}{\partial u}}$
denotes the well-known transitions matrices corresponding to the dynamic matrix

(70) $\frac{\partial S^i}{\partial y}$
and the input matrix

(71) $\frac{\partial S^i}{\partial u}$
of the local replacement subsystems. This equation is based on a constant (zero-order-hold) input variable 201, i.e. a constant input value over the next co-simulation time increment as illustrated in FIG. 2, but this can be generalized for other forms/kinds of the input 201 extrapolation, e.g. first-order-hold or an arbritary sequence of values, in general (herein omitted for sake of simplicity).

(72) Option 3: The input 201 for the next co-simulation time increment can be determined by solving the minimization problem (Equ. 27) as a Variational Approach in general.

(73) Option 4: In another aspect, the input 201 can be determined in an iterative manner, without solving the subsystems 308, by introducing an additional part 113 into the flowchart 100. After the determination of the input 109 and before the local simulation of the subsystems 103, the Error Differential Equation is iteratively solved for several times by using a redefined discretisation error:
∈.sup.i=ŷ.sup.i−u.sup.i,  (30)
for the next co-simulation time increment. On this basis new input variables 201 are iteratively determined at the coupling time instant prior (!) to the co-simulation increment.

(74) Option 5: The input (graph 201 in FIG. 2) for the next co-simulation time increment can be determined (truly) iteratively by utilizing the computed output variables 202 of the co-simulated subsystem over the next co-simulation time increment. By this approach the individual subsystems are (truly) iteratively solved for several times. In contrast to the other non-iterative co-simulation options mentioned above, this would lead to an iterative co-simulation approach. Iterations are performed until a predefined accuracy threshold is reached.

(75) The sections above describe the pre-step co-simulation algorithm. Depending on the available computing architecture two possibilities for implementation exist: a centralized-master 304 approach and a distributed-master 304, 307 approach, as shown in FIG. 3.

(76) Option 1: The centralized-master 304 approach is characterized by the implementation of the whole pre-step co-simulation algorithm 100 in a single computer-executable program running on a dedicated computing node 301. This approach can be considered as a single master approach, where all the necessary data, i.e. input and output variables of the subsystems and their time derivatives, respectively, and optionally the partial derivatives of the individual (replacement) subsystems, if provided by the subsystems itself, is communicated via the communication network 302 to the single master algorithm program 304. In this case, the centralized-master 304, optionally identify the replacement subsystems matrices (Eq. 10-11) of the individual subsystems 308 based input 201 and output 202 variables of the subsystems and their time derivatives, respectively, assembles the required matrices (Eq. 15-16, 101), evaluates the Error Differential Equation 106, determines the next co-simulation time increment 107, performs the prediction of the future progress of the approximated exact solution 108 as well as optimizes the input variables for the next co-simulation time increment 109 at the current coupling time instant.

(77) For realization of the pre-step co-simulation a dedicated order of execution of the centralized-master program 304, the computation of the co-simulation time increments of the individual subsystems 111 and the communication of data 112 has to be considered. This is ensured by the centralized-master program 304, which controls the individual subsystems 308, as described in the following. In a first instance, at each coupling time instant Tk, the master algorithm gathers the initial input 201 and output 202 variables of the parameters of the subsystems 308 and their time derivatives, respectively, as well as optionally partial derivatives from the individual (replacement) subsystems (Eq. 10-11) if determined at the individual computing nodes of the subsystems, optionally computes 104 required partial derivatives for the replacement subsystems (Eq. 10-11) of the individual subsystems 308 based on input 201 and output 202 variables and their time derivatives and assembles 105 the required matrices (Eq. 15-16, 101), evaluates the Error Differential Equation 106, determines the next co-simulation time increment 107, performs the prediction of the future progress of the approximated exact solution 108 as well as optimizes the input variables 109 for the next co-simulation time increment. In a second step, the determined inputs 201 and partial derivatives are stored to memory 305 sent to the corresponding subsystems via the communication network 110, which directly trigger the execution of the subsystem co-simulation step 111. At the end of the current subsystem time increment the subsystem simulation is paused and the resulting output variables 202 are sent back to the centralized-master node 301 via the communication network and are stored 110 by the centralized-master algorithm 304, which directly triggers the execution of the next co-simulation step. This procedure repeats for computation of subsequent co-simulation time increments until the end of the overall co-simulation is reached, where the subsystem simulations are forced to stop by the centralized-master 304.

(78) Option 2: The distributed-master 304, 307 approach is characterized by the separation of the pre-step co-simulation algorithm into several parts which are running on different computing nodes in order to improve the overall simulation performance in terms of speed. The first part is represented by a computer-executable program 304 running on a dedicated computing node, where all the necessary data, input 201 and output 202 variables of the subsystems 308 and their time derivatives, respectively, as well as optionally partial derivatives from the individual (replacement) subsystems (Eq. 10-11) if determined at the individual computing nodes of the subsystems, is communicated via the network 302 to this part of the master algorithm 304. This part of the algorithm optionally computes 104 required partial derivatives (Eq. 10-11) of the individual subsystems 308 based on input 201, input derivatives and output 202 variables and assembles 105 the required matrices (Eq. 15-16, 101), evaluates the Error Differential Equation 106, optionally determines the local next co-simulation time increments 107 and optionally performs a global prediction (model-based) of the future progress of the approximated exact solution 108 for the next co-simulation time increment at the current coupling time instant. Other parts are represented by a computer-executable program 307 running on computing nodes 301 where the subsystems 308 are executed too, where each of these parts locally and optionally computes 104 required partial derivatives (Eq. 10-11) of the individual subsystems 308 based on input 201, input derivatives and output 202 variables of the subsystems and their time derivatives, respectively, optionally determines the local next co-simulation time increments 107, optionally performs a global prediction (model-based or signal-based) of the future progress of the approximated exact solution 108 for the next co-simulation time increment and perform the optimization 109 of the input 201 variables for the next co-simulation time increment at the current coupling time instant.

(79) For realization of the pre-step co-simulation 100 a dedicated order of execution of the distributed-master program 304, 307, the computation of the co-simulation time increments 111 of the individual subsystems 308 and the communication of data via the communication network 302 has to be considered 103. This is ensured by the distributed-master program 304, 307) which controls the individual subsystems 308, as described in the following.

(80) In a first instance, at each coupling time instant Tk, the master algorithm 304 gathers the initial input 201 and output 202 variables of the subsystems 308 and their time derivatives, respectively, as well as optionally partial derivatives from the individual (replacement) subsystems (Eq. 10-11) if determined at the individual computing nodes of the subsystems, optionally computes 104 required partial derivatives for the replacement subsystems (Eq. 10-11) of the individual subsystems 308 based on input 201 and out variables and their time derivatives and assembles 105 the required matrices (Eq. 15-16, 101), evaluates the Error Differential Equation 106, optionally determines the local next co-simulation time increments 107 and, optionally, performs the prediction (model-based) of the future progress of the approximated exact solution 108. In contrast to the centralized-master approach, in a second step, the determined exact solution 204 or the prediction (model-based) future progress of the approximated exact solution 108 is optionally sent to the corresponding local masters 307 at subsystems computation nodes 301, where optionally the future progress of the approximated exact solution is locally predicted (model-based or signal-based), optionally the next co-simulation time increment 107 is determined and the input variables 201 for the next co-simulation time increment are optimized 109. After local optimization 109 the resulting input variables 201 are stored and communicated 110 to the other part of the master 304 via the communication network 302 and execution of the subsystem co-simulation step is triggered 111. At the end of the current subsystem time increment the subsystem simulation is paused, optionally the required subsystem partial derivatives are locally computed and stored into memory 305 and the resulting output variables 202 are sent back to the distributed master program 304, which directly triggers the execution of the next co-simulation step. This procedure repeats for computation of subsequent co-simulation time increments until the end of the overall co-simulation is reached, where the subsystem simulations are forced to stop by the distributed master.

(81) The following sections describe the related co-simulation device 300 for execution of the pre-step co-simulation method. FIG. 3 shows the architecture of the co-simulation device, which consists of different computing nodes 301 and a communication network 302 in general.

(82) A computing node 301 consists of at least one central processing unit, the CPU 303, a memory 305 and a network interface card, the NIC 306. The computing node can be instructed to carry out arbitrary sequences of arithmetic or logical operations automatically, which are defined within computer-executable programs. The memory 305 of the computing node is used, on the one hand, for storing information for immediate use in order to perform the arithmetic or logical operations (stored in a random-access memory, RAM storage) and used, on the other hand, for storing the programs and related data (stored in a hard-disc-drive, HDD storage). The network interface card 306 is a hardware component of a computing node that connects a computing node to a network of computing nodes 302, which consists of a network interface controller for data communication via local area (LAN) or a wide area networks (WAN) by utilizing a dedicated communication protocol, such as TCP/IP, UDP, CAN or for wireless communication Bluetooth or LTE, based on GSM/EDGE and UMTS/HSPA technologies, is used. For real-time applications typically EtherCAT or TT-Ethernet is implemented.

(83) For implementation of the pre-step co-simulation method the outlined algorithm 100, 200 is realized within computer-executable programs 304, 307 for the centralized- or distributed master approach. Subsystem simulation models are solved by tailored solvers within dedicated simulation tools (subsystems) 308. In case of a dedicated co-simulation scenario the different computer-executable algorithms 304, 307, 308 are distributed over an existing network of computing nodes 301, whereas also multiple computer-executable programs can run on a single computing node.

(84) In an exemplary application of the pre-step co-simulation method, the master program 304 is executed by the central processing unit 303 of the corresponding computing node 301 and, in a first step after initiation of the master program 304, the processor 303 conducts to compute the (replacement) subsystem matrices for all subsystems 101. The optionally stored partial derivatives from the original subsystems are read from the local HDD memory or are optionally determined by utilization of a system identification scheme based on the available input 201 and output 202 variables of the subsystems and their time derivatives, respectively, e.g. by application of the subspace identification method, to compute 104 the required partial derivatives (Eq. 10-11) of the individual (replacement) subsystems 308 and assembles 105 the required matrices (Eq. 15-16, 101). The processor conducts afterwards to write the resulting matrices into the local RAM memory of the computing node. These matrices are utilized by the subsequent instructions of the master program for evaluation of the Error Differenation Equation for approximating the exact output variables 106, subsequently for determining the length of the next co-simulation time increment 107, for the prediction of the future progress of the exact output variable 108 of the next co-simulation time increment, and for the optimization of the input variables 109 for the next co-simulation time increment. For these calculations the instructions of the master programme 304 causes the processor 303 to perform the related arithmetic and logical operations as well as the intermediate reading and writing of variables to the local RAM or HDD memory 305.

(85) In case of a distributed co-simulation scenario, implemented in a distributed computing network 300, where the different computer-executable programs (master 304, 307 and several subsystem simulation programs 308) are allocated to different computing nodes 301 within the network of computing nodes, the processor 303 of the master computing node 301 reads the corresponding values from the memory and forces the network interface controller of the network interface card 306 to transmit the input variables 201 and optionally their time derivatives to the individual subsystem simulation computation nodes 301, respectively, by utilizing the implemented communication network 302 and the corresponding communication protocol. As soon as the transmitted data is received at the subsystem computing nodes 301 the corresponding network interface controllers on the network interface cards 306 are indicating the arrival of new input variables to the subsystem 308 simulation, which stores these values to the local memory 305 and triggers the execution of the co-simulation time increment 111. The resulting output variables 202 of the individual subsystem simulation 308 at the end of the co-simulation time increment are stored into the local memory 305. The subsystem simulation program 308 executed by the local processor 303 conducts the transmission of the output variables by the local network interface controller via the network interface card 306 to the master computing node 301. As soon as the transmitted data is received by the master computing nodes the corresponding network interface controller on the network interface cards 306 are indicating the arrival of new output variables 201 to the master program 304, which writes these values to the local memory 305 and triggers the execution of the next co-simulation time increment.

(86) Another important case is represented by real-time applications. The herein disclosed pre-step co-simulation method is primarily developed for non-iterative co-simulation and thus, the method possesses determined runtime behavior and fulfills the main requirement for real-time application. On the one hand, the method can be applied for realization of a time-triggered co-simulation on a single computing node and, on the other hand, the method can be implemented in a distributed real-time co-simulation scenario, where a plurality of real-time computation nodes 301 is connected via an appropriate (real-time capable) communication network 302, e.g. EtherCAT, TT-Ethernet, UDP, CAN or LTE 5G. In this case, asynchronous communication between the individual subsystems is realized, i.e. the individual subsystems are time-triggered by the counter (clock) of the corresponding computation node. As the real-time behavior of the different computer-executable programs is ensured by the utilized computing hardware and the scheduler of the program hosting operating systems, the required input and output variables, and optionally the partial derivatives of the individual subsystems, are available in time for correct execution of the pre-step co-simulation method.

(87) Hence, a co-simulation computing device according to the present invention comprises a single or multiple computing nodes 301 connected via a communication network 302 to a network of computing nodes, a memory area 305 to store the input, input derivative, output and output derivative variables, the computed partial derivatives of the individual (replacement) subsystems (Eq. 10-11) and the assembled required matrices (Eq. 15-16) and processors having access to the memories 305 of the computing nodes 301 and configured to execute the inventive method.

(88) The described details represent the basis for the flowchart 100 depicted in FIG. 1. The flowchart is roughly separated into three main blocks: the computation of the required matrices 101, the determination of the input variables for the next co-simulation time increment 102 and the execution of the next co-simulation time increment 103 by utilizing the co-simulation device. Super block 101 comprises the computation 104 of the required partial derivatives (Eq. 10-11) of the individual (replacement) subsystems 308 based on input 201 and output 202 variables of the subsystems and their time derivatives, respectively, and assembles 105 the required matrices (Eq. 15-16, 101).

(89) Super block 102 comprises the evaluation of the global Error Differential Equation 106, using the computed global assembled matrices 101 for approximating the exact outputs variables for all subsystems, the, optional locally, determining of the adequate co-simulation time increment 107, the, optional locally, prediction of the future progress of the output variables 108. These predicted output variables 108 ensure the implicit character of the pre-step method and serve as fundamental basis for the, optionally local, optimization of the input variables 109 for the next co-simulation time increment.

(90) Super block 103 comprises the storage and communication of the determined input variables by the co-simulation device 110, which triggers the individual subsystems for execution of the next co-simulation time increment 111. After computation of the co-simulation time increment the resulting output variables are communicated and stored by the co-simulation device 112, which triggers the computation of the subsequent co-simulation time increment by repeating the process, starting with 101.

(91) Finally, an example is examined and the results are shown in FIG. 4. Three different simulations were performed: a monolithic simulation serving as reference for comparison, a co-simulation by applying the simplest possible but most often used coupling approach (called zero order hold, ZOH) and a pre-step co-simulation indicating the benefit and performance improvement of the herein outlined invention in terms of accuracy and stability. An overall stiff system was chosen consisting of two subsystems 308. For sake of simplicity a single output variable is illustrated. The vertical lines represent the communication points in time (coupling time points Tk) used with respect to a constant co-simulation time increment. FIG. 4 points out that the ZOH approach diverges due to the ZOH introduces co-simulation discretization error and the result from the pre-step co-simulation method is almost superimposable with the monolithic solution, indicating that co-simulation discretization errors are almost mitigated by the disclosed method.

(92) For demonstration purposes the pre-step co-simulation method is exemplarily applied to a distributed vehicle simulation and test environment, comprising a real vehicle mounted on a classical powertrain testbed and virtual wheels embedded within a virtual environment. High stiffness exists as the dynamics of the wheel slip significantly differs to those of the vehicle by several magnitudes. The mechanical link between the real vehicle and the virtual wheels is represented by the four half-shafts, where each one is connected to a torque-controlled e-drive, i.e. a so called Dynometer. The set torque to each e-drive controller is determined by the corresponding virtual wheel within the virtual environment (incl. modelled friction, road slope, etc.) based on the measured rotational speed (e.g. first parameter v) of the corresponding half-shaft of the real vehicle. Each half-shaft torque (e.g. second parameter F) and rotational speed (angular velocity) represent coupling signals, i.e. a power bond. The input signal into the wheel simulation is a rotational speed signal and the output signal is the resulting torque. As the overall simulation is performed in real-time, a parallel co-simulation is necessary and thus, at every coupling time instant the rotational speed signal and the torque signal has to be extrapolated. Due to extrapolation the herein addressed co-simulation discretization error is introduced and the energy-transfer over the half-shaft is distorted, which has a negative impact on the final simulation and test performance. In particular, a higher extrapolated wheel rotational speed will lead to a higher counter torque response over this co-simulation time increment. Hence, the feedback onto the source system (i.e. the vehicle on the powertrain testbed) is higher it can thus speed up more than it should have.

(93) Because rotational speed is a signal over the interface, the wheel will also speed up too much, compare to FIG. 4. This significantly affects test results, e.g. vehicle energy consumption, on short (over the current co-simulation time increment Tk) and long term (over the whole test run). The classical/typical procedure for mitigation is the utilization of very small sampling times (i.e. co-simulation time increment) leading to a massive over-sampling (10 kHz) and increased costs in terms of high-performance software and hardware components (e.g. EtherCAT real-time bus system).

(94) By application of the herein disclosed pre-step co-simulation method the co-simulation discretization error is compensated prior to the calculation of the co-simulation time increment, eliminating the problem described above. In this case the real vehicle mounted on the power-train testbed represents a computing node and the wheels simulated within the virtual environment are assigned to another computing node; the computing nodes are connected via a communication network as illustrated in FIG. 3. In particular, the speed signals as well as the torque signals at the individual half-shafts are modified according to FIG. 1 and FIG. 2 prior to the calculation of the wheel torque responses for the actual co-simulation time increment Tk, leading to a stabilizing effect as indicated in FIG. 4. This way the wheel accelerates as expected and the distortion is reduced to a minimum which enables the application of larger co-simulation time increments and, as a consequence, the usage of cheaper software and hardware components.

(95) It should be noted that the term “comprising” does not exclude other elements or steps and the article “a” or “an” does not exclude a plurality. Also elements described in association with different embodiments may be combined. It should also be noted that reference signs in the claims should not be construed as limiting the scope of the claims.