**1. Introduction**

The development of cleaner and more efficient engines requires the continuous measurement of in-cylinder pressure which is of fundamental importance for combustion optimization, air-fuel ratio control, noise and pollutant reduction. More specifically, cylinder pressure is a fundamental intermediate variable that indicates engine state and drives models for engine combustion control. High-bandwidth control of ignition and the air-fuel ratio may permit these engines to operate reliably near the lean-burn limit, with significant improvements in efficiency and decreased emission of CO and NOx. Fully utilizing this technology requires a key factor: real-time cylinder pressure [1]. The goal of Rizzoni's [2] method is to build a stochastic model for the combustion pressure process in the

spark-ignition engine. Control of torque balance can improve drivability and can suppress noise from light-duty diesel engines, but as a crucial part of the control algorithm, obtaining pressure data necessitates expensive pressure sensors and demands considerable computational time [3]. The advent of cylinder-pressure transducers seems promising, boasting a better approach to detect the in-cylinder pressure; however, its reliability, complex installation, high cost, and short working lifespan are full of controversy. As new technologies improve, such as computational methods, reconstructing cylinder pressure in a multi-cylinder internal combustion (IC) engines becomes more feasible. Bennett [4] proposed a robust adaptive-gradient descent-trained NARX neural network using both crank velocity and crank acceleration as inputs to reconstruct cylinder pressure in multi-cylinder IC engines. Another work has shown comprehensive results from an experimental assessment of a common rail diesel engine operated with neat diesel fuel and with ethanol substitutions [5]. With di fferent substitutions and di fferent substitution rates, pressure varies rapidly. Researchers have pursued other various approaches to develop a comparatively optimal way to obtain the desired cylinder pressure, directly or indirectly. Rizvi [6] proposed a novel method to detect engine misfire faults based on a hybrid model of the gasoline engine. As this method mainly employs cylinder pressure, closed-loop [7] combustion control [8] becomes feasible and misfire can be detected in multiple ways. However, these approaches vary in cost, reliability, robustness, accuracy, and convenience. Therefore, a low-cost [9] noninvasive soft-pressure sensor with a high level of reliability and accuracy is necessary. Payri [10] presented a step-by-step approach to optimize the signal processing both for o ffline and online applications based on the characteristics of the signal. Maurya [11] proposed that a method based on standard deviation of pressure, and pressure rise rate is used to find the minimum number of engine cycles to be recorded for averaging to ge<sup>t</sup> reasonably accurate pressure data independent of cyclic variability. Bhatti [12] proposed that two robust second-order sliding mode observers are employed that require two state mean value engine models based on inlet manifold pressure and rotational speed dynamics. The design of the second-order sliding mode observer depends largely on the model accuracy. And, more importantly, this method had the problems of chattering and large initial error. As a result, the tracking error of the steady-state speed was only about 4% and required a certain convergence time. Ferrari [13] proposed that pressure of nozzle opening and closing was identified by means of pressure data fitted by bench test. Then the time-frequency analysis was used to obtain the mean instantaneous frequency and adjust the control strategy of the injector. However, this method did not process data e ffectively, such as filter, before fast Fourier time-frequency transform, so this may a ffect the accuracy. Eriksson [14] calculated the corresponding relationship between pressure value and crankshaft rotation angle. The error of peak identification based on ion current in (Eriksson, 2003, Figure 4) was still large. Comparing with the corresponding relationship proposed by Eriksson [14], the cylinder pressure time curve is presented in this paper. It is based on burning theory and Law of energy conservation, and then modified by Extended-Kalman-Filter-optimized (EKF) engine model, which has a higher stability.

Observer-based method [12–15] is proven to be a good way to predict the output of SISO system. Engine cylinder pressure signal can be used to characterize the engine combustion state [16,17]. Due to the influence of multiple variables, it is di fficult to establish a precise physical model of cylinder pressure. However, by observing and analyzing the cylinder pressure signal, this signal is cyclical with variable frequency (associated with speed) and variable amplitude (associated with the peak combustion pressure). Meanwhile, the engine cylinder pressure signal has a delay characteristic, so it is necessary to predict the current cylinder pressure value by the input of the current moment.

In this paper, a novel method which can identify engine pressure data accurately will provide a spacious viewpoint and means to process signals, thereby achieving important information for practical engineering. This novel method involving Frequency-Amplitude Modulated Fourier Series (FAMFS) and Extended-Kalman-Filter-optimized (EKF) engine model for identifying the periodic signal with its variable frequency and amplitude is proposed. This cylinder pressure identification method is also called a virtual cylinder pressure sensor.

