Comparison of Physics-Based, Semi-Empirical and Neural Network-Based Models for Model-Based Combustion Control in a 3.0 L Diesel Engine

: A comparison of four di ﬀ erent control-oriented models has been carried out in this paper for the simulation of the main combustion metrics in diesel engines, i.e., combustion phasing, peak ﬁring pressure, and brake mean e ﬀ ective pressure. The aim of the investigation has been to understand the potential of each approach in view of their implementation in the engine control unit (ECU) for onboard combustion control applications. The four developed control-oriented models, namely the baseline physics-based model, the artiﬁcial neural network (ANN) physics-based model, the semi-empirical model, and direct ANN model, have been assessed and compared under steady-state conditions and over the Worldwide Harmonized Heavy-duty Transient Cycle (WHTC) for a Euro VI FPT F1C 3.0 L diesel engine. Moreover, a new procedure has been introduced for the selection of the input parameters. The direct ANN model has shown the best accuracy in the estimation of the combustion metrics under both steady-state / transient operating conditions, since the root mean square errors are of the order of 0.25 / 1.1 deg, 0.85 / 9.6 bar, and 0.071 / 0.7 bar for combustion phasing, peak ﬁring pressure, and brake mean e ﬀ ective pressure, respectively. Moreover, it requires the least computational time, that is, less than 50 µ s when the model is run on a rapid prototyping device. Therefore, it can be considered the best candidate for model-based combustion control applications.


Introduction
Nowadays, emission and fuel consumption reductions are the two main challenges for internal combustion engines, and in particular for diesel technology [1][2][3]. Many techniques have been proposed to achieve this aim, such as variable geometry turbocharger (VGT), high-pressure common rail systems [4][5][6][7][8], exhaust gas recirculation (EGR) [9], innovative combustion concepts such as HCCI and PCCI [10], and innovative combustion controls [11][12][13][14]. As far as these solutions are concerned, model-based combustion control is expected to make a significant contribution in the near future, thanks to the possibility of its integration with the emerging Vehicle-to-Everything (V2X) architectures. This conclusion is also supported by recent international projects, such as IMPERIUM H2020 [14], in which a significant reduction in CO 2 emissions will be achieved (−20%) for heavy-duty trucks, thanks to the use of model-based controllers coupled with V2X systems. The model-based control technology will also be boosted by the development of new sensors and by the increasing computational performance of the new multi-core processors that are now available for mobility applications.
Several advantages may be obtained from adopting a model-based controller instead of a conventional map-based one, such as the possibility of realizing a real-time optimization of the engine

•
Baseline physics-based model: this model was previously presented by the authors in [12] for control-oriented applications. However, in the present study, the model calibration procedure has been refined and the performance of the model has been improved with respect to the previously reported version. In this approach, the tuning parameters are modeled by means of power-law functions, in which the input parameters are the main engine operating variables. • ANN physics-based model: this is a new modeling approach, which is based on the use of ANNs to predict the tuning parameters of the aforementioned physics-based model.

•
Direct semi-empirical model: in this approach, semi-empirical correlations, which are constituted by power-law functions, are used to directly estimate MFB50, PFP, and BMEP. • Direct ANN model: this approach exploits feed-forward artificial neural networks to directly estimate MFB50, PFP, and BMEP. Details of the methodology used for the training and optimal selection of the number of neurons are also provided in this study.
The four different models have been assessed considering the same experimental dataset, for a 3.0 L FPT diesel engine for light-duty applications, and their performances have been compared under steady-state conditions and in transient operation over a WHTC.
In addition to the comparison of the performance of the previous models, another innovative contribution introduced in this paper concerns the methodology used for the optimal selection of the input parameters of the models. The proposed method is based on the sequential use of the Pearson correlation and partial correlation coefficients, which are used to identify the least and most robust correlation variables that should be excluded or included as model inputs. The use of the Pearson correlation and partial correlation analysis allows the computational effort required for the identification of the input parameters to be reduced significantly.
The paper is organized as follows. Section 2 summarizes the experimental setup and the engine specifications. Section 3 reports the description of the four investigated models. Section 4 is focused on the presentation of the methodology used for the selection of the model input variables. Finally, Section 5 reports the main results, which include: • A comparison of the performance of the physics-based model tuned using the methodology introduced in this paper (based on Pearson correlation and partial correlation analysis) with that of the previous version reported in [12]. This research is a first step towards the development of a model-based combustion controller, which will be developed and tested on the real engine, through rapid prototyping, in the near future.

Experimental Setup and Engine Conditions
The experimental tests were run on a Euro VI FPT F1C 3.0 L diesel engine (FPT Motorenforschung AG, Arbon, Switzerland), within a research project conducted in collaboration with FPT Industrial [11]. The activity was carried out in the dynamic test bench facility at the Politecnico di Torino. The main technical specifications of the engine are reported in Table 1. The experimental tests were conducted under steady-state and transient conditions. The steadystate tests have been divided into three categories ( Figure 2): • a complete engine map (123 tests, indicated with the blue circles in Figure 2).
• EGR-sweep tests (162 tests, carried out on the points indicated with the red diamonds in Figure 2). • sweep tests of the main injection timing (SOImain) and injection pressure (pf) (125 tests, carried out on the points indicated with the black circles in Figure 2). Validation of the methods was carried out over a WHTC. Details on these validations are given in Section 5.  An ETAS ES910 rapid prototyping device, which is equipped with a main Freescale PowerQUICCTM III MPC8548 processor with an 800 MHz clock, was used to verify the computational time required for the models investigated in this paper.

