A novel method for in-cylinder pressure prediction using the engine instantaneous crankshaft torque

Abstract

Current practices of condition assessment in large marine engines are largely based on the measurement of cylinder pressure using external kits, which poses challenges due to sensors synchronisation and durability issues, as well as the inability to perform continuous monitoring. For addressing these challenges, this study aims at developing a novel method to solve the inverse problem of predicting the pressure variations in all engine cylinders, by using the Instantaneous Crankshaft Torque (ICT) measurement for large internal combustion engines. This method is developed by considering the Initial Value Problem (IVP) technique along with the integration of a direct crankshaft dynamics model incorporating the sensitivity parameters and stability criteria calculation based on the Lyapunov Exponent (LE) as well as a state-of-the-art Nonmonotone Self-Adaptive Levenberg-Marquardt (NSALMN) optimisation algorithm. The method is tested for a number of case studies using different combustion models based on the Weibe and sigmoid functions, as well as for healthy, degraded and faulty engine conditions. The derived results demonstrate adequate accuracy exhibiting a maximum error of 0.3% in the prediction of the mean peak in-cylinder pressure. The analysis of the calculated sensitivity parameters resulted in the identification of the parameters that significantly impact the solution, thus providing improved insights for selecting the developed method settings. The developed method renders the continuous and non-intrusive in-cylinder pressures monitoring feasible, by using a permanently installed shaft power metre sensor with higher sample rates.

Introduction

Large two-stroke diesel engines are by design the most efficient internal combustion engine type, ¹ and as a result they are employed in a variety of ship propulsion and power generation applications.^2,3 However, under actual operating conditions, degradation effects and various faults significantly reduce the engine efficiency and availability, whilst leading to increased emissions and operational costs.^4–6 With the advent of digitalisation, monitoring and analysing various engine parameters, such as engine power, fuel consumption, working medium temperatures, scavenge air pressure and in-cylinder pressures, have become widespread among engine operators.^7,8 In this way, better insights for the engine’s condition can be obtained, which prompts maintenance and repair actions in order to attain the engine optimal efficiency and reliability as well as to reduce the associated operating cost and gaseous emissions, with tangible benefits on the ship operations sustainability. ⁹

Amongst the most prominent indicators for the engine condition is the measurement of the in-cylinder pressure, ¹⁰ as it provides information to analyse the combustion process quality,^11,12 blowby, ¹³ cylinder misfire and knock,^14,15 cylinder compression ratio and injection timing, ¹⁶ as well as injector degradation and failures. ¹⁷ However, the methods currently employed in determining all of the above, rely on direct measurements of the in-cylinder pressure by deploying pressure sensors to either one or multiple cylinders, which can be expensive, suffer from poor durability,^18–20 exhibit synchronisation issues ²¹ and most often do not allow for continuous monitoring. Nonetheless, recent technological improvements have led to the development of continuous pressure monitoring systems which are primarily utilised in high-value power plants, such as those with dual-fuel engines. ²² However, for the majority of typical large internal combustion engines such as those employed in the 62,100 ocean-going vessels world wide, ²³ the costs involved with such monitoring systems remain outside the reach of most shipowners and operators. ²⁴ Moreover, having measured the in-cylinder pressure, the non-trivial task of signal correction is required, in order to remove the noise and perform averaging over several cycles, and most importantly to remove offsets by calibrating the acquired signal with respect to the engine crank angle ²⁵ and the cylinder pressure. ²⁶

As a result, methodologies have been developed utilising alternative engine measurements to predict important features that belong to, or can be directly derived from the in-cylinder pressure diagrams in a reliable and low-cost manner. In specific, for automotive industry applications, the peak pressure position has been predicted using accelerometers mounted on every cylinder,^27,28 as well as an ICT sensor. ¹⁸ Furthermore, properties such as the burned mass fraction and combustion phasing were estimated by utilising measurements from an ICT sensor ²⁹ and an Instantaneous Crankshaft Speed (ICS) sensor. ³⁰ Regarding large diesel engines, limited applications of such methodologies can be found in the pertinent literature. In particular, for a slow speed two-stroke diesel engine, the indicated power in each cylinder was estimated by employing ICT sensor measurements with a black-box ANN model. ³¹ For high and medium speed four-stroke diesel engines, ICS sensor measurements were used to detect individual cylinders misfiring, ³² and diagnose fuel pump faults through in-cylinder pressure differences, respectively. ³³

In addition, further advances have been made by developing methodologies, which attempt to reconstruct the entire in-cylinder pressure diagram. Regarding automotive applications, engine block vibration signals were utilised,^34–36 however the signal quality is very sensitive to the sensors’ placement, which is often decided arbitrarily. ³⁷ Inverse physical modelling has also been utilised along with an ICS sensor,^19,20 however the inversion of a purely physical (white-box) crankshaft dynamics model presented challenges due to singularities occurring at the piston Top Dead Centre (TDC). To resolve such inversion problems, a grey-box model was developed utilising the ICS along with the measured in-cylinder pressure in one of the engine cylinders, in order to estimate the in-cylinder pressure in the remaining engine cylinders. ³⁸ Regarding larger diesel engines, only one application has been identified in the literature, which attempts to reconstruct the in-cylinder pressure diagrams of all cylinders of a large two-stroke diesel engine. ³⁹ For this case, the in-cylinder pressure function was parametrised by using two parameters; the polytropic expansion coefficient and a pressure amplitude factor. Hence, it remains uncertain whether by this method, the in-cylinder pressure can be predicted accurately in case of existing engine faults or degradation. Subsequently, the parametrised in-cylinder pressure was used along with the ICT sensor measurements as part of a grey-box model, that resulted in a maximum error of 12% for the in-cylinder pressure prediction. However, a change of 5% in peak in-cylinder pressure is sufficient to cause a substantial impact in engine performance, ⁴⁰ hence a method with smaller error tolerances is required.

It is deduced from the preceding discussion that three peripheral measurements (the engine block vibrations, ICS and ICT) have been used to establish a robust relationship with the in-cylinder pressure. In specific, the engine block vibration signals carry substantial non-pressure related information, which significantly complicates the analysis. ⁴¹ On the other hand, the ICS measurement (raw signal from the angle encoder) requires post processing, a comparatively higher sampling rate compared with the ICT measurements to capture the instantaneous shaft rotational speed with adequate accuracy, and exhibits a low signal to noise ratio. ³¹ Hence, the ICT forms a better choice compared to the engine block vibration and ICS measurements, due to its higher signal to noise ratio as well as its linear response with engine load for the case of large two-stroke engines. ⁴²

Therefore, from the preceding literature review, the following gap has been identified: there exists no methodology for large two-stroke diesel engines, which can predict the in-cylinder pressure using peripheral measurements, firstly for the cases of degraded or faulty engine operation, and secondly with adequate accuracy for the case of healthy engine conditions. In this respect, the aim of this study is to develop a method that utilises the ICT measurement from a single torque sensor, to predict the in-cylinder pressures pressure variations for every engine condition including healthy, and various faulty or degraded conditions, thus potentially offering a reliable and low cost solution for a continuous pressure monitoring system. To accomplish this aim, this study employs the extension of a direct crankshaft dynamics model, the parametrisation of the in-cylinder pressure function, the development of an appropriate inverse problem method and its implementation to predict the in-cylinder pressure for a number of realistic engine operating conditions.

Methodology

The methodological steps followed in this study are listed below:

1. Description of the reference system and the direct crankshaft dynamics model. This model is employed in step 3 as part of the inverse problem method, and in step 6 to simulate the engine ICT for healthy, faulty and degraded conditions.

2. Determination of the parametrised in-cylinder pressure functions, and selection of parameters to be estimated by the inverse problem method.

3. Derivation of sensitivity equations, to be utilised for the development of the inverse problem method in step 5.

4. Development of the inverse problem method, which estimates the in-cylinder pressure parameters by utilising ICT measurements.

5. Implementation of the inverse problem method for a number of case studies for validation and testing under various engine operating conditions. This is performed by using the direct crankshaft dynamics model for the reference system, to simulate the ICT measurements for healthy, faulty and degraded engine conditions.

6. Discussion and analysis of results.

Reference system and direct crankshaft dynamics model

The reference system used in this study is a large two-stroke diesel engine driving a generator, which operates at a nominal speed of 124 r/min. The engine has 10 cylinders with an inline configuration, and a Maximum Continuous Rating (MCR) of 15.5 MW. For developing the crankshaft dynamics model of this system, 16 degrees of freedom are considered, which represent the 3 masses within the engine torsional damper, the 10 engine cylinders, the engine camshaft drive, the engine flywheel and the electric generator rotor, as shown in Figure 1.

Figure 1.

Schematic of the direct crankshaft dynamics model.