Figure 1 shows a flow chart of method. Crankshaft speed [18–20] plays an important role in identifying the pressure as well as other applications, especially in misfire detection [21]. Instantaneous speed fluctuations [22], crankshaft segmen<sup>t</sup> acceleration, and transient rotational speed are the most widely utilized variables to locate misfire. Crank speed can also reveal and identify real-time engine combustion parameters [20], which is adopted in this paper. Speed data is affected by measurement noise from the crankshaft and process noise from the engine. These noise sources result in a huge accumulation error during pressure measurements, and EKF is chosen to optimize prediction and filtering [23,24].

**Figure 1.** Flow chart of method.

The crankshaft speed in one solution step is related to the speed of the last step, so the speed model is an iterative process. Firstly, a novel speed iteration model based on burning theory and Law of energy Conservation is proposed. This iteration model is associated with the throttle opening value and the load. The EKF is then used to estimate the optimal output of the speed iteration model, and this optimal output can then be used for computing the frequency and amplitude of the cylinder pressure cycle-to-cycle. Secondly, taking Otto cycle of standard four-stroke engine as an example, a 24th order Fourier series is chosen to fit the standard working cycle. Thirdly, associated with the variable frequency and amplitude, a frequency-amplitude modulation method is adopted to modulate the standard pressure model identified in step two. Finally, system performance is evaluated by an actual working engine. Therefore, the main contributions of this paper can be summarized as follows:


In this paper, on the one hand, as the key parameter of cylinder pressure identification, the speed iterative model using the EKF has been predicted with high accuracy. However, the overall deviation came from the calculation of pressure due to the unpredictable process of power stroke. On the other hand, due to only related to the crankshaft speed and the amount of air and fuel injection in engine, this method can be applied to multi-cylinder internal combustion (IC) engines including four-cylinder engines. More specifically, many key parameters of ignition engine such as spark advance angle (SAA) need to be considered in the process of building virtual in-cylinder pressure sensor model. Accordingly, at present, the application of the virtual sensor is only limited to spark ignition engines.

#### **2. Crankshaft Speed Model**

#### *2.1. Design of Speed Iteration Model*

As generally known, the angular acceleration of crankshaft in one cycle is produced by the output power of the engine, and the output energy per cycle is determined by the amount of air and fuel injected into engine. In internal combustion engines, the amount of air and fuel trapped in the cylinder is relatively influenced by engine speed through the volumetric e fficiency, and mainly depends on the engine load. The change of engine load will a ffect crankshaft speed; hence the amount of air and fuel is deemed relating to the speed. Therefore, it forms a closed-loop iteration process. The crankshaft speed model (CSM) is established in two parts: constant load (Part A) and variable load (Part B). In Part A, the amount of air and fuel is related to the speed of the previous cycle. The power produced through fuel burning compels the crankshaft to rotate; and in Part B, the engine overcomes the load with the price of crankshaft speed reduction. The complete model is then forged by the combination of these two parts.

#### 2.1.1. Constant Load Condition

In the following process, the load is constant. Cylinder pressure is directly related to the amount of fuel and air injection in engine. Thus, the air volume is equal to the engine displacement. The air mass flow of the intake manifold can be calculated as

$$
\dot{m}\_{\rm air} = \frac{S\_{eff} \mathbb{C}\_{\eta} \mathbb{C}\_{m} P\_{\rm im}}{\sqrt{T\_{\rm im}}} \tag{1}
$$

where, *Seff* is the active area of throttle valve; *Cq* is the flow coe fficient; *Cm* represents the quality factor, associating with the engine. Once the engine model is specified, these parameters such as *Seff* , *Cq* and *Cm* are determined. *Tim* and *Pim* are the temperature and pressure of intake manifold, respectively.

According to the actual working conditions, air mass is then obtained by Refs. [25–28]:

$$m\_{\rm air} = A(\varepsilon) \int\_0^t \dot{m}\_{\rm air} dt \approx \frac{\pi A(\varepsilon) \dot{m}\_{\rm air}}{N} = \frac{\pi A(\varepsilon)}{N} \frac{\text{SC}\_{\rm q} \text{C}\_{\rm m} P\_{\rm im}}{\sqrt{T\_{\rm im}}} (1 - \cos \theta) \tag{2}$$

where, the value for ε is the engine volumetric coe fficient and considered constant when the crankshaft speed does not dramatically change. The integration of air mass flow will inevitably produce a constant term. The sum of constant term and non-constant term is air mass. Since there must be a proportional relationship between the constant term and the non-constant term, in order to facilitate the subsequent calculation, the proportional relationship is defined as *<sup>A</sup>*(ε), so it is considered to be the engine parameter coe fficient related to ε. Once the engine model is specified, these two parameters are determined. In addition, *t* = π *N* , and *t* is opening time of intake valve, and *N* is the average value of speed in one working cycle.

