A Hybrid Physics-Based and Stochastic Neural Network Model Structure for Diesel Engine Combustion Events

Abstract

Estimation of combustion phasing and power production is essential to ensuring proper combustion and load control. However, archetypal control-oriented physics-based combustion models can become computationally expensive if highly accurate predictive capabilities are achieved. Artificial neural network (ANN) models, on the other hand, may provide superior predictive and computational capabilities. However, using classical ANNs for model-based prediction and control can be challenging, since their heuristic and deterministic black-box nature may make them intractable or create instabilities. In this paper, a hybridized modeling framework that leverages the advantages of both physics-based and stochastic neural network modeling approaches is utilized to capture CA50 (the timing when 50% of the fuel energy has been released) along with indicated mean effective pressure (IMEP). The performance of the hybridized framework is compared to a classical ANN and a physics-based-only framework in a stochastic environment. To ensure high robustness and low computational burden in the hybrid framework, the CA50 input parameters along with IMEP are captured with a Bayesian regularized ANN (BRANN) and then integrated into an overall physics-based 0D Wiebe model. The outputs of the hybridized CA50 and IMEP models are then successively fine-tuned with BRANN transfer learning models (TLMs). The study shows that in the presence of a Gaussian-distributed model uncertainty, the proposed hybridized model framework can achieve an RMSE of 1.3 × 10⁻⁵ CAD and 4.37 kPa with a 45.4 and 3.6 s total model runtime for CA50 and IMEP, respectively, for over 200 steady-state engine operating conditions. As such, this model framework may be a useful tool for real-time combustion control where in-cylinder feedback is limited.

1. Introduction

Stricter emissions regulations and renewed interest in reducing the carbon footprint of internal combustion engines have prompted the evolution of the fuel and air pathways of modern diesel engines, in order to maximize efficiency while minimizing exhaust emissions [1,2]. Some of the complex subsystems introduced in modern diesel engines include the turbocharger and high- and low-pressure exhaust gas recirculation loops [1]. High-fidelity models or combustion phasing feedback are essential to ensuring proper combustion management on these devices [3,4]. While some proposed modeling and control improvements require additional sensors [5,6], these changes would complicate the engine econometrics, and may not be economically viable [7]. Instead, powerful algorithms for online nonlinear model parameter estimation have made least-squares-based minimization methods feasible for real-time engine modeling and control [8,9,10,11]. These computationally efficient recursive algorithms may augment a model’s fidelity and improve the robustness of models so as to better handle system constraints and uncertainties. Recursive prediction approaches have also been shown to improve the predictive performance of ANN-based models [12]. Improvements in model accuracy could alleviate the need for expensive high-precision sensors in closed-loop control.

Over the years, a variety of control-oriented methods for modeling combustion phasing have been proposed, and there has been growing interest in utilizing data-driven techniques for modeling and control. Several researchers have demonstrated the use of online parameter optimization to capture or control combustion phasing in diesel engines [13,14]. In [13], the closed-loop control approach utilized online adjustment of the start-of-injection (SOI) timing and fuel rail pressure (RP) to achieve a desired CA50, and the authors showed that ANN-based control resulted in improved fuel consumption, nitrogen oxides (NOx), and soot emissions compared to a well calibrated look-up table. Tan et al. used extremum seeking and online optimization of engine model parameters to maximize IMEP [14]; an average error of 0.17 bar was achieved with this online-based strategy. The integration of other unique data-driven, empirical, and semi-empirical methods has also been explored in [15,16], where physics-based models were utilized to extend the range of machine learning predictions via the space transform. The method was effectively demonstrated for predictions of aerodynamic forces and hydraulic turbine efficiencies, as well as diesel engine calibration and combustion engine emissions. In [17], multiple ANN models were built with several input variables to capture an index of a complex transportation network. This approach was shown to improve the assessment of instabilities associated with the direct complex parameter capture encountered when using a single ANN model.

The use of ANNs has also been explored for the control of more complex engine architectures. In a recent study, Irdmousa led efforts to develop a data-driven model of combustion phasing for a reactivity-controlled compression ignition (RCCI) engine [18]. A linear parameter-varying support-vector machine (LPV-SVM) with multivariate inputs was used, and a maximum CA50 average tracking error of 0.7 CAD was achieved with an LPV-model-based controller. In [19], machine learning and combustion physics were integrated to predict CA50 and other combustion metrics for an RCCI engine. The random forests technique was implemented to develop correlative mapping between physics-based model coefficients and the engine’s operating condition, and an ANN was utilized for model coefficient determination. A standard deviation of 0.67 CAD, 1.19 CAD, and 0.223 bar was achieved for the predictions of the start of combustion (SOC), CA50, and IMEP, respectively. Similarly, Wang et al. [20] combined polynomial fitting functions and ANNs to predict diesel engine emissions and performance. A quadratic function was utilized to parameterize CA50, and the resulting model was used to optimize engine emissions and thermal efficiency. It was shown that with reasonable modifications to the models for wider engine operating conditions, the emissions prediction model could be used to achieve optimal performance in a combustion loop control. In addition, a semi-physical ANN-model-based prediction of the flame kernel development for a spark ignition (SI) engine was developed in [5], where critical physical parameters such as in-cylinder pressure, unburned gas density at spark timing, engine speed, fuel injected, and residual gas fraction were used as inputs to the ANN model. Model validation at 300 experimental steady-state operating points showed that a coefficient of determination (${\mathrm{R}}^2$) of 0.99 could be achieved. Recently, direct prediction of CA50 using an ANN has been utilized for combustion control in a heavy-duty natural gas engine [7], demonstrating that the need for an in-cylinder pressure sensor could be avoided. Excess air ratio, exhaust gas recirculation (EGR) ratio, intake manifold pressure, and spark timing were utilized as inputs to a 10-hidden-layer standard ANN model, which was trained with a Levenberg–Marquardt (LM) algorithm for a fixed speed test case of 1100 rpm with 0–30% EGR rate sweeps. An ${\mathrm{R}}^2$ of 0.998 was achieved in predicting CA50, and was deemed adequate for use in a sensorless combustion control.

While data-driven techniques can be utilized on their own, they can also be integrated with physics-based approaches to create various types of “grey box” models. For example, a direct classical ANN technique cascaded with a smooth variable structure filter (SVSF) for CA50 prediction was explored in conjunction with some physics-informed model inputs in [21]. The temperature, pressure, and oxygen fraction at intake valve closing (${\mathrm{T}}_{\mathrm{IVC}}$, ${\mathrm{P}}_{\mathrm{IVC}},{ \mathrm{F}}_{\mathrm{IVC}}$), SOI timing, and fuel mass were utilized as model inputs to the ANN model, and the CA50 estimation error was able to be reduced from 5 CAD to 1 CAD. In a previous work, the authors utilized an integrated physics-based and data-driven approach to capture the gas exchange processes of a modern diesel engine, and observed similar improvements in performance [1].

Further improvements in predictions have been observed with the incorporation of transfer learning approaches for combustion phasing modeling. In [22], a maximum mean discrepancy criterion was incorporated into an ANN to fine-tune the RCCI combustion phasing predictions from an LPV state–space combustion phasing model. A best fit rate (BFR) improvement in the CA50 was seen from 3.63%, without fine-tuning, to 85.52% BFR after fine-tuning. An improvement in IMEP BFR from 0% to 97.19% could be achieved with the transfer learning strategy when a new dataset with a different distribution was tested. Furthermore, in [23], an 11-input, 8-output Bayesian regularized (BR) 10-neuron hidden-layer ANN was coupled with a particle swarm optimization (PSO) algorithm to characterize the engine performance of an ethanol SI engine operating with a negative valve overlap. The performance indices, which were predicted with the joint BRANN–PSO network, were IMEP, NOx emissions, hydrocarbon emissions, carbon monoxide (CO) emissions, air flow, indicated efficiency, coefficient of variation of IMEP (${\mathrm{COV}}_{\mathrm{IMEP}}$), and indicated specific fuel consumption (ISFC). An ${\mathrm{R}}^{2 }$ above 0.93 was achieved for all of the output predictions. A final example can be seen in [24], where a multizone thermokinetic engine model was coupled with an ANN and a genetic algorithm to predict CA50, IMEP, CO, and total unburned hydrocarbons (THC); average errors of 1.2 CAD, 0.4 bar, 10 PPM, 0.8%, and 394 PPM were obtained for the predictions, respectively.

The literature demonstrates that leveraging physics to inform ANN models could enhance predictions of engine metrics and, as such, improve real-time combustion phasing control methods without the strict need for in-cylinder sensors. While prior work using either direct ANNs or integrated ANN and physics-based approaches to predict critical engine metrics has shown promise [1,5,6,7,13,15,16,18,19,20,21,22,23,24], some challenges with uncertainty and complexity still need to be addressed. Some existing methods assume linear relationships and, as such, are simpler to use for control, but have a more limited ability to extrapolate, and encounter greater inaccuracy due to the nonlinearity of engine processes [7,13,18,22]. Others have complex submodels of effects such as heat transfer loss and oxygen fraction that increase the computational burden [14,21]. In addition, methods that leverage ANNs can encounter instabilities and time-consuming heuristic model development [12,17].

In this study, these issues are addressed by segmenting the combustion phasing model into a lower model dimension via physical knowledge and then capturing the most influential and more deterministic combustion phasing input parameters with BRANNs [25]. The reduced model dimension allows for a hybridized physics-based CA50 prediction with a simple nonlinear model adaptation. The hybridized model outputs are also then successively fine-tuned with a BRANN ad hoc transfer learning model (TLM). This approach should allow model uncertainty to be better quantified in a deterministic and parametric way, since the ANN is integrated within a physical domain. In addition, the intelligent adaptation approach used in this work provides a simpler and more deterministic way to handle system constraints, since the input–output relationship of the combustion phasing can be explicitly identified in a physical domain. In this study, an intelligent adaptation is incorporated in a simple 0D Wiebe combustion model to capture CA50 using an operating-point-dependent optimization algorithm that circumvents excessive modeling and computational complexities.

