Ben Abir

Laurent Khaled

Mohamed El Duval

Ben Mongi

Daniel Gaïd

Simon

Abir Ben Khaled
email: abir.ben-khaled@ifpen.fr

Laurent Duval
email: laurent.duval@ifpen.fr

Mongi Ben Gaid
email: mongi.ben-gaid@ifpen.fr

Daniel Simon
email: daniel.simon@inria.fr

Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI

Keywords: FMI, parallel simulation, signal processing, polynomial extrapolation, real-time, contextbased decision

published or not. The documents may come

Introduction

During the design process of complex systems, such as in automotive, simulation is proven to be an indisputable step between concept design and prototype validation. Realistic simulations allow for the preliminary evaluation, tuning and possibly redesign of proposed solutions ahead of implementation, thus lowering the risks. To be confident in the result, building such simulations requires high-fidelity models both for the components and for their interaction.

Currently, building high-fidelity system-level models of cyber-physical systems in general and automo-tive cars in particular, is a challenging duty. One problem is the diversity of models, designed for different environments, provided by various multi-disciplinary teams. Distinctive environments are preferred for a specific use due to specific strengths (modeling language, libraries, solvers, cost. . . ). The FMI specification has been proposed to improve this issue [START_REF] Blochwitz | The functional mockup interface for tool independent exchange of simulation models[END_REF].

However, the simulation of high-fidelity models is time consuming, and reaching real-time constraints is out of the capabilities of single-threaded simulations running on single cores. Simulation speed-ups are needed, in particular by splitting the systems into submodels to be executed in parallel on currently available multi-core chips.

Unfortunately most of the existing simulation software are currently unable to exploit multi-core plateforms, as they rely on sequential Ordinary Differential Equations (ODE) and Differencial Algebraic Equations (DAE) solvers. The co-simulation approaches can provide significant improvements by allowing to simulate together models coming from different areas, and to validate both the individual behaviors and their interaction [START_REF] Valasek | Simulation Techniques for Applied Dynamics, chapter Modeling, simulation and control of mechatronical systems[END_REF]. The simulators may be exported from original authoring tools as Functional Mock-up Units (FMUs), and then imported in a co-simulation environment. Hence, they cooperate at run-time, thanks to the FMI definitions of their interfaces, and to the master algorithms of these environments.

The "FMI for model exchange" framework allows for solving independently the sub-models using custom solvers. In this context, several methods have been already proposed to perform real-time distributed simulation of complex physical models. For example, in [START_REF] Faure | Methods for real-time simulation of cyberphysical systems: application to automotive domain[END_REF], the study focused on the case of fixed-step solvers. Then, in [START_REF] Ben Khaled | Multicore simulation of powertrains using weakly synchronized model partitioning[END_REF], the study was extended to handle the case of variable time-step solvers.

However, accounting for the dependencies between the sub-models needs to synchronize them at some time intervals. Certainly, this synchronization avoids the propagation of numerical errors in the simulation results and guarantees their correctness. Unfortunately, these synchronization constraints also lead to waiting periods and idle time of some processors. Consequently it decreases the potential efficiency of the threaded parallelism existing in multi-core platforms.

To overcome this limitation, and to more efficiently exploit the available parallelism, the dependencies constraints should be relaxed as far as possible while keeping accumulated errors under control. In a first step, this can be performed by a well done system decomposition that minimizes the dependencies between the sub-models. For example, a method was proposed in [START_REF] Sjölund | Towards efficient distributed simulation in Modelica using transmission line modeling[END_REF] for distributed simulation using a technology based on bilateral delay lines called transmission line modeling (TLM), where the decoupling point is chosen when the variables change slowly and the timestep of the solver is relatively small. Unfortunately, most often perfect decoupling cannot be reached and data dependencies still exist between parallel blocks. Some synchronization between them must be kept tight through small communication steps between models, which prevents the variable time-step solvers to reach large integration steps.

It is proposed that the synchronization steps can be stretched out with limited deterioration of the simulation precision, thanks to a well-suited, albeit simple, context-based polynomial extrapolation of the exchanged data beyond the synchronization points between sub-models. This paper is organized as follows. First, a formal model of a hybrid dynamical system is given and a model of the integration errors due to slack synchronization is sketched in Section 2. The background on prediction and polynomial prediction algorithms are developed in Section 3. The principles for contextbased extrapolation, to cope with the hybrid nature of the models, are exposed in Section 4. The methodology is assessed in Section 5 using the model of an internal combustion engine.