Di fferent fuels and di fferent fuel ratios will strongly a ffect the combustion state of the engine [5]. The main elements in the fuel are oxygen, hydrogen, and carbon, with other elements being negligible. Thus, the relationship between mass fractions for the oxygen element *wo*, hydrogen element *wh*, and carbon element *wc* in 1 kg of fuel are shown below:

$$w\_{\mathcal{C}} + w\_{\mathcal{h}} + w\_{\mathcal{o}} = 1 \tag{3}$$

The theoretical air volume for 1 kg of fuel that is entirely burned is shown in Equation (4):

$$L\_O = \frac{22.4}{0.21} (\frac{w\_C}{12} + \frac{w\_H}{4} - \frac{w\_O}{32}) \tag{4}$$

The engine output power is determined by the calorific value of the gas mixture, which can be defined as (considering 1 kg of fuel)

$$Q\_{\rm mix} = \frac{h\_u}{22.4 \times \left(\lambda \frac{I\_Q}{22.4} + \frac{1}{M\_\odot}\right)}\tag{5}$$

where, λ is the excess-air factor; *hu* is the low calorific value of the fuel; and *M*τ is the relative molecular mass of fuel. Therefore, the power output related to the throttle valve is defined as

$$Q = \left[ (L\_{\text{O}} \rho\_{air} + 1) \frac{m\_{air}}{4 \rho\_{fuel} L\_{\text{O}} \rho\_{air}} \right] \frac{h\_u}{22.4 \times \left( \lambda \frac{L\_{\text{O}}}{22.4} + \frac{1}{M\_{\text{I}}} \right)} \tag{6}$$

where, ρ*air* and ρ*f uel* is the density of air and fuel, respectively. By substituting Equation (2) into Equation (6), Equation (6) then becomes

$$Q = \left| (L\_O \rho\_{\rm air} + 1) \frac{1}{4 \rho\_{\rm fuel} L\_O \rho\_{\rm air}} \frac{\pi A(\varepsilon) h\_u}{(\lambda L\_O + \frac{22.4}{M\_\pi})} \frac{\mathcal{SC}\_q \mathbb{C}\_m P\_{\rm im}}{\sqrt{T\_{\rm im}}} \right| \cdot \frac{(1 - \cos \theta)}{N} = B \frac{(1 - \cos \theta)}{N} \tag{7}$$

where,

$$B = (L\_{\bullet} \rho\_{air} + 1) \frac{1}{4 \rho\_{fuel} L\_{\bullet} \rho\_{air}} \frac{\pi A(\varepsilon) h\_{\rm u}}{(\lambda L\_{\bullet} + \frac{22.4}{M\_{\rm r}})} \frac{\rm SC\_{\rm q} C\_{\rm m} P\_{\rm im}}{\sqrt{T\_{\rm im}}} \tag{8}$$

As widely known, the total energy produced through burning cannot be completely converted to useful work. The process of energy transfer is shown in Figure 2.

**Figure 2.** Flowchart of Engine energy transfer.

η*t* is the theoretical thermal efficiency; η*j* is the loss coefficient; η*m* is the mechanical efficiency; and η*e* is the engine efficiency.

In a four-cylinder four-stroke engine, each cylinder burns once per cycle, and the crankshaft is rotated 720 degrees in one working cycle.

$$\mathcal{W}\_c = \mathcal{U}\_{kin}(q) + \mathcal{U}\_{pot}(q) = \frac{I}{2}\mathcal{N}^2 \tag{9}$$