This study presents several novel contributions to the area of combustion modeling. The first contribution is an intelligent adaptation of a simple 0D Wiebe combustion model to capture CA50, while the second is the introduction of a method of characterizing the more deterministic combustion phasing input parameters with ANNs within a physics-based structure. These additions provide modeling simplicity that has not been achieved in the literature, and these methods may enable reduced modeling effort and lower computational power to be required for combustion control. Moreover, as previously mentioned, Wang utilized an ANN to capture SOC for CA50 prediction in an SI engine [5], but the CA50 model needed independent modelling of the in-cylinder pressure as a model input, which can be challenging to capture accurately in the presence of uncertainties [1]. The combustion phasing model developed in this study avoids the need for an in-cylinder pressure model. This approach may provide a cost-effective way to implement combustion control, since fewer measurement devices would be required. A third contribution of this work is the exploration of the use of the BRANN model for combustion phasing modeling. Unlike in [22], the BRANN model is utilized in this work to develop the transfer learning model, which avoids unnecessary model validation. This reduces modeling effort and computation time, and provides high accuracy. In addition, the model fine-tuning can be done both offline and online. The handling of uncertainties by “Bayesifying” the ANNs and successively transferring the predictions from the hybridized framework to an ad hoc BRANN model also ensures robustness for real-time application. In summary, the hybridized model framework developed in this work provides robustness, along with the advantage of efficiently analyzing and accurately capturing the underlying physical and chemical interactions in the combustion phasing event, while utilizing fewer measurements from stock engine sensors.

The remainder of this paper is presented as follows: Section 2 details the experimental setup and testing conditions. Section 3 presents the individual model methodologies and validation. Section 4 compares the results of the various models, and Section 5 concludes the study and offers recommendations for future work.

2. Experimental Setup

Experimental studies for this work were conducted on a 2013 Volkswagen 2.0 L compression ignition engine equipped with a variable geometry turbocharger (VGT), cooled low-pressure exhaust gas recirculation (EGR), and high-pressure EGR. The engine system is shown in Figure 1. Flow through the air path on this system is dictated by the VGT vane nozzle position, low-pressure and high-pressure EGR valve positions, throttle valve position, and an exhaust flap that controls the system backpressure. Fresh air flow into the engine is measured by a Sierra Instruments Air-Trak system. The engine is also equipped with a common rail fuel injection system. Fuel flow is measured by a Sierra Instruments Fuel Trak system, and is also conditioned to keep a constant temperature and pressure.

The engine speed and torque are simultaneously measured by an eddy current dynamometer, and the data are monitored via the Sierra CADET software, which is installed on a computer as shown in Figure 2. Load and speed control are achieved by an intermediary engine modular box that interlinks the ECU and the Sierra CADET data monitoring software. In addition to the engine speed measurement, both a stock engine speed sensor—which is routed to the ECU—and a high-precision AVL crank-angle encoder are used to measure the engine speed. The engine speed data are monitored using the ECU software and AVL IndiCom. All data of other essential engine parameter sensors are obtained through the Sierra Data Acquisition System (DAQ) and monitored using the Sierra CADET software, as illustrated in Figure 2.

Data for the model creation and validation were collected at over 200 operating points that included variations in engine speed and torque, VGT and low-pressure EGR actuator settings, exhaust flap position, and charge air cooler (CAC) settings. During all tests, the throttle valve was kept fully open, and the high-pressure EGR valve was kept closed. Only low-pressure EGR use was explored in this work. The high-pressure EGR is used in a limited operating range and, as such, was not leveraged. At all operating points, a single fuel pulse was used, and the start of injection was altered to maintain a CA50 of 8 degrees after top dead center (aTDC). The actual CA50 and the IMEP were monitored in AVL IndiCom by using the conditioned pressure from a high-resolution Kistler pressure transducer mounted in cylinder three, along with the crank-angle signal obtained from an AVL crank-angle encoder. Figure 2 shows the interactions of the various measurement systems used to obtain data for this study.

The range of engine operating conditions explored in this work was scattered across the engine operation map, and included 212 operating conditions. The normal distribution plots of the data collected for engine speed, torque, VGT actuator settings, low-pressure EGR actuator settings, exhaust backflow flap valve settings, and the charge air cooler (CAC) settings are depicted in Figure 3. The normal distribution plots of the data collected for engine speed, torque, VGT actuator settings, low-pressure EGR actuator settings, and exhaust backflow flap valve settings show a quasi-normal distribution of the operating points. On the other hand, the charge air cooler (CAC) settings show a right-skewed distribution, as is consistent with the test protocol.

The engine is rated at a maximum speed of 4500 rpm and a maximum power output of 103 kW at 4000 rpm. The maximum torque is also rated at 320 Nm at 1750 rpm up to 2500 rpm.

Table 1. Range of engine testing conditions.

Parameter	Minimum	Maximum
Engine Speed	2000 rpm	4000 rpm
Engine Torque	50 Nm	150 Nm
Intake Manifold Pressure	13.6 psi	27.9 psi
Intake Manifold Temperature	306 K	395 K
Start of Injection	0° bTDC	10.2° bTDC

shows the range of conditions covered in the experiments.

3. Model Development

The data collected were used to train and validate three different modeling approaches to capture CA50 along with IMEP. These models are as follows:

• A multilayer perceptron (MLP) ANN for capturing CA50 with various engine input parameters, including the experimental start of combustion (${\mathrm{SOC}}_{\mathrm{Experm}}$) and burn duration (${\mathrm{BD}}_{\mathrm{Experm}}$). In this model framework, IMEP is also characterized by a classical MLP ANN using similar but fewer input parameters than the CA50 ANN;
• A physics-based model that utilizes parameterized physical SOC and BD submodels as the main model inputs to an adaptive simple Wiebe function for CA50 prediction. The Wiebe function is given by Equation (1), in which the SOC and BD are the main input parameters [4]:

  $$\mathrm{x}\left(\mathsf{\theta }\right)=1-{\mathrm{e}}^{-\mathrm{a}{\left(\frac{\mathsf{\theta }-\mathrm{SOC}}{\mathrm{BD}}\right)}^{\mathrm{n}}}$$

  (1)

  where the term $\mathrm{x}\left(\mathsf{\theta }\right)$ is the burn fraction over crankshaft positions $\mathsf{\theta }$, and the Wiebe shape factors $\mathrm{a}$ and $\mathrm{n}$ are parameterized via MATLAB’s $fmincon$ function to provide high-accuracy predictions across the engine’s operating range. A physics-based IMEP model was also developed. The physics-based IMEP model utilizes engine geometrical characteristics, fuel properties, in-cylinder temperature models, and measured fuel mass flow as the main model inputs;
• A hybrid model that utilizes ANN SOC and BD submodels with the same adaptive simple Wiebe function used in the physics-based approach. In this model framework, the SOC and BD are captured with a Bayesian regularized multilayer perceptron (BRMLP) artificial neural network and then fed into the adaptive Wiebe model, which is then transferred to an ad hoc BRANN transfer model to further fine-tune its predictions. In this hybrid framework, the IMEP is also characterized with a BRANN along with the SOC and BD models, but it is identified in parallel to the CA50 prediction, and also successively fine-tuned with its ad hoc BRANN TLM. Therefore, the hybrid framework features a decoupled multiple-input, single-output (MISO) strategy. While most inputs to the SOC BRANN and BD BRANN models are actuator signals, some input parameters to these submodels also need to be modeled or measured directly.

The subsections below describe the development of the ANN CA50 and IMEP models, the physics-based CA50 and IMEP models, and the hybridized CA50 and IMEP models. To ensure consistency in the flow of the models’ description, the direct ANN CA50 model’s development is discussed first, followed by that of the direct ANN IMEP model. The physics-based subsections describe the individual SOC and BD submodels for the physics-based CA50 prediction, followed by those for the physics-based IMEP. The subsection describing the hybridized CA50 and IMEP model framework discusses the SOC BRANN and BD BRANN models for the hybridized CA50 prediction as well as its BRANN transfer learning model. The IMEP-hybridized BRANN model is also discussed, followed by its transfer learning model.

3.1. ANN Models

As discussed previously, one promising way in which combustion phasing could be captured is with an ANN.

3.1.1. ANN CA50 Model

In this first approach, capturing the CA50 involved the use of direct measurements of the input parameters to train an MLP ANN. In this model framework, a standard approach was utilized for the ANN weight parameterization, using a Levenberg–Marquardt (LM) training algorithm. The mathematical implementation of the MLP and LM algorithm is described in the authors’ previous work [1]. The LM algorithm was utilized due to its known high computational power, low CPU memory usage, and high performance [26]. Figure 4 illustrates the CA50 and IMEP classical ANN model structure used in this work.

The direct CA50 ANN model was created using MATLAB’s deep learning toolbox. In the model setup, a large set of available engine variables were initially all used as inputs. These input parameters were then pruned to an optimal size using a sensitivity analysis. This was achieved by setting the mean square error (MSE) between the experimental values (CA50) and the predicted values ($\widehat{\mathrm{CA}50}$) as the optimization target. The model was then run 50 times to record the outcome. Next, while maintaining a fixed input, the model was again run 50 times but with the hyperparameters such as the hidden-layer size being varied with a one-step increment for every 50 runs until there was no significant reduction in the MSE. The optimal set of input variables and hidden-layer size were selected as those that gave the least MSE and lowest early-stopping epoch. The resulting hidden-layer size had 30 neurons. The variable set included nine variables, allowing CA50 to be expressed as follows:

$$\mathrm{CA}{50}_{\mathrm{ANN}}=\mathrm{f}\left({\mathsf{\delta }\mathsf{\Psi }}_{\mathrm{air}},{\mathsf{\delta }\mathsf{\Gamma }}_{\mathrm{fuel}},{ \mathsf{\delta }\mathsf{\Omega }}_{\mathrm{intake}}, \mathrm{N}, \mathrm{Torq},{ \mathrm{SOC}}_{\mathrm{Experm}},{ \mathrm{BD}}_{\mathrm{Experm}},{ \mathrm{T}}_{\mathrm{EM}},{\mathrm{P}}_{\mathrm{EM}}\right)$$

(2)

where ${\mathsf{\delta }\mathsf{\Psi }}_{\mathrm{air}}$ represents a particular combination of the air-handling actuator signals of the turbocharger’s vane nozzle position ($\mathrm{VNP})$, the low-pressure EGR valve position $\left(\mathrm{EVP}\right)$, and the exhaust flap valve position $\left(\mathrm{EFP}\right)$; the term ${\mathsf{\delta }\mathsf{\Gamma }}_{\mathrm{fuel}}$ represents the combination of the corresponding fuel-path actuator signals of the $\mathrm{SOI}$, injector electrical delay $({\mathsf{\tau }}_{\mathrm{elec}})$, and the fuel mass flow $({\mathrm{M}}_{\mathrm{fuel}})$; ${\mathsf{\delta }\mathsf{\Omega }}_{\mathrm{intake}}$ is the composite combination of the intake gas dynamics due to pressure, temperature, and mass changes; $\mathrm{N}$ is the engine speed; $\mathrm{Torq}$ is the engine torque; ${\mathrm{SOC}}_{\mathrm{Experm}}$ is the experimental start of combustion; ${\mathrm{BD}}_{\mathrm{Experm}}$ is the experimental burn duration; and ${\mathrm{T}}_{\mathrm{EM}}$ and ${\mathrm{P}}_{\mathrm{EM}}$ are the exhaust manifold temperature and pressure, respectively.