Motivation for extrapolation 2.1 Model formalization

Consider a hybrid dynamical system Σ described by a set of nonlinear differential equations:

Ẋ = f (t, X, D, U ) for t n ≤ t < t n+1 , Y = g(t, X, D, U ), where X ∈ R n X is the continuous state vector, D ∈ R n D is the discrete state vector, U ∈ R n U is the input vector, Y ∈ R n Y is the output vector and t ∈ R + is the time.
The sequence (t n ) n≥0 of strictly increasing time instants represents discontinuity points called "state events", which are the roots of the equation

h(t, X, D, U ) = 0.
The function h is usually called zero-crossing function or event indicator. It is used for event detection and location [START_REF] Zhang | Zerocrossing location and detection algorithms for hybrid system simulation[END_REF].

At each time instant t n , a new continuous state vector can be computed as a result of the event handler

X(t n ) = I(t n , X, D, U ),
and a new discrete state vector can be computed as a result of discrete state update

D(t n ) = J (t n-1 , X, D, U ).
If no discontinuity affects a component of X(t n ), the right limit of this component will be equal to its value at t n .

It is assumed that Σ is well posed in the sense that a unique solution exists for each admissible initial conditions X(t 0 ) and D(t 0 ) and that consequently X, D, U , and Y are piece-wise continuous functions in

[t n ,t n+1 ].
To execute the system in parallel, the model must be split into several sub-models. Assume for simplicity, that the system is decomposed into two subsystems as in Figure 1. Our approach generalizes to any decomposition into N blocks of system Σ. Therefore, the system can be written as:

Ẋ1 = f 1 (X 1 , X 2 , D 1 , U 1 ) Y 1 = g 1 (X 1 , X 2 , D 1 , U 1 ) and Ẋ2 = f 2 (X 1 , X 2 , D 2 , U 2 ) Y 2 = g 2 (X 1 , X 2 , D 2 , U 2 ) with X = [X 1 X 2 ] T and D = [D 1 D 2 ] T ,
Y 1 ∪ Y 2 = Y and Y 1 ∩ Y 2 = /
0. To perform the numerical integration of the whole multivariable system, each of these simulators needs to exchange, at communication points, the data needed by the others (see Figure 2). To speed up the integration, the parallel branches must be as independent as possible, so that they are synchronized at a rate P by far slower than their internal integration step h n (P ≫ h n ). Therefore, between communication points, each simulator integrates at its own rate (assuming a variable step solver), and considers that the data incoming from others simulators is hold as constant. It is likely that large communication intervals allow to speed up the numerical integration, but may result in integration errors and poor confidence in the final result. Modeling the errors induced by slack synchronization is a first step to find effective directions to improve the trade-offs between integration speed and accuracy.

Integration errors and parallelism

To compute the next state value X i (t n+1 ), i = 1, 2 (see Figure 3), the numerical solver needs at least the values of X i (t n ) and Ẋi (t n ) = f i (X(t n )) (e.g. for Euler integration). The inputs and discrete states are omitted for clarity. When computing Ẋ1 (t n ) = f 1 (X(t n )), the value of the local variable X 1 (t n ) is always available. This is not the case for X 2 (t n ), which is computed in a parallel branch. In fact, X 2 is only available in branch 1 at synchronization with interval P, which is larger than the integration step h n . In other words, X 2 (t n ) is available only when the time t n corresponds with a synchronization point t s (see Figure 2), otherwise its estimated value is the one transmitted at the previous synchronization point. Let us evaluate the evolution of integration errors due to slack synchronization between the parallel branches when computing

Ẋ1 (t n ) = f 1 (X(t n ))
. The analysis on Σ 1 remains valid for Σ 2 .

The influence of using a delayed value of X 2 in f 1 (.) (respectively X 1 in f 2 (.)) is due to the lack of updated data during a delay τ, represented by the difference between the current integration time t n and the last synchronization time t s as

τ = t n -t s ( 1
)
with t s = P t n P therefore

t s = lP when t n = l.P l ∈ N * (l -1)P when t n < l.P l ∈ N * leading to τ = 0 when t n = t s τ > 0 when t n > t s
Therefore, the induced error at t n+1 in the subsystem Σ 1 , denoted E 1 (t n+1 ), is the difference between X 1 (t n+1 ) for the unsplit model (2) and X1 (t n+1 ) for the split model (3):

Session 2A: FMI 2 DOI 10.3384/ECP14096225
Proceedings of the 10 th International ModelicaConference March 10-12, 2014, Lund, Sweden