Description of the Models
The experimental tests were conducted under steady-state and transient conditions. The steady-state tests have been divided into three categories ( Figure 2): • a complete engine map (123 tests, indicated with the blue circles in Figure 2). • EGR-sweep tests (162 tests, carried out on the points indicated with the red diamonds in Figure 2). • sweep tests of the main injection timing (SOI main ) and injection pressure (p f ) (125 tests, carried out on the points indicated with the black circles in Figure 2). The fuel is characterized by a density of 835 kg/m 3 at 14 °C, a viscosity of 2 mm 2 /s at 40 °C and a cetane number equal to 43. The considered injection pattern features two pilot shots and a main shot.
The engine is shown in Figure 1. It is equipped with a VGT turbocharger, a short-route cooled EGR system and a flap valve located on the exhaust side downstream from the turbine. Details about the layout of the engine, the test bench and the used sensors (including the main accuracy) have already been described in detail in [38] and are not repeated in this paper for the sake of brevity.
An ETAS ES910 rapid prototyping device, which is equipped with a main Freescale PowerQUICCTM III MPC8548 processor with an 800 MHz clock, was used to verify the computational time required for the models investigated in this paper. The experimental tests were conducted under steady-state and transient conditions. The steadystate tests have been divided into three categories ( Figure 2): • a complete engine map (123 tests, indicated with the blue circles in Figure 2).
• EGR-sweep tests (162 tests, carried out on the points indicated with the red diamonds in Figure 2). • sweep tests of the main injection timing (SOImain) and injection pressure (pf) (125 tests, carried out on the points indicated with the black circles in Figure 2). Validation of the methods was carried out over a WHTC. Details on these validations are given in Section 5.  Validation of the methods was carried out over a WHTC. Details on these validations are given in Section 5.

Physics-Based Model
Although the physics-based combustion model has already been presented in other studies, a summary of the approach and of the main equations is reported hereafter.
A scheme of the model is shown in Figure 3.

Physics-based Model
Although the physics-based combustion model has already been presented in other studies, a summary of the approach and of the main equations is reported hereafter.
A scheme of the model is shown in Figure 3. . Scheme of the physics-based combustion model [12]. N in Figure 3 indicates the engine speed, while pIMF and TIMF the intake manifold pressure and temperature, respectively, pEMF the exhaust manifold pressure, SOImain the start of injection of the main pulse, SOIpil,j the start of injection of the pilot pulse j, qpil,j the injected fuel quantity of pilot injection j, qtot the total injected fuel quantity, ValveEGR the opening position of the high pressure EGR valve and mair the fresh air trapped mass per cycle/cylinder. The first step consists in the evaluation of the chemical heat release, according to the following equations (AFM approach): where j indicates the generic injection pulse, HL indicates the lower heating value of the fuel, and K and  are the combustion rate coefficient and ignition delay coefficients, respectively.  Figure 3. Scheme of the physics-based combustion model [12]. N in Figure 3 indicates the engine speed, while p IMF and T IMF the intake manifold pressure and temperature, respectively, p EMF the exhaust manifold pressure, SOI main the start of injection of the main pulse, SOI pil,j the start of injection of the pilot pulse j, q pil,j the injected fuel quantity of pilot injection j, q tot the total injected fuel quantity, Valve EGR the opening position of the high pressure EGR valve and m air the fresh air trapped mass per cycle/cylinder. The first step consists in the evaluation of the chemical heat release, according to the following equations (AFM approach): The compression and expansion phases are modeled with polytropic processes: while the pressure at IVC (i.e., the starting condition) is correlated with the pressure in the intake manifold: The main combustion metrics, that is, MFB50, IMEP360 (gross Indicated Mean Effective Pressure) and PFP (Peak Firing Pressure), can now be evaluated. Finally, friction (FMEP) and pumping (PMEP) models are used, and the net IMEP (IMEP720) and BMEP (Brake Mean Effective Pressure) are estimated: The tuning parameters of the physics-based model are underlined in Equations (1)- (14). The Chen-Flynn approach [40] was adopted to estimate FMEP: where x 1-4 are fitting parameters. The intake oxygen concentration (O 2 ) is an operating variable that is closely correlated with the combustion process, and it is therefore used as an input variable in several sub-models. It is measured, on the experimental test bench, by means of a paramagnetic sensor, which is included in the test cell gas analyzer. However, if the developed combustion model is intended to be used for model-based control, this sensor is not available onboard, and O 2 therefore needs to be estimated by means of a specific sub-model.
The following function has been used to estimate O 2 [11,41]: where x 1-2 are fitting parameters, RAF is the relative air-to-fuel ratio and X r,EGR is the EGR rate. Therefore, it is necessary to estimate RAF and X r,EGR in order to evaluate the intake O 2 concentration. The relative air-to-fuel ratio is defined as follows: where m air is the trapped air mass, m f,inj is the injected fuel mass and AF st is the stoichiometric air-to-fuel ratio (taken as 14.4 in this work). The accuracy of the estimation of RAF depends to a great extent on the accuracy of the estimation of m air . If m air is obtained from the engine air flow sensor, its accuracy is generally not very high (the error is typically of the order of 5%).
An alternative way of estimating RAF (the method that was adopted in this study) is to develop a semi-empirical correlation, such as: where x 1-4 are fitted parameters. The correlation is fitted using a calibration dataset, in which the experimental values of RAF are obtained from steady-state measurements carried out using an accurate sensor ( Bosch Lambda sensor LSU 4.9 (Bosch, Gerlingen, Germany)). With reference to the EGR rate X r,EGR , it is defined as follows: X r,EGR = m EGR m EGR + m air (19) where m EGR is the trapped EGR mass. The EGR rate can in general be derived experimentally on the basis of the measured emissions and of the measured intake CO 2 concentration. However, if the combustion model is intended to be applied for control-oriented applications, a real-time estimation of the EGR rate is needed, through the use of a dedicated submodel. To this aim, several approaches can be used (e.g., see [41]).
The simplest method to estimate the EGR rate (which was adopted in this paper) is to resort to look-up tables that depend on the engine speed and pedal position. These look-up tables were derived on the basis of the analysis of steady-state measurements.
As far as the calibration procedure is concerned, as reported in the previous sections, the physics-based model is characterized by several tuning parameters, which are mentioned in Table 2. The calibration of these parameters is carried out in two steps. In the first step, the optimal values of all the parameters are identified, test by test, on the basis of the available experimental data. The process used for the identification of these values is reported in detail in [16], and is based on the best matching between the experimental and predicted heat release and in-cylinder pressure curves. At the end of this first step, a dataset of the optimal values of the model calibration parameters is available.
In the second step, on the basis of the dataset identified in step 1, the model calibration parameters are correlated with input variables, which are generally related to the engine operating conditions (e.g., intake manifold pressure/temperature, intake oxygen concentration, engine speed and load . . . ), using such mathematical methods as look-up tables, empirical or semi-empirical functions or machine learning tools. Semi-empirical correlations, based on power-law functions, are used in the baseline version of the physics-based model (see [12]), since they are capable of capturing the non-linear behavior of the correlations. However, in previous studies, the input parameters for these correlations were identified on the basis of a trial and error approach.
In the next sections, two aspects are investigated in order to improve the physics-based model with respect to previously reported versions:

1.
A new procedure for the identification of the optimal set of input parameters is developed. This procedure is based on the joint use of the Pearson correlation analysis, partial correlation analysis, and sensitivity analysis, and it allows the performance of the baseline physics-based model to be improved with respect to the previously developed versions.

2.
An alternative mathematical method (i.e., ANNs) is used to estimate the model calibration parameters: where x 0 , x 1 , . . . , x n are correlation variables. This approach has never been investigated before, and is denoted as "ANN physics-based model".

Direct Semi-Empirical Models
In addition to the baseline and ANN physics-based models, semi-empirical correlations were also developed in order to directly predict the MFB50, PFP, and BMEP metrics.
These correlations are based on power-law functions of the type: The identification of the x 0 , x 1 . . . x n correlation variables and the fitting of the model are analyzed in the next sections. The semi-empirical models are trained on the basis of the same experimental data used for the physics-based model calibration.

Direct Artificial Neural Networks
The final method investigated to predict MFB50, PFP, and BMEP is based on feed-forward artificial neural networks with one inner layer: where x 0 , x 1 , . . . , x n are the input variables of the neural networks. The networks are trained on the basis of the same experimental data used for the physics-based model calibration. The methodology used for the calibration of the model, including the selection of the optimal number of neurons, is discussed in detail in Section 5.

Selection of the Model Input Variables
Selecting the most appropriate input parameters as independent variables for all the previously mentioned approaches is a key-aspect to obtain accurate and robust models.
In this paper, optimal sets of input parameters have to be defined for: • The power-law correlations or the ANNs that evaluate the calibration parameters of the baseline physics-based model and of the ANN physics-based model, respectively.
In general, the input parameters are related to the engine operating conditions, and can be either directly measured quantities or derived quantities. The directly measured parameters (such as engine speed N, the start of injection of the pilot/main pulses SOI pil /SOI main , the total injection quantity of the pilot pulses q pil,tot , the injection quantity of the main pulse q main , the total injection quantity q tot , the intake/exhaust pressure and temperature p IMF ,T IMF /p EMF ,T EMF, the common rail pressure p f , the trapped air mass m air and the EGR valve opening position Valve EGR ) are the most commonly used input variables for control-oriented applications, since they are readily available for the ECU without the need of any additional sub-model. However, derived parameters (such as intake density ρ IMF , in-cylinder gas density and temperature at the start of each injection pulse ρ SOI,pil , ρ SOI,main, T SOI,pil , T SOI,main , the start of combustion of the pilot and main pulses SOC pil /SOC main , in-cylinder gas density and temperature at the start of combustion for each injection pulse ρ SOC,pil , ρ SOC,main, T SOC,pil , T SOC,main , the pressure difference between the exhaust and intake manifold ∆p exh-int ) or even outcomes of sub-models (such as PFP, MFB50) have also been used as input variables in models reported in the literature, since their correlation with the dependent variables is sometimes greater than that of the directly measured parameters. In this study, all these kinds of parameters were initially considered as candidate input variables to develop the models. Table 3 reports the list of the directly measured input parameters (1st column), of the generated parameters (from the directly measured ones or from the combustion model) which can be used as input parameters (2nd column), and of the dependent output variables which have been considered in the present investigation. A total of 46 parameters have been included. It should be noted that some of the dependent variables are the calibration parameters of the physics-based model, which needs to be estimated for the model application (3rd column), while MFB50, PFP, and BMEP are the combustion metrics that are adopted for combustion control applications, i.e., the final output variables of all the four modeling approaches that have been investigated in this study (4th column).
The procedure that was identified to select the optimal set of input parameters is based on the joint use of the Pearson correlation and partial correlation coefficients (see [42][43][44][45] for details concerning the use of these coefficients), and is reported in detail in Appendix A for the paper readability purposes. However, a synthetic description is reported hereafter.
Basically, the procedure involves 3 steps: 1. First, the Pearson correlation coefficient is calculated for all the possible combinations of the 46 variables reported in Table 3 (the results are reported in Appendix A in Table A1). This coefficient measures the linear dependence between two variables, x i and x j , and is defined as follows: where cov(x i , x j ) is the covariance between x i and x j , while var(x i ) and var(x j ) are the variances of x i and x j , respectively. The Pearson correlation coefficient varies between −1 and +1. After the analysis of the Pearson correlation coefficients, an initial candidate set of input variables is identified for each output quantity that has to be estimated (see Table A2 in Appendix A).