A total of 212 data points were utilized for the overall CA50 ANN model’s development, of which 80% of the data were used to train the model, 10% were used for validation, and the remaining 10% were used for model testing. Figure 5 shows the engine speed and torque combinations of the training, validation, and testing datasets used for the CA50 ANN model’s development.

Figure 5 shows the distribution of the data used for the training, validation, and testing subroutines of the CA50 ANN model. The data sectioning was randomized so as to enhance the model’s regularization. The desired model accuracy was reached with early stopping at 17 epoch iterations in ~0.28 s. The CA50 ANN model returned a Pearson product–moment correlation coefficient (PPMCC) of 81.8% and an RMSE of 0.27 CAD. Figure 6 shows the regression analysis of the CA50 ANN model’s results.

The first quadrant of the ANN CA50 model regression plots in Figure 6 shows the training performance, which has a reasonably high correlative measure (R = 0.95368). The second quadrant shows the validation of the trained model with data that it was not trained on, and shows reasonable generalization ability (R = 0.73327). The third quadrant of the plot shows the testing validation of the trained ANN CA50 model, with a different withheld dataset other than the training and validation datasets. The testing performance showed a lower but a statistically significant correlation measure (R = 0.15942). The last quadrant in Figure 6 provides the overall correlative measure (R = 0.81755) of the ANN CA50 model. The correlative metrics, along with the RMSE of 0.27 CAD, indicate a reasonably well-performing ANN CA50 model.

3.1.2. ANN IMEP Model

The IMEP ANN model followed a similar development routine as with the CA50 ANN model. The model was quite simple, utilizing only 6 main input parameters and 15 hidden-layer nodes. The mathematical function given in Equation (3) describes the general ANN IMEP structure. As with the CA50 ANN model, 80% of the data were used for training, 10% were used for model validation, and the remaining 10% were used for model testing.

$${\mathrm{IMEP}}_{\mathrm{ANN}}=\mathrm{f}\left({\mathsf{\delta }\mathsf{\Psi }}_{\mathrm{air}},{\mathsf{\delta }\mathsf{\Gamma }}_{\mathrm{fuel}},{ \mathsf{\delta }\mathsf{\Omega }}_{\mathrm{intake}}, \mathrm{N},{ \mathrm{T}}_{\mathrm{EM}},{\mathrm{P}}_{\mathrm{EM}}\right)$$

(3)

Figure 7 shows the classical IMEP ANN model validation and illustrates the excellent accuracy achieved by the ANN in these cases. The IMEP classical ANN model returned a testing set PPMCC of $99.37\%$ and an RMSE of $6.33 \mathrm{kPa}$ with a $0.20$ s computational runtime in 13 epoch iterations. The model was seen to be stable when random initializations were performed.

3.2. Physics-Based Models

A common way of predicting CA50 is with a physics-based approach. As such, physics-based models were developed to serve as a baseline for comparison. The CA50 was calculated using a Wiebe function that requires underlying SOC and BD models.

3.2.1. Physics-Based SOC Model

The physics-based SOC model used in this work was modified from that of [27], and is shown in (4). The SOC was taken to be the algebraic sum of the start of injection command ($\mathrm{SOI})$, delay due to the physical injector electrodynamics (${\mathsf{\tau }}_{\mathrm{elec}}$), and the ignition delay (${\mathsf{\tau }}_{\mathrm{ID}}$).

$${\mathrm{SOC}}_{\mathrm{phys}}= \mathrm{SOI}+{\mathsf{\tau }}_{\mathrm{elec}}+{ \mathsf{\tau }}_{\mathrm{ID}}$$

(4)

The electrical delay (${\mathsf{\tau }}_{\mathrm{elec}}$) takes into account the thermodynamic loss angle shift in the determination of the TDC by the AVL 365C encoder, as reported by the AVL IndiCom software, and was determined to be $1 \mathrm{CAD}$. The chemical ignition delay model (${\mathsf{\tau }}_{\mathrm{ID}}$) was modified from the form proposed in [28]. This model considers the ignition delay to be impacted by effects at three characteristic time periods representing the low (${\mathsf{\tau }}_{1\mathrm{r}}$), medium (${\mathsf{\tau }}_{2\mathrm{r}}$), and high (${\mathsf{\tau }}_{3\mathrm{r}}$)-temperature regions of ignition via Arrhenius-type equations. Due to the heterogeneous nature of diesel combustion [29], an additional term capturing the effects of turbulence (${\mathsf{\tau }}_{\mathsf{\psi }}$) was introduced to enhance the model’s accuracy [30]. This was achieved by capturing the influence of the intake dynamics and the piston movement via the mean piston speed $\left({\mathrm{S}}_{\mathrm{P}}\right)$ and the intake manifold pressure $\left({\mathrm{P}}_{\mathrm{IM}}\right)$. Equation (5) shows the general form for ignition delay, while Equation (6) shows the final model in which the underlying expressions are given:

$${\mathsf{\tau }}_{\mathrm{ID}}=\frac{1}{{\mathsf{\tau }}_{1\mathrm{r}}+{\mathsf{\tau }}_{2\mathrm{r}}}+\frac{1}{{\mathsf{\tau }}_{3\mathrm{r}}}+{\mathsf{\tau }}_{\mathsf{\psi }}$$

(5)

$${\mathsf{\tau }}_{\mathrm{ID}}=\frac{1}{{\mathrm{c}}_1{\mathsf{\varphi }}^{{\mathrm{c}}_3}{\mathrm{P}}_{\mathrm{IM}}^{{\mathrm{c}}_2}\exp \left(\frac{{\mathrm{c}}_{10}}{{\mathrm{T}}_{\mathrm{IM}}}\right)+{\mathrm{c}}_4{\mathrm{T}}_{\mathrm{IM}}^{{\mathrm{c}}_5}{\mathrm{P}}_{\mathrm{IM}}^{{\mathrm{c}}_6}\exp \left(\frac{{\mathrm{c}}_{11}}{{\mathrm{T}}_{\mathrm{IM}}}\right) }+\frac{1}{{\mathrm{c}}_7{\mathsf{\varphi }}^{{\mathrm{c}}_8}{\mathrm{P}}_{\mathrm{IM}}^{{\mathrm{c}}_9}\exp \left(\frac{{\mathrm{c}}_{12}}{{\mathrm{T}}_{\mathrm{IM}}}\right) }+({\mathrm{c}}_{13}{\mathrm{P}}_{\mathrm{IM}}+{\mathrm{c}}_{14}{\mathrm{S}}_{\mathrm{P}}+{\mathrm{c}}_{15})$$

(6)

In (6), ${\mathrm{T}}_{\mathrm{IM}}$ is the intake manifold temperature, and $\mathsf{\varphi }$ is the equivalence ratio. The ${\mathrm{c}}_{\mathrm{i}}$ terms are model constants that were found using the generalized reduced gradient (GRG) algorithm. The overall physics-based SOC model is found by combining (5) and (6), as follows:

$${\mathrm{SOC}}_{\mathrm{Phys}}= \mathrm{SOI}+\frac{1}{{\mathrm{c}}_1{\mathsf{\varphi }}^{{\mathrm{c}}_3}{\mathrm{P}}_{\mathrm{IM}}^{{\mathrm{c}}_2}\exp \left(\frac{{\mathrm{c}}_{10}}{{\mathrm{T}}_{\mathrm{IM}}}\right)+{\mathrm{c}}_4{\mathrm{T}}_{\mathrm{IM}}^{{\mathrm{c}}_5}{\mathrm{P}}_{\mathrm{IM}}^{{\mathrm{c}}_6}\exp \left(\frac{{\mathrm{c}}_{11}}{{\mathrm{T}}_{\mathrm{IM}}}\right) }+\frac{1}{{\mathrm{c}}_7{\mathsf{\varphi }}^{{\mathrm{c}}_8}{\mathrm{P}}_{\mathrm{IM}}^{{\mathrm{c}}_9}\exp \left(\frac{{\mathrm{c}}_{12}}{{\mathrm{T}}_{\mathrm{IM}}}\right) }+({\mathrm{c}}_{13}{\mathrm{P}}_{\mathrm{IM}}+{\mathrm{c}}_{14}{\mathrm{S}}_{\mathrm{P}}+{\mathrm{c}}_{15})+{\mathsf{\tau }}_{\mathrm{elec}}.$$

(7)

The ignition delay model is the main non-actuator component of the generalized parametric physics-based SOC model given by (6); as such, its analysis reflects the accuracy of the entire SOC model. The ignition delay model was observed to produce accurate predictions for all 212 engine operating points, with only a few exceptions. In fact, the model’s fidelity had a PPMCC of 97.4% and an RMSE of 1.02 CAD. A graphical representation of this model’s accuracy is shown in Figure 8.

This model is fairly accurate, but also complex. As such, the model’s 97.4% PPMCC was offset by the larger computational runtime of the parameter optimization. It can be seen in Equation (6) that the entire SOC model utilizes 15 decision variables, whose values are given in

Table 2. Tuning constants for the physics-based SOC model.

Model Constants	Value
${\mathrm{c}}_1$	$4.47\times {10}^{-9}$
${\mathrm{c}}_2$	$6.58$
${\mathrm{c}}_3$	$1.67$
${\mathrm{c}}_4$	$0.907$
${\mathrm{c}}_5$	$-0.0499$
${\mathrm{c}}_6$	$-0.972$
${\mathrm{c}}_7$	$3.53$
${\mathrm{c}}_8$	$-0.126$
${\mathrm{c}}_9$	$-9.21$
${\mathrm{c}}_{10}$	$-124$
${\mathrm{c}}_{11}$	$8.53$
${\mathrm{c}}_{12}$	$0.997$
${\mathrm{c}}_{13}$	$0.672$
${\mathrm{c}}_{14}$	$-1.24$
${\mathrm{c}}_{15}$	$-1.14$

. The SOC model development showed that the GRG algorithm utilized over 200 iterations for an approximately 230 s computational runtime in the global optimization subroutine.