A direct crankshaft dynamics model is utilised in this study. This model has been developed, validated and employed to investigate the reference system response considering various friction submodels and solver settings in a previous authors’ study. ⁴³ In this study, this model is employed to firstly develop the inverse problem method in Section Inverse problem method development, and secondly to provide the required input of ICT at the flywheel for the case studies presented in Section Case studies design. The direct model is of the lumped parameter type consisting of 16 Degrees of Freedom (DOFs), and uses as primary input the in-cylinder pressure diagrams for each cylinder, and as secondary input the engine load, the DOFs inertia, damping and stiffness coefficients (obtained from the engine torsional vibration analysis), the geometrical and other data obtained from the manufacturer documentation, as well as the settings (solver and submodels) defined by the user, as shown in Figure 1. Furthermore, the direct model is based on a system of 16 non-linear second order differential equations, the solution of which is used to calculate the ICT at the engine flywheel according to equations (1) and (2), respectively. These equations will also be employed for the development of the inverse problem method.

$$\begin{array}{c}\mathit{J}(\mathit{\Theta })\stackrel{··}{\mathit{\Theta }}+\mathit{D}(\mathit{\Theta }){\stackrel{·}{\mathit{\Theta }}}^2+\mathit{C}\stackrel{·}{\mathit{\Theta }}\hfill \\ +\mathit{K}\mathit{\Theta }-\mathit{T}({\mathit{p}}_{\mathrm{cy}}(\hat{\mathit{\Theta }}),\mathit{\Theta },\stackrel{·}{\mathit{\Theta }})=0\hfill \end{array}$$

(1)

$$T_w(t)=k_{15}[{\theta }_{16}(t)-{\theta }_{15}(t)]$$

(2)

In equation (1), $\mathit{\Theta }=[{\theta }_1,\cdots ,{\theta }_{16}{]}^{\mathrm{T}}$ , represents an $1\times 16$ crankshaft angular displacement vector, whereas $\stackrel{·}{\mathit{\Theta }}$ , $\stackrel{··}{\mathit{\Theta }}$ denote the angular velocity and acceleration vectors, respectively. In addition, $\mathit{C}$ and $\mathit{K}$ are $16\times 16$ tridiagonal structural damping and stiffness matrices respectively, whilst $\mathit{J}(\mathit{\Theta })$ and $\mathit{D}(\mathit{\Theta })$ are the variable inertia and variable damping matrices, respectively. In particular, these variable matrices map the vector $\mathit{\Theta }$ on to a $16\times 16$ diagonal matrix, such that $\mathit{J}(\mathit{\Theta })=\mathrm{diag}(J({\theta }_1),\cdots ,J({\theta }_{16}))$ and $\mathit{D}(\mathit{\Theta })=\mathrm{diag}(D({\theta }_1),\cdots ,D({\theta }_{16}))$ . Finally, $\mathit{T}({\mathit{p}}_{\mathrm{cy}}(\hat{\mathit{\Theta }}),\mathit{\Theta },\stackrel{·}{\mathit{\Theta }})$ is the $1\times 16$ forcing torque vector, which includes the torque induced by the in-cylinder pressure, as well as the engine friction and load torques. In specific, the forcing torque vector is a function of the angular displacement and its derivative, as well as the in-cylinder pressure ${\mathit{p}}_{\mathrm{cy}}(\hat{\mathit{\Theta }})$ , where $\hat{\mathit{\Theta }}$ is a subvector of $\mathit{\Theta }$ . This subvector contains the displacement of DOFs 4–8 and 10–14 which correspond to the engine cylinders, such that ${\mathit{p}}_{\mathrm{cy}}(\hat{\mathit{\Theta }})=[p_{\mathrm{cy}}({\theta }_4),\cdots ,p_{\mathrm{cy}}({\theta }_8),p_{\mathrm{cy}}({\theta }_{10}),\cdots ,p_{\mathrm{cy}}({\theta }_{14}){]}^{\mathrm{T}}$ . In equation (2), $k_{15}$ denotes the stiffness coefficient at the 15th DOF, whilst ${\theta }_{15}(t)$ and ${\theta }_{16}(t)$ represent the angular displacements as a function of time t of the 15th and 16th DOFs respectively, which are obtained from the solution of equation (1). A detailed derivation of the direct crankshaft dynamics model governing equations can be found in Tsitsilonis et al. ⁴³

In-cylinder pressure parametrisation

The parametrisation of the in-cylinder pressure curve for the closed part of the diesel cycle is a critical task. Therefore, the cylinder pressure curve should be adequately parametrised, such that faults or degradation is captured by changing the employed parameters. Moreover, the parametrised in-cylinder pressure function should employ as few parameters as possible, in order to reduce the mathematical and computational complexity. The engine in-cylinder pressure curve can be parametrised by using semi-empirical equations, which have demonstrated high accuracy and low computational complexity.^44–46 In specific, this is accomplished by using different equations to represent the three main in-cylinder processes (compression, combustion and expansion). The compression and expansion processes are modelled using polytropic relationships according to equations (3) and (4), whilst the combustion process is modelled by using a function that interpolates between the other two processes.^18,46 Hence, by considering the interpolating function, the analytic in-cylinder pressure function for the closed cycle is represented by equation (5).

$$p_{\mathrm{co}}(\theta )=p_{\mathrm{EVC}}{\left(\frac{V_{\mathrm{EVC}}}{V(\theta )}\right)}^{{\gamma }_{\mathrm{co}}}$$

(3)

$$p_{\mathrm{ex}}(\theta )=p_{\mathrm{me}}{\left(\frac{V_c}{V(\theta )}\right)}^{{\gamma }_{\mathrm{ex}}}$$

(4)

$$p_{\mathrm{cy}}(\theta )=(1-f_i(\theta ))p_{\mathrm{co}}(\theta )+f_i(\theta )p_{\mathrm{ex}}(\theta )$$

(5)

The parameter $p_{\mathrm{EVC}}$ denotes the pressure at the instant the exhaust valve closes, which is assumed equivalent to the scavenge air pressure, $p_{\mathrm{me}}$ denotes the maximum isentropic expansion pressure, ⁴⁷ whereas ${\gamma }_{\mathrm{co}}$ , ${\gamma }_{\mathrm{ex}}$ denote the compression and expansion polytropic coefficients. Moreover, $V_{\mathrm{EVC}}$ is the cylinder volume when the piston is at the exhaust valve close position, $V_c$ is the cylinder clearance volume and $V(\theta )$ is the cylinder volume as a function of the crank angle. Finally, $f_i(\theta )$ denotes the interpolating function.

For a diesel engine, there are two types of interpolating functions which can be used in equation (5). The most simplistic one is the sigmoid function, ⁴⁸ since it depends on two parameters according to equation (6). The second one is the single Weibe function, which is extensively used in combustion modelling applications,^46,47,49 and is represented by equation (7).

$$f_{i,s}(\theta )={(1+\exp \{k_1(\theta -k_2)\})}^{-1}$$

(6)