X 1 (t k+1 ) = X 1 (t k ) + h k f 1 (X 1 (t k ), X 2 (t k )), k ∈ {0, . . . , n} (2)
X1 (t k+1 ) = X 1 (t k+1 ) k = 0 X1 (t k ) + h k f 1 ( X1 (t k ), X2 (t k -τ)) k ≥ 1 (3)
In other words,

E 1 (t n+1 ) = n ∑ k=0 E 1 (t k ) + h n [f 1 (X 1 (t n ), X 2 (t n )) -f 1 ( X1 (t n ), X2 (t n -τ))] = E 1,p (t n ) + E 1,c (t n+1 )
where

E 1,c (t n+1 ) = h n [f 1 (X 1 (t n ), X 2 (t n )) -f 1 ( X1 (t n ), X2 (t n -τ))] E 1,p (t n ) = n ∑ k=0 E 1 (t k ) (4)
Here E 1,c (t n+1 ) is the current error generated at t n+1 whatever a synchronization or not. So, the global decoupling error E 1 (t n+1 ) is the result of the accumulation of past errors E 1,p (t n ) and the current error E 1,c (t n+1 ). As a conclusion, to achieve a correct result, two conditions must be met for the current (local) error and the global error:

• |E 1,c (t n+1 )| < ǫ loc : allowed local error • |E 1 (t n+1 )| < ǫ glo : allowed global error
These conditions can be satisfied by acting on some parameters. Indeed, in (4), the delay error depends on the integration steps h n and on the delay τ. The integration step h n is already adapted following the numerical solver strategy and the user-defined solver tolerance. The delay τ, however, depends on the last synchronization time t s , which is function of the synchronization period P.

The delay induced error tends to zero when the delay τ tends to zero, which means that the delay error can be eliminated with the synchronization interval set equal to the integration steps. In other words, all the parallel subsystems should be integrated at the same adaptive rate (in the case of adaptive synchronization period), or with same fixed time-step. These two assumptions are very restrictive, as they force to choose a global adequate time-step regardless the discontinuities and the stiffness of the sub-systems. Compared with the single-threaded simulation, the only possible speed-ups during a parallel execution would be brought by the brute force computation power of the multicore machine, reduced by the parallelization cost.

Therefore, considering a split model and a parallel execution, a trade-off must be found between acceptable simulation errors, thanks to tight enough synchronization, and simulation speed-ups thanks to decoupling between sub-models.

To add a degree of freedom to this trade-off achievement, we propose to extrapolate model inputs to compensate the stretching out of the communication steps between sub-models. Note that in this first approach of the polynomial extrapolation, the synchronization interval is considered as constant. Future enhancements will consider communication step size control, for which the error analysis and estimation can be inspired by [START_REF] Arnold | Error analysis and error estimates for co-simulation in FMI for model exchange and co-simulation V2.0[END_REF]. Extrapolation is sensitive for different reasons:

• prediction should be efficient: causal, sufficiently fast and reliable;

• there exist no universal prediction scheme, efficient with every signal;

• polynomial prediction may fail in stiff cases [START_REF] Arnold | Simulation Techniques for Applied Dynamics[END_REF] (cf. Section 4 for details).

We choose to base our extrapolation on polynomial prediction, which allows fast and causal calculations. The rationale is that, in this situation, the computing cost of a low-order polynomial predictor would be by far smaller than the extra model computations needed by shorten communication steps. Since such predictions would be accurate neither for any signal (for instance, blocky versus smooth signals) nor any signal behavior (slow variations versus steep onsets), we borrow a context-based approach , common with lossless image coders [START_REF] Weinberger | The LOCO-I lossless image compression algorithm: principles and standardization into JPEG-LS[END_REF], such as GIF or PNG formats. The general aim of these image coders is to predict a pixel value based on a pattern of causal neighboring pixels. Compression is obtained when the prediction residues possess smaller intensity values, and more generally a better distribution (concentrated around close-to-zero values) than the pixels in the original image. They may thus be coded on smaller "bytes", using entropy coding techniques. In images, one distinguishes basic "objects" such as smooth-intensity varying regions, or edges with different orientations. Based on simple calculation of the prediction pattern pixels, different contexts are inferred (e.g. flat, smooth, +45 o or -45 o edges, etc.). Look-up table predictors are then used, depending on the context.