Figure 9 shows a comparison of the resulting overall physics-based SOC model to the actual SOC. The prediction points with higher errors appear to have operating conditions with characteristics that were not well accounted for in the overall SOC parametric physics-based model.

Predictions that had errors outside the 90% confidence bounds were operating conditions at the extremes of the air-handling actuation signals for a certain engine speed and torque. For example, the case with 100% VNP, 100% EVP, 0% EFP at 3000 RPM, and 67 Nm torque had an absolute error of ~2.53 CAD. Furthermore, conditions that had inadequate air breathing characteristics and, thus, a tendency to misfire also had higher errors, such as 0% VNP, 60% EVP, 60% EFP at 2000 RPM, and 11.35 Nm, which had the maximum absolute error of ~4.01 CAD. It is also possible that, coupled with process noise, the GRG algorithm could not converge to a true global solution.

3.2.2. Physics-Based BD Model

The physics-based BD model was developed by parameterizing the physical contribution of the laminar flame speed and the rate of diffusion and late combustion in the burn profile [31,32,33]. A modified form of the van Tiggelen and Tabaczynski model capturing the laminar flame speed’s contribution was utilized, and is shown in Equation (8). The laminar flame speed model utilized the temperature at TDC (${\mathrm{T}}_{\mathrm{TDC}}$), fuel mass $\left({\mathrm{M}}_{\mathrm{fuel}}\right)$, engine charge flow (${\mathrm{W}}_{\mathrm{engine}}$) [1], and exhaust manifold temperature (${\mathrm{T}}_{\mathrm{EM}}$) as model inputs. The contribution of diffusion and late combustion, also depicted by Equation (9), was modeled empirically as a functional interaction of the SOC, SOI, equivalence ratio $\left(\mathsf{\varphi }\right)$, excess air ratio ($\mathsf{\lambda }$), intake manifold temperature $\left({\mathrm{T}}_{\mathrm{IM}}\right)$, and the air-handling actuator positions (VNP, EVF, and EVP) expressed in percentages. To improve the robustness of the model, the measured excess air ratio was used, and the equivalence ratio was calculated from the air and fuel flow measurements.

$${\mathsf{\Delta }\mathsf{\theta }}_{\mathrm{flame}}={\mathrm{d}}_1+{\left[{\mathrm{d}}_2{\mathrm{T}}_{\mathrm{TDC}}\sqrt{\frac{1}{{\mathrm{T}}_{\mathrm{EM}}}{\left[{\mathrm{M}}_{\mathrm{fuel}}\right]}^{{\mathrm{d}}_3}{\left[{\mathrm{W}}_{\mathrm{engine}}\right]}^{{\mathrm{d}}_4}\exp \left(\frac{-{\mathrm{d}}_5}{{\mathrm{d}}_6{\mathrm{T}}_{\mathrm{EM}}}\right) }\right]}^{\mathrm{d}7}$$

(8)

$$\begin{gathered} {\mathsf{\Delta }\mathsf{\theta }}_{\mathrm{diff}+\mathrm{lc} }={\left({\mathrm{d}}_8{\mathsf{\theta }}_{\mathrm{SOC}}\mathsf{\varphi }\right)}^2+\left({\mathrm{d}}_9{\mathsf{\theta }}_{\mathrm{SOI}}\right)+{\left({\mathrm{d}}_{10}{\mathsf{\theta }}_{\mathrm{SOC}}{\mathsf{\theta }}_{\mathrm{SOI}}\right)}^2+\left({\mathrm{d}}_{11}+{\mathrm{d}}_{12}\mathsf{\varphi }+{\mathrm{d}}_{13}{\mathsf{\varphi }}^{{\mathrm{d}}_{14}}\right) \\ +\left({\mathrm{d}}_{15}{\mathrm{T}}_{\mathrm{IM}}{\mathsf{\lambda }}^{{\mathrm{d}}_{16}}\right)+\left({\mathrm{d}}_{17}\mathrm{VNP}·\mathrm{EVF}·\mathrm{EVP}\right) \end{gathered}$$

(9)

The final burn duration model, given by Equation (10), was evaluated as the algebraic sum of the laminar flame speed’s contribution and the contribution of diffusion and late combustion. The ${\mathrm{d}}_{\mathrm{i}}$ terms are model constants that were also found using the GRG algorithm.

Table 3. Tuning constants for the physics-based BD model.

Model Constants	Value
${\mathrm{d}}_1$	$3$
${\mathrm{d}}_2$	$5.03\times {10}^{95}$
${\mathrm{d}}_3$	$-12.1$
${\mathrm{d}}_4$	$48.31$
${\mathrm{d}}_5$	$1.91\times {10}^5$
${\mathrm{d}}_6$	$3.48$
${\mathrm{d}}_7$	$0.0223$
${\mathrm{d}}_8$	$-0.000221$
${\mathrm{d}}_9$	$-0.254$
${\mathrm{d}}_{10}$	$3.73\times {10}^{-5}$
${\mathrm{d}}_{11}$	$-234$
${\mathrm{d}}_{12}$	$-5.50$
${\mathrm{d}}_{13}$	$2355$
${\mathrm{d}}_{14}$	$0.02795$
${\mathrm{d}}_{15}$	$0.0008125$
${\mathrm{d}}_{16}$	$2.555$
${\mathrm{d}}_{17}$	$4.47\times {10}^{-7}$

shows the values of the ${\mathrm{d}}_{\mathrm{i}}$ terms in the optimized physics-based BD model.

$$\begin{gathered} {\mathrm{BD}}_{\mathrm{Phys}}={ \mathrm{d}}_1+{\left[{ \mathrm{d}}_2{\mathrm{T}}_{\mathrm{TDC}}\sqrt{\frac{1}{{\mathrm{T}}_{\mathrm{EM}}}{\left[{\mathrm{M}}_{\mathrm{fuel}}\right]}^{{\mathrm{d}}_3}{\left[{\mathrm{W}}_{\mathrm{engine}}\right]}^{{\mathrm{d}}_4}\exp \left(\frac{-{\mathrm{d}}_5}{{\mathrm{d}}_6{\mathrm{T}}_{\mathrm{EM}}}\right) }\right]}^{\mathrm{d}7}+{\left({\mathrm{d}}_8{\mathsf{\theta }}_{\mathrm{SOC}}\mathsf{\varphi }\right)}^2+\left({\mathrm{d}}_9{\mathsf{\theta }}_{\mathrm{SOI}}\right)+ \\ {\left({\mathrm{d}}_{10}{\mathsf{\theta }}_{\mathrm{SOC}}{\mathsf{\theta }}_{\mathrm{SOI}}\right)}^2+\left({\mathrm{d}}_{11}+{\mathrm{d}}_{12}\mathsf{\varphi }+{\mathrm{d}}_{13}{\mathsf{\varphi }}^{{\mathrm{d}}_{14}}\right)+\left({\mathrm{d}}_{15}{\mathrm{T}}_{\mathrm{IM}}{\mathsf{\lambda }}^{{\mathrm{d}}_{16}}\right)+\left({\mathrm{d}}_{17} \mathrm{VNP}·\mathrm{EVF}·\mathrm{EVP}\right) \end{gathered}$$

(10)

The burn duration model has higher errors than the SOC model, as illustrated in Figure 10. The parametric physics-based BD model returns a PPMCC of $97.0\%$ and an RMSE of $1.77$ CAD. The GRG algorithm utilized over 906 iterations for a 42.1 s computational runtime in the parametric optimization subroutine.

3.2.3. Physics-Based CA50 Model

After the physics-based SOC and BD models were developed, they were then used as inputs to the CA50 model, which is derived from Equation (1) and can be expressed as follows:

$$\mathrm{CA}{50}_{\mathrm{Phys}}\left(\mathrm{k}\right)={\mathrm{SOC}}_{\mathrm{Phys}}\left(\mathrm{k}\right)+{\left[\frac{\ln 2 }{{\mathrm{a}}_{\mathrm{wieb}}{}_{\mathrm{k}}}\right]}^{\frac{1}{{\mathrm{n}}_{{\mathrm{wieb}}_{\mathrm{k}}}}}{ \mathrm{BD}}_{\mathrm{Phys}}\left(\mathrm{k}\right).$$

(11)

Figure 11 shows the model structure of the physics-based CA50 and IMEP models.

The Wiebe function is a simple and widely accepted 0D combustion model, which is highly amenable for control; however, it generally lacks high predictive capability [31], and cannot adapt to changes that may occur as an engine ages [8]. Consequently, the Wiebe function becomes less effective if it is globally parameterized for the ECU [14,31]. In order to improve its predictive robustness for wider engine operating conditions, the Wiebe function was made adaptable by the use of a recursive least-squares model parameter optimization. In Equation (11), $\mathrm{k}$ represents a particular engine operating condition. For each k operating point, ${\mathrm{a}}_{\mathrm{wieb}}$ and ${\mathrm{n}}_{\mathrm{wieb}}$ (which are the Wiebe shape factors) are updated by minimizing the cost function given by (12), in which $\mathrm{CA}{50}_{\mathrm{target}} \left(\mathrm{k}\right)$ is the desired optimal trajectory. The minimization objective is realized through MATLAB’s $fmincon$ nonlinear model parameter optimization function.

$${\mathsf{\Theta }}_{\mathrm{Phys}}\left(\mathrm{k}\right)=\parallel \mathrm{CA}{50}_{\mathrm{target}}\left(\mathrm{k}\right)-\mathrm{CA}{50}_{\mathrm{Phys}}\left(\mathrm{k}\right){\parallel }^2$$

(12)

The validation of the resulting model predictions of the physics-based CA50 model is shown in Figure 12.

The physics-based CA50 model returned a PPMCC of $~95.0\%$ and an RMSE of $0.17$ CAD with a $49.4$ s computational runtime using an iteration convergence tolerance of $1\times {10}^{-99}$ for the online parameter optimization subroutine in an unperturbed environment.

3.2.4. Physics-Based IMEP Model

The IMEP quantifies the work output of an engine [29]. Due to the computational complexities involved in in-cylinder pressure modeling for a larger engine operating range [1,34], the physics-based IMEP model (${\mathrm{IMEP}}_{\mathrm{Phys}})$ was similarly parametrized as in [34]. As shown in Equation (13), it was expressed as a function of displaced volume (${\mathrm{V}}_{\mathrm{dis}})$, fuel mass flow rate $({\mathrm{M}}_{\mathrm{fuel}})$, intake valve closing temperature (${\mathrm{T}}_{\mathrm{IVC}})$, exhaust valve closing temperature $({\mathrm{T}}_{\mathrm{EVC}}),$ start-of-combustion temperature $({\mathrm{T}}_{\mathrm{SOC}})$, end-of-combustion temperature $\left({\mathrm{T}}_{\mathrm{EOC}}\right)$, and the fuel-specific heat at constant volume $({\mathrm{c}}_{\mathrm{v}}$).