$$f_{i,w}(\theta )=\{\begin{array}{c}1-\exp \{-a{\left(\frac{\theta -{\theta }_0}{\mathrm{\Delta }\theta }\right)}^{m+1}\}\mathrm{for}\theta >{\theta }_0\hfill \\ \multicolumn{1}{c}{0\mathrm{for}\theta \le {\theta }_0}\end{array}$$

(7)

For the sigmoid function, $k_1$ is a dimensionless parameter, which affects the rate of heat release as shown in Figure 2(a), whereas $k_2$ denotes the crank angle where 50% of heat release is achieved. For the Weibe function, the parameters ${\theta }_0$ and $\mathrm{\Delta }\theta$ represent the start of combustion and combustion duration respectively, whilst m is the Weibe parameter which controls the overall shape of the function as shown in Figure 2(b). Since the Weibe function is an over-parametrised function, the parameter a is always kept constant with a value of 6.9078 as this represents a combustion efficiency of 99.9%. ⁴⁹

Figure 2.

(a) Sigmoid functions with different $k_1$ parameter values and (b) Weibe functions with different m parameter values.

The performance of the inverse problem method is evaluated for both interpolating functions as described in Section Case studies design.

Having obtained the equations (3)–(7), each of their included variables can be classified as known or unknown parameters, as summarised in Table 1. This information will be employed in Sections Sensitivity and stability analysis and Inverse problem method development, so that the sensitivity equations and subsequently the inverse problem method can be formulated accordingly, in order to estimate the values of the unknown parameters. The parameters classified as known are the ones, which are either known a priori, can be obtained from basic manufacturer documentation, or can be easily measured on site using existing permanently installed sensors. In specific, these parameters belong to equations (3) and (4) and include the in-cylinder pressure at the exhaust valve close position ( $p_{\mathrm{EVC}}$ ), approximated by the scavenge air receiver pressure, which is always monitored in larger power plants with sufficient accuracy.^39,50 In addition, parameters $V_{\mathrm{EVC}}$ and $V_c$ can be calculated using information from the manufacturer documentation. ⁵¹ Finally, the polytropic coefficients ${\gamma }_{\mathrm{co}}$ and ${\gamma }_{\mathrm{ex}}$ can be estimated with sufficient accuracy a priori, by using air ideal gas tables, ⁵² as they do not exhibit significant changes under faulty or degraded engine conditions. ⁴⁵

Table 1.

Pressure function parameters classification.

Equation	Parameter	Classification
Compression (equation (3))	$p_{\mathrm{EVC}}$	K
	$V_{\mathrm{EVC}}$	K
	${\gamma }_{\mathrm{co}}$	K
Expansion (equation (4))	$p_{\mathrm{ex}}$	U
	$V_c$	K
	${\gamma }_{\mathrm{ex}}$	K
Sigmoid (equation (6))	$k_1$	U
	$k_2$	U
Weibe (equation (7))	${\theta }_0$	U
	$\mathrm{\Delta }\theta$	U
	m	U

K: known parameter; U: unknown parameter.

The parameters classified as unknown are the ones that their values will be estimated in this study, by using the ICT measurement. These parameters greatly effect the shape of the in-cylinder pressure diagram and can only be estimated using the in-cylinder pressure measurements, which are associated with inherit drawbacks, as explained in Section Introduction. In specific, the unknown parameters are, the maximum expansion pressure ( $p_{\mathrm{me}}$ ) used in equation (4), the parameters $k_1$ and $k_2$ of the sigmoid function in equation (6), and the start of combustion ( ${\theta }_0$ ), combustion duration ( $\mathrm{\Delta }\theta$ ) and Weibe parameter (m) of the Weibe function in equation (7).

The parametrised in-cylinder pressure function (equation (5)) can be incorporated to the direct model governing equations (equation (1)). This is firstly performed by grouping the unknown parameters for the n th cylinder into a parameter vector ${\mathit{\beta }}^{(n)}$ , and subsequently concatenating the ${\mathit{\beta }}^{(n)}$ vectors into a larger unknown parameters vector $\mathit{B}$ . These parameter vectors can then be included in the in-cylinder pressure vector found in the forcing torque term of equation (1). This is demonstrated in equations (8)–(10) for the case of the sigmoid function parameters, given the reference system is a 10-cylinder engine and the engine cylinders are represented by DOFs 4–8 and 10–14.

$${\mathit{\beta }}^{(n)}={\left[\begin{array}{ccc}{p}_{\mathrm{me}}^{(n)}\hfill & k_1^{(n)}\hfill & k_2^{(n)}\hfill \end{array}\right]}^{\mathrm{T}}$$

(8)

$$\mathit{B}={\left[\begin{array}{ccccc}{\mathit{\beta }}^{(1)}\hfill & \cdots \hfill & {\mathit{\beta }}^{(n)}\hfill & \cdots \hfill & {\mathit{\beta }}^{(10)}\hfill \end{array}\right]}^{\mathrm{T}}$$

(9)

$${\mathit{p}}_{\mathrm{cy}}(\hat{\mathit{\Theta }},\mathit{B})={\left[\begin{array}{ccc}{p}_{\mathrm{cy}}({\theta }_4,{\mathit{\beta }}^{(1)})\hfill & \cdots \hfill & p_{\mathrm{cy}}({\theta }_{14},{\mathit{\beta }}^{(10)})\hfill \end{array}\right]}^{\mathrm{T}}$$

(10)

Subsequently, the direct model governing equations have become functions of the unknown parameters vector $\mathit{B}$ according to equations (11) and (12).

$$\begin{array}{c}\mathit{J}(\mathit{\Theta })\stackrel{··}{\mathit{\Theta }}+\mathit{D}(\mathit{\Theta }){\stackrel{·}{\mathit{\Theta }}}^2+\mathit{C}\stackrel{·}{\mathit{\Theta }}\hfill \\ +\mathit{K}\mathit{\Theta }-\mathit{T}({\mathit{p}}_{\mathrm{cy}}(\hat{\mathit{\Theta }},\mathit{B}),\mathit{\Theta },\stackrel{·}{\mathit{\Theta }})=0\hfill \end{array}$$

(11)

$$T_w(t,\mathit{B})=k_{15}({\theta }_{16}(t,\mathit{B})-{\theta }_{15}(t,\mathit{B}))$$

(12)

It should be noted that the angular displacement, velocity and acceleration vectors are functions of time and $\mathit{B}$ , such that $\mathit{\Theta }=\mathit{\Theta }(t,\mathit{B})$ , $\stackrel{·}{\mathit{\Theta }}=\stackrel{·}{\mathit{\Theta }}(t,\mathit{B})$ , $\stackrel{··}{\mathit{\Theta }}=\stackrel{··}{\mathit{\Theta }}(t,\mathit{B})$ , since the model response depends on the values of the unknown parameters vector. However, for notational convenience the functional dependencies of the above are not explicitly stated in equation (11).

Sensitivity and stability analysis

The sensitivity equations are employed in order to determine the rate of change of a differential equation’s solution, with respect to specific parameters of that differential equation. ⁵³ In specific, the sensitivity equations are useful when, as it is often the case, a closed-form expression of the differential equation response is not available, hence symbolic differentiation cannot be directly performed. In this study, the sensitivity equations are used in order to calculate the objective function gradient, and to perform the stability analysis, which warns the user of possible chaotic behaviour when implemented within the inverse problem method, as described in Section Inverse problem method development .

Thus, in order to aid the derivation of the aforementioned sensitivities, the generic form of the sensitivity equations should be obtained first given an arbitrary differential equation, in accordance to equations (13)–(15): ⁵³

$$\stackrel{·}{y}=h(x,y,\mathit{c})$$

(13)

$$\frac{\partial \stackrel{·}{y}}{\partial \mathit{c}}=\frac{\partial }{\partial \mathit{c}}h(x,y,\mathit{c})$$

(14)

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial y}{\partial \mathit{c}}={\mathbf{J}}_y\frac{\partial y}{\partial \mathit{c}}+\frac{\partial h}{\partial \mathit{c}}\mathrm{with}\frac{\partial y}{\partial \mathit{c}}{|}_{t=0}=0$$

(15)

Equation (15) is a first order linear differential equation with variable coefficients, which is solved for the dependent variable $\partial y/\partial \mathit{c}$ . The dependent variable therefore denotes the sensitivity of the arbitrary differential equation $\stackrel{·}{y}=h(x,y,\mathit{c})$ , with respect to a parameter vector $\mathit{c}$ . Moreover, the term ${\mathbf{J}}_y$ denotes the Jacobian of equation h with respect to its dependent variable y. Similarly to equation (15), the sensitivity with respect to the initial conditions $y_0$ can also be obtained by treating the initial conditions as another parameters vector, according to the following equation: ⁵³

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial y}{\partial y_0}={\mathbf{J}}_y\frac{\partial y}{\partial y_0}\mathrm{with}\frac{\partial y}{\partial y_0}{|}_{t=0}=\mathit{I}$$

(16)

Following the above, the sensitivity equations for this study will be firstly derived with respect to the unknown parameters vector $\mathit{B}$ . The derived sensitivity equations will be used to calculate the gradient of the objective function as described in Section Inverse problem method implementation. By appropriately differentiating equation (11), the following equations are derived:

$$\begin{array}{c}\mathit{J}\frac{\partial \stackrel{··}{\mathit{\Theta }}}{\partial \mathit{B}}+(2\mathit{D}+\mathit{C}-\frac{\partial \mathit{T}}{\partial \stackrel{·}{\mathit{\Theta }}})\frac{\partial \stackrel{·}{\mathit{\Theta }}}{\partial \mathit{B}}-\frac{\partial \mathit{T}}{\partial \mathit{B}}\hfill \\ +(\frac{\partial \mathit{J}}{\partial \mathit{\Theta }}\stackrel{··}{\mathit{\Theta }}+\frac{\partial \mathit{D}}{\partial \mathit{\Theta }}{\stackrel{·}{\mathit{\Theta }}}^2-\frac{\partial \mathit{T}}{\partial \mathit{\Theta }}+\mathit{K})\frac{\partial \mathit{\Theta }}{\partial \mathit{B}}=0\hfill \end{array}$$

(17)

$$\begin{array}{c}\mathit{J}\frac{\mathrm{d}}{\mathrm{d}{t}^2}\frac{\partial \mathit{\Theta }}{\partial \mathit{B}}+{\mathbf{J}}_{\stackrel{·}{\mathit{\Theta }}}\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial \mathit{\Theta }}{\partial \mathit{B}}+{\mathbf{J}}_{\mathit{\Theta }}\frac{\partial \mathit{\Theta }}{\partial \mathit{B}}-\frac{\partial \mathit{T}}{\partial \mathit{B}}=0\hfill \\ \mathrm{with}\frac{\partial \mathit{\Theta }}{\partial \mathit{B}}{|}_{t=0}=0\hfill \end{array}$$

(18)

Equation (18) is a linear second order differential equation with variable coefficients, where ${\mathbf{J}}_{\mathit{\Theta }}$ and ${\mathbf{J}}_{\stackrel{·}{\mathit{\Theta }}}$ are the Jacobian matrices of equation (11) with respect to $\mathit{\Theta }$ and $\stackrel{·}{\mathit{\Theta }}$ , respectively. This differential equation is integrated for the dependent variable $\partial \mathit{\Theta }/\partial \mathit{B}$ , which represents the angular displacement sensitivity with respect to the unknown parameters, and is defined according to the following equation:

$$\frac{\partial \mathit{\Theta }}{\partial \mathit{B}}=\left[\begin{array}{c}\frac{\partial {\theta }_1}{\partial \mathit{B}}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\frac{\partial {\theta }_{16}}{\partial \mathit{B}}}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial {\theta }_1}{\partial {\beta }_1}\hfill & \cdots \hfill & \frac{\partial {\theta }_1}{\partial {\beta }_N}\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\frac{\partial {\theta }_{16}}{\partial {\beta }_1}}& \cdots \hfill & \frac{\partial {\theta }_{16}}{\partial {\beta }_N}\hfill \end{array}\right]$$

As a result, the sensitivity of the flywheel torque $T_w$ with respect to the unknown parameters vector $\mathit{B}$ , can be obtained by employing the appropriate elements from the $\partial \mathit{\Theta }/\partial \mathit{B}$ vector, and by differentiating equation (12), according to the following equation:

$$\frac{\partial T_w}{\partial \mathit{B}}=k_{15}(\frac{\partial {\theta }_{16}}{\partial \mathit{B}}-\frac{\partial {\theta }_{15}}{\partial \mathit{B}})$$

(19)

To analyse the stability of the solution from equation (11), the local sensitivity with respect to the direct model initial conditions ${\mathit{\Theta }}_0$ and ${\stackrel{·}{\mathit{\Theta }}}_0$ needs to be derived in a similar way to equation (16). In specific, this is performed by treating the initial conditions as another parameter vector and differentiating appropriately. Subsequently, the initial conditions have to be combined in one vector ${\mathbf{\Phi }}_0=[{\mathit{\Theta }}_0^{\mathrm{T}},{\stackrel{·}{\mathit{\Theta }}}_0^{\mathrm{T}}{]}^{\mathrm{T}}$ as shown in equation (20), and the same differentiation procedure is employed as with equation (18).

$$\begin{array}{c}\mathit{J}\frac{\mathrm{d}}{\mathrm{d}{t}^2}\frac{\partial \mathit{\Theta }}{\partial {\mathbf{\Phi }}_0}+{\mathbf{J}}_{\stackrel{·}{\mathit{\Theta }}}\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial \mathit{\Theta }}{\partial {\mathbf{\Phi }}_0}+{\mathbf{J}}_{1}Q\frac{\partial \mathit{\Theta }}{\partial {\mathbf{\Phi }}_0}=0\hfill \\ \mathrm{with}\frac{\partial \mathit{\Theta }}{\partial {\mathbf{\Phi }}_0}{|}_{t=0}=\mathit{I}\hfill \end{array}$$

(20)

By utilising the solution of equation (20), the stability of the system of equations can be quantified by using the Lyapunov Exponents (LE). ⁵⁴ In specific, the LE is provided in accordance to equation (21), and it quantifies the average exponential rate of divergence of nearby solutions from a given current solution, which is obtained provided a set of initial conditions.

$$\psi ({\mathbf{\Phi }}_0)=\underset{n\to \infty }{lim}\frac{1}{n}\sum _{i=0}^{n-1}\ln \left|\frac{\partial \mathit{\Theta }(t_n)}{\partial {\mathbf{\Phi }}_0}\right|$$

(21)

Therefore, the case where the LE exists and is positive implies that nearby solutions will diverge from the current solution if a small change in the initial conditions is introduced. Consequently, there exists a very sensitive dependence of the differential equations solution with respect to the initial conditions. In contrast, negative values of the LE imply that changes in initial conditions will eventually converge to the current solution given a sufficiently large length of time. Hence, the solution is deemed insensitive to changes in the initial conditions. In this case, a smaller value of the LE denotes a faster convergence to the current solution.

Extended direct crankshaft dynamics model

Having derived the sensitivity equations in the previous section, it is observed that the values for their variable coefficients depend on the solution of the direct model governing equations. Therefore, it is common practice when employing the sensitivity equations to solve them simultaneously with the original model equations, in order to enhance the computational efficiency of the solution procedure. ⁵⁵ As a result, the governing and the sensitivity equations can be simultaneously solved as part of an extended version of the direct crankshaft dynamics model, hereinafter referred to as the extended model.

In particular, the direct model shown in Figure 1 is extended such that it can provide as output the following two sensitivities: the flywheel torque with respect to the unknown parameter vector (equation (19)), and the solution of equation (11) with respect to the initial conditions (equation (20)). The former sensitivity is used for the objective function gradient calculation as described in Section Inverse problem method development, whereas the latter is used to perform the stability analysis as described in Section Sensitivity stability analysis. The extended model is therefore configured to calculate both sensitivities as shown in Table 2, so either one or both can be derived as output by the user as required. Moreover, it should be noted that the model input remains the same for the extended model as it was for the direct model.

Table 2.

Direct and extended model output; extended model input remains the same as the direct model.

Model output	Direct model	Extended model
$T_w$ (equation (2))	×	×
$\partial T_w/\partial \mathit{B}$ (equation (19))		×
$\partial \mathit{\Theta }/\partial {\mathbf{\Phi }}_0$ (equation (20))		×

The extension of the direct model is accomplished by firstly transforming the governing equations and the two sets of sensitivity equations into a pseudo-explicit form. This implies that the aforementioned equations are transformed, such that their left hand side is a vector that consists of the first and second dependent variables’ time derivatives. For this, the vectorisation operator $\mathrm{vec}(·)$ needs to be utilised, which transforms the matrices of the dependent variables from the sensitivity differential equations (equations (18)–(20)) to vectors. Therefore, by using the subscript notation to denote the partial derivatives with respect to $\mathit{\Theta }$ and by vectorising aforementioned matrices, it follows that:

$$\mathrm{vec}({\mathit{\Theta }}_{\mathit{B}})={\left[\begin{array}{c}{\left(\frac{\partial {\theta }_1}{\partial \mathit{B}}\right)}^{\mathrm{T}},\cdots ,{\left(\frac{\partial {\theta }_{16}}{\partial \mathit{B}}\right)}^{\mathrm{T}}\hfill \end{array}\right]}^{\mathrm{T}},$$

$$\mathrm{vec}({\mathit{\Theta }}_{{\mathbf{\Phi }}_0})={\left[\begin{array}{c}{\left(\frac{\partial {\theta }_1}{\partial {\mathbf{\Phi }}_0}\right)}^{\mathrm{T}},\cdots ,{\left(\frac{\partial {\theta }_{16}}{\partial {\mathbf{\Phi }}_0}\right)}^{\mathrm{T}}\hfill \end{array}\right]}^{\mathrm{T}}$$

Since the time derivative is a linear operator where $\frac{d}{\mathrm{dt}}\mathrm{vec}(·)=\mathrm{vec}(d(·)/\mathrm{dt})$ , the governing equations and the two sets of sensitivity equations can be reformulated as follows:

$$\left[\begin{array}{c}\stackrel{·}{\mathit{\Theta }}\hfill \\ \multicolumn{1}{c}{\stackrel{··}{\mathit{\Theta }}}\end{array}\right]=\mathbf{F}(\mathit{\Theta },\stackrel{·}{\mathit{\Theta }},\mathit{B}),$$

$$\left[\begin{array}{c}\mathrm{vec}({\stackrel{·}{\mathit{\Theta }}}_{\mathit{B}})\hfill \\ \multicolumn{1}{c}{\mathrm{vec}({\stackrel{··}{\mathit{\Theta }}}_{\mathit{B}})}\end{array}\right]={\mathbf{S}}_1({\mathit{\Theta }}_{\mathit{B}},{\stackrel{·}{\mathit{\Theta }}}_{\mathit{B}},\mathit{B}),$$

$$\left[\begin{array}{c}\mathrm{vec}({\stackrel{·}{\mathit{\Theta }}}_{{\mathbf{\Phi }}_0})\hfill \\ \multicolumn{1}{c}{\mathrm{vec}({\stackrel{··}{\mathit{\Theta }}}_{{\mathbf{\Phi }}_0})}\end{array}\right]={\mathbf{S}}_2({\mathit{\Theta }}_{{\mathbf{\Phi }}_0},{\stackrel{·}{\mathit{\Theta }}}_{{\mathbf{\Phi }}_0},{\mathbf{\Phi }}_0)$$

where system of equations $\mathbf{F}$ , ${\mathbf{S}}_1$ and ${\mathbf{S}}_2$ refer to the pseudo-explicit form of equations (11), (18) and (20) respectively. Hence, the aforementioned systems of equations can be combined to one larger system according to equation (22). This larger system of equations is therefore solved using the a piecewise Linear Time-Invariant (LTI) system solver, which is the default solver for the direct model. ⁴³

