1 Introduction
Even though the electrification of vehicle powertrain has picked up an aggressive speed in the automotive industry, the majority will still be Hybrid Electric Vehicles (HEVs) at least for the next half century[1]. In the HEVs, the power is provided by two sources: a motor and an internal combustion engine[2]. Thus, to realize the economic and green transportation system of the future society, it is of great importance to increase efficiency and reduce pollutant emissions of the engine by the combustion control[3]. A big issue in the combustion control is how to make a trade off between efficiency and abnormal combustion such as the engine knock[4]. The combustion process is dramatically influenced by the Fuel Injection Timing (FIT) in the diesel engines or the Spark Advance Timing (SAT) in the gasoline engines[5]. Generally, the FIT or SAT with the Maximal Brake Torque (MBT) is chosen to realize the highest efficiency of the energy transformation[6]. However, the FIT or SAT with the MBT might cause a selfignition in an unburned mixture, called end gas. The selfignition often brings the pressure oscillations in the cylinder chamber of the engine. The phenomenon of the pressure oscillation is coupled with a metallic sound from the wall of the cylinder chamber. Thus, it is named the knock event[3], or knock for abbreviation. Small amount of slight knock is good for the combustion efficiency while frequent knock can cause serious cylinder damages and a increase of pollutant emission[6]. The trade off is to operate the engine at a boundary or borderline with a specific knock probability[3]. To realize the boundary control or borderline control, many researches have converged to a statistical feedback control framework[5]
. The statistical feedback controller does not response to the raw measurements directly. Firstly, the statistical properties of the measurements such as the probability density function or probability mass function are estimated. Then, the feedback controller adjusts the control input to realized the desired probability distribution. To evaluate the performance of the statistical knock controller, it needs a batch of experiments since one good case does not mean that the controller would succeed statistically
[3]. The repeats of the experiment tests causes a big expense. To decrease the cost on the experimental tests, it is useful to develop a knock simulator to test the statistical knock controller instead of the experimental tests.The stateoftheart methodologies of designing the knock simulator can be classified into 2 main streams. The first stream is the combustion physical modelbased simulation
[7]. In the combustion physical modelbased simulator, the determinacy of the combustion is described by the Wiebe function and LivengoodWu integration. The stochastic part is modeled by adding some noises on the deterministic part. Then, the cylinder pressure obtained from the ’stochastic’ heat release profile exhibits cycletocycle variations. The knock signal is obtained by setting a threshold for the peak cylinder pressure or the integration of the cylinder pressure. If the threshold is surpassed, the cycle is identified as a knock cycle, otherwise, it is a cycle without knock. The other stream of the knock simulation methodology is recognized as stochastic processbased or Markov processbased simulation[8]. The stochastic processbased simulator focuses particularly on the knock signal itself. It investigates the statistical properties of the knock intensity, which can be calculated from the pressure oscillation[9], or measured by vibration sensor[10]. The probability distribution of the knock intensity varies as the input of the engine or the operating condition of the engine changes. The common method is to calibrate the map from the input and the operating condition to the parameters of the probability distribution. Then, statistical simulation is done according to the parameters of the probability distribution[8]. The map adopts the linear model or exponential model. The probability distribution of the knock intensity adopts Exponentially modified Gaussian distribution
[11, 12]. However, the exist method still have drawbacks:
The conventional simulators only approximate a conditional mean and variance of the data assuming that the data is single Gaussian. Under the single Gaussian assumption, only a very limited statistics can be represented by the conventional simulators. If the mixture distribution model is adopted, it is able to obtain a more complete description of the probability distribution;

The relation from input and operating condition to the parameters of the knock intensity distribution is not well addressed with a simple model, polynomial models or exponential models. The error of the model would bring bias between the distribution of the actual data and the simulated data.

Except for the establishment of the knock simulator, the question that how to validate the simulator quantitatively has not been answered yet in the exist researches.
To address above drawbacks, the knowledge of statistics is necessary. This paper proposes a statistical knock simulator based on the Mixture Density Network (MDN) and the acceptreject method. At first, the statistical analysis is conducted with the experimental data. Regarding the knock intensity as a stochastic process, it is statistical independent and identically distributed under identical input. Moreover, the distribution, or the parameters of the distribution, can be regarded as a function of the input. The MDN is applied to approximate the function from input signal to the knock intensity distribution. The acceptreject algorithm is used to generate knock intensity according to the MDNbased knock intensity distribution model. The proposed method is evaluated in experimental databased validation.
The rest of the paper is organized as followings: Section 2 presents the data analysis results, assumptions about the knock intensity and the problem formulation. Then, section 3 describes the proposed statistical simulator. Section 4 gives the results of experimental databased validations. Finally, section 5 concludes the whole paper.
2 Experimental Data Analysis and Problem Formulation
2.1 Experimental Details and Data Preprocessing
Engine system  Detail 