$${\mathrm{IMEP}}_{\mathrm{Phys}}=\mathrm{f}\left({\mathrm{M}}_{\mathrm{fuel}},{\mathrm{c}}_{\mathrm{v}},{\mathrm{V}}_{\mathrm{dis}},{\mathrm{T}}_{\mathrm{IVC}},{\mathrm{T}}_{\mathrm{SOC}},{\mathrm{T}}_{\mathrm{EOC}},{\mathrm{T}}_{\mathrm{EVC}}\right)$$

(13)

${\mathrm{T}}_{\mathrm{IVC}}$ was assumed to be the same as the measured intake manifold temperature $\left({\mathrm{T}}_{\mathrm{IM}}\right)$, and the other model arguments ${\mathrm{T}}_{\mathrm{SOC}}$, ${\mathrm{T}}_{\mathrm{EOC}}$, and ${\mathrm{T}}_{\mathrm{EVC}}$ were modeled as given in Equations (14)–(17); they are only briefly discussed in this work. More details of their modeling development can be found in [34]. ${\mathrm{T}}_{\mathrm{SOC}}$ is found by assuming isentropic compression via Equation (14).

$${\mathrm{T}}_{\mathrm{SOC}}={\mathrm{T}}_{\mathrm{IVC}}·{\left(\frac{{\mathrm{V}}_{\mathrm{IVC}}}{{\mathrm{V}}_{\mathrm{SOC}}}\right)}^{{\mathrm{n}}_{\mathrm{comp}}-1}$$

(14)

where ${\mathrm{V}}_{\mathrm{IVC}}$ and ${\mathrm{V}}_{\mathrm{SOC}}$ are the in-cylinder volume at intake valve closing and the start of combustion, respectively; these were determined using the slider–crank theory volume model [29]. The term ${\mathrm{n}}_{\mathrm{comp}}$ is the polytropic compression index obtained from a simple polytropic pressure fitting of the predicted and measured compression pressure phase.

${\mathrm{T}}_{\mathrm{EOC}}$ is described as follows:

$${\mathrm{T}}_{\mathrm{EOC}}={\mathrm{T}}_{\mathrm{SOC}}+\left({\mathrm{a}}_0+{\mathrm{a}}_1·{\mathsf{\theta }}_{\mathrm{SOC}}+{\mathrm{a}}_2·{\mathsf{\theta }}_{\mathrm{SOC}}^2\right)·\mathsf{\Delta }\mathrm{T}$$

(15)

in which $\mathsf{\Delta }\mathrm{T}$ accounts for heat transfer and is given by:

$$\mathsf{\Delta }\mathrm{T}=\frac{\mathrm{LHV}·\left(\frac{{\mathrm{M}}_{\mathrm{fuel}}}{{\mathrm{W}}_{\mathrm{engine}}}\right)·\mathsf{\varphi }}{{\mathrm{c}}_{\mathrm{v}}·\left(\frac{{\mathrm{M}}_{\mathrm{fuel}}}{{\mathrm{W}}_{\mathrm{engine}}}·\mathsf{\varphi }+1\right)}.$$

(16)

where $\mathrm{LHV}$ is the fuel’s lower heating value. The ${\mathrm{T}}_{\mathrm{EOC}}$ parameter is then used to compute ${\mathrm{T}}_{\mathrm{EVC}}$ via:

$${\mathrm{T}}_{\mathrm{EVC}}=\left[{\mathrm{T}}_{\mathrm{EOC}}·{\left(\frac{{\mathrm{V}}_{\mathrm{EOC}}}{{\mathrm{V}}_{\mathrm{EVO}}}\right)}^{{\mathrm{n}}_{\exp }-1}\right] ·{\left(\frac{{\mathrm{V}}_{\mathrm{EVO}}}{{\mathrm{V}}_{\mathrm{exh}}}\right)}^{{\mathrm{n}}_{\exp }-1 }.$$

(17)

The model parameters ${\mathrm{a}}_0,{ \mathrm{a}}_1$, and ${\mathrm{a}}_2$ are tuning constants found by model calibration, and whose values are discussed later in this subsection, while ${\mathrm{n}}_{\exp }$ is the expansion polytropic index found similarly as for ${\mathrm{n}}_{\mathrm{comp}}$. The ${\mathrm{V}}_{\mathrm{EOC}}$, ${\mathrm{V}}_{\mathrm{EVO}}$, and ${\mathrm{V}}_{\mathrm{exh}}$ parameters are the volume at the end of combustion, exhaust valve opening, and exhaust manifold, respectively. The resulting model arguments are then fed into the model in Equation (18) to obtain the physics-based IMEP model.

$${\mathrm{IMEP}}_{\mathrm{Phys}}={\mathrm{C}}_1·{\mathrm{M}}_{\mathrm{fuel}}\frac{{\mathrm{c}}_{\mathrm{v}}}{{\mathrm{V}}_{\mathrm{dis}}}\left({\mathrm{T}}_{\mathrm{IVC}}-{\mathrm{T}}_{\mathrm{SOC}}+{\mathrm{T}}_{\mathrm{EOC}}-{\mathrm{T}}_{\mathrm{EVC}}\right)$$

(18)

where ${\mathrm{C}}_1$ is a model tuning constant that is also found by model optimization.

In this study, two modeling optimization approaches were explored for the physics-based IMEP model: First, an online parameterization approach that utilized a recursive least-squares algorithm to update each model constant for each operating condition using MATLAB’s fmincon was explored. Figure 13 shows the variations of the various modeling constants obtained in the adaptive physics-based IMEP model.

The second method was a static parameterization approach that utilized the GRG algorithm to find global tuning constants for the whole operating point set. For this second approach, a regressor (based on similar argumentation to that given in [35]) given by Equation (19) was developed to globally fit the constants obtained from the online parameterization. The ${\mathrm{G}}_{\mathrm{i}}$ terms are the resulting global model constants obtained for the physics-based IMEP model.

$${\mathrm{M}}_1={\mathrm{G}}_1·{\left[{\mathrm{G}}_3{\mathrm{N}}^{{\mathrm{G}}_4}+{\mathrm{G}}_5{\mathsf{\varphi }}^{{\mathrm{G}}_6}+{\mathrm{G}}_7{\mathrm{W}}_{\mathrm{engine}}^{{\mathrm{G}}_8}+{\mathrm{G}}_9{\left(\frac{{\mathrm{M}}_{\mathrm{fuel}}}{{\mathrm{W}}_{\mathrm{engine}}}\right)}^{{\mathrm{G}}_{10}}+{\mathrm{G}}_{11}{\mathrm{M}}_{\mathrm{fuel}}^{{\mathrm{G}}_{12}}+{\mathrm{G}}_{13}{\left({\mathrm{W}}_{\mathrm{engine}}^{{\mathrm{G}}_{15}}{\mathsf{\varphi }}^{{\mathrm{G}}_{16}}{\mathrm{N}}^{{\mathrm{G}}_{17}}{\left(\frac{{\mathrm{M}}_{\mathrm{fuel}}}{{\mathrm{W}}_{\mathrm{engine}}}\right)}^{{\mathrm{G}}_{18}}·{\mathrm{M}}_{\mathrm{fuel}}^{{\mathrm{G}}_{19}}\right)}^{\mathrm{G}14}\right]}^{{\mathrm{G}}_2}$$

(19)

The ${\mathrm{G}}_{\mathrm{i}}$ terms and their resulting optimized values are shown in

Table 4. Global tuning constants for the physics-based IMEP model.

Global Constants	${\mathit{C}}_1 \mathbf{Constants}$	${\mathit{a}}_0 \mathbf{Constants}$	${\mathit{a}}_1 \mathbf{Constants}$	${\mathit{a}}_2 \mathbf{Constants}$
${\mathrm{G}}_1$	$-1.12022$	$8.18030$	$1.98573$	$191433.3$
${\mathrm{G}}_2$	$0.11789$	$-2.07501$	$-0.01369$	$541.240$
${\mathrm{G}}_3$	$-0.27385$	$0.76698$	$-0.16477$	$0.10354$
${\mathrm{G}}_4$	$-9.24762$	$-0.64384$	$-2.38677$	$-6.12314$
${\mathrm{G}}_5$	$0.31595$	$659.915$	$0.46292$	$0.24445$
${\mathrm{G}}_6$	$0.00388$	$4.02354$	$-2.10743$	$0.09418$
${\mathrm{G}}_7$	$0.08581$	$2407.80$	$6.00093$	$-128.950$
${\mathrm{G}}_8$	$-0.01172$	$0.65975$	$-0.06542$	$-217.089$
${\mathrm{G}}_9$	$-0.41197$	$1.01783$	$3.8952\times {10}^{-5}$	$0.74714$
${\mathrm{G}}_{10}$	$0.00637$	$-0.00201$	$7.7920\times {10}^{-9}$	$-0.00163$
${\mathrm{G}}_{11}$	$4098.27$	$1$	$-0.78875$	$0.56064$
${\mathrm{G}}_{12}$	$2.15874$	$1$	$20.12791$	$5.49783$
${\mathrm{G}}_{13}$	$0.01003$	$1.26320$	$25.18099$	$-0.18746$
${\mathrm{G}}_{14}$	$2.8246\times {10}^{-4}$	$-0.00950$	$21.49473$	$-0.01356$
${\mathrm{G}}_{15}$	$-0.36759$	$1\times {10}^{-99}$	$-2.00973$	$17.99333$
${\mathrm{G}}_{16}$	$-1.09780$	$1\times {10}^{-99}$	$1.15165$	$-122.763$
${\mathrm{G}}_{17}$	$-7.57761$	$1\times {10}^{-99}$	$-3.66506$	$35.90610$
${\mathrm{G}}_{18}$	$1.72289$	$1\times {10}^{-99}$	$2.41961$	$-7.61119$
${\mathrm{G}}_{19}$	$0.94137$	$1\times {10}^{-99}$	$7.51263$	$17.95118$

.

Figure 14 compares the performance of the adaptive and static physics-based IMEP models with the experimental results, and demonstrates that similar results are achieved with both approaches in most cases.