$$\left[\begin{array}{c}\stackrel{·}{\mathit{\Theta }}\hfill \\ \multicolumn{1}{c}{\stackrel{··}{\mathit{\Theta }}}\\ \multicolumn{1}{c}{\mathrm{vec}({\stackrel{·}{\mathit{\Theta }}}_{\mathit{B}})}\\ \multicolumn{1}{c}{\mathrm{vec}({\stackrel{··}{\mathit{\Theta }}}_{\mathit{B}})}\\ \multicolumn{1}{c}{\mathrm{vec}({\stackrel{·}{\mathit{\Theta }}}_{{\mathbf{\Phi }}_0})}\\ \multicolumn{1}{c}{\mathrm{vec}({\stackrel{··}{\mathit{\Theta }}}_{{\mathbf{\Phi }}_0})}\end{array}\right]=\left[\begin{array}{c}\mathbf{F}(\mathit{\Theta },\stackrel{·}{\mathit{\Theta }},\mathit{B})\hfill \\ \multicolumn{1}{c}{{\mathbf{S}}_1({\mathit{\Theta }}_{\mathit{B}},{\stackrel{·}{\mathit{\Theta }}}_{\mathit{B}},\mathit{B})}\\ \multicolumn{1}{c}{{\mathbf{S}}_2({\mathit{\Theta }}_{{\mathbf{\Phi }}_0},{\stackrel{·}{\mathit{\Theta }}}_{{\mathbf{\Phi }}_0},{\mathbf{\Phi }}_0)}\end{array}\right]$$

(22)

Inverse problem method development

Optimisation algorithm selection

The inverse problem method is the means of accomplishing the aim of this study, that is predicting the in-cylinder pressures based on the ICT measurement. The development of this method is based on the selection of an optimisation algorithm, and subsequently, it is implemented by combining this algorithm with the extended crankshaft dynamics model, as described further below in this Section. In specific, as explained in Section Reference system the direct model, the in-cylinder pressure is used as the primary input to the direct crankshaft dynamics model to predict the engine ICT. In the inverse problem method, the measured ICT from a torque metre installed at the engine flywheel will be utilised as the input, and the calculated unknown parameters of the in-cylinder pressure function listed in Table 1 will be provided as output. In principle, this is a parameter estimation problem of the governing differential equations employed in the direct crankshaft dynamics model.

Therefore, the in-cylinder pressure parameters estimation can be treated as an IVP, where the differential equations are repeatedly solved by employing an optimisation algorithm, which varies the values of their parameters in order to match the measured data. ⁵³ The advantages of this IVP technique include its computational simplicity,^56,57 and that the highly complex inverse structure of the problem (which may often be ill-posed) is completely bypassed. ⁵⁸ However, in order for this technique to work effectively, a very accurate estimate of the initial parameters needs to be provided, and the differential equations should have small sensitivity to changes in their initial conditions. ⁵⁹ The latter will be examined by conducting the direct model’s stability analysis as described in Section Sensitivity and stability analysis. However, the former does not pose a challenge for this study, as a relatively accurate estimate for the initial values of the unknown parameters vector can be obtained. This can be done by means of least squares fitting of equation (5) to the in-cylinder pressure diagrams, ¹⁵ which can be obtained from the manufacturer documentation (e.g. engine shop tests) and are available for a number of engine operating points.

The optimisation algorithm is a basic part of the IVP technique, hence its selection is of critical importance. Overall, optimisation algorithms can be categorised as derivative-free, stochastic and gradient-based. ⁶⁰ Derivative-free algorithms are simplistic, yet their convergence is time consuming, and they lack robustness for more complex nonlinear objective functions. ⁶¹ Stochastic algorithms can handle extremely complex objective functions with multiple minima and discontinuities. However they are computationally complex and require a very large number of function evaluations, which results to an increased execution time. ⁶¹

Gradient based algorithms and particularly Levenberg-Marquardt (LM) algorithms, are extensively used for nonlinear parameter identification problems due to their robustness and quick convergence near the minimum. ⁶⁰ Their application would be especially advantageous for the case of this study, as the analytic gradient of the direct model can be obtained, hence improving the execution time and the solution robustness compared to the use of a numerically estimated gradient. In this respect, this study employs the NSALMN algorithm proposed by Kimiaei, ⁶² which is an especially robust version of the LM algorithm.

Inverse problem method implementation

In order to provide the necessary foundation to describe the implementation of the NSALMN optimisation algorithm, a description of the objective function and the algorithm steps should be provided first. In particular, the NSALMN algorithm employs an iterative procedure to minimise the least squares objective function, which depends on the unknown parameters vector $\mathit{B}$ , in accordance to the following equation:

$$O(\mathit{B})=\frac{1}{2}||{\mathit{T}}_w(\mathit{B})-{\tilde{\mathit{T}}}_w|{|}_2^2$$

(23)

where $||·|{|}_2$ is the Euclidean norm, ${\mathit{T}}_w(\mathit{B})$ is the time vector which consists of the predicted torque at the flywheel obtained from the extended model such that ${\mathit{T}}_w(\mathit{B})=[T_w(t_1,\mathit{B}),\cdots T_w(t_N,\mathit{B}){]}^{\mathrm{T}}$ , and ${\tilde{\mathit{T}}}_w$ is the time vector which consists of the measured ICT data recorded by a torque metre at the flywheel within one engine revolution. The basis of the algorithm is the LM step, which makes a linear approximation of the objective function at the k th iteration in accordance to the following equations:

$$\nabla O_k={\mathbf{J}}_{{\mathit{B}}_k}^{\mathrm{T}}({\mathit{T}}_w({\mathit{B}}_k)-{\tilde{\mathit{T}}}_w)$$

(24)

$$\mathrm{\Delta }{\mathit{B}}_k={({\mathbf{J}}_{{\mathit{B}}_k}^{\mathrm{T}}{\mathbf{J}}_{{\mathit{B}}_k}+{\lambda }_k\mathbf{I})}^{-1}\nabla O_k$$

(25)

where $\nabla O_k$ is the objective function gradient, $\mathbf{I}$ is the identity matrix and ${\lambda }_k$ is a the LM damping parameter. This parameter indirectly controls the LM step size and direction, the calculation of which is specific to the NSALMN algorithm and is further described below in this Section. Furthermore, ${\mathbf{J}}_{{\mathit{B}}_k}$ is the sensitivity matrix of the predicted flywheel torque with respect to the unknown parameters vector, which is obtained from the extended model for every time t within one engine revolution such that:

$${\mathbf{J}}_{{\mathit{B}}_k}={\left[\begin{array}{c}\frac{\partial T_w(t_1,{\mathit{B}}_k)}{\partial {\mathit{B}}_k},\cdots ,\frac{\partial T_w(t_N,{\mathit{B}}_k)}{\partial {\mathit{B}}_k}\hfill \end{array}\right]}^{\mathrm{T}}$$

The developed method for the solution of the inverse crankshaft dynamics problem employs the NSALMN algorithm. The flowchart of this method that includes the method steps and order of operations is depicted in Figure 3. The method is initiated by providing the following input: (a) an accurate estimate for the initial values of the unknown parameters vector ${\mathit{B}}_0$ , which can be found as mentioned in the second paragraph of Section Inverse problem method implementation; (b) the the initial LM damping parameter value ${\lambda }_0$ ; (c) the iteration counter k is set to zero; (d) the measured ICT data ${\tilde{\mathit{T}}}_w$ ; (e) the secondary input as indicated in Figure 1 (used to calculate the objective function, and as input to the extended model respectively); (f) the optimisation algorithm stopping condition $\epsilon$ , which pertains to the tolerance in the objective function gradient where the minimum is identified; (g) the optimisation algorithm step acceptance criterion $\mu$ .

Figure 3.

Flowchart of the developed method to solve the inverse problem.

The method then proceeds to the calculation of the in-cylinders pressure variations (for one engine revolution) by employing equation (5), based on the unknown parameters vector $\mathit{B}$ (its initial value is used in the first iteration whereas the value calculated in the last successful iteration is used in every other iteration). Subsequently, by employing the in-cylinder pressure curves and the secondary input, the extended model provides as output the flywheel ICT ( ${\mathit{T}}_w$ ), and the sensitivity matrix of the flywheel torque with respect to the unknown parameters vector ( ${\mathbf{J}}_{{\mathit{B}}_k}$ ). Having available ${\tilde{\mathit{T}}}_w$ , ${\mathit{T}}_w$ and ${\mathbf{J}}_{{\mathit{B}}_k}$ , the optimisation algorithm objective function and its gradient can be evaluated in accordance to equations (23) and (24), respectively. The values of the objective function and its gradient are then compared against the optimisation algorithm stopping conditions. In case where the algorithm stopping criteria are accomplished (true conditions), the optimised values of the unknown parameters ${\mathit{B}}^{*}$ are provided as output. Otherwise (false conditions), the optimisation algorithm proceeds by calculating the trial step size using the ${\lambda }_k$ values provided from the last successful iteration, or set by the user if its the first iteration.