No. of cylinders  4cylinder 
Arrangement  In line 
Displacement  2499 cc 
Compression ratio  18:1 
Maximum output  87kW 
Maximum torque  280 Nm 
Idle speed  800 rpm 
Maximum engine speed  4200 rpm 
Fuel system  Direct injection 
Ignition system  Compressed ignition 
The data used in this study was collected from the 4JK1TC, a commonrail direct injection intercooled turbodiesel engine. The specifications of the 4JK1TC is listed in Table. 1. During the collection of the database of the engine tests, the 4JK1TC was operated at Wide Open Throttle (WOT) under a variety of speed, manifold pressure, and FIT conditions. These variables affect an engine’s propensity to knock, and their values were chosen to span both a range of knock intensities as well as the normal operating envelope of the engine. The database essentially defines a multidimensional grid in which engine speed is varied in increments of 400 rpm in the range 8002000 rpm, and a variety of manifold pressures are applied, typically within the range 3 bar to 8 bar. Besides, the air/fuel ratio is constantly fixed at 14.6. At each point in the grid, the FIT corresponding to Borderline knock (BL) was identified, and the FIT was then varied typically within the range BL4 to BL+2 degrees.
In total, a series of 118 different operating points were considered. Three data sets were recorded at each point, containing data for 300 consecutive combustion events in each individual record. Each record has the highspeed incycle recordings of cylinder pressure for the cylinder . The knock intensity metric can be obtained by calculating the power spectral density of the cylinder pressure data over the resonant frequency from the incylinder pressure. The details of the calculation process can be referred to [13].
2.2 Statistical Analysis of Data Sets
2.2.1 Statistical independence
Fig. a) shows the cycletocycle knock intensity series obtained from the 300 cycles taken from cylinder of the engine operated at 800 rpm, wide open throttle, and one degrees of fuel injection timing relative to the borderline knock condition. The figure shows that the knock intensity has random cycletocycle variation even under steady operating condition. The scatter plot of Fig. b) gives additional evidence to that the knock intensity has little or no prior cycle correlation. A more quantitative analysis is obtained by calculating the autocorrelation function of the knock intensity for different cycle lags as
(2.1) 
The results are summarized in Fig. 2. As shown in Fig. 2 a), most of the correlations for cycle lags fall within the probability interval which means there is not statistical significance in most cases. However, there are very rare cases in which the correlation falls outside these limits. Although this means that the prior cycle correlation is statistically significant, it should also be noticed that the absolute magnitude of this correlation is still very small, less than 0.2 in all cases. Fig. 2 b) shows the results of autocorrelation function at Lag 1 for the entire database. In most cases, the autocorrelation function at Lag 1 is below 0.2. The rare cases with value larger than 0.2 are all the operating conditions with the increased instability of the combustion process. According to the above statistical independence analysis results and discussions, the assumption about the knock intensity under a given operating condition can be summarized as:
Assumption 2.1.
The knock intensity is independently and identically distributed random variables for a given operating condition.
2.2.2 Probability distributions
The results of the distribution fit performance are summarized in Fig. . The comparisons of the relative frequency and cumulative density function of the knock intensity are shown which exhibits that Gaussian mixture models have better performance on fitting the relative frequency and cumulative probability. Besides, a more concise measure of the goodness of fitting is to compare the fitting errors of each model calculated from historical data for all operating conditions. The fitting error is defined as:
(2.2) 
where denotes the th measured output variable in a data record of length , represents the modelestimated output[14]. Small values indicate good model fits. Here, defining as the cumulative probability or relative frequency, and , where denotes the measured value for the knock intensity. The approach was replicated here, giving the results shown in Fig. 4. As the number of kernels increases, both fitting errors decrease, sharply from one to two and then gently from three. Thus, for the tradeoff of accuracy and complexity, the number of kernels can be chosen as values from 25.
According to the above statistical analysis of the knock intensity distribution, the following assumption can be summarized:
Assumption 2.2.
Gaussian mixture models have better accuracy on fitting the distribution of the knock intensity than single Gaussian model. As numbers of the used kernels increases, the fitting accuracy increase.
Fig. 5 a) shows a example of the knock intensity relative frequency curves for a range of different fuel injection timing taken from cylinder 1 of the engine operated at 1200 rpm and manifold pressure of 7 bar. The corresponding cumulative probability curves are plotted in Fig. 5 b). As the fuel injection timing is advanced, the curves of relative frequency and cumulative probability shifts to the right side which means the means and modes of the distribution clearly increase. This reflects the strong influence of fuel injection timing on the engine’s propensity to knock intensity. Similar results are obtained at a range of different engine speed from cylinder 1 of the engine operate at manifold pressure of 7 bar and borderline +2 or at a range of different manifold pressure from cylinder 1 of the engine operated at 1200 rpm and borderline +2. These observations highlight the fact that the fuel injection timing, manifold pressure and engine speed all influence the knock intensity distribution strongly. More details are shown in Fig. 6. The characteristic parameters of knock intensity’s distribution, such as mean, variance and manifold pressure, vary as the fuel injection timing, manifold pressure and engine speed change. According to the above discussion, the assumption on the relationship between the operating condition and the knock intensity distribution can be summarized:
Assumption 2.3.
The parameters of the knock intensity distribution are affected by the input and operating condition of engine, such as fuel injection timing, manifold pressure and engine speed. The parameter vector can be regarded as a function of the input and operating condition of engine.
2.3 Problem Formulation
Denote the parameter vector of the engine operating conditions as and the knock intensity as . According to Assumption 2.1, 2.2 and 2.3, probability density of depends on and expressed as . The addressed problem is to obtain an approximation probability density of the probability density using a training sample set . The problem is expressed as
(2.3) 
3 Proposed Statistical Knock Simulator
The schematic of the proposed simulator is illustrated in Fig. 7. Firstly, a Mixture Density Network is trained using the training sample set. The obtained Mixture Density Network outputs the parameter vector of the approximated probability density function of the knock intensity. The parameter vector is conditioned on the input and operating conditions of the engine such as FIT, engine speed and manifold pressure. Then, an AcceptReject algorithm is applied to generate the knock intensity data corresponding to the conditioned probability density function. The generated data will be checked whether their distribution are consistent with the data in testing sample set.
3.1 Mixture Density Networks
In the framework of MDN, a mixture model is used to approximate the probability density of the target data instead of using the Gaussian distribution model. The mixture model has the flexibility to model completely general distribution functions [15]. The approximate probability density of the target data is then represented as a linear combination of kernel functions in the form
(3.1) 
where is the number of components in the mixture model. The parameters are called mixing coefficients which is a prior conditional probabilities (conditioned on ) of the knock intensity having been generated from the th component of the mixture. Note that the mixing coefficients are also taken to be functions of the input vector and satisfies the following constraints
(3.2) 
(3.3) 
The functions denote the conditional density of the knock intensity for the th kernel. Various choices for the kernel functions are possible. In this research, the chosen kernel functions are Gaussian of the form
(3.4) 
where and represent the centre and variance of the th kernel respectively. In principle, a Gaussian mixture model with kernels given by (3.4), can approximate any given density function to arbitrary accuracy if the mixing coefficients and the Gaussian parameters (means and variances) are correctly chosen[15]. Therefore, the representation given by (3.1) and (3.4) is completely general. Formally, as given in [16], the above discussions can be summarized in Theorem 3.1.
Theorem 3.1.
For a given , , let and let for be such that and . Assume that is a continuous squareintegrable function. The integrated squared bias and integrated variance of the Gaussian mixture density model expressed by (3.1) and (3.4) have asymptotic behavior
(3.5) 
and
(3.6) 
respectively. The firstorder asymptotic Approximation of Mean Integrated Squared Error(AMISE), is thus given by
(3.7) 
The asymptotically optimal value of is minimizer of the AMISE
(3.8) 
giving the minimum value
(3.9) 
Remark 3.1.
Theorem 3.1 gives the upper boundary of AMISE when using the Gaussian mixture density model expressed by (3.1) and (3.4) since can be optimized to get smaller AMISE instead of fixing them as two constants. Namely, if are well estimated, the AMISE of optimal satisfies
(3.10) 
Even the upper boundary converges to 0 if and increase to , thus the Gaussian mixture density model is able to converge to the actual probability density theoretically.
On the other hand, for any given , the mixture model expressed by (3.1) gives a general formalism for modelling an arbitrary conditional density function . Denote the parameter vector of the mixture model as
(3.11) 
The parameter vector is taken to be a continuous function of
. The continuous function can be approximated by using a conventional neural network which takes
as input and as output. According to [18], conventional neural network can approximate any given continuous function to arbitrary accuracy. Formally, the approximate property is summarized in Theorem 3.2.Theorem 3.2.
Given a bounded nonlinear activation function
for which there exists , then for any arbitrary distinct samples , there exist , , and , , such that(3.12) 
Theorem 3.1 and 3.2 are the theoretic basis of MDN. The basic structure is illustrated in Fig. 8. By choosing a mixture model with a sufficient number of kernel functions, and a neural network with a sufficient number of kernel functions, the MDN can approximate any conditional density as closely as desired. Define the parameter of the neural network as , problem given by (2.3) can be then equivalently replaced by the problem defined as
(3.13) 
Therefore, the probability density can be approximated by a Mixture Density Network. The parameter vector of the Mixture Density Network is obtained through solving problem expressed by (3.13) with a training sample set . The gradient decent algorithm can be used for training the parameter vector [15]. The obtained Mixture Density Network has a property that it has the maximal likelihood to generate the data with the identical distribution of the training sample set.
3.2 AcceptReject algorithm for statistical simulation
It is difficult to directly generate data since the approximate probability density is too complex. AcceptReject method can be applied to simulate data of arbitrary probability density. The philosophy of AcceptReject is summarized as:
Theorem 3.3.
Let and let
be a simple density function (for example, uniform distribution) that satisfies
for some constant . Then, to simulate , it is sufficient to generate and , until . means uniform distribution conditioned on .The proof of Theorem 3.3 can refer to Chapter 2 of [19]. Based on Theorem 3.3, the AcceptReject algorithm is formed which is summarized in Algorithm 1.
Remark 3.2.
If well approximates the probability density , , Algorithm 1 can generate the data which has the same conditional probability density of .
4 Validation Result
This section gives the validation results using experimental data set. Both steady and transient cases are concerned.
4.1 Validation methodology
In the validation, the data for testing is different from the data used for training the parameters of the MDN models.
For the steady case, the data set mentioned in section 2 is used. The data set is firstly separated into 118 different subsets according to the operating conditions. Each subsets is given a label . Then, 118 different validation simulations were implemented. In each simulation, the subset was chosen as the test set and the rest subsets were training set. The MDNbased simulators with kernel numbers from 1 to 20 were validated in each simulation. Each simulator generated 500 groups of data set in each test. In every group, the number of simulated data is 900.
For the transient case, the training set is the whole data set mentioned in section 2. The data of the test set was obtained in another experiment under a transient case.
4.2 Steady case
Fig. 9 shows the comparison between the data from the test set and the data generated by the simulators. Fig. 9 a) gives the cycletocycle plots of real data, simulated data using MDN with single Gaussian, 2 Gaussian kernels, 10 Gaussian kernels respectively. The engine was operated at speed of 2000 rpm. Besides, the manifold pressure was 7 bar and the relative fuel injection timing was operated at the borderline knock condition. Fig. 9 shows results of more relative fuel injection timings, from 3 to 1. Intuitively, the simulated data is distributed more closely to the real data if more kernels are used in the model.
A more clear demonstration of the characteristics of data distribution shown in Fig. 9 is given in Fig. 10. The simulator using 10 Gaussian kernels exhibits better fitting performance. The quantitative results of the fitting performance evaluation are plotted in Fig. 11. The statistics of the fitting error for each simulator are plotted in Fig. 11 a). Instead of using modelestimated output as , the cumulative probability calculated from the simulated data is used. Since each simulator generated 500 groups in each validation simulation, the mean and variance of the fitting error were calculated. The solid line with circles represents the mean error and the dashed lines is the bounds of 99 confidence interval. The resulted plots show that mixture models exhibit much better fitting performance than the single Gaussian model. Besides, Fig. 11 b) gives the KolmogorovSmirnov statistics which is expressed as[20]
(4.1) 
The results of the KolmogorovSmirnov statistics have the same trend as the fitting error. Fitting error represents the mean square error and the KolmogorovSmirnov statistics show the maximal absolute value of error. Increasing the kernel number can decrease both.
4.3 Transient case
In the transient case, the engine was operated at 1200 rpm. The manifold pressure is 7 bar. The relative fuel injection timing was firstly kept as constant and then changed to another constant. Fig. 12 shows the relative fuel injection timing evolution of the transient case. The corresponding test set data and simulated data by several different simulators are shown in Fig. 13. Mixture models exhibit better performance than the single Gaussian model.
5 Conclusion
This paper proposes a statistical knock simulator based on the mixture density network and the acceptreject algorithm. The knock intensity is a stochastic process. With a identical input, it is independent and identically distributed. Moreover, the distribution is a function of the input. The mixture density network is applied to approximate the function from input signal to the knock intensity distribution. The acceptreject algorithm is used to generate knock intensity according to the mixture density networkbased knock intensity distribution model. The proposed method is evaluated in experimental databased validation. Using mixture models can improve the performance of the simulator. With more kernels, the simulator is able to output the knock intensity that has closer distribution with the real data. However, there is still future work should be done. Especially, how to quantitatively evaluate the simulator performance for the transient case.
References
 [1] Lewis, A. M., Kelly, J. C., and Keoleian, G. A. (2014). Vehicle lightweighting vs. electrification: life cycle energy and GHG emissions results for diverse powertrain vehicles. Applied Energy 126(1) 13–20.
 [2] Chau, K. T., and Wong, Y. S. (2002). Overview of power management in hybrid electric vehicles. Energy Conversion and Management 43(15) 1953–1968.
 [3] Eriksson, L., and Nielsen, L. (2014). Modeling and control of engines and drivelines, 1st ed. John Wiley and Sons Ltd.
 [4] Bares, P., Selmanaj, D., Guardiola, C., and Onder, C. (2018). A new knock event definition for knock detection and control optimization. Applied Thermal Engineering 131 80–88.
 [5] Payton, J. J., Spelina, J., and Frey, J. (2013). Likelihoodbased control of engine knock. IEEE Transactions on Control System Technolodge 21(6) 2169–2180.
 [6] Kiencke, U., and Nielsen, L. (2005). Automotive Control Systems for Engine, Driveline, and Vehicle, 2nd ed. Berlin, Germany: Springer.
 [7] Di, D., and Shen, T. (2018). Simulation of knock probability in an internal combustion engine. Phys. Rev. E 98(1) 2102–2011.
 [8] Payton, J. J., Spelina, J., and Frey, J. (2016). Controloriented knock simulation. SAE Int. J. Engines 9(2) 1201–1209.
 [9] Panzani, G., Ostman, F., and Onder, C. H. (2018). Engine knock margin estimation using incylinder pressure measurements. IEEE/ASME transactions on Mechatronics 22(1) 301–311.
 [10] Boubai, O. (2000). Knock detection in automobile engines. IEEE Instru. Meas. Mag. 3(3) 2428.
 [11] Spelina, J., Payton, J. J., and Frey, J. (2014). Characterization of knock intensity distributions: Part 1: statistical independence and scalar measures. P. I. Mech. Eng. DJ. Aut., 228(2) 117–128.