In the proposed approach, we build a heuristic table of contexts (in Section 4) based on a short frame of past samples, and affect a pre-determined polynomial predictor to obtain a context-dependent extrapolated value. We now review the principles of extrapolation.

Causal polynomial prediction 3.1 Background on prediction

This section is dedicated to a peculiar instance of discrete time series, or signal, forecasting. The neighboring topics of prediction or extrapolation represent a large body of knowledge in signal processing [START_REF] Wiener | Extrapolation, interpolation, and smoothing of stationary time series[END_REF], econometrics [START_REF] Brown | Smoothing, Forecasting and Prediction of Discrete Time Series[END_REF] or control [START_REF] Box | Time Series Analysis: Forecasting and Control. Probability and Statistics[END_REF].

In the present case, we consider a real-valued, regularly sampled signal u, with period P, known at synchronization or communication intervals. Prediction in general assumes the knowledge of signal formation models. Since very little is assumed on the signal's dynamics (no behavioral/explicit model is available, periodicity and regularity are unknown), and as we operate under real-time conditions, implying strong causality, only a tiny fraction of time series methods are practically applicable. Zeroth-order hold or nearestneighbor extrapolation is probably the most natural, the less hypothetical, and the less computationally expensive forecasting method. It consists in using the latest known sample as the predicted value. It possesses small (cumulative) errors when the time series is relatively flat or its sampling rate is sufficiently high, with respect to the signal's dynamics. In other words, it is efficient when the time series is sampled fast enough to ensure small variations between two consecutive sampling times. However, it indirectly leads to undersampling related disturbances, that affect the signal content. They appear as quantization-like noise, offset or peak flattening.

In our co-simulation framework, communication intervals are not chosen arbitrarily small for computational efficiency. Thus, the slow variation of inputs and outputs cannot be ensured in practice. Hence, borrowing additional samples from the past known data and using higher-order extrapolation methods could be beneficial, provided a trade-off of cost and error is met. Different forecast methods of various fidelity and complexity may be efficiently evaluated. We focus here on polynomial methods, for their simplicity and ease of implementation, following initial works in [13, Chapter 16].

Notations

We denote by P (δ ,λ ) the least-squares polynomial predictor of degree δ ∈ N and prediction length λ ∈ N * .

The prediction length λ represents the number of past samples required for each prediction, performed in the least-squares sense [14, p. 227 sq.]. For convenience, we use a 0-last-sample-index convention: we re-index the frame of the λ past samples such that the last known sample is indexed by 0. Computations for the prediction at relative time τ (loosely denoted by u(τ)), defined in (1), thus require past samples {u 1-λ , u 2-λ , . . . , u 0 }. We first recall principles and formulas for a standard least-squares, degree-two or parabolic prediction. The general equations are derived next.

Polynomial prediction of degree δ = 2

We look for the best fitting parabola, i.e. with degree δ = 2, u(t) = a δ + a δ -1 t + a δ -2 t 2 to approximate the set of discrete samples {u 1-λ , u 2-λ , . . . , u 0 }. The prediction polynomial P (2,λ ) is defined by the vector of polynomial coefficients a = [a 2 , a 1 , a 0 ] T . They are determined, in the least-squares sense [START_REF] Stigler | Gauss and the invention of least squares[END_REF], by minimizing the squared or quadratic, error:

e(a) = 0 ∑ l=1-λ u l -(a 2 + a 1 l + a 0 l 2 ) 2 .
Note that the l indices here are non-positive, between 1-λ and 0. The minimum error is obtained by solving the following system of equations (zeroing the derivatives with respect to each of the free variables a i ):

∀i ∈ {0, 1, 2}, ∂ e(a) ∂ a i = 0 namely:                    0 ∑ l=1-λ l 0 u l -(a 2 l 0 + a 1 l 1 + a 0 l 2 ) = 0, 0 ∑ l=1-λ l 1 u l -(a 2 l 0 + a 1 i l + a 0 l 2 ) = 0, 0 ∑ l=1-λ l 2 u l -(a 2 l 0 + a 1 l 1 + a 0 l 2 ) = 0.
(5)