2.
In the second step, the partial correlation coefficient is evaluated between the dependent variables and each input variable of the initial candidate set, when removing the effect of the remaining input variables of the same set. It should be recalled that the partial correlation coefficient is able to measure the linear dependence between two variables, x i and x j , when the effect of the other input variables variable (z 1, z 2 , . . . , z i ) is removed: The partial correlation coefficient is more robust than the Pearson coefficient. However, the number of parameters should not be too high, otherwise, the consistency of the method decreases. Also, the partial correlation coefficient varies between −1 and +1. As reported in Appendix A (see Table A3), the simultaneous analysis of the Pearson and of the partial correlation coefficients has allowed the least and most robust correlation variables of the initial candidate set to be identified for each dependent output variable. Therefore, some variables were excluded from the initial candidate set identified in step 1 (since they showed low values of the Pearson and partial correlation coefficients), while other variables were kept since they showed high values of both the Pearson and partial correlation coefficients.

3.
The remaining correlation variables of the initial candidate set which showed intermediate values of the Pearson and partial correlation coefficients were analyzed on the basis of a sensitivity analysis, in order to identify those that had to be kept and those that had to be excluded. A power law function, such as that reported in Equation (21), was used to model each output dependent parameter, and all the possible combinations of input variables were identified by evaluating, for each combination, the related model fitting precision, which was quantified by a correlation coefficient R adj : where n exp is the total number of experimental data and k is the total number of free coefficients. The adjusted correlation coefficient R adj takes into account the number of available experimental tests used for the model fitting, in relation to the number of free parameters of the model. In the end, the combinations which led to the best trade-off between the prediction accuracy and the number of input parameters were selected. Details on the sensitivity analysis have not been reported in this paper for the sake of brevity.
It should be noted that the proposed approach, i.e., the preliminary selection of the input variables through Pearson and partial correlation coefficients, and the subsequent sensitivity analysis carried out on a reduced set of variables, is more computationally efficient than a pure sensitivity analysis approach which includes the entire set of input variables.
In fact, the number of possible combinations that have to be investigated in the sensitivity analysis can be calculated through the following formula: (26) where N c is the number of possible combinations for the sensitivity analysis and j is the number of independent variables. According to Equation (26), N c increases to a great extent as j increases. Therefore, the combined use of the Pearson correlation and partial correlation analysis is a very useful and efficient way of decreasing the computational effort of the sensitivity analysis, since it allows the poorly correlated variables to be excluded a-priori, or the highly correlated variables to be selected, so that the sensitivity analysis (which is highly time consuming, since it requires model fitting of each possible combination), is only carried out for a reduced set of variables.

Results and Discussion
The improvement in the prediction performance of the physics-based model (compared to the previous version reported in [12]), when adopting the new methodology for the input parameter selection described in Section 4, is reported first. Subsequently, the performances of the four investigated models described in Section 3 are compared, in terms of accuracy at steady-state conditions and over WHTC, and in terms of required computational time. The main advantages and drawbacks of each kind of approach are also discussed.

Improvement to the Baseline Physics-Based Combustion Model
The performance of the baseline physics-based combustion model has been improved, with respect to the previous version reported in other papers [12]. The set of input variables for the model correlations has been revised, on the basis of the Pearson correlation and partial correlation analysis, according to the procedure described in the previous section.
The final list of selected input variables for each dependent variable and the corresponding fitting precision are shown in Table 4. It should be noted that some correlations are improved with respect to pure power-law functions. Table 4. Final list of the selected input variables, the related correlations, for each tuning parameter of the physics-based model and the corresponding fitting precision in terms of correlation coefficient (R 2 ), root mean squared error (RMSE) and relative RMSE (RMSE rel ).

SOI,main
Improved correlation :   [12]) and the revised model which embeds the improved correlations, concerning the estimation of BMEP, PFP, and MFB50, for the steady-state tests reported in Figure 2. Figures 4 and 5 report the accuracy of the original physics-based model with the baseline calibration (presented in [12]) and the revised model which embeds the improved correlations, concerning the estimation of BMEP, PFP, and MFB50, for the steady-state tests reported in Figure 2.
It can be seen that, although the accuracy in the estimation of BMEP is somewhat similar, an improvement in the estimation of PFP, and MFB50 can be obtained, since the related RMSE values decrease from 2.47 bar to 1.82 bar and from 0.86 deg to 0.63 deg, respectively.

Direct Semi-empirical Models of MFB50, PFP, and BMEP
The procedure used to identify the input parameters of the physics-based model, which was described in Section 4, was also adopted to select the input variables of the direct semi-empirical models. The correlations identified for MFB50, PFP, and BMEP are reported in Table 5.  It can be seen that, although the accuracy in the estimation of BMEP is somewhat similar, an improvement in the estimation of PFP, and MFB50 can be obtained, since the related RMSE values decrease from 2.47 bar to 1.82 bar and from 0.86 deg to 0.63 deg, respectively.

Direct Semi-Empirical Models of MFB50, PFP, and BMEP
The procedure used to identify the input parameters of the physics-based model, which was described in Section 4, was also adopted to select the input variables of the direct semi-empirical models. The correlations identified for MFB50, PFP, and BMEP are reported in Table 5.

ANN-Based Models
Artificial neural networks have a great capacity to capture complex nonlinear relationships between variables through relatively simple mathematical operations, and require very little computational effort. Therefore, they are good candidates for control-oriented applications. In this section, ANNs have been used to: • Replace the power law-based correlations of the physics-based models shown in Table 4. The resulting new physics-based model is referred to as "ANN physics-based model".