where, *Ukin*(ϕ) is the kinetic energy of the system; *Upot*(ϕ) is the potential energy of the system; *J* is the rotational inertia parameter (in this paper, it is equal to a flywheel's rotational inertia parameter [24]); and ϕ is the crankshaft angle. By combining Equations (7)–(9), the CSM for Part A is shown in the following equation, and the relationship between *Nk*|*k* and *Nk*|*k*+<sup>1</sup> is shown in Figure 3.

$$
\eta\_{\epsilon} B \frac{(1 - \cos \theta)}{N\_{k|k}} = \frac{I}{2} N\_{k|k+1}^2 \tag{10}
$$

**Figure 3.** Relationship of speed in the iteration model.

#### 2.1.2. Variable Load Condition

The change of engine load will a ffect crankshaft speed. This signal is discretized by computing the mean value of crank speed in power stroke. Since the time step of power stroke is enough small, Equations (11) and (12) are approximately satisfied:

$$T - T\_L = J \frac{dN}{dt} = J \frac{N\_{k+1|k+1} - N\_{k|k+1}}{t\_{\text{cyc}}} \quad \Rightarrow \quad N\_{k+1|k+1} = \frac{T - T\_L}{J} t\_{\text{cyc}} + N\_{k|k+1} \tag{11}$$

$$t\_{cyc} = \frac{1}{\frac{2N\_{k|k+1}}{60}} = \frac{30}{N\_{k|k+1}}\tag{12}$$

where, *J* is the rotational inertia; *TL* is the load torque; and *tcyc* is the working cycle time. *T* is the engine output torque and can be obtained by the empirical model [25], and its coe fficients are obtained by calibrating the actual engine.

$$\begin{aligned} T &= -181.3 + 379.36m\_{\text{kif}} + 21.91R\_{\text{A/F}} - 0.85R\_{\text{A/F}}^2 + 0.26\sigma + 0.0028\sigma^2 + 0.027N\_{\text{k}\parallel\text{k}+1} - 0.000107N\_{\text{k}\parallel\text{k}+1}^2 + \dots \\ &0.00048N\_{\text{k}\parallel\text{k}+1}\sigma + 2.55m\_{\text{k}\parallel\text{}} - 0.05\sigma^2m\_{\text{k}\parallel\text{}} \end{aligned} \tag{13}$$

where, *mair* is the mass of air in cylinder for combustion (g); σ is the spark advance angle (SAA).

In this case, only the load *TL* is variable, thereby satisfying Equation (14).

$$T = C + DB\frac{1 - \cos\theta}{N\_{k\bar{k}}} + 0.027N\_{k\bar{k}+1} - 0.000107N\_{k\bar{k}+1}^2 + 0.00048N\_{k\bar{k}+1}\sigma \tag{14}$$

where *C* = −181.3 + 21.91 *RA*/*F* − 0.85 *R*<sup>2</sup> *A*/*F* + 0.26σ + 0.0028σ2, *D* = 379.36 + 2.55σ − 0.05σ2. From this analysis, the model for Part B can be described by:

$$N\_{k+1\vert k+1} = \frac{\mathcal{C} + \text{DB} \frac{1-\text{q}\,\text{m}\Omega}{N\_{\text{k}\vert k}} + 0.027\text{N}\_{\text{k}\vert k+1} - 0.000107\text{N}\_{\text{k}\vert k+1}^2 + 0.00048\text{N}\_{\text{k}\vert k+1}\sigma - T\_L}{f} \cdot \frac{\text{50}}{N\_{\text{k}\vert k+1}} + N\_{\text{k}\vert k+1} \tag{15}$$

By combining Equations (12) and (15), the speed iteration model can then be established:

*Nk*+1|*k*+<sup>1</sup> = *C*+*DB* 1−cos θ *Nk*|*k* +0.027 ) 2η*eB* (<sup>1</sup>−cos θ) *Nk*|*k J* −0.000107 2η*eB* (<sup>1</sup>−cos θ) *Nk*|*k J* +0.00048 2η*eB* (<sup>1</sup>−cos θ) *Nk*|*k J* <sup>σ</sup>−*TL J* · 30 2η*eB* (<sup>1</sup>−cos θ) *Nk*|*k J*+ 2η*eB* (<sup>1</sup>−cos θ) *Nk*|*k J* (16)

Equation (16) can be expressed as

$$N\_{k+1|k+1} = \mathcal{g}(N\_{k|k}, \theta) \Rightarrow N\_{k+1} = \mathcal{g}(N\_k, \theta) \tag{17}$$

where, *g*(·) is a nonlinear function.

> The relationship among *Nk*|*k*, *Nk*|*k*+1, and *Nk*+1|*k*+<sup>1</sup> is depicted in Figure 3.

#### *2.2. Design of Extended Kalman Filter*

Under actual working conditions, the throttle opening value can be decomposed into two parts: an ideal opening value and throttle position signal noise mainly due to EMI. In measuring the crankshaft speed process, speed data is a ffected by measurement noise from the crankshaft and process noise from the engine. Engine vibration caused by engine knock, etc., is transmitted to the crankshaft. In addition, the disturbance caused by road surface irregularity is also gradually transmitted to crankshaft. According to Kalman filter theory, these fluctuations are considered as process noise and need to be processed. These noise sources result in a huge accumulation error during pressure measurements, and the Kalman filter (KF) [29,30] is the most preferable filtering method for measurement and process noise. The KF is used to predict an optimal output for the current step based on the optimal output of the previous step and the observed value of the current step. However, for nonlinear systems, the KF cannot achieve optimal prediction. Therefore, the Extended Kalman filter (EKF) [31–33] is proposed by using the local linear property of nonlinear systems. State prediction is accomplished by using EKF gain to update state and error covariance. After linearization by computing the nonlinear function's first order Taylor series and Jacobian matrix at the operating point, the EKF algorithm matches the KF and it comprises two parts: time update and state update. A discrete-time state-space model with a zero mean and a *Qk* variance in the white process noise (*wk*−<sup>1</sup>) is shown as

$$X\_k = \mathcal{g}(X\_{k-1}) + w\_{k-1} \tag{18}$$

