Machine Learning Approaches to Learn HyChem Models

Arrhenius.jl is built using the programming language of Julia to leverage the rich ecosystems of auto-differentiation and differential equation solvers. Arrhenius.jl does two types of differentiable programming: (i) it can differentiate elemental computational blocks. For example, it can differentiate the reaction source term with respect to kinetic and thermodynamic parameters as well as species concentrations; (ii) it can differentiate the entire simulator in various ways, such as solving the continuous sensitivity equations [17] as did in CHEMKIN [18] and Cantera [19] and in adjoint methods [20, 21]. The first type of differentiation is usually the basis of the second type of high-level differentiation. Arrhenius.jl offers the core functionality of combustion simulations in native Julia programming, such that users can readily build the applications on top of Arrhenius.jl and exploit various approaches to do high-level differentiations.

Figure 1 shows a schematic diagram of the structure of the Arrhenius.jl package. The Arrhenius.jl reads in the chemical mechanism files in YAML format maintained by the Cantera community, and the chemical mechanism files contain the kinetic models, thermodynamic, and transport databases. The core functionality of Arrhenius.jl is to compute the reaction source terms and mixture properties, such as heat capacities, enthalpies, entropies, Gibbs free energies, etc. In addition, Arrhenius.jl offers flexible interfaces for users to define neural network models as submodels and augment them with existing physical models. For example, one can use a neural network submodel to represent the unknown reaction pathways and exploit various scientific machine learning methods to train the neural network models, such as neural ordinary differential equations [9, 22] and physics-informed neural network models [23, 24]. One can then implement the governing equations for different applications with these core functionalities and solve the governing equations using classical numerical methods or neural-network-based solvers, such as physics-informed neural networks [23, 24]. Arrhenius.jl provides solvers for canonical combustion problems, such as simulating the auto-ignition in constant volume/pressure reactors and oxidation in jet stirred reactors.

Compared to the legacy combustion simulation packages, Arrhenius.jl can not only provide predictions given the physical models, but also optimize model parameters given experimental measurements. By efficiently and accurately evaluating the gradient of the solution outputs to the model parameters, experimental data can be incorporated into the simulation pipeline to enable data-driven modeling with deep learning algorithms.

Following the work of [7], we employ the IDT dataset compiled in [7], including the IDTs of the lean, stoichiometric, and rich F-24/air mixtures measured at 20 bars in a rapid compression machine (RCM) and a shock tube (ST) over the temperature range of  625 K to  1250 K. The RCM is driven pneumatically with a hydraulic stop apparatus, and the initial chamber temperature and pressure are achieved by adjusting the compression ratio [25]. Low-temperature (up to 700 K) IDTs were measured in the RCM. The shock tube uses a polycarbonate diaphragm to separate the driver and driven sections, the lengths of which are 2.744 m and 4.227 m, respectively. A tailored mixture of He and N₂ was used as the driver gas, and premixed F-24/air mixtures were introduced to the driven section.

For the base model, we adopt a HyChem version of the Jet-A skeletal mechanism [5] that covers the NTC regime. It has 48 species and 254 reactions. In addition, details of the formulations and stoichiometric coefficients in the F-24 model are the same as in [7], and they are listed in Table 1. For simplicity and consistency, all of the simulations are carried out under constant volume conditions using Arrhenius.jl.

As the only available dataset is the IDT, the key step in utilizing SGD for the optimization of the HyChem model is to compute the gradient of the IDT with respect to the kinetic parameters. Recent work [26, 17, 27, 28] have substantially advance the algorithms for computing the gradient of IDT to model parameters, this work employs the sensBVP method [26]. The approach converts the initial value problem (IVP) to a boundary value problem (BVP) by treating the temperature at ignition as a boundary condition and the IDT as a free variable to solve. The algorithm is implemented with Arrhenius.jl and further enhanced with the auto-differentiation, multi-threading, and sparse matrix techniques. For the size of the skeletal HyChem model considered in this work, the computational cost of the sensitivity of IDT is negligible comparing to computing the IDT.