•
Directly simulate the MFB50, PFP, and BMEP combustion metrics, analogously to the semi-empirical models shown in Table 5.
In both cases, the sets of input variables that were used for the ANNs are the same as those shown in Table 4 (physics-based model) and Table 5 (semi-empirical models). The experimental data used for ANN training is constituted by the 410 tests shown in Figure 2.
In general, the implementation of ANNs involves three steps: the selection of the input variables, the training phase, and the testing phase. The available experimental dataset is usually split randomly for training, testing and validation purposes. The main parameters that have to be selected are the net type, the number of inner layers and the number of neurons for each layer, the activation function, and the training algorithm.
As far as the network type is concerned, feed-forward networks are typically used in the literature, and have also been used in this paper. These types of networks contain an input layer, some hidden layers, and an output layer. The back-propagation training algorithm is generally used as a training algorithm for feed-forward neural networks [46]. The Bayesian regularization training function 'trainbr' is one of the best-known BP training algorithms; it can lead to very precise ANNs [46,47] and is able to provide good generalization for difficult, small or noisy datasets [48], even though its number of epochs is usually higher. In this work, the training dataset is not large and the network structure is quite simple, therefore a relatively large epoch number should not lead to an unacceptably long computation time. For these reasons, 'trainbr' was selected as the training algorithm.
In order to speed up the learning process for a high number of networks, the inputs are usually normalized before sending them to the input layer. In this work, the 'mapminmax' process function [49] has been used to normalize the inputs, so that all the values fall into the [−1, 1] interval. Similarly, network outputs can also be associated with processing functions [49].
The activation function that was used in this work is the 'tansig' function. It was found, with reference to the number of hidden layers, that one layer was enough to simulate the parameters of the physics-based model and the combustion metrics (MFB50, PFP, BMEP). The use of one hidden layer for feed-forward neural networks is also quite common in the literature [47,50,51], especially when the number of training data is not very high.
The MSE minimization criterion was used for the training process, as suggested in [52]. The ANN performance is evaluated by comparing the network outputs with the experimental data through regression analysis [53]. The values of the mean-squared-error MSE and the determination coefficient R 2 are used to measure the performance of ANNs. Details of the ANN parameters developed on the MATLAB (R2016b, MathWorks, Natick, MA, USA) platform are shown in Table 6.
The optimal number of neurons in the hidden layer depends on the specific application, and it cannot, therefore, be identified a-priori. The ANN predictive performance is influenced to a great extent by the number of neurons. According to the literature, a too large number of neurons may lead to overfitting, while a too small number of neurons may make the ANN too simple to be able to capture the characteristics of a complex system (that is, underfitting [37]). Therefore, a trade-off must be made between the number of hidden layer neurons and the predictive performance of the ANN, and a sensitivity analysis is always needed in order to identify what the optimal number of neurons is for a specific application. In this section, a suitable number of neurons is selected for each of the dependent parameters that have to be estimated (i.e., those reported in Tables 4 and 5). All the available steady-state experimental data (410 tests) were used to generate the ANNs. Overall, 80% of the data was used as training data, and the other 20% was used to test the trained ANNs. Considering the smallness of the training dataset, the training process features some uncertainties, because the initial value of the weights of the neural networks is generated randomly. Thus, when the training process of a given neural network is repeated, its performance may vary to a great extent, even when the same input training dataset is used. Moreover, the experimental dataset for each training process is randomly divided into training and testing datasets with a certain ratio. This can also influence the performance of the trained neural networks when the training process is repeated. However, it has been shown in the literature that these uncertainties decrease considerably, and may even be ignored, when the size of the training dataset is sufficiently large [54], but this was not the case in this study.
Therefore, 200 training and testing repetitions were carried out for each ANN, with random initial weights and random splitting of the experimental dataset into training and testing datasets. The determination coefficient R 2 and the mean square error MSE were used to measure the performance of the ANNs. The mean values of R 2 and MSE for 200 trials were calculated for the training and testing of each ANN.
The trade-off between the number of hidden layer neurons and the performance of the ANN was explored for each ANN that simulated the parameters of the physics-based model (in the "ANN physics-based model") and for each ANN that simulated the MFB50, PFP, and BMEP combustion metrics (in the "direct ANN model"). The results of the sensitivity analysis are shown in Appendix B.
The sensitivity analysis allowed the optimal number of hidden layer neurons to be identified for each parameter. The most suitable number of neurons for each parameter is listed in Table 7, along with the mean values of R 2 and MSE for the training and testing phases. In general, it can be seen that the most suitable neuron number is no higher than 7 in all cases.
Once the sensitivity analysis had been carried out, the best ANN out of the 200 training trials was selected on the basis of the criteria reported in Appendix B for each parameter.
After this stage, the following investigations were conducted 1.
With reference to the physics-based model, the power law-based correlations of the calibration parameters were replaced by the corresponding ANNs, and the performance of the resulting model (i.e., the "ANN physics-based model") was evaluated. Table 8 shows a comparison between the accuracy of the ANN-based correlations and of the power-law based ones in the estimation of the parameters of the physics-based model. It can be observed that all the precisions based on ANNs are higher than those based on the power law correlations. ANNs have in fact been proved to have a powerful ability to catch the non-linear characteristics between input and output parameters 2.
The performance of the ANNs that were used to directly simulate MFB50, PFP, and BMEP was evaluated A detailed comparison of the performance of the four developed approaches, for the estimation of the MFB50, PFP, and BMEP combustion metrics, is reported in the next section.

