Model-Based Catalyst Screening and Optimal Experimental Design for the Oxidative Coupling of Methane

1. Introduction

Synthesis of platform chemicals via catalyzed reactions lead to a wide product distribution because every catalyst has different active sites, composition and response to operating conditions, making it complex to identify their role in reaction mechanisms. This challenges the selective and economical manufacturing of target products at scale, as evidenced by studies on the oxidative coupling of methane (OCM) where the selective formation of ethane/ethylene (C2) for the polymer manufacturing chain is limited by the thermodynamically favored over-oxidation [1]. The Edisonian approach to material design has been replaced by information-driven platforms that seamlessly integrate digitized database with modeling and optimization for hypothesis-driven design decisions [2]. Central to these platforms are machine learning (ML) surrogates that map properties of candidate materials in the database to process performance outcomes, so that appropriate materials can be recommended for experimental synthesis [3]. This work seeks to use the high-throughput experimental database for a variety of mixed metal oxide catalysts ($M1-M2-M3O_4/\mathrm{Support}$) to train ML surrogates for catalyst screening and to device future experimental strategies that meet the selectivity-conversion targets for the extensively studied OCM chemistry.

The digitized data to build catalyst informatics platforms curated exhaustively from literature [4] are associated with inconsistencies in data, methods and reproducibility because of which there has been a shift towards high-throughput experimentation (HTE) [5], and high-throughput theoretical calculations (HTCs) [6] to reliably record catalyst performance across scales from the level of reaction energetics to process operations. Web-based visualization tools to deploy exploratory data analysis on HTE data using co-ordinated multiple views (CMVs) to discover apparent trends in the reaction performance across a variety of catalysts and operating conditions can provide insights for future experimentation [7]. Sophisticated ML tools to uncover the not so apparent insights require quantitative descriptors of a catalyst from elemental properties (atomic numbers, electron affinity, ionization energy, density) of constituent metal atoms from the periodic table to characterize its activity [8], or HTC-based reaction energetics descriptors from density functional theory [9,10]. Once the catalyst design space has been quantified by descriptors, unsupervised clustering can be used to identify catalyst groupings based on how they impact reaction performance, across different experimental conditions [11]. ML has been used to develop supervised descriptor-based reaction performance prediction models, and to minimize the time and cost in strategizing recommendations for physical experiments or theoretical calculations to guide exploration of the design space for materials discovery [12]. Descriptor-based ML models have been used to screen electrocatalysts for carbon-dioxide reduction [13], and also for the adaptive electrocatalyst and photocatalyst discovery either by human-in-the-loop learning, where the ML model is updated once the outcome has been observed via experimental/theoretical runs at algorithmically sampled points of the design space [14,15,16]. Alternatively, descriptor-based ML models have also been used for goal-driven exploration via Bayesian optimization or evolutionary genetic algorithms [17], to create self-driving laboratories that integrate databases (literature, HTE, HTC), ML and automated experimentation [18,19].

Most of the aforementioned approaches are yet to reveal catalyst candidates for OCM chemistry with a C2 yield> 30%, a threshold considered practical for industrial applications that are limited by the maximum achievable C2 yields because the reactant methane, is much less reactive to oxygen than the target C2 products leading to selectivity-conversion tradeoffs. Analysis of 1868 literature reported OCM catalysts, reveals that most of them barely meet 20% C2 yields, with just ∼ 12 of them surpassing the thresholds for feasible industrial production [4]. The inconsistencies of literature-reported data (missing data, mass balance errors), not only pose an obstacle to reproducibility but are shown to result in poorly trained regression models to predict reaction performance that register prediction outliers on literature data with C2 yields greater than 30%. For instance, the support vector regression trained on literature-mined data for OCM chemistry to predict C2 yields has $R^2\sim 05-0.6$, which is not impressive, because of which catalyst candidates discovered by it when used as a surrogate in Bayesian optimization lacks diversity in predicted materials, with a narrow field around $L{a}_2{O}_3$ derivatives, and a maximum C2 yield of ∼ 15-16% [20]. To ensure reliability of the database used to propose catalyst candidates for OCM chemistry, HTE data has been used with informatics tools for visualization, supervised ML and catalyst networks to uncover patterns among dynamically evolving factors like catalyst synthesis, composition and operating conditions on reaction performance [21]. However, going beyond the interpolation filling abilities of ML in multi-dimensional data to predict rare targets with C2 yields > 30% when the HTE dataset used to train the ML surrogates covers yields capped as much lesser values ∼ 20 %, is still an ongoing research effort. In that spirit, this manuscript represents an effort to create informed serendipity using ML surrogates to enhance discovery of catalyst candidates by avoiding a narrow field focus.