where, *Xk* is the state vector of the *kth* step, and *g*(·) is a nonlinear function.

The state observation can be written as

$$Z\_k = h(\mathfrak{x}\_k) + v\_k \tag{19}$$

where *vk* is the white measurement noise with a zero mean; *Rk* is the variance; and *h*(·)is a nonlinear function.

The prediction is shown below:

$$
\hat{X}\_{k|k-1} = \mathcal{g}(X\_{k-1}) \tag{20}
$$

The first-order Taylor series for *g*(·) at *X*ˆ *k*−1|*k*−1 and *h*(·) at *X*ˆ *k*|*k*−1 are obtained as shown in Equations (21) and (22), respectively.

$$\mathbf{G}\_{k} = \left. \frac{\partial \mathbf{g}}{\partial \mathbf{X}} \right|\_{\hat{X}\_{k-1|k-1}} = \mathbf{g} \left( \hat{\mathbf{X}}\_{k-1|k-1} \right) + \frac{\delta \mathbf{g}(\mathbf{X}\_{k-1})}{\delta \mathbf{X}\_{k-1}} |\_{\hat{X}\_{k-1} = \hat{\mathbf{X}}\_{k-1|k-1}} \left( \mathbf{X}\_{k-1} - \hat{\mathbf{X}}\_{k-1|k-1} \right) \tag{21}$$

$$H\_k = \left. \frac{\partial h}{\partial \mathbf{X}} \right|\_{\mathbf{X}\_{k|k-1}} = h \Big( \hat{\mathbf{X}}\_{k|k-1} \Big) + \frac{\delta h(\mathbf{X}\_k)}{\delta \mathbf{X}\_k}|\_{\mathbf{X}\_k = \hat{\mathbf{X}}\_{k|k-1}} \Big( \mathbf{X}\_k - \hat{\mathbf{X}}\_{k|k-1} \Big) + \upsilon\_k \tag{22}$$

Let:

$$
\Delta\_k = X\_k - \hat{X}\_k \tag{23}
$$

Defining *Pk*+<sup>1</sup> as the covariance of Δ*k*+1, the transcendental error covariance can be calculated as

$$P\_{k+1|k} = G\_k P\_k G\_k^T + Q\_k \tag{24}$$

The major purpose of the EKF is to obtain a minimum-variance state estimate. The gain matrix of the EKF is computed subject to minimization of the estimation error covariance.

$$\operatorname{Tr}\Big[P\_{k+1}^{\mathrm{x}}\big] = \min\_{K\_{k+1}} \mathbb{1} \Rightarrow K\_{k+1} = P\_{k+1\mid\mathbb{k}} H\_k^T \left(H\_k^T P\_{k+1\mid\mathbb{k}} H\_k^T + R\right)^{-1} \tag{25}$$

State prediction is accomplished by using EKF gain to update state and error covariance.

$$
\hat{X}\_{k+1} = \hat{X}\_{k+1|k} + K\_{k+1}(Z\_k - H\_k \hat{X}\_{k+1|k}) \tag{26}
$$

$$P\_{k+1} = (E - \mathcal{K}\_{k+1}H)P\_{k+1\mid k} \tag{27}$$

The observed vector is written as follows, where the observed matrix *H* = 1:

$$Z\_k = H\_k \mathbf{x}\_k + \mathbf{v}\_k \tag{28}$$

#### **3. Calculation of Cylinder Pressure**

In this paper, due to only related to the crankshaft speed and the amount of air and fuel injection in engine, this method can be applied to multi-cylinder internal combustion (IC) engines including four-cylinder engines. More specifically, many key parameters of ignition engine such as spark advance angle (SAA) need to be considered in the process of building virtual in-cylinder pressure sensor model. Accordingly, at present the application of the virtual sensor is only limited to spark ignition engines. In addition, now there are many ways to describe the status of engine working cycle such as Otto-cycle, Diesel-cycle, Sabtache-cycle, and Atkinson-cycle. Among the many cycles, the Otto-cycle optimally depicts the four-stroke ignition combustion engine. Taking Otto cycle of standard four-stroke engine as an example, but the virtual sensor method is not limited to the engine described by Otto cycle.

In this cycle, the crankshaft rotates 720◦ as one working period, where the 720◦ can be divided into four sub-cycles each with a crankshaft rotation of 180◦ [34]. A diagram of the standard working cycle is shown in Figure 4. These four smaller cycles are known as: intake stroke, compression stroke, power stroke, and exhaust stroke. It is the power stroke that involves the burning and expanding process. According to the Otto-cycle, the cylinder pressure during the intake stroke and the exhaust stroke is equal to atmosphere. The compression stroke and the expanding process are considered as isentropic, and the burning process is regarded as constant volume combustion [35].

**Figure 4.** Diagram of the standard working cycle.