In tandem with the calculation of the objective function and its gradient, the stability analysis is also performed. It should be noted that this is not an original part of the NSALMN algorithm, however it is incorporated as part of the developed inverse problem method. In specific, by using the sensitivity equations calculated from the extended model and equation (21), all of the Lyapunov Exponents (LEs) are calculated corresponding to each one of the sensitivities with respect to the initial conditions $\partial \mathit{\Theta }/\partial {\mathbf{\Phi }}_0$ . If they are positive, a warning message is displayed to indicate the possibility of chaotic behaviour in the extended model differential equations, in line with the theory described in Section Sensitivity and stability analysis.

A modified LM ratio ${\hat{\rho }}_k$ specific to the NSALMN algorithm is then calculated, which denotes the agreement of the linear approximation made by the LM step as compared to the exact value of the objective function in the current iteration. Subsequently, the step acceptance criterion is defined by comparing ${\hat{\rho }}_k$ with a user-defined lower bound $\mu$ , below which the agreement of the linear approximation to the exact value of the objective function is deemed unacceptable, and the step is rejected. Hence, if the step is accepted the algorithm proceeds to the next iteration; otherwise the step is recalculated. In either of the above cases, the NSALMN algorithm uses a self-adaptive method to update the LM damping parameter, by considering it as a function of the modified LM ratio and the objective function value. This method renders the LM step behaving similar to the more robust gradient descent step when further from the minimum, or the faster converging Newton-Gauss step when closer to the minimum. The specific calculation procedures for the modified LM ratio and damping parameter can be found in Kimiaei. ⁶²

Case studies design

A series of case studies are designed in this study to allow for the testing of the developed method to solve the crankshaft dynamics inverse problem under ideal and realistic operating conditions. However, prior to listing all the case studies, it should be noted that the direct model is employed to provide the mock up ICT data, which are used as input to the inverse problem method. In specific, the direct model simulates the measured ICT by using as input the in-cylinder pressure for each engine cylinder, in accordance to the conditions dictated from each case study, hereinafter referred to as the measured in-cylinder pressures. Thus, the simulated ICT measurement is utilised as the input of the inverse problem method, which is in turn validated by using the method’s output (i.e. the optimised unknown in-cylinder pressure parameters) to reconstruct the in-cylinder pressure diagrams and compare them with the in-cylinder pressures used as input to the direct model. In addition, the inverse problem method performance is evaluated by utilising the following Key Performance Indicators (KPIs):

1. Number of successful iterations.

2. Ratio of number of failed to successful iterations.

3. Method execution time.

4. Average error and standard deviation between the reconstructed in-cylinder pressure diagrams and the direct model in-cylinder pressure input, as a percentage of the average peak pressure.

To test the performance of the developed inverse problem method, the designed case studies are classified under the following two main categories: (a) modelling inaccuracies and (b) distant initial conditions. For each of these two main categories, the inverse problem method is tested for the three subcategories of healthy, degraded and faulty engine conditions. Finally, for each subcategory, the extended model is employed using either the sigmoid or Weibe functions to facilitate a comparison in the performance of the inverse problem method. All the case studies are performed for the 100% engine load and are summarised in Table 3.

Table 3.

List of case studies and their respective identification codes.

Engine condition	Distant initial conditions		Modelling inaccuracies
	Sigmoid	Weibe	Sigmoid	Weibe
Healthy	S1	W1	S4	W4
Degraded	S2	W2	S5	W5
Faulty	S3	W3	S6	W6

For the case studies involved in the distant initial conditions category, the initial conditions of the unknown parameters vector ${\mathit{B}}_0$ that initialise the inverse problem method, are chosen to correspond to a different engine load than the load considered for the designed case studies. In specific, the values of ${\mathit{B}}_0$ are chosen to correspond to 70% MCR, whilst the measured ICT is simulated for 100% MCR. This represents the realistic case where in-cylinder pressure measurements are not available for a wide range of engine loads, hence the inverse problem method will be forced to employ initial conditions different from the actual parameter values.

For the case studies corresponding to the modelling inaccuracies category, the mock up ICT is simulated using the direct model, which employs an accurate coefficients matrix engine friction submodel. The latter models the engine friction by using manufacturer-tuned coefficients proportional to the DOF angular speed, instead of a less accurate constant engine friction submodel which is employed by the extended model in the inverse problem method (the specific details for both friction submodels are reported in Tsitsilonis et al. ⁴³ ). By employing different friction submodels, discrepancies in the modelling are purposely introduced. This emulates a realistic scenario, where the physical system (simulated by the direct model with the coefficients matrix engine friction submodel) is not perfectly represented by the extended model (employed by the inverse problem method), which utilises the simplistic constant engine friction submodel. The latter introduces a non-zero residual problem, which can potentially pose additional challenges for the accurate prediction of the in-cylinder pressure. Finally, the modelling inaccuracy is chosen to be introduced in the engine friction rather than in other parts of the engine model, since this part is a very frequent source of errors in dynamic engine modelling. ⁶³

Regarding the first subcategory of healthy engine conditions, the direct model simulation runs are performed by using as input the same in-cylinder pressure for all cylinders. The initial in-cylinder pressure function parameters are determined by performing a best fit of the in-cylinder pressure equation with the measured in-cylinder pressure data at 100% MCR reported in Guerrero and Jiménez-Espadafor ³¹ for the case of the reference system. This indicates that there are no discrepancies in the condition of each cylinder, and thus the engine is assumed to be in as-healthy-as-new condition.

The degraded engine conditions are simulated by using a random number generator to change the parameters of the healthy in-cylinder pressure functions for each cylinder. In specific, changing the in-cylinder pressure function parameters with respect to the healthy conditions results in subtle alterations in the shape of each individual in-cylinder pressure diagram, thus simulating the degraded conditions for all engine cylinders. As a result, the derived in-cylinder pressure function parameters are employed to generate the degraded in-cylinder pressure diagrams, which are used as input to the direct model for deriving the mock up ICT measurement. To derive the model parameters for engine degraded conditions, a uniform distribution was employed to generate values varying within a range of $\pm 10\%$ from the healthy parameters, which is typical of degraded engine conditions.^64,65 Subsequently, the unknown parameters vector $\mathit{B}$ is adjusted according to the following equation:

$$\mathit{x}=[x_1,x_2,\cdots x_N{]}^{\mathrm{T}}~(1-a)+2a\mathcal{U}(0,1)$$

(26)

$${\mathit{B}}_d=\mathrm{diag}(\mathit{x})\mathit{B}$$

(27)

where $a=0.1$ representing the desired percentage variation from the healthy conditions, wheras $\mathcal{U}(0,1)$ denotes the standard uniform distribution from 0 to 1, and therefore $\mathit{x}$ is a vector of the same size as $\mathit{B}$ that contains uniformly distributed numbers from 0.95 to 1.05 in accordance to equation (26). The unknown parameters vector $\mathit{B}$ is multiplied by a diagonal matrix containing the vector $\mathit{x}$ in accordance by equation (27), which results to the degraded engine parameters vector ${\mathit{B}}_d$ .

Finally, the faulty engine conditions are simulated by introducing a complete misfire on the first engine cylinder, whilst the other cylinders have the in-cylinder pressures corresponding to the degraded engine conditions at 100% MCR. It should be noted that misfire events can be identified easily by monitoring the exhaust gas temperatures of each cylinder. However, in this study, the misfire is used as an extreme scenario to test the capabilities and verify the robustness of the developed in-cylinder pressure prediction method.

Results and discussion

Case studies KPIs

Having tested the inverse problem method for all the case studies listed in Table 3, it was deduced that the method was successfully converged for all case studies, and the involved stability analysis did not indicate at any instance the potential occurrence of a chaotic behaviour. Therefore, the performance of the inverse problem method can be analysed holistically for every case study in accordance to the KPIs introduced in Section Case studies design, the obtained results of which are presented in Figure 4. Specific representative case studies were selected to be analysed in detail as presented below in this Section.

Figure 4.

(a) Successful iterations of inverse problem method, (b) failed iterations as a percentage of the successful iterations of the inverse problem method, (c) execution time of inverse problem method and (d) mean pressure error of predicted in-cylinder pressure by the inverse problem method as a percentage of the mean peak pressure (the error bars refer to the error standard deviation).

Firstly, the performance of the inverse problem method employing the Weibe and sigmoid functions is compared. As shown in Figure 4(a), when the Weibe function is employed as compared to the sigmoid function, approximately three additional successful iterations are required on average to achieve convergence, except for cases S3 and W3, where the same number of successful iterations takes place for both functions. Furthermore, as shown in Figure 4(b), failed iterations are observed only for cases W3–W6, which correspond the Weibe function, and range from 9% to 38% of the successful iterations number. In contrast, for cases corresponding to the sigmoid function, no failed iterations take place. Hence, it is deduced that the Weibe function is more computationally demanding, due to the additional parameter that it introduces for each cylinder, which is also a highly nonlinear one. This is confirmed from the results presented in Figure 4(c), where the inverse problem method execution time is consistently greater for all case studies employing the Weibe function, compared to the respective results for the case studies employing the sigmoid function. Regarding the accuracy of the inverse problem method, as shown in Figure 4(d), the mean pressure error does not indicate a specific trend, however the error standard deviations are slightly smaller for the case of the sigmoid function by an average of 0.04% for all case studies. This implies that the error does not depend on whether a Weibe or a sigmoid function is used as interpolating function. In addition, the maximum observed error is 0.3%, which indicated that the developed method is quite accurate in predicting the in-cylinder pressure diagrams.