The system in (5) may be rewritten as: Let m d = ∑ λ -1 l=0 l δ -d u -l (here the indices l are positive) denote the (δ -d)-th moment of the frame u i , and m the vector of moments [m 2 , -m 1 , m 0 ] T . We express the sums of integer powers by Σ d λ = ∑ λ -1 i=0 i d . Closedform expressions exist for Σ d λ , involving Bernoulli sequences [START_REF] De Bruyn | Formulas for 1 + 2 p + 3 p + . . . + n p[END_REF]. For instance:

                   0 ∑ 1-λ u l = a 2 0 ∑ 1-λ l 0 + a 1 0 ∑ 1-λ l 1 + a 0 0 ∑ 1-λ l 2 , 0 ∑ 1-λ lu l = a 2 0 ∑ 1-λ l 1 + a 1 0 ∑ 1-λ l 2 + a 0 0 ∑ 1-λ l 3 , 0 ∑ 1-λ l 2 u l = a 2 0 ∑ 1-λ l 2 + a 1 0 ∑ 1-λ l 3 + a 0 0 ∑ 1-λ l 4 .
• Σ 0 λ = λ , • Σ 1 λ = (λ -1)λ /2, • Σ 2 λ = (λ -1)λ (2λ -1)/6, • Σ 3 λ = (λ -1) 2 λ 2 /4, • Σ 4 λ = (λ -1)λ (2λ -1)(3λ 2 -3λ -1)/30.
We now form the matrix Z (2,λ ) of sums of powers (depending on δ = 2 and λ ):

Z (2,λ ) =   Σ 0 λ -Σ 1 λ Σ 2 λ -Σ 1 λ Σ 2 λ -Σ 3 λ Σ 2 λ -Σ 3 λ Σ 4 λ   .
The system in (5) rewrites:

  m 2 -m 1 m 0   =   Σ 0 λ -Σ 1 λ Σ 2 λ -Σ 1 λ Σ 2 λ -Σ 3 λ Σ 2 λ -Σ 3 λ Σ 4 λ   ×   a 2 a 1 a 0   or m = Z (2,λ )
× a. Now we want to find the value predicted by P (2,λ ) at time τ. Let τ 2 = [1, τ, τ 2 ] T be a vector of τ powers. Then u(τ) is equal to a 2 + a 1 τ + a 0 τ 2 = τ T 2 × a. Finally, Z (2,λ ) is always invertible, provided that λ > δ . Its inverse is denoted Z (-2,λ ) . It thus does not need to be updated in real-time. It may be computed off-line, numerically or even symbolically. Hence:

u(τ) = τ T 2 × Z (-2,λ ) × m .
The vector τ 2 and Z (-2,λ ) are fixed, and the product τ T 2 × Z (-2,λ ) may be stored at once. Thus, for each prediction, the only computations are the update of the vector m and his product with the aforementioned stored matrix. It thus enables look-up-table-based predictions, which helps to reduce propagation errors in matrix computations.

General formulas

Inferring from the previous example, we easily get a more generic extrapolation pattern in its matrix form:

1 τ • • • τ δ u(τ) = ×       Σ 0 λ -Σ 1 λ • • • (-1) δ Σ δ λ -Σ 1 λ . . . . . . . . . . . . . . . . . . . . . (-1) δ Σ δ λ • • • • • • Σ 2δ λ       -1 ×      m δ -m δ -1
. . .

(-1) δ m 0      . Note τ δ = [1, τ, • • • , τ δ ] T , then: u(τ) = τ T δ Z (-δ ,λ ) m .
As in the previous case, only m and one matrix product need be computed in real-time. When δ = 0, one easily sees that:

u(τ) = m 0 Σ λ 0 = u 1-λ + • • • + u 0 λ ,
that is, the running average of past frame values, reducing to the zeroth-order hold when λ = 1. Although the matrix formulation is convenient, actual computation does not require true matrix calculus, especially for small degrees δ . For instance, P (1,3) yields the simple estimator form:

u(τ) = τ 2 (u 0 -u -2 ) + 1 6 (5u 0 + 2u -1 -u -2 ).

Context-based extrapolation

Actual complex systems usually present nonlinearities and discontinuities, so that it is hard to predict their future behavior from past observations. Moreover the considered models are generated using the FMI for model exchange framework, which does not provide the inputs' derivatives (conversely with the FMI for co-simulation architecture). Hence the previously described polynomial prediction cannot correctly extrapolate along all the system trajectories.

For example, [START_REF] Arnold | Simulation Techniques for Applied Dynamics, chapter Numerical methods for simulation in applied dynamics[END_REF] studies a method based on a sequential implementation of continuous dynamical systems that uses a constant, linear or quadratic extrapolation and a linear interpolation to improve the accuracy of the modular time integration. The study shows that the method is successful for non-stiff systems and it fails for the stiff case. The context-based extrapolation is then performed to account for steps, stiffness, discontinuities or weird behavior, and use adapted extrapolation to limit excessively wrong prediction.