The pressure at the end of compression stroke is also the initial pressure of the burning process. The peak occurs at the end of the burning process and is known as the peak pressure of the cylinder. This pressure is an important variate to understand the combustion e fficiency of the burning process and to monitor the state of the cylinder.

#### *3.1. Cylinder Pressure Peak Time Confirmation*

It's generally accepted that when the spark advance angle (SAA) is 10–15◦ after top-dead-center (TDC), combustion e fficiency becomes maximum, while engine vibration reaches a minimum [35]. In order to obtain the maximum power of engine output at Optimal-Spark-Advance-Angle (OSAA), the OSAA need to be adjusted continuously.

In the case of high speed and low load of engine, the OSAA should be increased. This is due to the delay period of fuel combustion in the cylinder. The faster the speed is, the bigger the OSSA. On the contrary, with the speed increased, the stronger the turbulence formed by combustion gas in the cylinder is, the faster the combustion speed and the lower the advance angle is. combustion turbulence, which positively a ffects the burning speed, will increase with faster crank speeds. In this case, the corresponding OSAA delay period will be small and because of the delay caused by fuel burning in the power stroke, the OSAA should be slightly smaller. Additionally, in this paper, the pressure identification based on the power flow is adopted. The comparison between identified pressure signal and actual pressure signal can be realized by computing each power. Hence the SAA is deemed as constant. Due to the interaction and method selected, it's deemed that the OSAA is constant at 15◦, and the maximum power output will occur when the crankshaft angle is 10◦ after the TDC where pressure in the cylinder reaches its peak value.

#### *3.2. Computation of Cylinder Pressure Peak Value*

After the burning process, the exhaust gases in cylinders is primarily composed of *CO*2, *H*2*O* and *<sup>N</sup>*2.According to Equations (7) and (8), the specific heat capacity at constant volume for the burned gas is

$$\mathcal{C}\_{\text{mix}} \approx \left( \frac{\frac{w\_{\text{C}}}{12}}{\frac{w\_{\text{C}}}{12} + \frac{w\_{\text{H}}}{8}} \mathcal{C}\_{\text{CO}\_2} + \frac{1}{2} \cdot \frac{\frac{w\_{\text{H}}}{8}}{\frac{w\_{\text{C}}}{12} + \frac{w\_{\text{H}}}{8}} \mathcal{C}\_{\text{H}\_2\text{O}} \right) \cdot \frac{2}{9} + \mathcal{C}\_{\text{N}\_2} \cdot \frac{7}{9} \tag{29}$$

where, *CCO*2 is the specific heat capacity at a constant volume for *CO*2, and *C H*2*O* is that for *H*2*O*.

As computed in Equation (30), Cylinder temperature varies with output power.

$$
\Delta T = \frac{\mathcal{Q}}{\mathcal{C}\_{\text{mix}} m} \tag{30}
$$

where, *m* is the mass of the burned gas, and according to conservation of mass theory:

$$m = m\_{air} + m\_{fuel} \tag{31}$$

According to the ideal gas state equation,

$$pV = nRT\tag{32}$$

Since the burning process is deemed as isovolumetric, the peak pressure in cylinder is written as:

$$p\_{peak} = \frac{nR\left(\frac{Q}{C\_{mix}\left(m\_{air} + m\_{fuel}\right)} + T\_{comp}\right)}{V} \tag{33}$$

where, *n* is the amount of substance which can be gained from the coe fficient of the fuel burning equation; *V* is the volume of cylinder at top dead center; *Tcomp* is the temperature at the end of the compression stroke; and *R* is the ideal gas constant.

#### **4. Modeling of in-Cylinder Pressure**

As previously mentioned, the cylinder pressure signal is periodic with variating frequency and amplitude, FAMFS is proposed to modulate the cylinder pressure [36]. A high-factorial Fourier series with variating frequency and amplitude is adopted to modulate the cylinder pressure.

## *4.1. Frequency-Modulated Fourier-Series*

If there is a periodic signal *f*(*x*), then for every *x* there exists a positive value *L* which makes *f*(*x*) = *f*(*x* + *L*) correct, where *L* is known as the fundamental period and ω0 = 2π *L* is termed the fundamental frequency. The Fourier-series is written as [37]:

$$f(\mathbf{x}) = a\_0 + \sum\_{n=1}^{\infty} \left( a\_n \cos nx \frac{2\pi}{L} + b\_n \sin nx \frac{2\pi}{L} \right) \tag{34}$$

where, *a*0 is the intercept and *an*, *bn* are coe fficients of the Fourier-series.

The frequency-modulated method [38] allows the instantaneous frequency of the carrier signal to vary with the change law of the delivery signal. This modulation can be classified as primary or secondary according to the category of the modulation e ffects.