Most works outlined herein develop descriptor-based ML to predict C2 yields and $C{H}_4$ conversion using random forest regression, or neural network formalisms with mass balance reconciliation for the same [22]. However, using these ML models to develop catalyst performance curves by tuning operating conditions that maximize both C2 selectivity ($S_{C_2}$) and methane conversion ($X_{C{H}_4}$) for each of the catalysts, followed by using these performance curves to screen $M1-{\left(M2\right)}_{1-2}-M3O_x/\mathrm{Support}$ type catalysts with respect to the popularly used $Mn-N{a}_2-W{O}_4/Si{O}_2$ for OCM chemistry, is yet to be investigated. Most ML models have been rationalized in terms of feature importances of the descriptors in predicting reaction outcomes, however, ML model validation to deduce activation barriers via lumped reaction kinetics models for methane conversion, and subsequent overoxidation of C2, is also pending investigation. Additionally, the desciptor-based ML models have also been used to test an evolutionary framework of exploring combinations of catalyst descriptors and process conditions in tandem, that maximize the two-fold selectivity-conversion targets, in an attempt to propose new candidates for synthesis. A multi-objective optimization routine using the NSGA-2 genetic algorithm has been contrasted against a Bayesian optimization routine to propose candidates and operating conditions, in an attempt to analyse field focus in proposing candidates, and the number of generations/ sampling iterations that are required to arrive at proposed entities with the highest achievable C2 yields. Assessing the synthesis feasibility of the proposed candidates and their experimental validation is out of scope of this manuscript.

2. Methodology

2.1. Dataset

The HTE database for OCM chemistry using 40 types of $M1-{\left(M2\right)}_{1-2}-M3O_x/\mathrm{Support}$ catalyst, and 19 references across 216 experimental conditions leading to ∼ 12700 data points [5], has been used in this work. The dataset is hosted on a web-based informatics platform called the Catalyst Acquistion by Data Science (CADS), and records reaction outcomes at the end of each run in sequentially programmed experimental campaigns across combinations of temperature (900, 850, 800, 775, 750 and 700 °C), total reagent flowrate (10, 15, 20 ml/min), $C{H}_4/O_2$ ratio (2,3,4 and 6 mol/mol), and contact times (0.75, 0.50, or 0.38 s) [23]. The conversion of methane and yields/selectivity of the target products, ethane and etylene (C2 products), as well as the undesired products (carbon dioxide, carbon monoxide) as a result of over-oxidation have been recorded. The design space is defined by quantitative descriptors of the catalyst and the operating conditions. Catalyst descriptors encompass atomic numbers of its constituent metal atoms (M1, M2, M3), the composition of these metal elements in mol%, and finally the nature of the support indicated by its unique index identification, all of which have been indicated in the CADS data repository. Descriptors of the reaction conditions encompass temperature (T), contact time (t), total flow rate (${\dot{V}}^{\circ }$), methane flow rate (${\dot{V}}_{CH4}^{\circ }$) and $C{H}_4:O_2$ , and have also been furnished in the dataset. Hence, there are 12 descriptors in all, when it comes to defining the design space that has been investigated via HTE to find optimal combinations that maximize methane conversion and C2 selectivity for OCM chemistry.

2.2. Random Forest Regression

Data collection under controlled experimental conditions eliminates most inconsistencies due to variability across experimental platforms. Yet, one cannot avoid uncertainties in the recorded measurements either because of material balance violation or instrumental errors [24]. The use of an ensemble model that learns aggregate predictions in this work, runs a lower risk of overfitting such biases [25], placing lesser emphasis on the need to eliminate data points with higher mass balance uncertainties. Also, random forest regressors learn via decision thresholds on descriptors to segment the design space in which aggregate predictions are made, making it agnostic to scale of the features and eliminating the need for much data pre-processing. When highly parametrized set of ML models like neural networks that run the risk of overfitting to the artefacts in data is used, good pre-processing, efforts to embed the training with mass balances, and even truncation of datapoints with mass balance violation beyond a fixed threshold have been widely considered [22,26].

Random forests are an ensemble model comprising many decision trees. Each decision tree has binary nodes, where the number of samples n at each node splits into $n_1$ and $n_2$ samples based on the optimal segmentation of the $j^{th}$ descriptor into 2 subnodes, basd on a decision threshold $\theta$, as outlined in Eqn 1