Keeping with the previous 0-last-sample-index convention, and for the sake of simplicity, we first define a measure of variation based on the last three samples:

d 0 = u 0 -u -1 and d 1 = u -1 -u -2
, the last and previous differences. Their absolute values are compared with two thresholds, γ 0 and γ -1 , respectively. We then define three complementary conditions:

• O if |d i | = 0; • C i if 0 < |d i | ≤ γ i ; • C i if |d i | > γ i ;
We can now define the six-context Table 1, and examples for their associated heuristic polynomial predictors. The six contexts form a partition, i.e. they are mutually exclusive, and cover all possible options for a hybrid dynamical system. They are illustrated in Figure 4. Their names represent their behavior. For instance, the flat context addresses steady signals, for which a mere zeroth-order hold suffices, hence P (0,1) . The calm context represents a sufficiently sampled situation, where value increments over time remain below fixed thresholds. In this case, the signal is relatively regular, and could be approximated by a quadratic polynomial, for instance P [START_REF] Valasek | Simulation Techniques for Applied Dynamics, chapter Modeling, simulation and control of mechatronical systems[END_REF][START_REF] Sjölund | Towards efficient distributed simulation in Modelica using transmission line modeling[END_REF] . For the "flat" and "jump" contexts, there is an additional procedure which consists in resetting the extrapolation to prevent inaccurate prediction. For example, when context 1 is chosen just after context 5, the quadratic extrapolation P (2,5) requires 5 valid samples, whereas the last 3 only are relevant.

Our two-threshold is relatively simple. Hence, the choice of the thresholds γ 0 and γ -1 , is potentially crucial. For instance, fixed values may reveal inefficient under important amplitude or scale variation of signal. Hence, we have chosen here to compute them, in a running manner, on the past frame {u 1-ω , . . . , u -3 }. With excessively low thresholds, high-order extrapolations would be rarely chosen, loosing the benefits of predictions. Too high thresholds would in contrast suffer from any unexpected jump or noise. As the contexts are based on backward derivatives, we have used in the simulations presented here the mid-range statistical estimator of their absolute values. This amounts to set: 1.

γ 0 = γ -1 = 1 2 max i∈[1-ω,...,-3] (|u i -u i+1 |). Table 1: Summary of the six-context Table. n(ame) # |d -2 | |d -1 | d -2 .d -1 (δ , λ ) f(lat) 0 O O O (0, 1) c(alm) 1 C 1 C 2 any (2, 5) m(ove) 2 C 1 C 2 any (0, 1) r(est) 3 C 1 C 2 any (0, 2) t(ake) 4 C 1 C 2 > 0 (1, 3) j(ump) 5 C 1 C 2 < 0 (0,
5 Case study

Engine simulator

In this study, a Spark Ignition (SI) RENAULT F4RT engine has been modeled. It is a four-cylinder in-line Port Fuel Injector (PFI) engine in which the engine displacement is 2000 cm 3 . The air path (AP) consists in a turbocharger with a mono-scroll turbine controlled by a waste-gate, an intake throttle and a downstreamcompressor heat exchanger. The engine model was developed using the Mod-Engine library [START_REF] Benjelloun-Touimi | From physical modeling to real-time simulation: Feedback on the use of modelica in the engine control development toolchain[END_REF]. ModEngine is a Modelica [START_REF] Fritzson | Principles of object-oriented modeling and simulation with Modelica 2[END_REF] library that allows the modeling of a complete engine with diesel and gasoline combustion models. Requirements for the ModEngine library were derived from the existing IFP-Engine AMESim1 library. Mod-Engine contains more than 250 sub-models. It has been developed to allow the simulation of a complete virtual engine using a time-scale related to fractions of the crankshaft angle. A variety of elements are available to build representative models for engine components, such as turbocharger, wastegate, gasoline or Diesel injectors, valve, air path, EGR loop etc. Mod-Engine is currently functional in the Dymola tool 2 .

The engine model and the split parts were imported into xMOD model integration and virtual experimen-tation tool [START_REF] Ben Gaïd | Heterogeneous model integration and virtual experimentation using xMOD : Application to hybrid powertrain design and validation[END_REF], using the FMI export features of Dymola. The engine model has 118 state variables and 312 event indicators (of discontinuities).

Decomposition approach