Comparison of the Four Different Models under Steady-State and Transient Conditions
In this section, the performance of the four different models has been tested under steady-state and transient operating conditions.
First, the performance of the developed models was compared under steady-state conditions. The results are shown in Figure 6. It can be seen, from the previous charts, that the accuracy of the ANN physics-based model is higher than that of the baseline physics-based model concerning MFB50 and PFP. This indicates that the ANNs are more powerful than the power-function-based empirical functions in catching the nonlinear characteristics between the input and dependent parameters of the model for these two metrics under steady-state conditions. Overall, the accuracy of the direct ANN model is the best,  It can be seen, from the previous charts, that the accuracy of the ANN physics-based model is higher than that of the baseline physics-based model concerning MFB50 and PFP. This indicates that the ANNs are more powerful than the power-function-based empirical functions in catching the nonlinear characteristics between the input and dependent parameters of the model for these two metrics under steady-state conditions. Overall, the accuracy of the direct ANN model is the best, while the accuracy of the semi-empirical model is the lowest. However, the absolute values of the error of the semi-empirical model are still acceptable, considering that this kind of model ignores the detailed simulation of the in-cylinder physical process and requires a limited number of free parameters to be tuned (as a consequence, it also requires a limited number of experimental calibration tests). This suggests that the semi-empirical model, based on power functions, is also suitable for control-oriented combustion control applications, since its computational effort is much less than that of the physics-based model.
The performance of the developed models was then evaluated and compared under transient operating conditions over WHTC. The results are shown in Figures 7-9. Table 9 reports a summary of the accuracy of the models over WHTC and under steady-state operating conditions.  Table 9 that the accuracy of the direct ANN model is the best among all the different approaches, not only for steady-state operation (RMSE equal to 0.25 deg, 0.85 bar and 0.071 bar for MFB50, PFP, and BMEP, respectively), but also for transient operation over WHTC (RMSE of 1.1 deg, 9.6 bar and 0.7 bar for MFB50, PFP, and BMEP, respectively).
The A remarkable result is the deterioration of the accuracy of the ANN physics-based model (compared to the baseline physics-based model) which occurs for transient operation, while this approach leads to more accurate results for steady-state conditions. Therefore, this approach lacks robustness, compared to the baseline physics-based model.       Finally, the semi-empirical model shows a lower level of accuracy than the baseline physics-based model and the direct ANN model. However, the accuracy deterioration can still be considered acceptable, considering the simple structure of this model and the low calibration effort that is required, as discussed above.
Overall, the direct ANN model can be considered as the best candidate for control-oriented applications, since: • it is much simpler in structure and theory and easier to build than the physics-based model; • it does not require any detailed knowledge or modeling of the physical and chemical processes; • it is robust not only under steady-state conditions, but also for transient operation; • its accuracy is high even when a not so large dataset of experimental tests (410) is used for training; • it features the best accuracy under steady-state and transient operating conditions; • the required computational time is less than that of the physics-based model (see the next section).
The main drawback of ANNs is that they are a black-box type models, so they do not provide detailed information about the physical and chemical processes which occur inside the combustion chamber (e.g., heat release trend, in-cylinder pressure, temperatures . . . ).

Analysis of the Computational Time
The computational time of all the models (excluding the ANN physics-based model, which was found to be the least robust approach) was estimated when they were run on the ETAS ES910 device. The results are reported in Table 10. It can be seen that the direct ANN direct model and the semi-empirical model feature the lowest computational times, so that they can be considered the most suitable approaches to be implemented on an ECU for onboard control. The physics-based model requires the highest computational time. However, it is also compatible with real-time combustion control applications, since it is much lower than the typical engine cycle time. On the basis of the results reported in this study, a model-based combustion controller will be developed in the near future, and will be tested on the engine through the ETAS ES910 rapid prototyping device.