$$\begin{array}{cc}\hfill R_1(j,\theta )& =\left\{x_{ij}\right|x_{ij}\le \theta \}\forall i=1,2,\cdots n_1\hfill \\ \hfill R_2(j,\theta )& =\left\{x_{ij}\right|x_{ij}>\theta \}\forall i=1,2,\cdots n_2\hfill \end{array}$$

(1)

The descriptor j and its value is chosen such that the residual sum of squares is minimized for the binary split at that tree node, achieved by minimizing the residual sum of sqaures as in Eqn 2, where ${\overline{y}}_1$ and ${\overline{y}}_2$ are the average target response of the samples in each of the subnodes $R_1(j,\theta )$ and $R_2(j,\theta )$.

$$\underset{j,\theta }{min}\left[\underset{{\overline{y}}_1}{min}\sum _{x_{ij}\in R_1(j,\theta )}{(y_i-{\overline{y}}_1)}^2+\underset{{\overline{y}}_2}{min}\sum _{x_{ij}\in R_2(j,\theta )}{(y_i-{\overline{y}}_2)}^2\right]$$

(2)

Repeating the process, stratifies the d dimensional design space into W regions $R_1,R_2,\cdots R_W$ to generate a decision tree given in Eqn 3, where $\mathbb{I}(.)$ is the indicator function, i.e. $\mathbb{I}=1$ if $x\in R_w$ else it is 0, w is the number of decision threshold splits, and ${\overline{y}}_w$ is the average target response of the samples in the region $R_w$.

$$f\left(x\right)=\sum _{w=1}^W{\overline{y}}_w\mathbb{I}(x\in R_w)$$

(3)

This procedure is repeated on an ensemble of decision trees, $f_t\left(x\right)$, where $t=1,2,3\cdots Ntrees$ to generate a random forest model that aggregates the predictions across the learners in order to map the descriptors to the target response, as $\overline{f}:X\to Y$ given in Eqn 4

$$\overline{f}\left(x\right)=\frac{1}{N_{trees}}\sum _{t=1}^{N_{trees}}{f}_t\left(x\right)$$

(4)

2.3. Model Validation by Power-Law Reaction Kinetics

This work demonstrates the ability of ML models to connect descriptors to parameters of lumped kinetic models to facilitate interpretability when used as a catalyst informatics tool. Evidence of ML being rationalized either by incorporating domain knowledge as constraints during its training or postfacto by the use of metrics like SHAP, feature importance and locally interpretable model explanations are quite popular. Here, in addition to just feature importance, we seek to interpret the relationship between descriptors and kinetic parameters using a lumped kinetic model for the otherwise complex reaction mechanism for OCM chemistry that involves both the gas and surface reactions [27]. The reactions considered in this work are as follows:

$$\begin{array}{cc}\hfill 2C{H}_4+0.5O_2& \to C_2{H}_6+H_{2}O\hfill \\ \hfill C_2{H}_6+0.5O_2& \to C_2{H}_4+H_{2}O\hfill \\ \hfill C{H}_4+1.5O_2& \to CO+2H_{2}O\hfill \\ \hfill C{H}_4+2O_2& \to C{O}_2+2H_{2}O\hfill \end{array}$$

Once the random forest regression has been trained using the descriptors, the stoichiometry of the global reaction scheme by combining the above equations can be used to regress the lumped power-law kinetic model given in Eqn 5.

$$\begin{array}{cc}\hfill r_{CH4}& =k_{10}exp\left[\frac{-E{a}_{CH4}}{RT}\right]{\overline{P}}_{CH4}^a{\overline{P}}_{O2}^b\hfill \\ \hfill r_{C2}& =k_{20}exp\left[\frac{-E{a}_{C2}}{RT}\right]{\overline{P}}_{CH4}^{a*}{\overline{P}}_{O2}^{b*}\hfill \end{array}$$

(5)

Latin hypercube sampling (LHS) is used to randomly sample the process operation descriptors comprising temperature (T), contact time (t), total inlet volumetric flow rate (${\dot{V}}^{\circ }$), inlet volumetric flowrate of methane (${\dot{V}}_{CH4}^{\circ }$) and $C{H}_4:O_2$ molar ratio, for a given set of catalyst descriptors comprising details of the elemental metals, their molar percentages and support ID. The reaction is performed in a continuous flow reactor, and is operated at a pressure P of 1 bar under isothermal conditions. Using this information the partial pressures of the reactants at the start of the reaction is calculated. The random forest regressor is then used to predict the conversion and selectivity for each set of descriptors across all samples, to obtain ${\widehat{X}}_{CH4}\%$ and ${\widehat{S}}_{C2}\%$, using which the reaction rates are expressed as given in Eqn 6, where ${\dot{n}}_{CH4}^{\circ }$ is the molar flowrate of methane entering the reactor. The HTE data has been reported for 1g mass of catalyst support impregnated with the mixed metal oxides [5], and hence the kinetic models implicitly fit the specific reaction rates.