The partitioning of the engine model is performed by separating the four-cylinder from the air path (AP), then by isolating the cylinders (C i , for i ∈ [1, 2, 3, 4]) from each other. This kind of splitting allows for the reduction of the number of events acting on each subsystem. In fact, the combustion phase raises most of the events, which are located in the firing cylinder. The solver can process them locally during the combustion cycle of the isolated cylinder, and then enlarge its integration time-step until the next cycle.

From a thermodynamic point of view, the cylinders are loosely coupled, but a mutual data exchange does still exist between them and the air path. The model is split into 5 components and governed by a basic controller denoted CTRL. It gathers 91 inputs and 98 outputs.

Tests and results

Tests are performed on a platform with 16 GB RAM and 2 "Intel Xeon" processors, each running 8 cores at 3.1 GHz.

Reference simulations

The model validation is based on the observation of some quantities of interest as the intake and exhaust manifold pressures, air-fuel equivalence ratio and torque. These outputs are computed using LSO-DAR which is a variable time-step solver with a rootfinding capability to detect the events occurring during the simulation. It has also the ability to adapt the integration method depending on the observed system stiffness.

The simulation reference Y ref is built from the integration of the entire engine model, the solver tolerance (tol) being decreased until reaching stable results, which is reached for tol = 10 -7 (at the cost of an unacceptable slow simulation speed).

Then, to explore the trade-offs between the simulation speed and precision, simulations are run with increasing values of the solver tolerance until reaching a desired relative integration error Er, defined by ( 6)

Er(%) = 100 N . N-1 ∑ i=0 Y ref (i) -Y (i) Y ref (i) (6)
with N the number of saved points during 1 s of simulation. Iterative runs showed that the relative error converge to a desired error (Er ≤ 1 %) for tol = 10 -4 . The single thread simulation of the whole engine with LSODAR and tol = 10 -4 provides the simulation execution time reference, to which the parallel versions are compared. When using the split model, each of its 5 components is assigned to a dedicated core and integrated by LSODAR with tol = 10 -4 .

Effect of the context-based extrapolation on accuracy

To explore the effect of extrapolation on accuracy, the communication step has been set to 250 µs in a first set of experiments. This value has been chosen to provide acceptable results for the accuracy (Er ≈ 1 %), while being large enough to make extrapolation useful.

The tests show that performing only a fixed polynomial prediction (conventional first and second order extrapolation) on the engine model fails, with integration errors larger than for the reference simulation. This is due to the hybrid nature of the model, for which the extrapolation failures are caused by discontinuities, and also by sharp variations of some variables at specific instants. These cases totally waste the gain in precision due to successful extrapolation in the other parts of the state trajectories.

In contrast, using the context-based polynomial predictor, the outputs of the simulation are almost always closer to the reference trajectory than those computed when considering the inputs hold as constant (Figure 5). Figure 6 shows that using context-based extrapola-

Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI tion, the prediction step is discarded when there is a discontinuous behavior in the signal, and that the degree of the predictor is adapted according to the signal slope (Figure 6). The cumulative relative integration error on a long simulation run is computed in Table 2. It shows that the context-based extrapolation efficiently decreases this error for the chosen variables, for example by 63 % for the temperature and by 72.5 % for the fuel density.

Effect of the context-based extrapolation on simulation time

The ultimate objective of extrapolation is to decrease the simulation by stretching out the synchronization interval, while keeping the relative integration error Er inside predefined bounds. Indeed, widening the communication step from 100 µs to 250 µs without extrapolation (Figure 7) saves time but increases the error (e.g. 6.97 % for the burned gas density and 340.5 % for the fuel density). Using the extrapolation for the 250 µs step fortunately decreases the relative error to values close to, or below, those measured for the 100 µs step with frozen inputs.

Pressure

Temperature Air density Fuel density Burned gas density 0 0.5 Table 3 shows the simulation speed-up compared with the single-threaded reference. First note that when splitting the model into 5 threads integrated in parallel on 5 cores, the speed-up is supra-linear w.r.t. the number of cores. Indeed, the containment of events detection and handling inside small subsystems allows for solvers accelerations, enough to over-compensate the multi-threading costs. Secondly, it appears that combining the enlarged communication step and the context-based extrapolation, the 10 % extra speed-up is reached without loss for the relative error. Even more surprising, using the extrapolation slightly speeds-up the simulation, possibly because the inputs shaped by the predictor enables a faster convergence of the solver step.

Conclusion and future work