In an unperturbed environment, the adaptive physics-based IMEP model returned a PPMCC of 100% and an RMSE of $1.79\times {10}^{-6}$ kPa in a 50.9 s computational runtime, while the more computationally expensive static approach returned a $95.16\%$ PPMCC and an RMSE of 17.7 kPa. Table 4 shows that a significant number of decision variables (76 parameters) had to be tracked for the static approach. This approach could be very useful for use by the ECU during engine operating conditions that may exhibit saturation [36], or when excessive model adaptation may be detrimental for engine control.

3.3. Hybridized Models

The hybridized model framework for the CA50 and IMEP features an integration of the ANN and physics-based models.

3.3.1. Hybridized CA50 Model

The hybridized CA50 model utilized ANN submodels of the input parameters (SOC and BD) to the simple 0D Wiebe combustion model. The development of the ANN SOC and BD models will first be discussed, followed by an examination of the hybrid Wiebe model.

ANN SOC Model

An ANN model was created to characterize the SOC. To ensure model robustness while integrating the SOC ANN model into the hybridized framework, a stochastic neural network was utilized. Stochastic neural networks involve the development of an ANN with probabilistic transfer functions or weights [37]. This technique ensures a model’s high robustness to epistemic uncertainties, especially when there are limited data. The SOC ANN was developed using MATLAB’s deep learning toolbox, and followed a similar heuristic development to the direct classical ANN models developed in Section 3.1. However, the SOC ANN model featured a Bayesian regularized multilayer perceptron (BRMLP) artificial neural network with nine input nodes, five hidden-layer nodes, and one output node. The BRMLP architecture was chosen because it provides excellent regularization performance coupled with simplicity. As such, the MLP should be able to capture complex dynamics, and the Bayesian regularization (BR) will adequately handle epistemic uncertainties [1,10,38]. A brief mathematical description of the BRANN is given below.

The Bayesian regularized ANN is based on the premise that some prior distributions can be imposed on the ANN model parameters in order to shrink large weights in the anticipation of achieving a more generalized infinitesimal mapping error between the estimated and target parameters [38]. Therefore, given the training dataset $\mathrm{D}=\left\{\mathrm{SOC}, {\left({\mathrm{p}}_{\mathrm{i}}\right)}_{\mathrm{i}=1,\dots ,\mathrm{n}}\right\}$ of the inputs $({\mathrm{p}}_{\mathrm{i}})$ as a vector of $\mathrm{i}$ samples of the target variable $(\mathrm{SOC})$, a sum of squared estimation errors ${\mathrm{E}}_{\mathrm{D}}\left(\mathrm{D}|\mathrm{w},\mathrm{M}\right)$ that maps $\mathrm{n}$ inputs to the target variable with some model coefficients ($\mathrm{w}$) for a particular ANN architecture $(\mathrm{M})$ can be defined as follows:

$${\mathrm{E}}_{\mathrm{D}}\left(\mathrm{D}|\mathrm{w},\mathrm{M}\right)={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{SOC}}_{\mathrm{i}}-{\widehat{\mathrm{SOC}}}_{\mathrm{i}}\right)}^2$$

(20)

where ${\widehat{\mathrm{SOC}}}_{\mathrm{i}}$ is the estimated variable. Equation (20) is minimized with early stopping regularization, and the BRANN technique adds a sum of squares of the architecture weights term ${\mathrm{E}}_{\mathrm{w}}\left(\mathrm{w}|\mathrm{M}\right)$, which penalizes large weights in expectation of an optimal input–output functional mapping for unseen datasets. A gradient-based optimization technique is used to minimize the resulting function given in Equation (21), which equates to penalizing a log-likelihood [38].

$$\mathrm{F}={ \mathsf{\beta }\mathrm{E}}_{\mathrm{D}}\left(\mathrm{D}|\mathrm{w},\mathrm{M}\right)+{ \mathsf{\alpha }\mathrm{E}}_{\mathrm{w}}\left(\mathrm{w}|\mathrm{M}\right)$$

(21)

In Equation (21), $\mathsf{\beta }$ and $\mathsf{\alpha }$ are hyperparameters, which are estimated adaptively. The weight decay coefficient $(\mathsf{\alpha })$ in the weight decay term ${\mathsf{\alpha }\mathrm{E}}_{\mathrm{w}}\left(\mathrm{w}|\mathrm{M}\right)$ has been shown to favor small values of $\mathrm{w}$, which then maximizes the generalization ability of the model [38,39]. The posterior densities of the weights are near zero for large Bayes estimates of $\mathsf{\alpha }$, so the weights effectively disappear and the model slashes connections in the network. Therefore, $\mathsf{\alpha }$ and $\mathsf{\beta }$ are adaptively predicted to balance model complexity and goodness of fit [38]. The strength connections ($\mathrm{w}$) of the network are randomly initialized for training. After the data are taken, the density function for the weights can be updated according to Bayes’ rule. The Bayes’ posterior distribution of $\mathrm{w}$ given $\mathsf{\alpha },\mathsf{\beta },\mathrm{D}$, and $\mathrm{M}$ is depicted by Equation (22).

$$\mathrm{P}\left(\mathrm{w}|\mathrm{D},\mathsf{\alpha },\mathsf{\beta },\mathrm{M}\right)=\frac{\mathrm{P}\left(\mathrm{D}|\mathrm{w},\mathsf{\beta },\mathrm{M}\right)\mathrm{P}(\mathrm{w}|\mathsf{\alpha },\mathrm{M})}{\mathrm{P}(\mathrm{D}|\mathsf{\alpha },\mathsf{\beta },\mathrm{M})}$$

(22)

$\mathrm{P}\left(\mathrm{D}|\mathrm{w},\mathsf{\beta },\mathrm{M}\right)$ is the likelihood function of $\mathrm{w}$, and $\mathrm{P}\left(\mathrm{w}|\mathsf{\alpha },\mathrm{M}\right)$ is the prior distribution of weights under $\mathrm{M}$, which is the expectation of achieving the dataset D given $\mathrm{w}$. $\mathrm{P}\left(\mathrm{D}|\mathsf{\alpha },\mathsf{\beta },\mathrm{M}\right)$ expressed by Equation (23) is the normalization factor or evidence for hyperparameters $\mathsf{\alpha }$ and $\mathsf{\beta }$, and is independent of $\mathrm{w}$ [38].

$$\mathrm{P}\left(\mathrm{D}|\mathsf{\alpha },\mathsf{\beta },\mathrm{M}\right)=\int \mathrm{P}(\mathrm{D}|\mathrm{w},\mathsf{\beta },\mathrm{M})\mathrm{P}\left(\mathrm{w}|\mathsf{\alpha },\mathrm{M}\right)\mathrm{dw}$$

(23)

In order to maximize the marginal likelihood of achieving the dataset D with respect to the hyperparameters $\mathsf{\alpha }$ and $\mathsf{\beta }$, the normalization factor needs to be maximized. In this study, the weights ($\mathrm{w})$ were assumed to be identically distributed, with each being admissible in a Gaussian distribution $\left(\mathsf{\alpha },\mathrm{M}\right)~ \mathrm{N}\left(0,{\mathsf{\alpha }}^{-1}\right)$. The integral given in Equation (23) is intractable; therefore, a Laplacian approximation to the marginal density can be computed and utilized as follows:

$$\log \left[\mathrm{P}\left(\mathrm{D}|\mathsf{\alpha },\mathsf{\beta },\mathrm{M}\right)\right] \approx \mathrm{K}+\frac{\mathrm{n}}{2}\log \left(\mathsf{\beta }\right)+\frac{\mathrm{m}}{2}\log \left(\mathsf{\alpha }\right)-{\left|\mathrm{F}\left(\mathsf{\alpha },\mathsf{\beta }\right)\right|}_{\mathrm{w}={\mathrm{w}}_{\left(\mathsf{\alpha },\mathsf{\beta }\right)}^{\mathrm{MAP}}}-\frac{1}{2}\log \Vert \mathsf{H}\Vert { }_{\mathrm{w}={\mathrm{w}}_{\left(\mathsf{\alpha },\mathsf{\beta }\right)}^{\mathrm{MAP}}}$$

(24)

where K is a constant and $\mathsf{H}=\frac{{\partial }^2}{\partial \mathrm{w}\partial \mathrm{w}\prime }\mathrm{F}\left(\mathsf{\alpha },\mathsf{\beta }\right)$ is the Hessian matrix. An iterative update of the hyperparameters is evaluated at the previous-step values of the tuning parameters given in Equations (25) and (26) with (${\mathrm{w}}^{\mathrm{MAP}}$)—the maximum a posteriori estimate of the early-stopping objective function, F. Utilizing a Gauss–Newton approximation to the Hessian matrix, an effective number of parameters, $\mathsf{\gamma }=\mathrm{n}-2{\mathsf{\alpha }}^{\mathrm{MAP}}\mathrm{tr}{\left({\mathsf{H}}^{\mathrm{MAP}}\right)}^{-1}$, can be computed and passed on to compute the new estimates of $\mathsf{\alpha }$ and $\mathsf{\beta }$.

$${\mathsf{\alpha }}_{\mathrm{new}}=\frac{\mathsf{\gamma }}{2{\mathrm{E}}_{\mathrm{w}}\left({\mathrm{w}}^{\mathrm{new}}\right)}$$

(25)

$${\mathsf{\beta }}_{\mathrm{new}}=\frac{\mathrm{n}-\mathsf{\gamma }}{2{\mathrm{E}}_{\mathrm{D}}\left({\mathrm{w}}^{\mathrm{new}}\right)}$$

(26)

More details on the BRANN mathematical derivation can be found in [38]. BRANN therefore uses the LM algorithm for computing the posterior modes in BR. As mentioned earlier, the LM algorithm provides high robustness, and has the ability to converge at a local minimum even when there is a high absolute initial deviation from the optimal value trajectory [40,41].

The BRANN model’s development was carried out in a similar manner as that of the direct CA50 ANN model. While there is no single best approach to any ANN’s development [25], the heuristic method used required relatively lower computational complexity due to the localized analysis of a less complex output parameter resulting from the reduction in the overall combustion phasing model’s dimension. The BRANN SOC model was run 20 times to achieve input parameter pruning, and it was observed to be stable when trained with random initializations. The desired model accuracy was reached with early stopping at only nine epoch iterations. In summary, a mathematical description of the ANN SOC model is given by (27) as follows:

$${\mathrm{SOC}}_{\mathrm{ANN}}=\mathrm{f}\left({\mathsf{\delta }\mathsf{\Psi }}_{\mathrm{air}},{\mathsf{\delta }\mathsf{\Gamma }}_{\mathrm{fuel}},{ \mathrm{W}}_{\mathrm{engine}}, \mathrm{N},{\mathrm{T}}_{\mathrm{EM}},{\mathrm{P}}_{\mathrm{EM}}\right)$$

(27)

A total of 212 data points were utilized for the overall SOC ANN model’s development, of which 80% of the data were used to train the model, and 20% were used for model testing. The BRANN SOC model returned a PPMCC of $99.3\%$ and an RMSE of $0.45$ CAD. The BRANN SOC model had a $~0.77$ s computational runtime. Figure 15 shows the model validation of the BRANN SOC model.

ANN BD Model

Similar to the ANN SOC model, the ANN BD model was developed using the BRMLP ANN in MATLAB’s deep learning toolbox. The ANN BD model featured 10 input nodes, 30 hidden-layer nodes, and 1 output node. The model’s development followed a similar process to that of the ANN SOC model, but the ANN SOC model’s output was also utilized as an input to the ANN BD model. Again, 212 data points were utilized for the model’s development, of which 80% were used for training and 20% for testing. The resulting condensed mathematical representation of the BD ANN model is given below in Equation (28):

$${\mathrm{BD}}_{\mathrm{ANN}}=\mathrm{f}\left({\mathsf{\delta }\mathsf{\Psi }}_{\mathrm{air}},{\mathsf{\delta }\mathsf{\Gamma }}_{\mathrm{fuel}},{ \mathrm{SOC}}_{\mathrm{ANN}},{ \mathrm{W}}_{\mathrm{engine}}, \mathrm{N},{ \mathrm{T}}_{\mathrm{EM}},{\mathrm{P}}_{\mathrm{EM}}\right)$$

(28)

The BRANN BD model had a $~0.32$ s runtime and accuracy of $~99.0\%$ PPMCC, as illustrated in Figure 16.

After the BRANN SOC and BD models were developed, they were then passed on to the adaptive physics-based Wiebe model to predict the CA50 given by Equation (29):

$$\mathrm{CA}{50}_{\mathrm{Hyb}-\mathrm{Wiebe}} \left(\mathrm{k}\right)={\mathrm{SOC}}_{\mathrm{ANN}}\left(\mathrm{k}\right)+{\left[\frac{\ln 2 }{{\mathrm{a}}_{\mathrm{wieb}}{}_{\mathrm{k}}}\right]}^{\frac{1}{{\mathrm{n}}_{{\mathrm{wieb}}_{\mathrm{k}}}}}{ \mathrm{BD}}_{\mathrm{ANN} }\left(\mathrm{k}\right)$$

(29)

The interconnection between these submodels is shown in Figure 17 for the hybrid model, and the performance of the CA50 hybrid Wiebe model in an unperturbed environment is shown in Figure 18. The hybridized CA50 Wiebe model returned a PPMCC of 100% and an RMSE of $1.13\times {10}^{-7}$ CAD at a 48.4 s computational runtime for all 212 steady states’ operating conditions.

3.3.2. Hybridized IMEP Model

A hybridized IMEP model was also developed and compared to the ANN and physics-based IMEP models described in Section 3.1. The hybridized IMEP model followed a similar development routine to that of the ANN IMEP in Section 3.1, but the hybrid model featured the use of a BRANN. The charge flow rate $({\mathrm{W}}_{\mathrm{engine}})$ obtained via the charge flow model was also used as a model input instead of the direct input of the intake gas dynamics measurements. The resulting functional representation of the hybrid IMEP model is given in Equation (30), and the model structure is shown in Figure 19.

$${\mathrm{IMEP}}_{\mathrm{Hyb}-\mathrm{BRANN}}=\mathrm{f}\left({\mathsf{\delta }\mathsf{\Psi }}_{\mathrm{air}},{\mathsf{\delta }\mathsf{\Gamma }}_{\mathrm{fuel}},{ \mathrm{W}}_{\mathrm{engine}}, \mathrm{N},{\mathrm{T}}_{\mathrm{EM}},{\mathrm{P}}_{\mathrm{EM}}\right)$$

(30)

As demonstrated in Figure 20, the hybridized IMEP model is very accurate, with a testing set PPMCC of $99.5\%$ and an RMSE $4.94 \mathrm{kPa}$. The model had a $0.71$ s computational runtime in 17 epoch iterations in an unperturbed environment.

4. Impact of Variations and Uncertainty

While each of these model options works, their performance can also be affected by variations and uncertainty in a real-world environment. As such, the performance of the three CA50 and IMEP model frameworks was also evaluated in a stochastic environment. The models were perturbed with a Gaussian distributed disturbance (as shown in the block diagrams in Figure 4, Figure 11, Figure 17 and Figure 19) to assess their robustness.

4.1. CA50 Perturbation

The ANN CA50 model, the physics-based CA50 model, and the hybrid CA50 model were perturbed with additive Gaussian noise. The perturbation signal was applied as an external additive noise to the SOC and BD models in the physics-based and hybrid CA50 models, and also as an internal additive noise to the model arguments of the SOC and BD models, in order to assess the robustness of the physics-based and hybrid CA50 models. The effects of the perturbations were assessed based on the output of the CA50 models. Figure 21 shows the resulting CA50 model performances with the application of the perturbation signal, both externally and internally.

When perturbing the model states with Gaussian noise, it was observed that the direct classical ANN CA50 model did not deviate significantly from the target values. On the other hand, excessive model deterioration took place when either internal or external perturbation was included in the parametric physics-based SOC or BD models. This could result in egregious CA50 physics prediction errors as they propagate into the online parameterization, as shown in Figure 21. The internal perturbation of the physics-based SOC and BD models resulted in a CA50 model PPMCC of 6.94% and an RMSE of 15.18 CAD. The external perturbation of the physics-based SOC and BD models resulted in a CA50 model PPMCC of 70.65% and an RMSE of 0.48 CAD. As such, the perturbation exercise clearly showed the degradation of the physics-based CA50 model framework due to the global parameterization characteristics of the SOC and BD physics-based submodels.

In the internal perturbation of the BRANN SOC and BD models with Gaussian noise for the hybridized CA50 model’s framework, a CA50 model PPMCC of 100% and an RMSE of $8.8\times {10}^{-8}$ CAD were observed. This exceptional model performance (as depicted in Figure 21) in the presence of input uncertainties could be explained by the high extrapolative capabilities of the highly regularized BRANN SOC and BD models in the hybridized CA50 model framework. It could be seen that the additive Gaussian noise causing drifts in the input parameters for the BRANN SOC and BD models did not result in any change in the 100% PPMCC that was found for the unperturbed case. However, a decrement in RMSE from the unperturbed case was observed for the internally perturbed case. This is due to the ability of the BRANN models to better predict numerically precise values for the online parameter optimization algorithm with drifting help from the Gaussian noise. This is consistent with the prior Gaussian distribution assumption made on the hyperparameters in the development of the BRANN. On the other hand, the PPMCC for the hybridized CA50 model framework during external perturbation of the BRANN SOC and BD models was observed to have dropped to ~92.1% with an increase in the RMSE to 0.25 CAD. Although a 7.88% decrease in accuracy was observed in the external perturbed case for the hybridized CA50 model, the resulting correlation coefficient was still within a 90% confidence interval, which may still make the hybridized CA50 model framework an acceptable modeling strategy for real-time control, regardless of how the system is perturbed. To ensure consistency in the robustness, an ad hoc model was created to correct the model performance, and is discussed below.

CA50 BRANN Transfer Learning Model (TLM)

Due to possible model deterioration from uncertainties in the hybridized framework, an ad hoc transfer learning model (TLM) was developed with a BRANN to ensure the hybridized model’s robustness in a stochastic environment. The hybridized model framework’s outputs from the adaptive Wiebe model were treated as transitional dummy output steps and then passed on to the TLM as a subroutine for a final model fine-tuning step. This could be achieved by minimizing the error between the Wiebe model predictions and the target optimal trajectory, as shown in Equation (31). This realization was based on the neural network transfer learning model’s premise that, given stochastic variations in the distributions in the oracle of an ANN, the stored knowledge from previous distributions can be utilized to enhance the predictions in the presence of stochastic perturbations [22]. To implement this, an ad hoc BRANN model was developed by utilizing the transitional $\mathrm{CA}{50}_{\mathrm{Hyb}-\mathrm{Wiebe}}$ perturbed output as a model argument and the optimal CA50 trajectory as the model target. As such, given changes in the stochastic environment, a model re-substitution is carried out via less extensive ad hoc retraining. The transfer learning method utilized in this work is simple, since the already robust output of the hybridized model framework—the $\mathrm{CA}{50}_{\mathrm{Hyb}-\mathrm{Wiebe}}$—is used as a model input, and is highly likely to be very close to the optimal trajectory $\mathrm{CA}{50}_{\mathrm{target}}$. Therefore, the mapping between these two variables would be close to affine, which then becomes computationally highly tractable to optimize. In the event that high nonlinearity in the mapping occurs from excessive nonlinear model perturbations, the utilization of the BRANN is still useful and, as such, a high prediction accuracy is invariably guaranteed.

$${\mathsf{\Theta }}_{\mathrm{final}}\left(\mathrm{k}\right)=\parallel \mathrm{CA}{50}_{\mathrm{target}}\left(\mathrm{k}\right)-\left\{\mathrm{CA}{50}_{\mathrm{Hyb}- \mathrm{TLM}}{}^{\left(\mathrm{k}\right)}+{\vartheta }_{\mathrm{Gaussian}}\right\}{\parallel }^2$$

(31)

In Equation (31), ${\vartheta }_{\mathrm{Gaussian}}$ represents the attributed external stochastic changes in the environment, and ${\mathsf{\Theta }}_{\mathrm{final}}$ becomes the new overall objective function that accounts for stochastic variations in the modeling environment. Computational tractability for real-time implementation was ensured in this step by utilizing the BRANN, which circumvents unnecessary model validation steps. In addition, only 50% of the $\mathrm{CA}50{ }_{\mathrm{Hyb}-\mathrm{Wiebe}}$ outputs were used for the ad hoc training, and the remaining 50% for model testing. The final overall CA50 model was then condensed to a simple relation, given as follows:

$$\mathrm{CA}{50}_{\mathrm{Hyb}- \mathrm{TLM}} \left(\mathrm{k}\right)=\mathrm{f}\left(\mathrm{CA}{50}_{\mathrm{Hyb}-\mathrm{Wiebe}}{}^{\left(\mathrm{k}\right)}\right).$$

(32)

Figure 22 shows the model regression of the CA50 BRANN TLM. The regression analysis indicates that R = 1 for the training and testing subroutines. This shows that the TLM had achieved its model fidelity augmentation objective. The high R = 1 correlative measure seen in the plots is consistent with the linearly separable or near-affine premise invoked for the output of the hybridized Wiebe model and the target CA50.

4.2. IMEP Perturbation

A similar perturbation exercise was carried out for the IMEP models. Due to the nature of the IMEP model frameworks, only internal model perturbation was feasible. Figure 23 shows the model performances in a stochastic environment due to perturbations with Gaussian noise. The adaptive physics-based IMEP model had the least error when the respective IMEP models were perturbed with Gaussian noise—with a PPMCC of 100% and an RMSE of $1.79\times {10}^{-6}$ kPa utilizing a 75.7 s computational runtime—followed by the hybrid BRANN IMEP, which had a PPMCC of 99.69% and an RMSE of $4.47$ kPa with a 2.75 s computational runtime. The ANN IMEP model had a PPMCC of −7.33% and an RMSE of 82.63 kPa with a 0.22 s computational runtime. This model deterioration affirms the possible degrading performance and instabilities associated with ANN models that are not physics-based. It may be the case that due to the lack of physical deterministic input structures for the ANN IMEP, some input parameters that could compensate for perturbations were excluded in the pruning process. The static physics-based IMEP model showed the highest model deviation in the stochastic environment, returning a PPMCC of $88.2\%$ and an RMSE of 87.89 kPa with a 50.25 s computational runtime. Similarly, to ensure high model robustness, the hybrid IMEP model was fine-tuned with a TLM, and is discussed below.

IMEP BRANN Transfer Learning Model (TLM)

The IMEP hybrid TLM development followed a similar routine to that of the CA50 hybrid TLM. The perturbed ${\mathrm{IMEP}}_{\mathrm{Hyb}-\mathrm{BRANN}}$ was used as an input to train a BRANN TLM, as depicted in Equation (33). In this case, 80% of the data were used for the transfer learning and the remaining 20% to test the model’s predictive performance. The model utilized 15 hidden-layer neurons and 1000 epoch iterations.

$${\mathrm{IMEP}}_{\mathrm{Hyb}- \mathrm{TLM}} \left(\mathrm{k}\right)=\mathrm{f}\left({\mathrm{IMEP}}_{\mathrm{Hyb}-\mathrm{BRANN}}{}^{\left(\mathrm{k}\right)}\right)$$

(33)

The BRANN IMEP TLM returned a PPMCC of 99.8% for the testing and an RMSE of 4.37 kPa in a 3.63 s computational runtime. There was a 2.2% RMSE improvement in the performance of the hybrid IMEP model in the stochastic environment.

5. Discussion

After calibration of each of the models and various submodels, the results were compared. In this section, the results of the parametrized physics-based SOC model and the ANN SOC model are first compared, followed by a comparison of the respective BD models. The results of the CA50 models along with IMEP are then explored.

5.1. SOC Model Results

In contrast to the physics-based SOC prediction, the BRANN SOC model had a PPMCC of $99.3\%$ and an RMSE of $0.45$ CAD. A high PPMCC of $99\%$ was observed for the SOC BRANN testing set, showing the model’s high predictive capability. Figure 24 compares the performance of the parametric physics-based SOC model and the BRANN SOC model with the experimental SOC. The BRANN SOC model had better model performance than the physics-based model. The operating conditions that may have exhibited complex dynamics could be more accurately captured with the BRANN model without incurring a high computational cost.

The maximum absolute error for the BRANN SOC model occurred at an operating condition that also had a significant model error in the parametric physics-based model. This was the operating point with 20% VNP, 20% EVP, 20% EFP air-handling actuator combinations at 2500 RPM, and 150.4 Nm, which had an absolute error of $3.35$ CAD for the parametric physics-based model and an absolute error of $2.95$ CAD for the ANN model.

5.2. BD Model Results

The parametric physics-based BD model returned a PPMCC of $97.0\%$ and an RMSE of $1.77$ CAD, while the BD BRANN model had a PPMCC of $99.0\%$ and an RMSE of $1.12$ CAD. A comparison between the performances of the physics-based BD model, the BD BRANN model, and the experimental values is shown in Figure 25. The high model accuracy of $~97.0\%$ for the parametric physics-based BD model could again be explained by the model’s complexity, as depicted in Equation (18). Seventeen decision variables were utilized to capture the burn duration. The GRG algorithm utilized over $900$ iterations for a $~42.1$ s computational runtime in the parametric optimization subroutine. The BRANN BD model, on the other hand, utilized 17 epoch iterations, and had a $~0.32$ s runtime and accuracy of $~99.0\%$ PPMCC.

5.3. CA50 Model Results

The BRANN CA50 TLM development’s training subroutine was executed with 281 epoch iterations in 2.20 s, and returned a PPMCC of 100% for the testing set. When the BRANN CA50 TLM was cascaded in an unperturbed adaptive hybrid Wiebe CA50 model framework, a PPMCC of a 100% and an RMSE of $3.83\times {10}^{-8}$ CAD at a $43.2$ s computational runtime were recorded. This was a significant improvement on the results from the classical direct ANN model, which had returned an overall PPMCC of $82\%$ and a testing PPMCC of $16\%$. Similarly, when the BRANN CA50 TLM was cascaded in the hybridized framework in a Gaussian stochastic environment, an RMSE of $1.34\times {10}^{-5}$ CAD with a $45.4$ s computational runtime was observed. Compared to the direct ANN CA50 model, the physics-based CA50 model framework had a better performance in an unperturbed environment, as it returned a PPMCC of $~95.0\%$ and an RMSE of $0.17$ CAD with a $49.4$ s computational runtime. The CA50 model comparison is therefore shown for the unperturbed physics-based CA50, direct ANN CA50 models, and the hybridized CA50 TLM in order to effectively show the improvements made in the TLM modeling step in a normalized environment. Figure 26 shows the comparison of the various CA50 models’ frameworks.

The hybrid BRANN CA50 TLM clearly has superior model performance, both in stochastic and in unperturbed environments. Hence, the hybrid BRANN CA50 TLM provides a solid framework for use in model predictive combustion control. The high model fidelity would make it feasible for sensorless combustion control. Although the engine may experience drifts in inputs due to aging and environmental factors, the results indicate that a consistently robust extrapolative probability could be guaranteed.

5.4. IMEP Model Results

The IMEP results show less variation between the different options, as illustrated in Figure 27. The adaptive physics-based IMEP model was very accurate (a PPMCC of $100\%$ and an RMSE of $1.79\times {10}^{-6} \mathrm{kPa}$), while the more computationally expensive static physics-based IMEP model had more significant errors (a PPMCC of $95.2\%$ and an RMSE of $17.7 \mathrm{kPa}$). The direct classical ANN IMEP model returned a PPMCC of $99.4\%$ and an RMSE of $6.33$ kPa, while the hybrid BRANN IMEP model had a PPMCC of $99.6\%$ and an RMSE $4.94 \mathrm{kPa}$. A fair improvement of the IMEP model’s accuracy was therefore seen in the hybrid BRANN IMEP model compared to that of the ANN IMEP model, but there was an offset in computational time, with the ANN taking 0.2 s in contrast to the 0.71 s runtime of the BRANN.

The higher model accuracy of the adaptive physics-based IMEP model could be explained by the complexity of the parameterized physics-based model and the recursive nature of the operating-point-dependent algorithm that was utilized for the model optimization. A relatively high computational cost was incurred (~50.85 s) for the 212 steady-state operating conditions in the process of achieving the high accuracy of the adaptive physics-based IMEP model, while the hybrid BRANN IMEP model required a much smaller computational runtime of 0.71 s. For the adaptive physics-based IMEP model, some prior constraints on the upper and lower limits of the tuning constants were required in order to ensure the high model performance for this approach. This was determined by a trial-and-error procedure, thereby increasing the computational complexity.

As mentioned earlier in Section 4, only a 2.2% improvement in accuracy was recorded for the implementation of the hybrid IMEP BRANN TLM. Hence, the IMEP TLM model step could be bypassed by the hybrid BRANN IMEP model and still be able to achieve reasonable model accuracy in a stochastic environment. However, the IMEP TLM may be useful for control development when there are hard constraints, such as the availability of only a low-bandwidth microcontroller, which may require a step down in computational runtime.

6. Conclusions and Future Work

We have demonstrated the use of an integrated approach using combustion physics and ANNs for combustion phasing modeling of a modern diesel engine for use in real-time combustion and engine load control. Steady-state parametrized physics-based models for the SOC and BD, along with IMEP, were developed along with simple Bayesian regularized artificial neural network models for these parameters. In addition, an adaptive nonlinear parameter optimization algorithm was used for the prediction of CA50 using the developed SOC and BD models as inputs to a 0D simple Wiebe combustion model, whose outputs are fine-tuned with an ad hoc BRANN model in the event of stochastic environmental changes. The hybrid model frameworks for CA50 and IMEP showed superior performance, and will be leveraged in future work for uncertainty and disturbance estimator (UDE)-based combustion control. Other advanced methods will also be explored in future work, such as the use of particle swarm optimization for in absentia training of the BRANN transfer model, as well as the use of parsimonious neural networks (PNNs) to learn interpretable physical laws whilst directly capturing CA50 and IMEP. The main highlights of this study are as follows:

• Due to the recursive nature of the proposed hybrid frameworks, in addition to the augmentation provided by the CA50 and IMEP TLMs, the fidelity of the submodels may not matter much. This could reduce modelling effort and avoid the need to track all uncertainties pertaining to the submodels. As such, lower quality data may be able to be utilized for the SOC and BD models. Drift in the various measurements due to engine aging and other environmental factors for the IMEP model may also not be much of a concern;
• TLMs for the hybrid IMEP developed in the same fashion as in this study may not be necessary. However, it may be useful for control in the long term, where physical constraints of the control processors may be inevitable;
• The robustness of the models of less complex parameters, such as IMEP, may be feasible in a physics-only model framework, but there would be the need to develop complex physics-based models and, thus, incur some computational costs and increased computational complexity.

The hybridized models augmented with TLMs that were developed in this study showed superior steady-state performance as compared to the physics-based-only and ANN-only frameworks.