$$\begin{array}{cc}\hfill r_{CH4}& ={\dot{n}}_{CH4}^{\circ }\frac{{\widehat{X}}_{CH4}}{100}\hfill \\ \hfill r_{C2}& ={\dot{n}}_{CH4}^{\circ }\frac{{\widehat{X}}_{CH4}{\widehat{S}}_{C2}}{100}\hfill \end{array}$$

(6)

The stoichiometry of the global reaction scheme and the predicted methane conversion from the random forest regression is then used to calculate the partial pressure of the species at the end of the reaction. The average of the initial and the final partial pressures for oxygen and methane, along with the reaction rate expressions in Eqn 6 is substituted into the power-law kinetic expression in Eqn 5 for all the LHS sample points to estimate via regression the kinetic parameters viz. the Arrhenius pre-exponential factors for methane conversion ($k_{10}$) and C2 formation ($k_{20}$), their corresponding apparent activation energies ($E{a}_{CH4},E{a}_{C2}$), and the orders of the species in each of the reactions (a,b,a*, b*). It must be noted that the HTE datasets report reaction performances only at the end of the contact time, owing to which fitting kinetic models to species concentration profiles is approximated by the average of the initial and final partial pressures. The estimated kinetic parameters are therefore treated as coarse estimates to characterize the impact of different catalysts on the reaction performance across varying operating conditions sampled by LHS.

2.4. Genetic Algorithm for Multi-Objective Optimization

The reaction performance for OCM chemistry in terms of methane conversion and C2 selectivity are known to exhibit a tradeoff, and depend on the operating conditions and the type of catalyst used. The combination of descriptor values $x\in R^{d\times 1}$, that maximizes both reaction performance indicators is posed as a multiobjective optimization problem (Eqn 7), where $f_{X_{CH4}}(.)$ and $f_{S_{C2}}(.)$ are the trained random forest regressors to predict the corresponding indicators.

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\displaystyle \underset{x}{max}}{f}_{X_{CH4}}\left(x\right)\hfill \\ {\displaystyle \underset{x}{max}}{f}_{S_{C2}}\left(x\right)\hfill \end{array}\right.\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathrm{S}.\mathrm{T}.& lb\le x\le ub\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & M1\%+M2\%+M3\%=100\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\dot{V}}^{\circ }-{\dot{V}}_{CH4}^{\circ }-{\dot{V}}_{CH4}^{\circ }\frac{1}{C{H}_4:O_2}>0\hfill \end{array}$$

(10)

The multiobjective optimization is constrained by limits of the descriptors (Eqn 8), given by the their range bounds in the HTE dataset, and Eqn 10 to ensure that the inlet volumetric flow of the inert Ar gas is non-negative. When solving the above optimization to develop performance curves for catalyst formulations outlined in the HTE datasets, the elemental metals (M1, M2, M3), their corresponding molar percentages and their supports are fixed descriptor values, with the decision variables comprising just the operating conditions. However, when it comes to proposing different catalysts, their descriptors along with the operating conditions are treated as variables in the decision space in solving the multiobjective optimization, wherein the additional constraint in Eqn 9 enforces closure in the molar percentages of the elemental metals of the catalyst. The elemental metals and their supports are treated as categorical descriptor values but the rest of the descriptors are continuous. The multiobjective optimization is solved using an evolutionary appraoch via the NSGA-2 genetic algorithm [28], using 50 individuals, over 200 iterations with a mutation probability of 0.2 and crossover probability of 0.8. The constraints are implemented via the Delta penalty approach [29], where the fitness is penalized for invalid individuals by constant factor delta that is subtracted from the objectives we seek to maximize.

2.5. Bayesian Optimization for Adaptive Experimentation

Bayesian optimization is a sequential global optimization approach that iteratively samples the design space of decision variables using a probabilistic surrogate model [30], like a Gaussian process regressor that captures the distribution of target predictions, P(y|x) for a given $x\in R^{d\times 1}$, and an acquisition function like expected improvement (EI) to guide sampling as given in Eqns 11-13, subject to the constraints in Eqns 8-10.