The two main categories of the distant initial conditions and the modelling inaccuracies are compared (corresponding to the first and second columns in Table 3, i.e. case studies S1–S3 and W1–W3 vs S4–S6 and W4–W6). As shown in Figure 4(a), the distant initial conditions as compared to the modelling inaccuracies case studies take on average three additional iterations for the Weibe, or one less successful iteration for the sigmoid function. In addition, as shown in Figure 4(b), with the exception of case study W3, the failed iterations correspond entirely to the modelling inaccuracies case studies with the Weibe function. Furthermore, as expected the method execution time directly correlates with the total number of successful and failed iterations as shown in Figure 4(c). Overall, this indicates that for the distant initial conditions as compared to the modelling inaccuracies case studies, the inverse problem method performs faster by an average of 9.7 s for the Weibe function, or slower by an average of 0.5 s for the sigmoid function.

Most importantly, as shown in Figure 4(d), for the case of the distant initial conditions, despite case studies S3 and W3, which are discussed in detail below, the in-cylinder pressure mean error is almost negligible. This is expected as a zero-residual objective function has been formulated. In contrast, for all the case studies corresponding to the modelling inaccuracies, the in-cylinder pressure mean error only slightly increases (by an average of 0.07%) and a standard deviation of 0.25% of the mean maximum pressure, which although small, is significantly greater as compared to the distant initial conditions case studies. This indicates that from all the case studies, the ones corresponding to the modelling inaccuracies can result to the largest in-cylinder pressure mean errors. However, the error is still sufficiently small, and this is a favourable characteristic of the inverse problem, since a small effect on the in-cylinder pressure error is caused even for considerable modelling inaccuracies of the engine friction model.

The three subcategories of the healthy, degraded and faulty engine conditions are compared. Regarding the Weibe and sigmoid functions for both the modelling inaccuracies and distant initial conditions case studies, as shown in Figure 4(a) and (b), the faulty engine conditions (S3, W3 and S6, W6) seem to require the greatest number of iterations to achieve convergence; 13 and 8, respectively. Special notice is taken for the case study W3 due to the large number of failed iterations, which approaches 38% of the total successful iterations, as shown in Figure 4(b), as well as the largest mean in-cylinder pressure error standard deviation shown in Figure 4(d). In addition, case study S3 also exhibits considerably larger mean in-cylinder pressure error than the other two case studies in its respective category, which correspond to the distant initial conditions using the sigmoid function. This implies that the faulty engine condition represented with a misfire on cylinder no. 1, is the most computationally intensive case study. This is attributed to the fact that the initial conditions, which are always provided as those of a healthy engine with identical in-cylinder pressures for all cylinders, lie further from the objective function minimum. However, by examining the healthy and degraded engine conditions for both the Weibe and sigmoid functions, small or no differences are observed in the total number of successful or failed iterations respectively. Moreover, the inverse problem method appears to be the fastest as shown in Figure 4(c), when compared to the faulty engine conditions, by an average of 8.9 and 8.2 s for the Weibe and sigmoid functions, respectively. This implies that the developed method behaves very similarly, between the healthy and degraded engine conditions as the objective function minima between these two engine conditions are fairly similar despite the changes of the in-cylinder pressure parameters from cylinder to cylinder for the degraded engine conditions.

Individual case studies analysis

Following the KPIs analysis presented in the preceding section, a more specific discussion is provided for two representative case studies as explained below. Case study W3 is considered for detailed analysis due to the following reasons: (a) it exhibits the largest in-cylinder pressure error standard deviation of 0.25% as presented in Figure 4(d)); (b) it has by far the largest percentage of failed iterations, which contributes in exhibiting the third slowest execution time of 24.5 s (Figure 4(b) and (c), respectively); and (c) it is amongst the most complex case studies, as it employs the highly non-linear Weibe function and the engine experiences faulty conditions (deviations in the in-cylinder pressures of each individual cylinder and a misfire on cylinder no. 1).

Furthermore, case study S4 will also be analysed as justified in the following. It is amongst the fastest performing case studies with an execution time 10.9 s with no failed iterations (Figure 4(b) and (c), respectively). Moreover, it employs the sigmoid function, it corresponds to the engine healthy engine conditions and belongs to the case studies category of the modelling inaccuracies. Therefore, case study S4 is by far the one with the greatest difference in setup from case study W3. Both case studies S4 and W3 demonstrate interesting characteristics among the entire spectrum of the investigated case studies.

Starting with the case study S4, the engine ICT variations at the flywheel for one engine revolution are shown in Figure 5(a) (this figure shows both the ICT provided as input to and the ICT predicted by the inverse problem method). The ICT counter-intuitively exhibits five peaks instead of the expected 10 peaks, which would correspond to 10 engine cylinders. This is due to the non-linearities introduced primarily by the large reciprocating and rotating masses of the engine’s piston assembly.^43,66 It is observed that the measured ICT data sufficiently match the ICT from the extended model deployed within the inverse problem method.

Figure 5.

(a) Measured ICT data simulated by the direct model, and predicted ICT by the inverse problem method at the final iteration, for case study S4, (b) flywheel ICT residuals between measured and predicted flywheel torque at final iteration of the inverse problem method, for case study S4, (c) measured and predicted in-cylinder pressure by the inverse problem method, for case study S4 and (d) error between the simulated and predicted in-cylinder pressure by the inverse problem method, for case study S4.

As shown in Figure 5(b), the residuals between these ICT variations fluctuate at approximately two orders of magnitude less than the measured ICT data, which demonstrates a very small difference. In specific, this difference is also positive, which signifies that the measured ICT data is greater than the predicted ICT, which is derived by using the extended direct model (called by the inverse problem method) employing the simplified constant engine friction model. This results from how the inverse problem optimisation algorithm approaches the objective function minimum, which partially depends on the provided initial conditions. Furthermore, the torque residuals exhibit a high frequency fluctuation, which indicates that there are intrinsic physical differences between the measured and predicted ICT. However, this frequency is higher than the 40th engine order, which for the case where a torque metre is used to obtain the measured ICT data (in contrast to this study where a simulation is performed) can be considered as noise that needs to be filtered.

Figure 5(c) presents the predicted pressure diagrams by the inverse problem method (by employing the extended direct model) along with the pressure diagrams (fitted from the measured pressure reported in Guerrero and Jiménez-Espadafor ³¹ for the engine operation at 100% load), which were used as input by the direct model. As this case study corresponds to healthy engine conditions, the measured pressure diagrams of each engine cylinder are identical and shifted in the crank angle axis by each cylinder phase angle. Most importantly as shown in Figure 5(d), the predicted in-cylinder pressure diagrams exhibit a very small error compared to the measured in-cylinder pressure diagrams (which were provided as input in the direct model to simulate the ICT data), reaching a maximum of 0.8 bar for cylinder no. 9. As also mentioned in the preceding paragraphs, these obtained results confirm that the modelling inaccuracies in the engine friction sumbodel do not result in significant errors in the predicted in-cylinder pressure. Similarly to the ICT residuals, the in-cylinder pressure error is primarily positive, which indicates that the simulated in-cylinder pressure is, for most crank angle positions, greater than the predicted in-cylinder pressure. Furthermore, the pressure error reaches its greatest values for each individual cylinder after the respective cylinder’s TDC, near the in-cylinder peak pressure, with the maximum error incurred in cylinder no. 9. This takes place because changes in the in-cylinder pressure function unknown parameters affect the shape of the in-cylinder pressure curve mostly near its peak, which occurs close to the TDC location. Hence the error can become more pronounced at the TDC location.

For the case study W3 that corresponds to faulty engine conditions with a misfire on cylinder No. 1, degraded conditions of all the other cylinders and use of Weibe function, the ICT predicted by the inverse problem method as well as the ICT provided as input to the inverse problem method are shown in Figure 6(a). It is observed that the ICT exhibits four peaks (instead of five as discussed for the case study S4). In specific, the ICT peak that is located after the cylinder no. 9 TDC is not present when a misfire on cylinder no. 1 occurs, which demonstrates that the effect of the misfire on the flywheel ICT lags the actual misfire event by approximately 36° crank angle. This occurs as a result of the crankshaft’s flexibility, hence the changes in torque due to the different forces exerted in each piston do not get transmitted instantly to the flywheel.^31,43 In this case study as shown in Figure 6(b), since no modelling inaccuracies are introduced, the ICT residuals are three orders of magnitude smaller than the measured ICT, which is also one order of magnitude less than those for case study S4 discussed above. Furthermore, a similar high frequency fluctuation is not visible as it was found for the case study S4, which is attributed to the fact that no modelling inaccuracies are introduced for this case study. In addition, the average residuals for case studies W3 and S4 have an average value of approximately 0 Nm and $2.5\times {10}^4$  Nm respectively, which demonstrates that in the case study W3, a zero-residual problem has been formulated and no inaccuracies in the modelling occurred.