It is worthy to note that the gradient of laminar flame speed to kinetic parameters is also usually computed using the boundary value problem approach as in the sensBVP approach for IDTs. Thus, the computational cost for the gradient computation of flame speed is also negligible comparing to computing the flame speed. Therefore, the gradient-based optimization approach proposed in this work can be readily extended to the experimental dataset of flame speed as well.

Following the HyChem approach, the 16 steps of the pyrolysis submodel, listed in Table 1, are subject to optimization, and the rest of the core model is fixed. Although in the work of [7], only the pre-factor $A$ and activation energy $Ea$ are optimized, the SGD optimizer can easily scale up to high-dimension, such that all three parameters, including the non-exponential temperature dependence factor $b$, are optimized. Consequently, the optimization involves 48 parameters.

As discussed in the introduction, we would like to regularize the model to achieve an optimized model that can fit the experimental datasets well but only moderately change the kinetic parameters from the base model. To facilitate regularization, the model parameters are scaled and centered at zero, such that the relative changes of Arrhenius parameters is to be optimized. Specifically, we define the scaled parameters as

where the subscript ₀ refers to the base model. The dimension of $Ea$ is specified as $kcal/mol$ rather than $cal/mol$ as a change of 2 $kcal/mol$ in $Ea$ is roughly equivalent to changes of $e$ times in $A$, such that the changes in $A$ and $Ea$ will be balanced and avoid stiffness in the parameter space.

The loss function is defined as the mean square error (MSE) between the predicted IDTs in the logarithmic scale using the optimized model and the measured IDTs, as shown in Eq. 2:

The choice of logarithmic scale is similar to the choice of relative error between measured and predicted IDT as in [7]. The Adam [29] optimizer with the default learning rate of 0.001 is adopted. Weight decaying and early stopping [30] are employed to regularize the parameters, such that the optimization will prefer kinetic parameters that are close to their original values.

Figure 2 show the comparisons between the experimentally measured IDTs and the predictions using the base HyChem model, the optimized model using genetic algorithms (GA), and the current model. As pointed out by [7], the base model substantially underpredicts the IDTs in the NTC regimes for all three equivalence ratios, while it performs reasonably well in the high-temperature regimes. Conversely, both of the optimized models substantially improve the performance in the NTC regime.

Since the slope of the response curve is indicative of the global activation energy, it is interesting to note that the current model shows a smaller global activation energy at lower temperature regimes, where the temperature is below 700 K, than the one suggested by the experimental data and the GA model. Such global activation energy is related to the key reactions that govern the low-temperature auto-ignition [31, 32], and thus it is important to accurately predict the global activation energy to gain physical insight of the ignition behavior of F-24. The discrepancies in the predicted global activation energy can be attributed to several plausible reasons:

1. 1.

   The current model indeed underpredicts the activation energy due to large uncertainties in the training data, since the IDT measurements in the NTC regime are usually associated with large uncertainties. In addition, the simulations for training and predictions are carried out under constant volume assumptions for simplicity and consistency, which might not be appropriate for two-stage ignition [33]. These uncertainties in the training dataset and modeling may result in a biased model. In such case, the data re-balance procedure in [7] can be employed to increase the weights of low-temperature data on the model training.

2. 2.

   The reported measurements are subject to heat loss, and therefore overestimate the global activation energy. The heat loss for experiments with long testing time in RCM usually has a substantial effect on the IDTs. The reported experimental data using the temperature at the end of compression gives higher global activation energy than the actual one. In such case, the effective temperature [34, 35] is usually a better representative temperature than the temperature at the end of compression.

On the other hand, the SGD optimization tends not to over-fit the data at low temperatures compared to the GA optimization, since multiple regularization techniques are employed to prevent over-fitting. From this point of view, the experimental data are likely giving an overly large global activation energy than the actual one. Nevertheless, good data-driven modeling not only can help us develop a predictive model from experimental data but also help reveal potential experimental uncertainties and guide further experimental designs to advance our physical understandings.

To assess the generalization capability of the learn model, we further evaluate the performance of the optimized models in predicting laminar flame speeds, although the model is not optimized against flame speed data. To the best of our knowledge, there are no experimental data of flame speed for F-24 in the open literature. Thus, we apply the model to predict flame speed for Jet-A as F-24 is derived from Jet-A and might share similar flame propagation behaviors. The predicted and measured flame speeds for Jet-A are presented in Fig. 3.