$$\begin{array}{cc}\hfill f\left(x\right)& =\frac{f_{X_{CH4}}\left(x\right)f_{S_{C2}}\left(x\right)}{100}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill x^{*}& \leftarrow \underset{x\in R^{d\times 1}}{\arg \max }f\left(x\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill EI\left(x\right)& =\sum _y(y-f\left(x^{*}\right))P\left(y\right|x)\hfill \end{array}$$

(13)

The idea is to start with an initial number of LHS samples, say 10, from the space of decision variables ($x\in R^{d\times 1}$), and use the above method to sample by exploitation to find the most likely optimal solutions based on the posterior distribution, while also resorting to exploration by sampling from points in areas with low probability density in order to be able to find the combination of descriptor values that maximize the yield of the desired C2 products. This approach has been widely used to encourage serendipity while navigating the combinatorial explosion of the decision space of design descriptors for the goal-driven enumeration of candidates in material science [31].

3. Results and Discussion

Section 3.1 presents the results from estimating kinetic parameters for different catalysts in the HTE dataset for OCM chemistry using power-law models based on the predictions of the descriptor-based random forest regressors. Section 3.2 discusses the performance curves of the given catalyst by tuning the decision space of operating conditions to maximize reaction performance indicators. The catalysts are screened with respect to a reference with the aid of the performance curves. Section 3.3 is an attempt to use two different techniques viz. multiobjective and Bayesian optimization to navigate both the catalyst and operating condition descriptors to propose new candidates.

3.1. Assessment and Validation of Random Forest Regression via Kinetic Parameter Estimation

A random forest (RF) regressor model is fitted to map the descriptors to the reaction performances by way of methane conversion and C2 selectivity using a 5-fold cross validation for model hyperparameter tuning to prevent it from overfitting to the training data. An 85% train-test split is used for the same and predictions are assessed on completely new test data. The parity plots shown in Figure A1 indicates that the RF model adequately captures trends in the training data, and generalizes well on the test data too. Error from instruments or intrinsic phenomena like sintering [24] can potentially lead to mass balance inconsistencies in HTE datasets quantified in terms of the total carbon balance error based on the difference between the methane conversion and the yields of the products formed, as follows:

$$\mathrm{Total}\mathrm{carbon}\mathrm{balance}\%=\frac{X_{CH4}-Y_{C2H6}-Y_{C2H4}-Y_{CO}-Y_{CO2}}{X_{CH4}}\times 100$$

The ensemble nature of the RF regressors make them robust to fitting such biases as can be seen in Figure A2, where the total carbon balance is uncorrelated with the prediction errors of the RF models for both conversion and selectivity. RF regressors also have the inherent procedure of calculating feature/descriptor importance [32], based on the optimal choice of descriptor that most reduces the residual sum of squares (Eqn 2).

It can be seen from Figure 1 that the operational descriptors like temperature, inlet volumetric flowrate of methane and the $C{H}_4:O_2$ ratio are more important than any of the catalyst descriptors to the prediction model for methane conversion. While for the C2 selectivity model, the catalyst descriptors encompassing atomic numbers of metal elements M1 and M2, and the support ID are important in addition to the operational descriptors of temperature and inlet volumetric flowrate of methane and $C{H}_4:O_2$ ratio. Although the catalyst surface active sites at OCM conditions are still unknown, the rate determining step involves hydrogen abstraction by C-H bond cleavage of methane either via surface-active oxygen (Langmuir-Hinshelwood kinetics) or via oxygen from the lattice sites (Mars-Van Krevelen), and requires high temperatures [33]. Coupling of methyl radicals to form ethane that dehydrogenates to ethylene is thermodynamically less favored than its further oxidation to $C{O}_x$ because of which, although high flowrates of reagents are known to increase methane conversion, lower proportions of oxygen is used [34]. Also, the use of M2 type promoters in mixed metal oxide catalysts of $M_1{M}_3{O}_4$ type are found to suppress further oxidation by hindering the exposure of the tetrahedral $M_3{O}_4^{2-}$ active site, thereby increasing C2 selectivity [35]. Hence, prediction of methane conversion is dominated by the identified operating condition descriptors, and that of C2 selectivity is governed by the said catalyst descriptors as seen in Figure 1.

Once the RF models have been duly fit and assessed, it is important to validate them via fitting power-law kinetic parameters for each of the catalysts in the HTE dataset. As specified in Section 2.3, the points sampled via LHS in the space of the descriptors for the operating conditions are outlined in Table A1, and are used to regress the rate expressions to estimate the kinetic parameters for each catalyst. The $MnN{a}_{2}W{O}_4/Si{O}_2$ is a popular OCM catalyst that has registered high experimental C2 yields ∼ 14-27% and stability [36] and hence has been chosen as reference with respect to which the other catalysts are screened in this work. Figure 2b,c presents the regression to estimate the kinetic parameters for methane conversion and C2 formation, using the RF model predicted conversion and selectivity values at the LHS points as shown in Figure 2a.

The same procedure is followed for the other catalysts, where Figure A3 reports the regression fits, Figure A4 reports the species orders, and Figure 3c reports the apparent activation energies. The apparent activation energies for $MnN{a}_{2}W{O}_4/Si{O}_2$ in Figure 3a is within the literature ballpark ∼ 200-270 kJ/mol [27]. Clear groupings among the catalysts when the supports are varied for the $MnN{a}_{2}W{O}_4$ catalyst in Figure 3a, and the metal atoms are varied for a fixed $Si{O}_2$ support in Figure 3b have been observed. Transition metal oxide supports are reported to have better C2 yields than $Si{O}_2$ at similar process conditions [6], and even when porous aluminosilicates/zeolite-like supports or SiC are used as it leads to the formation of highly dispersed active sites after calcination [11], thereby having lower apparent activation energies for OCM reactions, as seen in Figure 3a. It must be noted that the regression calculated activation energies for zeolite-like supported $MnN{a}_{2}W{O}_4$ catalysts is very low, and is not the case in reality. This is primarily an artefact of approximately estimating kinetic parameters using RF predicted reaction KPIs at the end of contact time, and at the given LHS sample points by neglecting concentration gradients and potentially different chemistries. However, it is a reasonable approach to screen for trends and groupings among different catalysts. Similarly, in Figure 3b shows 3 groupings among alkali/alkaline earth metals, transition metals, and lanthanides/actinide metals on a fixed $Si{O}_2$ support with increasing apparent activation energies is seen. This can largely be attributed to the electronic properties of the mixed metal oxides by way of their electronegativity and ionization energies that impact the activation of methane and gas phase oxygen [37].

3.2. Performance Curves for Catalyst Screening

A good approach to compare how different catalysts i.e. ones with different supports or different mixed metal oxides impact the reaction KPIs of OCM chemistry is to use kinetic models to ascertain the best performance which can be achieved by a given catalyst over a range of operating conditions [38]. Using RF models as kinetic surrogates to maximize both methane conversion and C2 selectivity for a given catalyst with fixed loading of metal atoms as outlined in the HTE dataset, in the decision space of just the 5 operating condition descriptors has resulted in the S-X performance curves shown in Figure 4.

The operating conditions corresponding to the best C2 yields obtained from these S-X performance curves have been tabulated in Table A2 and Table A3, and shows on average across all catalysts ∼ 15% improvement over the experimental values of the HTE dataset within the limits of the total carbon balance (TCB%). The C2 yield is defined as the product of the methane conversion and C2 selectivity that the RF models have been trained to predict, and combines the tradeoff between 2 of the reaction KPIs for the OCM reactions. Comparatively screening the catalysts based on the C2 yields in the locus of their S-X performance curves from multiobjective optimization, with respect to the reference catalyst reveals distinct groupings as seen in Figure 5. The figures indicate the RMS distance of the points on the performance curves of the catalysts from that of the reference, also the standard deviation of the performance curve C2 yields of the reference within which those of the other catalysts lie, and also the correlation of the same with the reference. The points are annotated with the maximum C2 yields of the performance curves that have been tabulated with their associated operating conditions in Appendix B.

These RF model-based screening plots reveal insights into interactions between the mixed metal oxides and the supports that contribute to surface properties linked to catalyst activity and C2 selectivity. Na-Mn-W oxides supported on $Si{O}_2$ is subject to severe restructuring owing to Na-induced phase change crystallization, but its superior thermal stability at high temperatures required for OCM reactions has made it a popular baseline [39]. The use of supports with different acid, basic or amphoteric nature to explore its synergies with Mn and W for redox cycles for $O_2$ activation and $C{H}_4$ activation, respectively can easily be deduced from Figure 5a. The use of acidic metal oxide supports like $Ti{O}_2$ acts as a sink for alkali metal dopants, and stabilizes active species ($Mn{O}_x$, $W{O}_x$) to facilitate lower temperatures to activate gas phase oxygen, thereby suppressing undesirable further oxidation to form $C{O}_x$[40]. The non-selective oxidation during high temperature exothermic reactions when acidic metal oxide supports are used is mitigated by using less acidic supports like $Si{O}_2$ instead of $A{l}_2{O}_3$, alkali metal promoters ($M_2$) to neutralize acidic sites in the support, altogether using basic metal oxides for supports (MgO, BaO, CaO) despite their inability to stabilize $W{O}_x$ required for $C{H}_4$ activation, and even by using an inert gas stream like Ar to dissipate hotspots [37]. $Si{O}_2$ and alumino silicate-based materials like zeolite supports are known to phase transform by crystallization due to the alkali metal dopant (Na), which is conducive for the dispersion and stabilization of the $W{O}_x$ active species [41]. The phase change of the support is seen to cause a drastic decrease in the surface area for $A{l}_2{O}_3$, $Zr{O}_2$ and $Si{O}_2$, but a sharp increase for $SiC$ supports that have an added benefit of thermal stability [38]. The surface area and porosity of the supports correlates positively with C2 yield owing to better dispersion of active sites, as with zeolite supports [42]. However, if the supports are highly porous and also acidic in nature, alkali metal dopants are used to poison surface sites to limit excessive unselective oxidation [37]. It must be noted that neither the role of active sites nor individual supports is clear in OCM chemistry because even blank tests have shown good activity and C2 selectivity leading to questions about the actual contribution of the support to the chemistry [36]. However, it must be pointed that surface area, porosity, acidity/basicity, thermal stability and phase change amenability of the different supports for Na-Mn-W oxides justify the insights deduced from the patterns of the screening plot in Figure 5a.

The screening plot in Figure 5b presents a comparative assessment of the C2 yields from the S-X performance loci for different mixed metal oxide catalysts on $Si{O}_2$ supports with respect to the reference catalyst. The patterns reveals a coupling between the electronic properties of the metal elements and the OCM reaction performance. Host oxides doped with alkali/alkaline earth metals are reported to increase C2 selectivity due to their low electronegativity and ionization energy [11], with mainly Na known to induce crystallization of $Si{O}_2$ supports, while most of these dopants are known to distort the active site ($W{O}_x^{2-}$) for methane activation [37]. But the stability of these dopants at the harsh reaction conditions is challenged, as evidenced by catalyst degradation when the highly volatile (low melting point metal dopant) Li is used because of which either promoted lanthanide oxides, or alkali dopants with higher melting points are generally preferred for OCM reactions. The positive correlation between the conductivity of the M1 metals with catalytic performance has been discussed, with Manganese (Mn) having the high electrical conductivity and hence preferred [43]. The type of M1 oxide and its oxidation state is known be impacted by the nature of the oxo anion, and the choice of tungsten (W) and Molybdenum (Mo) as the M3 elements are discussed to have reasonable catalytic performance [43]. In line with this discussion, the three major groups identified in Figure 5b are seen to comprise alkali/alkaline earths, transition metals, and lanthanides.

3.3. Proposed Candidates across Combinations of Catalysts and Operating Conditions

Different metal elements in the mixed metal oxides, their molar proportions, interactions with the support they are impregnated on, and even the operating conditions are known to exhibit synergies impacting reaction performance, from the discussions thusfar. Most literature, until recently had deemed to achieve the best performance only if Mn, Na or K, and W were present [43]. The combinatorial explosion in navigating the design space of these catalyst and operating condition descriptors to create serendipity in proposing new candidates, is demonstrated by goal-driven sampling discovery via Bayesian optimization for C2 yield maximization as shown in Figure 6a. Alternatively, using the NSGA-2 genetic algorithm to solve multi-objective optimization in the decision space of catalyst and operating condition descriptors to maximize methane conversion and C2 selectivity, was seen to take longer to converge to values of maximum C2 yields, as seen in Figure A5.

The descriptor values of the 50 individuals (annotated with black numbers in increasing order of the C2 yields) at the end of 3000 solution generations have been visualized in a 2d t-SNE plot of Figure 6b. Additionally, the descriptors corresponding to the best C2 yield among the samples picked by the acquisition function in each iteration of the Bayesian optimization have also been visualized in the 2d t-SNE plot, and are annotated in pink by their sampling indices in increasing order of the C2 yields. The descriptors corresponding to the experimental dataset are also visualized on the t-SNE plot, and all the points are shaded by their associated C2 yields predicted by the RF models. The catalysts proposed by the multiobjective optimization was found to have a narrow field focus in proposing candidates, as most of the catalysts proposed were $MnW{O}_4$ doped with alkali/alkaline earth metals on either SiC, SiCnf or $Si{O}_2$ supports at their respective optimal operating conditions. Hence, only the descriptor candidates corresponding to the best C2 yields among the sampling iterations from Bayesian optimization have been presented in Table 1. It can be seen that a number of new candidates like mixed metal oxides of transition metals on a wide variety of supports have been proposed to have comparable performance as the best performing Mn-Na-W family of baseline catalysts. However, no new mixed metal lanthanide oxides have been proposed. Lanthanide group elements have been reported to hinder the exposure of $W{O}_x^{2-}$, and thereby lower catalytic activity [11]. The feasibility of synthesizing the tabulated catalysts, characterizing them and designing experiments around the specified optimal operating conditions, encourages goal-driven approaches to experimentation in the future.

4. Conclusions

This work develops a descriptor-based random forest regression model that maps to the reaction KPIs of methane conversion and C2 selectivity furnished in the HTE dataset for OCM chemistry, spanning a wide variety of catalysts and operating conditions. The synergies among the catalyst constituents and operating conditions impact the conversion of reactant methane to selectively form ethane and ethylene by suppressing their undesirable total oxidation to form $C{O}_x$. Ethane and ethylene are viable platform chemicals in polymer processing. The feasibility of their industrial production necessitates the OCM reactions to be designed with catalysts and operating conditions resulting in C2 yields greater than 30%. Most literature and trial-and-error experimental efforts have fallen short. Also, the maximum reported C2 yield in the HTE dataset in this work is $20$%, within the limits of the total carbon balance. This limits the models trained herein to make higher predictions. However, the focus is to deploy the models in (i) screening catalysts to identify relationships between catalyst properties and reaction KPIs, (ii) optimizing the operating conditions for catalyst formulations in the HTE datasets to maximize both methane conversion and C2 selectivity, and (iii) proposing new catalysts and their optimal operating conditions by a goal-driven Bayesian optimization for C2 yield maximization to guide future experimentation. The ability of the RF models to capture lumped kinetics has been validated and is shown to reveal patterns among $Si{O}_2$ supported catalysts, and among $MnW{O}_4$ catalysts across different supports, in the space of the estimated kinetic parameters. The RF models were found to improve the C2 yield by ∼ 15% on average when used to optimize operating conditions for catalysts in the HTE dataset to meet both methane conversion and C2 selectivity targets. Screening the catalysts in the space of the best performance achieved across a range of operating conditions along the S-X curves are found to reveal similar patterns, as in the space of the kinetic parameters estimated by the RF model. A number of transition metal oxides on different supports have been proposed by the Bayesian optimization routine, but lanthanide metal oxides were not sampled.

Appendix A.1. Assessing the Random Forest Regression Model

Figure A1. Parity plots between the model predictions for methane conversion and C2 selectivity against the experimental data shaded by the total carbon balance.

Figure A1. Parity plots between the model predictions for methane conversion and C2 selectivity against the experimental data shaded by the total carbon balance.

Figure A2. Random forest regressor predictions are shown to be uncorrelated to the total carbon balance

Figure A2. Random forest regressor predictions are shown to be uncorrelated to the total carbon balance

Appendix A.2. Tabulation of LHS Sample Points

Table A1. Points sampled by LHS from the space of operating condition descriptors

T(°C)	time(s)	$\dot{\mathit{V}}$(ml/min)	${\dot{\mathit{V}}}_{\mathit{C}\mathit{H}4}$ (ml/min)	CH4:O2 (mol:mol)
756.17	0.60	21.62	10.92	15.41
747.83	0.41	21.04	10.42	16.67
751.17	0.53	20.88	10.08	12.74
752.83	0.47	21.63	10.75	12.29
749.50	0.72	21.13	10.58	19.54
754.50	0.66	21.21	10.25	10.70

Table A1. Points sampled by LHS from the space of operating condition descriptors

T(°C)	time(s)	$\dot{\mathit{V}}$(ml/min)	${\dot{\mathit{V}}}_{\mathit{C}\mathit{H}4}$ (ml/min)	CH4:O2 (mol:mol)
756.17	0.60	21.62	10.92	15.41
747.83	0.41	21.04	10.42	16.67
751.17	0.53	20.88	10.08	12.74
752.83	0.47	21.63	10.75	12.29
749.50	0.72	21.13	10.58	19.54
754.50	0.66	21.21	10.25	10.70

Appendix A.3. Fits and Orders of Power-Law Kinetic Parameter Estimation

Figure A3. Regression fits for the estimation of power-law kinetic parameters

Figure A3. Regression fits for the estimation of power-law kinetic parameters

Figure A4. Estimated orders for the consumption of methane (a),(b), and those for the C2 formation (c),(d) for different catalysts.

Figure A4. Estimated orders for the consumption of methane (a),(b), and those for the C2 formation (c),(d) for different catalysts.