Conclusions
A comparison of the performance of physics-based, semi-empirical and artificial neural network (ANN)-based models has been carried out in this paper to estimate the main combustion metrics for an FPT Euro VI 3.0 L F1C diesel engine. The models were developed for combustion control-oriented applications. The combustion metrics that are predicted by the models are MFB50, PFP, and BMEP. These metrics are interesting for control-oriented applications, since they are closely related to the thermal efficiency and to the pollutant formation process of the engine. Four different types of modeling approaches have been investigated. The first approach is based on a previously developed low-throughput mean-value physics-based model, which is capable of simulating the heat release rate, the in-cylinder pressure, and the related metrics on the basis of the accumulated fuel mass principle. In this approach, the tuning parameters of the model are estimated on the basis of physically-consistent correlations, which are modeled using power-law functions. The second investigated approach is based on the joint use of artificial neural networks (ANNs) and physical modeling: instead of using the power-law functions, ANNs have been used to identify the parameters of the physics-based model. The third investigated approach is based on a direct estimation of MFB50, PFP, and BMEP on the basis of semi-empirical correlations. The fourth approach is based on the use of ANNs to directly estimate MFB50, PFP, and BMEP. The performance of the four models was evaluated under steady-state conditions and transient operation over WHTC.
Moreover, a new methodology, which is based on the sequential use of the Pearson correlation and partial correlation analysis, has been developed to identify the optimal set of input model parameters in order to a priori identify the correlation variables that should be excluded and the most robust correlation variables.
The main results can be summarized as follows: ( Overall, the direct ANN model can be considered the best candidate for control-oriented applications, since: • it is much simpler in structure and theory and easier to build than the physics-based model; • it does not require any detailed knowledge or modeling of the physical and chemical processes; • it is robust not only under steady-state conditions, but also in transient operation; • its accuracy is high, even when a not so large dataset of experimental tests (410) is used for training; • it features the best accuracy under steady-state and transient operating conditions • the required computational time is lower than that of the physics-based model.
The main drawback of this kind of model is that it is of the black-box type. Thus, it does not provide detailed information about the physical and chemical processes that occur inside the combustion chamber, such as the heat release trend, in-cylinder pressure, or temperatures. The authors would like to thank Omar Marello for his efforts in discussing and carrying out this research.

Conflicts of Interest:
The authors declare no conflict of interest. Greek Symbols ∆p exh-int difference between the exhaust and intake manifold pressure γ isentropic coefficient ρ density ρ IMF density in the inlet manifold ρ SOC charge density at the start of combustion ρ SOI charge density at the start of injection τ ignition delay coefficient

Appendix A Correlation Analysis for the Selection of the Model Input Variables
In order to make the correlation analysis more synthetic, the considered variables have been associated with indexing numbers (IN), as shown in Table A1.
For the sake of clarity, Figure A1 reports the temporal sequence of the calculation of the different variables used in the physics-based model. The Pearson correlation coefficient was calculated, using Matlab R2016b, in order to identify the correlations between the variables identified in Table A1. The results are reported in Figure A2. The numbers in the rows and columns in Figure A2 are the indexing numbers of the variables shown in Table A1. The background color changes from green to white and from white to red when the Pearson coefficient value varies in the (−1, 0) range or in the (0, 1) range, respectively. The linear relationship between two variables is considered strong when their Pearson coefficient value is in the (−1.0, −0.5) or (0.5, 1) ranges, moderate when it is in the (−0.5, −0.3) or (0.3, 0.5) ranges, weak when it is in the (−0.3, −0.1) or (0.1, 0.3) ranges, and non-existent or very weak when from it is in the (−0.1, 0.1) range [43]. The Pearson correlation coefficient was calculated, using Matlab R2016b, in order to identify the correlations between the variables identified in Table A1. The results are reported in Figure A2.
The numbers in the rows and columns in Figure A2 are the indexing numbers of the variables shown in Table A1. The background color changes from green to white and from white to red when the Pearson coefficient value varies in the (−1, 0) range or in the (0, 1) range, respectively. The linear relationship between two variables is considered strong when their Pearson coefficient value is in the • the input parameters of the correlations which are used to estimate the calibration parameters of the physics-based models (i.e., those reported in Table 3), or the input parameters of the ANNs that are used to estimate the same calibration parameters; • the input parameters of the semi-empirical correlations that are used to directly estimate MFB50, PFP, BMEP; • the input parameters of the ANNs that are used to directly estimate MFB50, PFP, BMEP.  Table A1.
As can be seen in Figure A2,  correlated to each other, since their crossed Pearson coefficient is almost equal to 1. Therefore, during the identification of the input variables of the models, it is sufficient to select only one of them as representative of each variable group. Moreover, it should be noted that the variables which are calculated at an earlier stage are expected to be closely correlated to the variables that are calculated in a subsequent stage (see Figure A1).
The following rules were identified to select the set of input variables to use in the correlations of the physics-based model or in the semi-empirical and ANN-based models, on the basis of the information derived from the Pearson correlation analysis: • the temporal sequences of the variables in the model have to be taken into account, i.e., variables cannot be used as the independent variables if they are calculated or occur after the dependent Generated Variables

Variables in Pressure Model
Pearson correlation coefficient  Table A1.
As can be seen in Figure A2,  ], are very closely correlated to each other, since their crossed Pearson coefficient is almost equal to 1. Therefore, during the identification of the input variables of the models, it is sufficient to select only one of them as representative of each variable group. Moreover, it should be noted that the variables which are calculated at an earlier stage are expected to be closely correlated to the variables that are calculated in a subsequent stage (see Figure A1). The following rules were identified to select the set of input variables to use in the correlations of the physics-based model or in the semi-empirical and ANN-based models, on the basis of the information derived from the Pearson correlation analysis: • the temporal sequences of the variables in the model have to be taken into account, i.e., variables cannot be used as the independent variables if they are calculated or occur after the dependent variable, according to the scheme in Figure A1; • the directly measured engine operation variables (namely N, p f , SOI, fuel injection quantity, p IMF , T IMF , Valve EGR ) should be considered as independent variable candidates with top priority, because they are the primary cause of the changes in engine operation. Moreover, these variables are usually very precise because they are measured directly in real-time. The generated parameters (ρ SOI , T SOI ), based on SOI, should be considered with higher priority than those based on SOC (ρ SOC , T SOC ), because of the uncertainty chain for the estimation of SOC on the basis of SOI; • the dependent variables (such as τ pil , τ main , SOC, K pil , K 1,main , K 2,main , Q f,evap , Q ht,glob , n, n', ∆p IMF , PMEP, FMEP) should be considered with a lower priority as correlation parameters, because they are usually less precise than the directly measured parameters. The X r,EGR , RAF, O 2 , ρ SOC and T SOC generated parameters should also be considered with a lower priority because of their relatively low precision; the total number of independent variable candidates should not be too high (in this study a maximum of 15 variables was chosen for the estimation of each dependent one), in order to ensure that the partial coefficient order is not too large and to make the model more robust to input uncertainty.
Some considerations should also been made concerning the m air and p EMF variables. Although m air is measured directly, it is usually not a precise or stable variable, since the accuracy of the engine air mass flow sensor is not very high (~5%), and its use, therefore, needs to be evaluated carefully. Moreover, many medium and heavy-duty engines are not generally equipped with air mass flow sensors and the exhaust manifold pressure is generally not measured in production engines, so its use as an input variable for the model should also be evaluated carefully (in this study it was avoided).
On the basis of the rules explained above, an analysis of the Pearson correlation coefficients was carried out to select independent variable candidates for each dependent variable. A summary is reported in Table A2.
It should be recalled that some of the dependent variables are the calibration parameters of the physics-based model (see Table 3), and they need to be modeled through correlations or ANNs, while other dependent variables (i.e., MFB50, PFP, BMEP) are the direct outcome of the semi-empirical model or of the direct ANN model, i.e., the metrics to be used for combustion control applications.
After the previous step, the partial correlation coefficient between each input variable of the initial candidate set reported in Table A2 and the corresponding dependent variable was calculated. The results are reported in Table A3, where the previously estimated Pearson correlation coefficients are also reported. The simultaneous analysis of the Pearson correlation and of the partial correlation coefficients is very useful, since it allows the least and most robust correlation variables to be identified. Table A2. Initial candidate set of the input variables for each output parameter, identified on the basis of the Pearson correlation analysis.

Dependent Variable Initial Candidate Set of the Independent Variables Selected Using Pearson Correlation Analysis
RAF N, p f , SOI pil , q tot,pil , q tot , p IMF , T IMF , X rEGR , ρ SOI,pil , ρ SOI,main τ pil N, p f , SOI pil , q tot,pil , p IMF , T IMF , Valve EGR , RAF, ρ SOI,pil , T SOI,pil τ main N, p f , SOI pil , q tot,pil , q tot , T IMF , m air , X rEGR , RAF, ρ SOI,pil , ρ SOI,main , ρ SOC,pil , T SOI,pil , ∆p exh−int , τ pil K pil N, p f , SOI pil , q tot,pil , p IMF , T IMF , Valve EGR , ρ SOI,pil , T SOI,pil , τ pil , ∆p exh−int K 1,main N, p f , SOI main , q tot,pil , q tot , p IMF , T IMF , Valve EGR , RAF, ρ SOI,main , ρ SOC,pil K 2,main N, p f , SOI pil , q tot,pil , q tot , T IMF , m air , X rEGR , RAF, ρ SOI,pil , ρ SOI,main , T SOI,main , K 1,main Q ht,glob N, p f , SOI pil , q tot,pil , q tot , T IMF , m air , X rEGR , RAF, ρ SOI,pil , ρ SOC,pil , τ main , In fact, input variables that are characterized by low Pearson and partial correlation coefficients can be removed from the candidate list, while input variables which are characterized by a high value of the Pearson correlation and the partial correlation coefficient at the same time can be selected as very robust correlation variables. Particular attention should also be paid to the input variables which, although showing a strong correlation with the dependent variables, are characterized by a limited precision (since they are derived parameters or measured parameters with low accuracy sensors). On the basis of the previous considerations, the least robust correlation variables are highlighted in Table A3 with a gray background, while the highly correlated variables are marked with an orange background, and finally the highly correlated variables characterized by a low precision are marked with a blue background.
It was found that some input variables (i.e., those marked with a white background in Table A3) are characterized by an intermediate value of the correlation coefficients. In order to understand whether they should have been included or excluded from the input candidate list, a sensitivity analysis was carried out, considering all the possible combinations of these variables and estimating the accuracy of the related models, as reported in Section 4.
As reported in the same section, it should be noted that the proposed approach, i.e., the preliminary selection of the input variables through Pearson and partial correlation coefficients, and a subsequent sensitivity analysis carried out on a reduced set of variables, is more computationally efficient than a pure sensitivity analysis approach which includes the entire set of input variables.   (Figures A3-A5) and the combustion metrics ( Figure A6). analysis was carried out, considering all the possible combinations of these variables and estimating the accuracy of the related models, as reported in Section 4.
As reported in the same section, it should be noted that the proposed approach, i.e., the preliminary selection of the input variables through Pearson and partial correlation coefficients, and a subsequent sensitivity analysis carried out on a reduced set of variables, is more computationally efficient than a pure sensitivity analysis approach which includes the entire set of input variables.

Appendix B-Training and Selection of the ANNs
Figures A3-A6 report the mean values of R 2 and MSE for 200 training trials of the ANNs that simulate the parameters of the physics-based model ( Figures A3, A4, A5) and the combustion metrics ( Figure A6). It can be noted, from the charts, that the training precision is generally higher than the testing precision, for both R 2 and MSE, in all cases. It is obvious that the training precision will always increase when the number of neurons increases, since the neural network is tested over the same dataset as that used for training. Moreover, it can be seen that the training precision increases slightly as the number of neurons increases when the number of neurons is above a certain value. As far as the testing precision is concerned, although some fluctuations occur, a rising trend of R 2 and a decreasing trend of MSE are also in general observed at the beginning when the number of neurons increases, and the trends then tend to become stable or even reverse. Reversing of the trend indicates overfitting. The precision of both the training and testing is quite low when the number of neurons is low, and this means that a too small number of neurons cannot capture the data characteristics accurately (underfitting). An interesting phenomenon can be seen with reference to the trends of the training precision of MFB50, PFP, and BMEP ( Figure A6), since they fluctuate less than the other parameters of the physics-based model (Figures A3-A5). This can be ascribed to the fact that the calibration precision of the ANNs that simulate MFB50, PFP, and BMEP is higher than that of the other parameters.    It can be noted, from the charts, that the training precision is generally higher than the testing precision, for both R 2 and MSE, in all cases. It is obvious that the training precision will always increase when the number of neurons increases, since the neural network is tested over the same dataset as that used for training. Moreover, it can be seen that the training precision increases slightly as the number of neurons increases when the number of neurons is above a certain value. As far as the testing precision is concerned, although some fluctuations occur, a rising trend of R 2 and a decreasing trend of MSE are also in general observed at the beginning when the number of neurons increases, and the trends then tend to become stable or even reverse. Reversing of the trend indicates  It can be noted, from the charts, that the training precision is generally higher than the testing precision, for both R 2 and MSE, in all cases. It is obvious that the training precision will always increase when the number of neurons increases, since the neural network is tested over the same dataset as that used for training. Moreover, it can be seen that the training precision increases slightly as the number of neurons increases when the number of neurons is above a certain value. As far as the testing precision is concerned, although some fluctuations occur, a rising trend of R 2 and a decreasing trend of MSE are also in general observed at the beginning when the number of neurons increases, and the trends then tend to become stable or even reverse. Reversing of the trend indicates Once the sensitivity analysis had been carried out, the best ANN out of the 200 trials was selected for each parameter. The automatic selection algorithm is shown hereafter (Algorithm A1): Where, i is the trial number; ANN_best is the best ANN selected by this algorithm; ANN_i is the trained ANN for i th trial; R2(i,1), R2 (i,2) and R2 (i,3) are the R 2 precisions for training, testing and total data validation, respectively; MSE (i,1), MSE (i,2) and MSE (i,3) are the MSE precisions for training, testing and total data validation, respectively; abs is the absolute value function.
Finally the best ANN out of the 200 trials is selected for each parameter, and the fitting results are compared with experimental data, as shown in Figures A7-A10. It can be observed that both the training and testing precisions are quite high (all the R 2 values are higher than 0.85) and there is not much difference between them. This indicates that all the selected ANNs present quite good fitting and prediction performances.