While it is difficult to assess which model performs better given the uncertainties in the experimental data, we can use the base model as a reference, for the base model is calibrated for high-temperature pyrolysis. The base model is expected to perform well for flame propagation that is dominated by high-temperature chemistry. Overall, both SGD and GA models predict the similar flame speed with the base model, with the differences within 2 cm/s. On the other hand, the GA model systematically predicts slower laminar flame speeds than the base model. The predictions by the current model agree with those by the base model at fuel-lean conditions and show similar trends as those by the GA model under fuel-rich conditions. This could be potentially attributed to the regularization in the SGD optimization. As shown in Tables 2 and 3, the kinetic parameters in the current model are closer to the base model compared to those in the GA model. For the GA model, the activation energy for the direct decomposition reaction of R1 is changed a lot, and such changes may lead to systematic changes in flame speed under both fuel-lean and rich conditions. For the SGD model, it is possible that the flame speed is only sensitive to the changes of kinetic parameters under fuel-rich conditions, as the fuel-lean and rich conditions correspond to the different radical pools. Future work could exploit advanced pathway analysis tools, such as global pathway analysis [36], to further investigate the different trends under lean and rich conditions. Nevertheless, the SGD model tends to less affect the model parameters for high-temperature chemistry while improving the low-temperature chemistry from experimental data. Thus, comparing to the GA model, the SGD model is more likely to be able to extrapolate to conditions beyond the range of the training dataset.

While the experimental dataset of laminar flame speed can be readily incorporated in the SGD approach as discussed in the methodology section, the results in Fig. 3 suggest that learning the HyChem model from the ignition delay time dataset might be sufficient to develop a model that works well for both ignition and flame propagation with the help of regularization. Such generalization capability is especially useful recognizing that experimental measurements of laminar flame speed are usually limited to near-atmospheric pressures due to the onset of turbulence at elevated pressures [37, 38]. Nevertheless, the SGD approach can potentially enable us to develop a data-driven kinetic model by providing only the high-pressure ignition delay time dataset but can work well for predicting both ignition and flame propagation at engine relevant elevated pressures.

We then present the details of the learning algorithms that produce the SGD-optimized model. Figure 4 shows the loss history in the training and validation of the model. 76 out of 96 samples were randomly chosen as the training dataset, and the rest were for the validation. Within each epoch, the training dataset was randomly shuffled to introduce randomness to the training process to prevent over-fitting. Mini-batch optimization was adopted such that the model parameters were updated for each sample and 76 times for each epoch. Cross-validation and early-stopping were exploited to prevent over-fitting. For instance, the training is stopped when the loss in the validation reaches a plateau, although the loss in the training is still decreasing. The validation loss reached a plateau at around 60 epochs, although we let the training continues to 110 epochs to confirm the plateau. We thus select the model at 58 epochs, which is the closest checkpoint near 60 epochs, and the checkpoint is saved every time the validation loss reaches a new minimal. In summary, we adopted multiple techniques, including mini-batching, cross-validation and early-stopping, to regularize the model and prevent over-fitting. The regularization in SGD in general leads to better generalization performance, which is one of the most important advantages comparing to GA optimization, although both SGD and GA can give a similar fitness in a small dataset with less than one hundred parameters.

As far as the computational cost is concerned, SGD optimization is also much more efficient than genetic algorithms. To estimate the computational costs for both algorithms, we denote the computation of the IDT for all 96 samples as one computational unit. Since the computational cost for gradient evaluation is negligible, the current SGD optimization for F-24 costs about 100 units. For comparison, the GA optimization for the same datasets and base model costs about $10^{5}$ units, with 32 particles and 300 generations [7]. Such direct comparison reveals 1000 times of acceleration, which also results in a reduced need for computational power. For example, the current training was performed on a desktop workstation for hours, while the GA optimization in [7] was performed on a high-performance cluster. In addition, the acceleration factor scales with the size of the dataset and the dimension of the model parameter. Therefore, such computational saving is expected to be more significant when the kinetic system is more complex and more data are available.