The aim of this work is to speed up the numerical integration of hybrid dynamical systems, eventually until reaching a real-time execution, while keeping the integration errors inside controlled bounds. The basic approach consists in splitting the system into submodels, which are integrated in parallel. Using large synchronization intervals between the branches allows for numerical integration speed-ups. However, slack In this paper, the errors caused by the slack synchronization are modeled, giving directions to find effective trade-offs between integration speed and accuracy. Then, an approach of stretching out the communication steps while keeping a predefined integration precision is proposed. Rather than using costly small integration and communication steps, it uses extrapolations of the behavior of the models over the synchronization intervals. Test results on a hybrid dynamical engine model, show that well chosen context-based extrapolation allows for an effective speed-up of the simulation with negligible computing overheads.

This work shows that properly-chosen contextbased extrapolation, combined with model splitting and parallel integration, can potentially improve the speed/precision trade-off needed to reach real-time simulation. However, the accuracy could be widely improved by accessing on the input derivatives of the models. This is the case for the FMI for co-simulation, and it would be highly useful to also integrate this feature in the FMI for model exchange. Future works intend to improve the context-based extrapolation algorithm, to make it more subtly aware of data freshness and even more decrease the prediction induced integration errors. Another possibility is to process the input signals to separate them into simpler components, easier to predict with different predictors, and to cope with noise. When it comes to polynomials, wavelet pre-processors [START_REF] Chaux | Noise covariance properties in dual-tree wavelet decompositions[END_REF] could be useful, as they play an important role in polynomial model fitting.

Figure 2 :

 2 Figure 2: Σ split into Σ 1 and Σ 2 for parallel simulation.

Figure 3 :

 3 Figure 3: System's internal composition.

  Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI 228 Proceedings of the 10 th International ModelicaConference March 10-12, 2014, Lund, Sweden DOI 10.3384/ECP14096225

Session 2A: FMI 2 DOI 10 .

 210 3384/ECP14096225 Proceedings of the 10 th International ModelicaConference March 10-12, 2014, Lund, Sweden

Figure 4 :

 4 Figure 4: Illustration for contexttable in Table 1.

Figure 5 :

 5 Figure 5: Airpath output: pressure.

Figure 6 :

 6 Figure 6: Context behavior during simulation.

Figure 7 :

 7 Figure 7: Cumulative relative error using different communication steps.

Session 2A: FMI 2 DOI 10 .

 210 3384/ECP14096225 Proceedings of the 10 th International ModelicaConference March 10-12, 2014, Lund, Sweden synchronization intervals may generate integration errors in the final result.

  where T denotes the matrix transpose.Here, U 1 are the inputs needed for Σ 1 and U 2 are the inputs needed for Σ 2 . In other words, U 1 ∪ U 2 = U

	Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using
	FMI
	U	Σ	Y	U	U1 U2	X1	Σ 1 Σ 2	X2	Y1 Y2	Y
	Figure 1: System splitting for parallelization.
	and U 1 ∩ U 2 can be an empty set or not according to
	the achieved decoupling.
	In the same way, Y 1 are the outputs produced by Σ 1
	and Y 2 are the outputs produced by Σ 2 . In other words,
	226			Proceedings of the 10 th International ModelicaConference	DOI
						March 10-12, 2014, Lund, Sweden	10.3384/ECP14096225

Σ1 Σ2 Initialization Exchange 1 Exchange 2 Integration step h n Communication step P Communication step P Integration step h n Σ t s t s+1 t n t n+1 t s+2

table in Table

 in

Table 2 :

 2 Relative integration error.

	Outputs	Er(%) w/o extrapolation w/ extrapolation Er(%)
	Pressure	0.499	0.304
	Temperature	0.511	0.19
	Air density	0.784	0.31
	Fuel density	3.55	0.978
	Burned gas density	4.99	3.47

Table 3 :

 3 Simulation speed-up.

	Communication time 100 µs	250 µs
	Extrapolation	No	No	Yes
	Speed-up	8.9 10.01 10.07

f 1 (X(t n )) f 2 (X(t n )) X 2 (t n ) X 2 (t n ) X 1 (t n ) X 1 (t n ) Solver Solver X 1 (t n+1 ) X 2 (t n+1 ) X 1 (t n ) . X 2 (t n ) . Σ1 Σ2 X=[X 1 X 2 ] T

www.lmsintl.com/imagine-amesim-1-d-multi-domainsystem-simulation

2 www.3ds.com/products/catia/portfolio/dymola Session 2A: FMI 2 DOI 10.3384/ECP14096225 Proceedings of the 10 th International ModelicaConference March 10-12, 2014, Lund, Sweden