Figure 6.

(a) Measured ICT data simulated by the direct model, and predicted ICT by the inverse problem method at the final iteration, for case study W3, (b) flywheel ICT residuals between measured and predicted flywheel torque at final iteration of the inverse problem method, for case study W3, (c) measured and predicted in-cylinder pressure by the inverse problem method, for case study W3 and (d) error between the simulated and predicted in-cylinder pressure by the inverse problem method, for case study W3.

Moreover, as shown in Figure 6(c), the measured in-cylinder pressure variations almost coincide with the predicted ones. However, as shown in Figure 6(d), the in-cylinder pressure error is significantly greater by one order of magnitude compared to the respective errors observed in the case study S4. This error is still acceptable, however, since the peak pressure values reach an average of approximately 100 bar for the cylinders no. 2–10, the maximum pressure error reaches a value of around 2.2 bar, which occurs for cylinder no. 8 close to its TDC. The much larger magnitudes for the in-cylinder pressure errors obtained in the case study W3 are attributed to the significantly greater non-linearities of the Weibe function as compared to the sigmoid function.

These non-linearities introduce radical topological changes in the shape of the objective function, which especially when providing initial conditions further from its minimum, can pose a challenge to the optimisation algorithm’s convergence. This is demonstrated by the visualised two-dimensional subdomains of the objective function for the sigmoid and Wiebe functions shown in Figure 7. In specific, for the case of the sigmoid function (Figure 7(a)), the objective function contours exhibit a mildly elliptical shape, which greatly simplifies the task of the optimisation algorithm. In contrast, for the case of the Weibe function (Figure 7(b)), the objective function contours form a narrow canyon shape that contains a number of local minima, which can potentially pose challenges for the optimisation algorithm. As observed from the presented path of the optimisation algorithm, it gets trapped close to the local minimum within the narrow canyon. The algorithm subsequently resets and the optimisation process recommences from the closest point outside the canyon, where the objective function slope is more pronounced and thus, convergence can be achieved. Hence, when the Weibe function is employed, the optimisation algorithm has to follow a more complex path to reach its minimum, and as a result, its efficiency and performance is not as consistent. However despite the above, as shown from both Figure 6(a) and (b), the cylinder no. 1 with misfiring was identified correctly based on the predicted in-cylinder pressure by the inverse problem method and with the best accuracy, as the prediction exhibits the smallest maximum pressure error compared to all other cylinders.

Figure 7.

(a) Two-dimensional subdomain of the objective function by varying the unknown parameters $k_1$ and $k_2$ (in normalised form) for cylinder no. 2, for the case study S4 and (b) two-Dimensional subdomain of the objective function by varying the unknown parameters ${\theta }_0$ and m (in normalised form) for cylinder no. 2, for the case study W3.

The inverse problem method was also tested for two other engine loads of 55% and 70%. These tests were performed by using the same case studies as above, with the corresponding changes made in their set up in order to accommodate for the different engine loads being tested. The results obtained were similar to the ones described above for 100% load.

In addition, from the preceding analysis, it was observed that the developed inverse model succeeded in identifying the sigmoid and Weibe functions parameters (three and four parameters per cylinder, respectively), and as a result the in-cylinder pressure diagrams for all the investigated case studies. In this study, the ICT was provided with a resolution of one data sample per degree crank angle. Therefore, it is concluded that the inverse model is characterised as identifiable ⁶⁷ for both the sigmoid and Weibe functions, provided that the ICT signal resolution is at least one sample per degree crank angle.

Sensitivity analysis

The results from the S4 and W3 case studies can be correlated to specific features of the flywheel torque sensitivities with respect to the unknown in-cylinder pressure parameters for improving the understanding of the inverse problem. As the sensitivities exhibit a consistent oscillatory behaviour throughout a complete crankshaft revolution, the sensitivity analysis herein was conducted by using the Euclidean norm of all the sensitivities. The respective results are presented in Figure 8.

Figure 8.

(a) Normalised sensitivities Euclidean norm of the unknown parameters with respect to the flywheel torque, for case study S4 and (b) normalised sensitivities Euclidean norm of the unknown parameters with respect to the flywheel torque, for case study W3.

By considering the sensitivities corresponding to parameters $k_1$ and m presented in Figure 8(a) and (b) respectively, it is observed that their Euclidean norms exhibit on average the smallest values across all cylinders, as compared to the sensitivities of the other unknown parameters. This indicates that large changes in these two parameters result in small changes in the engine ICT. Therefore, when considering the inverse problem, small variation of the engine ICT could result in large changes in the values of parameters $k_1$ and m for all cylinders. As a result, it is expected to be more challenging for the developed method to accurately determine these parameters, since small deviations from the actual values of the ICT measurement (e.g. due to measurement errors or inaccuracies in the engine dynamics model) can result in large errors in the estimations of these parameters.

Furthermore, by comparing the parameter sensitivities across individual cylinders, it is observed from Figure 8(a) and (b) that the cylinder further away from the flywheel exhibits larger sensitivity norm for all parameters. This implies that the unknown parameters of the cylinders further from the flywheel exhibit a greater impact on the flywheel torque calculation. As a result, the engine ICT is less sensitive to changes in the parameters of the cylinders closest to the flywheel.

In addition, as observed by the results presented in Figure 7, due to the misfire on cylinder no. 1, the sensitivities of all the unknown parameters are zero. This is due to the fact that combustion does not take place when a misfire occurs. As a result, the unknown parameters, which relate to the respective cylinder’s combustion, do not contribute to the engine’s ICT. Finally, the sensitivity norms corresponding to parameter $p_{\mathrm{me}}$ for cylinders no. 2–10, are identical in both case studies, as shown in Figure 8(a) and (b). This occurs since this parameter is employed in the polytropic expansion relationship, which is modelled by equation (4) in both case studies. As a result, the $p_{\mathrm{me}}$ parameter sensitivity remains the same regardless of whether the Weibe or sigmoid interpolating functions are used.

Conclusions

This study developed a novel method to solve the inverse problem of determining the in-cylinder pressure, based on the measurement of the crankshaft ICT. This method was based on the IVP technique and employed an optimisation algorithm along with a lumped parameter crankshaft dynamics model. The applicability of the method was tested for various cases including healthy, degraded and faulty engine conditions, as well as considering the lack of in-cylinder pressure data, and modelling inaccuracies. It was found that this method accurately predicted the in-cylinder pressure with a maximum error of less than 0.5% from the respective measured values.

In specific, from the two interpolating functions that were tested, the Weibe function performed the slowest response by an average of 11.2 s as compared to the sigmoid. However, there was no significant effect to the in-cylinder pressure error. As a result, it is recommended to use the Weibe interpolating function, since despite the increased computational burden, it is able to model the combustion more accurately, due to the employed additional parameter, compared to the sigmoid function. Furthermore, the inverse problem method demonstrated its robustness when tested for both distant initial conditions and modelling inaccuracies, with a slightly errors on the predicted in-cylinder pressure by 0.07% on average for the latter case. Testing at healthy, degraded and faulty engine conditions provided lower execution time was slower for the degraded and healthy engine conditions by an average of 8.5 s as compared to the healthy conditions, whereas no significant changes were observed for the in-cylinder pressure error.

Further insight into the inverse problem of estimating the in-cylinder pressures by utilising the engine ICT at the flywheel was also obtained, by considering the individual sensitivities of the unknown in-cylinder pressure parameters. In specific, the cylinders further from the flywheel exhibited smaller sensitivities and thus have a smaller impact on the engine ICT. In addition, the parameters $k_1$ and m of the sigmoid and Weibe functions respectively, exhibited the smallest sensitivities. As a result, small deviations from the actual values of the ICT measurement (e.g. due to measurement error or modelling inaccuracies) could result in large errors in the estimation of these parameters.

Therefore, it can be concluded that the developed method for predicting the in-cylinder pressure based on the engine ICT measurement exhibited adequate performance. This method implementation has the potential for a significant change in the maritime industry norm by employing a permanently installed shaft torque metre sensor with the capability of obtaining higher sample rates, thus replacing the use of external in-cylinder pressure measurement kits, or expensive multi-sensor continuous pressure monitoring systems. This will allow the continuous and non-intrusive monitoring of the in-cylinder pressure without requiring extensive retrofitting to deploy multiple sensor networks.

Future studies will include testing of the developed method by using the ICT data scheduled to be obtained from large marine four-stroke and two-stroke engines. Further development will include modifying the method to be implemented in on-line applications, as well as expanding the method to integrate decision support capabilities.