With a defined carrier signal *xc*(*t*) = *A* cos(<sup>2</sup>π *fct*) and a transmitted signal *y*(*t*), the modulating signal can be written as:

$$\alpha\_{\varepsilon}(t) = A \cos(2\pi \int\_{0}^{t} [f\_{\varepsilon} + f\_{\Delta}y(\tau)]d\tau) \tag{35}$$

where, *fc* is the center frequency of the carrier signal; *A* is the amplitude; *fc* + *f*Δ*y*(τ)is the instantaneous frequency; and *f*Δ is the frequency deviation gain.

Based on the engineering practice, the 24-order FAMFS is selected. It should be noted that the identification process is completed o ffline, and the only online requirement is adjustment of the frequency and amplitude of the FAMFS. Thus, the instantaneity of method is well guaranteed. The 24th order FAMFS is shown below:

$$f(\mathbf{x}) = A a\_0 + A \sum\_{n=1}^{24} \left( a\_n \cos n(2\pi \sum\_{t\_i=0}^t \frac{t}{N\_n} y(t\_i)) + b\_n \sin n(2\pi \sum\_{t\_i=0}^t \frac{t}{N\_n} y(t\_i)) \right) \tag{36}$$

where, *t* is the total running time; *Nn* is the number of working cycles; and *ti* is time sequence.

#### *4.2. Pressure Model*

The optimal output of the Kalman observer is the speed of each power stroke. For a four-cylinder engine, the cylinders work in a logical sequence, known as firing order. Generally, the firing order is 1-3-4-2. By sampling the optimal output of the EKF, according to the working sequence, the optimal power-stroke speed for each cylinder is obtained:

$$\begin{array}{l} N\_{\text{man1}} = N\_{1+k\text{u}} \\ N\_{\text{man2}} = N\_{2+k\text{u}} \\ \cdot \\ \cdot \\ N\_{\text{numi}} = N\_{i+k\text{u}} \end{array} \tag{37}$$

where, *Nnumi* is the speed of the *ith* cylinder; *Ni*+*kn* is the optimal speed of (*i* + *kn*)*th* cycle; and *k* is a natural number.

According to actual working conditions, the cylinder pressure signal can be deemed as a FAMFS [39,40] signal with a center frequency at zero and a unit gain of frequency deviation. Thus, the pressure identification is ultimately written as:

$$F(\mathbf{x}) = A\_k a\_0 + P\_{atm} + A\_k \sum\_{n=1}^{24} \left( a\_n \cos n \left( 2\pi \sum\_{t\_i=0}^t \frac{t}{N\_n} y(t\_i) \right) + b\_n \sin n \left( 2\pi \sum\_{t\_i=0}^t \frac{t}{N\_n} y(t\_i) \right) \right) \tag{38}$$

$$y(t\_i) = \frac{\pi N\_k}{15} \qquad \sum\_{n=1}^{k} \frac{2\pi}{\omega\_n} \le t\_i \le \sum\_{n=1}^{k+1} \frac{2\pi}{\omega\_n} \tag{39}$$

$$A\_k = \frac{nR\left(\frac{Q\_k}{C\_{\text{mix}\parallel}(m\_{\text{air}\parallel k} + m\_{\text{f}\text{au}\parallel k})} + T\_{\text{comp}}\right)}{V} \quad \sum\_{n=1}^{k} \frac{2\pi}{\alpha\_n} \le t\_l \le \sum\_{n=1}^{k+1} \frac{2\pi}{\alpha\_n} \tag{40}$$

where, *t* is time, and *Ak* is the theoretical pressure peak at the *kth* working cycle.

#### **5. Validation and Results**

When the engine is in normal combustion state, its air-fuel ratio usually exceeds 0.85. A ratio of 0.7 or 0.8 will make the burning process of inside cylinder unstable and unpredictable, and this condition only happens in a very special condition like startup in an extremely cold environment or engine malfunctioned. No matter in which condition, it's not a normal working status so the prediction is not that useful. The validity of the proposed method is verified by comparing the identification value of the virtual cylinder pressure sensor with the measured value through the cylinder pressure sensor.

In the following section, data collected from a genuine 2.0-L four-stroke four-cylinder engine was provided by the FAW Group Corporation R&D Center (Changchun, China). The engine model is EA211 and produced by VW (Volkswagen) (Changchun, China). Data was collected over the course of 100 working cycles adopting Kistler Model 6052A Pressure Sensor to test the performance of the proposed method.

#### *5.1. Set Air-Fuel Ratio Coe*ffi*cient at 0.85*

The parameters and setup of the engine are shown in Figure 5 and Table 1.

**Figure 5.** Parameters and setup of Engine at A/F ratio 0.85.

**Table 1.** The parameters and setup of the engine.