[12]
Spelina, J., Payton, J. J., and Frey, J.
(2013). Characterization of knock intensity distributions: Part 2: parametric models.
P. I. Mech. Eng. DJ. Aut., 227(12) 1650–1660.  [13] Shen, X., Zhang, Y., and Shen, T. (2019). Cylinder pressure resonant frequency cyclic estimationbased knock intensity metric in combustion engines. Applied Thermal Engineering, 158 113756.
 [14] Devore, J. L. (2009). Probability and statistics for engineering and the sciences. Brooks/Cole, Boston.
 [15] McLachlan, G. J., and Basford, K. E. (1988). Mixture Models: Inference and Applications to Clustering. New York: Marcel Dekker.

[16]
Botev, Z. I., Grotowski, J. F., and Kroese, D. P.
(2010). Kernel density estimation via diffusion.
The Annuals of Statistics, 38(5) 2916–2957. 2722460.  [17] Wand, M. P., and Jones, M. C. (1995). Kernel Smoothing. Chapman and Hall, London. 1319818.

[18]
Huang, G, and Babri, H. A.
(1998). Upper Bounds on the Number of Hidden Neurons in Feedforward Networks with Arbitrary Bounded Nonlinear Activation Functions.
IEEE transactions on Neural Networks 9(1) 224–229.  [19] Robert C. P., and Casella, G. (2005). Monte Carlo Statistical Methods, 2nd ed. SpringerVerlagBerlin, Heidelberg.
 [20] Stephens M. A. (1973). EDF Statistics for Goodness of Fit and Some Comparisons. Journal of the American Statistical Association, 69(347) 730–737.
Comments
There are no comments yet.