Figure 6 compares the crankshaft speed outputs of the EKF with the measured results acquired from the genuine engine. As seen in Figures 6 and 7, from 30 cycles to 70 cycles the engine speed measured showed signs of rising and began to decline after reaching the maximum at 70 cycles. The EKF method accurately tracks this trend and its maximum percentage of error is only 9.6%. This shows acceptable agreemen<sup>t</sup> between identified and actual values.

0 10 20 30 40 50 60 70 80 90 100 -10 Cycle/times


**Figure 7.** Error percentage of speed tracking at A/F ratio 0.85.

Once the optimal speed from the EKF engine model is predicted, the peak pressure of cylinder can be calculated. Furthermore, the 24th order FAMFS is adopted for identification purposes, and its parameters are shown is Table 2:


**Table 2.** Parameters of the 24th order Fourier series.

A comparison of the actual cylinder pressure with the model-identified values is shown in Figure 8. It shows that in early stage of identification process, the identified values are basically in agreemen<sup>t</sup> with the measured values, whether in phase or amplitude. However, the phase deviation between measured and identified values begins to emerge and become larger with time. This is due to the cumulative error caused by the fitting error of the 24-order Fourier series and the speed tracking error.

For Figure 8, the area under the curve is multiplied by the piston area to represent the piston impulse. Hence, the cumulative error percentage of piston impulse is the cumulative error percentage of cylinder pressure. Figure 9 shows the curves of cumulative piston impulse including the actual and identification situation within 100 cycles. The two coincide basically. And Figure 10 shows that its cumulative error percentage is less than 1.8%. This agreemen<sup>t</sup> as well as the permitted errors, represents good performance from the proposed method.

**Figure 8.** Pressure identification at A/F ratio 0.85.

**Figure 9.** Piston impulse accumulated at A/F ratio 0.85.

**Figure 10.** Error percentage of piston impulse accumulated.

In general, the engine is a complex and sophisticated system, and small changes in valve or atmospheric temperature, fuel rail pressure, vibration from irregularity, etc., can severely impact the state of the engine; thus, besides a direct comparison of pressure and the cumulative piston impulse, comparison of piston impulse per cycle is also needed. As depicted in Figure 11, the actual and identified IMEP of each cycle one by one are contrasted. Figure 12 shows that these two error percentages are below 5.4%. And then the proposed method is proven to possesses a grea<sup>t</sup> ability to track the IMEP calculated by actual pressure.

**Figure 11.** IMEP at A/F ratio 0.85 per cycle.

**Figure 12.** Error percentage of IMEP at A/F ratio 0.85 per cycle.

#### *5.2. Set Air-Fuel Ratio Coe*ffi*cient at 0.95*

The parameters and setup of the engine are shown in Figure 13.

**Figure 13.** Parameters and setup of Engine at A/F ratio 0.95.

Figure 14 shows the comparing result of speed from EKF and a genuine engine at A/F ratio 0.95. Error percentages of the speed in Figure 15 is below 1.7%. Calculated and actual values have a grea<sup>t</sup> uniformity. Results show that proposed method possessed a grea<sup>t</sup> performance on tracking speed.

The peak pressure of cylinder can be calculated by using the predicted speed. Furthermore, the 24th order FAMFS is adopted for identification purposes, and its parameters are the same as Table 2. A comparison of the actual cylinder pressure with the model-identified values is shown in Figure 16.

**Figure 14.** Speed tracking at A/F ratio 0.95.

 **Figure 16.** Pressure identification at A/F ratio 0.95.

Figure 16 shows that the comparing result of actual and identified cylinder pressure. Figure 17 also shows the curves of cumulative piston impulse including the actual and identification situation but at A/F ratio 0.95 within 100 cycles. The two coincide also basically. As depicted in Figure 18, the actual and identified IMEP of each cycle one by one are contrasted. Compared with Figures 8 and 16, the proposed method possesses a better ability to track the cylinder pressure at A/F ratio 0.95 than at A/F ratio 0.85. The closer the air-fuel ratio coefficient is to the optimal A/F ratio, the more predictable the combustion state of the engine is. This is also proved by the cumulative error percentage of piston impulse is less than 1.4% in Figure 19, and the error percentage of IMEP in each cycle is less than 4.9% in Figure 20. On the contrary, the phase difference of identified cylinder pressure began to emerge and become larger as the air-fuel ratio coefficient getting far away from optimal Air-Fuel ratio coefficient. But no matter how the A/F ratio coefficient changed, if it is within reasonable range, error percentage of cylinder pressure is within tolerance.

**Figure 17.** Piston impulse accumulated at A/F ratio 0.95.

**Figure 18.** IMEP at A/F ratio 0.95 per cycle.

**Figure 19.** Error percentage of piston impulse accumulated.

**Figure 20.** Error percentage of IMEP at A/F ratio 0.95 per cycle.
