Multiscale modeling accelerates materials engineering and discovery through in silico inverse design and optimization. At the microscopic scale, thermodynamic and kinetic parameters like enthalpy of formation, entropy, and activation energy of single atoms or molecules are estimated via ab initio calculations using density functional theory (DFT) (Parr, 1983) and molecular dynamics (MD) coupled with statistical mechanics. Other semi-empirical parameters such as pre-exponential factors are calculated using transition state theory (TST) (Pechukas, 1981). These parameters are then embedded into ordinary differential equations (ODEs) to represent elementary reaction rates. Using domain knowledge or Monte Carlo (MC) techniques, the ODEs are then combined to form MK models (Hansen et al., 2007) that describe elementary, surface-level chemistries of catalytic reactions (Prasad et al., 2010; Stamatakis and Vlachos, 2012; Stamatakis, 2014). At the macroscopic scale, tools such as continuum equations (Hadjiconstantinou and Patera, 1997), computational fluid dynamics (CFD) (Dixon and Nijemeisland, 2001), and mathematical reactor models (Rawlings and Ekerdt, 2002) are used to evaluate scale-appropriate mass and heat transfer correlations. Using these tools, scientists and engineers set material design targets for catalyst development and determine the value of additional data to reduce uncertainty. Detailed reactor and process design frameworks further enable advanced control (Murase et al., 1970; Elnashaie et al., 1988; Gentric et al., 1999; Abel et al., 2000; Hagh, 2003), reactor arrangement and configuration optimization (Hillestad, 2010; Rafiee and Hillestad, 2012, 2013; Manenti et al., 2014; Fischer and Freund, 2021), and reactor network optimization through superstructure optimization (Balakrishna and Biegler, 1992; Kokossis and Floudas, 1994; Lakshmanan and Biegler, 1996; Gong and You, 2018).

Several modeling frameworks have been proposed over the last 20 years for the multiscale development of catalytic processes. Molecular-scale tools such as quantum mechanics, molecular dynamics, and kinetic Monte Carlo (KMC) for property estimation, coupled with data science tools including parameter estimation and identifiability analysis form the basis for a “multiscale simulation ladder” (Vlachos et al., 2006). Raimondeau and Vlachos (2002) demonstrated the computational feasibility of multiscale simulations by generating data at the molecular level (MC simulations and other surface simulators) to fit models for macroscopic reactor modeling. Sutton et al. (2018) successfully established feedback loops between first-principles kinetic Monte Carlo (KMC) simulations and computational fluid dynamics (CFD) to extract chemical insights for rate-determining steps at the molecular-scale and optimal reactor operation mode at the process-scale. Their framework accurately predicts reactor-scale selectivity and conversion under unexplored conditions for the catalytic oxidation of carbon monoxide. Partopour and Dixon (2019) developed a multiscale approach for transient simulation of fixed bed reactors using CFD to study the overall transient behavior of the reactor as well as local effects such as species distribution inside the particles under dynamic conditions. Quiceno et al. (2006) incorporated the MK model for catalytic partial oxidation of methane over platinum catalyst into CFD models and were able to perfectly match experimental results from literature while Maestri et al. (2008) proposed a hierarchy of validated kinetic models for steam and dry reforming of methane starting from MK models to two-step reduced models fit for implementation in reactor design. More recently, Karst et al. (2015) and Partopour and Dixon (2016) used MK model reduction strategies to mitigate the computational demand associated with the incorporation of MK models in process-scale reactor design. These reduction strategies are reviewed in detail below. While these frameworks enable some multiscale design, they were demonstrated using relatively simple reaction networks such as methane reforming and water gas shift reaction that involve O ( 10 ) elementary steps; how to scale these approaches to more complex reaction systems with multiple products or O ( 1 0 3 ) steps remains an open question. Moreover, in-depth knowledge of both micro- and macro-scopic phenomena is required to successfully implement these modeling frameworks.

We now review three categories of approaches—MK model reduction, machine learning emulators, and ROK models—to incorporate MK information in macroscale models (Solle et al., 2017).

The goal of MK model reduction is to create a simplified reaction network that emulates the predominant model trends (e.g., conversion) by combining and removing elementary reactions from the original (large) network. Sensitivity analysis followed by principal component analysis (PCA) (Wold et al., 1987) is the most widely-used strategy for MK model reduction in literature (Mhadeshwar and Vlachos, 2005; Salciccioli et al., 2011a; Karst et al., 2015; Partopour and Dixon, 2016; de Carvalho et al., 2018). In this method, the significance/importance of each elementary reaction step in the MK model is quantified using sensitivity matrices that record the change in MK model performance (e.g., conversion, product yield, selectivity) resulting from the one-at-a-time removal of the elementary steps from the model. PCA is then performed using the sensitivity matrix to eliminate the least sensitive elementary steps based on eigenvalue analysis from the network until a threshold is met. Recent successful applications of PCA-based MK model reduction include analysis of the MK model for water gas shift reaction on platinum catalyst consisting of 46 elementary steps (Mhadeshwar and Vlachos, 2005) and on nickel catalyst consisting of 19 reversible elementary steps (de Carvalho et al., 2018), ethylene hydrogenation and ethane hydrogenolysis MK model consisting of 32 reversible elementary steps (Salciccioli et al., 2011a), and partial ethylene oxidation MK model with 17 reversible steps (Partopour and Dixon, 2016). Karst et al. (2015) used the PCA-based approach followed by reaction pathway analysis to reduce the methane chemistry model from Maestri et al. (2009) from 41 elementary steps to 21 steps with an overall average error of 0.05% in predictions. The reduced-order model was then used to maximize the yield of hydrogen by optimizing the feed composition and reactor temperature profile and the results were validated to predict hydrogen production within 0.13% of the predictions obtained by simulating the MK model at the optimal conditions. Our work scales these approaches to multi-product systems with two orders of magnitude more elementary reaction steps.

Backward feature elimination (BFE) is another popular approach for MK model reduction. To start, BFE creates a library of candidate reduced MK models by sequentially removing a single reaction from the model in each step and recalibrating the reduced MK model using experimental data. The reduced MK models with the largest prediction errors are removed from the library and the process repeats until an acceptance criterion such as minimum model size or maximum prediction error is met. Often, simulated annealing or other heuristic search algorithms are then used to reparameterize the reduced MK models and minimize prediction error compared to the original MK model. Avşar (2017) used this approach to study carbon monoxide conversion from the water gas shift reaction on platinum consisting of 46 elementary steps. Although BFE yielded better-predicting reduced models with 5.8% lower root mean squared error (RMSE) compared to PCA in the study by Avşar (2017), the BFE approach cannot guarantee optimal reduction of MK models because it is a heuristic search.

Despite their potential advantages, PCA- and BFE-enabled MK model reduction strategies, so far, have only been applied to reaction networks consisting of 40–50 reaction rates, providing up to 50% reduction in model size without losing considerable model fidelity. In this work, we consider the HZSM-5-catalyzed propylene oligomerization reaction model fully described by Vernuccio et al. (2019), which contains 4,243 reactions, 2,345 net reaction rates, and 909 state variables (ODEs). These ODEs, defined as the partial pressure of each species with respect to the number of catalyst active sites along the reactor, were integrated using the double precision differential/algebraic sensitivity analysis code (DDASAC) solver (Warren, 1995) as an initial value problem. Each MK model integration takes a few minutes to complete. Additionally, large MK models are often sensitive to the initial operating conditions (Becerra et al., 2021), i.e., small perturbations in initial conditions can cause large changes in the output. Therefore, we anticipate a MK model embedded in process flowsheet with recycle loops would need to be integrated tens to hundreds of times to converge the flowsheet in sequential modular simulation environments (e.g., AspenPlus) which would make process simulation or optimization cumbersome and computationally intensive. Finally, the above-described MK model reduction methods are intrusive by nature and require nontrivial modifications to MK model/code to implement. For this reason, practitioners often prefer using modeling strategies, if possible, which require less software modifications and in-depth MK modeling expertise.

Alternately, machine learning methods can be used to train surrogates such as polynomial regression models (Wang et al., 2018), artificial neural networks (ANNs) (Teske, 2014; Miriyala et al., 2016; Kotidis and Kontoravdi, 2020), or distribution models (Hutter et al., 2021) to emulate the behavior of MK models. For example, Wang et al. (2018) correlated combustion reaction data at different operating conditions using neural network surrogate models. They introduced constraints based on the proximity of model response to available data to avoid model extrapolation for data selection, experimental design, and property estimation. Miriyala et al. (2016) replaced the detailed kinetic model for polymerization with neural networks in multi-objective reactor optimization problems which reduced the number of function calls by up to 90% and, in turn, reduced the solution time. Recently, Bradley and Boukouvala (2021) successfully demonstrated up to 67% reduction in computation time without loss of prediction accuracy using a two-step approach for reaction kinetics parameter estimation. First, they trained neural ordinary differential equations (NODEs) to represent state derivative information and then used the NODEs for parameter estimation. These and similar results highlight two key benefits of surrogate models: they require limited knowledge of reaction kinetics and have the potential to improve computation time. However, model extrapolation using purely machine learning-based (surrogate) models remains challenging due to the absence of model physics (Von Stosch et al., 2014). Trust region methods are one way to overcome this limitation (Eason and Biegler, 2016). For example, Franzoi et al. (2021) use a trust region based adaptive sampling strategy to develop highly accurate multiple-input single-output (MISO) surrogate models to replace Arrhenius-form kinetics in the Williams-Otto optimization problem (Williams and Otto, 1960). Using this approach, they iteratively refined the surrogate model and performed optimization, ultimately obtaining solutions within 0.016% of the best known objective value.

Hybrid models are emerging to overcome the limitations of both first-principles and purely data-driven models. The seminal work of Psichogios and Ungar (1992) combines a first-principles model and ANNs to model the operation of a fedbatch reactor; the first-principles model emulated known reaction kinetics while the black-box ANN was used to model dynamic behavior arising from unknown interactions between living organisms. Since this seminal work, neural networks have remained the preferred form of surrogate that are coupled with some form of first-principles model in hybrid modeling (Qi et al., 1999; Zahedi et al., 2005; Mosavi et al., 2019; Bangi and Kwon, 2020). Von Stosch et al. (2014) provide a review of recent advances in coupling non-parametric surrogates and first-principles models. More recently, Yang et al. (2020) used a hybrid approach to combine a first-principles, rate-based model with neural network surrogate models that emulates fluid catalytic cracking reactions and improves prediction accuracy by 9% compared to stand-alone neural network models. Tsopanoglou and del Val (2021) combined mechanistic models for material balances and kinetics with neural network surrogate models for in situ bioprocess monitoring. De Jaegher et al. (2021) use neural differential equations to model colloid-membrane collisions, coupled with mechanistic models, to study fouling in electrodialysis membranes. The neural differential equations elucidate missing physics for ion-exchange membranes and motivate the need for further research towards developing fully-mechanistic models for membrane fouling. Foad et al. (2022) trained hybrid models—reduced-order model derived using PCA combined with a deep neural network (DNN)—to emulate the transient behavior of a nuclear reactor and were able to emulate the output power trajectory of the reactor with only 2% deviation. While hybrid models are promising, they are often demonstrated using reaction networks with only tens of reactions (Potočnik et al., 2003; Bhutani et al., 2006; Bayer et al., 2020).

ROK modeling strategies, a proven practice in literature (Jacob et al., 1976; Coxson, 1995; de Andrade Lima and Hodouin, 2005; Radmanesh et al., 2006; Oliveira et al., 2010), use a priori knowledge of the reaction mechanism such as the rate form (Laidler, 1984) to postulate and train first-principles based lumped-rate models. These models, once validated, can be extrapolated to make conversion and product distribution predictions. For example, Zong et al. (2010) proposed an ROK model for the heavy oil catalytic cracking and maximizing of iso-paraffins (MIP) process and validated the model using industrial data. In addition to predicting product distribution and product quality with 7.5% error, the validated model was able to predict expected kinetic behavior for the MIP process. Similarly, Ying et al. (2015) proposed an ROK for the catalytic methanol-to-olefins (MTO) process which, after calibration, predicted conversion of major olefins within 5% of experiments. More recently, Tsu et al. (2019) developed a library of candidate rate equations and used Monte Carlo algorithms to sample from the candidate equations to create kinetic models. Wang et al. (2021) demonstrated the tractability of ROK models in reactor design by embedding an ROK model for CaO-CO₂ carbonation in a bubbling bed reactor model; simulations indicated that the ROK model successfully emulate temperature, concentration, and particle size effects observed in laboratory experiments.

ROK models are especially useful for modeling multi-product reaction networks. Jacob et al. (1976) first presented the idea of developing ROK models for complex reaction networks such as fluid catalytic cracking (FCC) to obtain a more “fundamental basis” to relate feed composition and process variable effects to product distribution. Similarly, Tabak et al. (1986) proposed the use of ROK models to model and study catalytic oligomerization. Based on the simplified mechanism proposed by Tabak et al. (1986), Borges et al. (2007) and Oliveira et al. (2010) developed ROK models using rate expressions lumped by the carbon number of the species that are capable of emulating low conversion product distributions for ethylene, propylene, and butene oligomerization. Owing to their practical applicability towards modeling complex, multi-component chemical reaction systems, ROK models continue to be the preferred modeling strategy for systems such as catalytic cracking (Ebrahimi et al., 2018; Sani et al., 2018; Sun et al., 2018; John et al., 2019) and pyrolysis (Yan et al., 2019; Schubert et al., 2020; Zhang and Sarathy, 2021). Beyond modeling complex petrochemical processes, ROK models have also been developed for a variety of chemical processes including gold ore cyanidation (de Andrade Lima and Hodouin, 2005), biomass pyrolysis (Radmanesh et al., 2006), and solid fuel combustion (Moríñigo and Hermida-Quesada, 2016).

We focus on the HZSM-5 catalyzed propylene oligomerization MK model developed by Vernuccio et al. (2019). The elementary steps in this model involve physisorption of a gas-phase olefin into the pores of the zeolite followed by protonation by a Brønsted acid site to form a reaction intermediate (carbenium ion or covalent alkoxide). The addition of a physisorbed olefin to a protonated intermediate results in the formation of a higher carbon number intermediate which finally deprotonates to form a heavier product olefin (Vernuccio et al., 2019). Following the formation of oligomers resulting from dimerization and subsequent oligomerization of the olefin monomer reactant (Sarazen et al., 2016), skeletal isomerization and cracking via β-scission occur, resulting in a mixture of olefins from ethylene to nonene. The reaction network for this model is developed using an automated network generator, NetGen (Broadbelt et al., 1994). The automated generation process for an oligomerization network is in principle infinite because of the nature of the oligomerization steps which lead to the formation of heavier species. However, the experimental data collected in Vernuccio et al. (2019) for model validation showed a negligible amount of C_>9 species in the low pressure regime. For this reason, the reaction network was limited to molecular species and reaction intermediates with carbon number lower or equal to 9. Vernuccio et al. (2019) track 909 species through 4,243 reactions and 2,345 net rates to model the formation of species with two to nine carbon atoms with a pure propylene feed (75% propylene and 25% inert). The kinetic parameters for each elementary step were estimated using transition state theory, Evans–Polanyi relationships, and thermodynamic data. Originally developed for low-pressure and low-conversions operations, the MK model assumes the rate of hydride transfer to be negligible and does not track conversions to aromatics and paraffins. This assumption is justified when the reaction system primarily includes small olefins with only few abstractable hydrogens. This prevents hydride transfer steps from occurring which in turn limits the production of paraffins. In order to verify the applicability of the MK model to higher pressure operations, we developed an extended version of the model which includes hydride transfer steps between molecular species and carbenium ions/alkoxides as well as cyclization steps which lead to the formation of C₅ and C₆ ring structures. This extended MK model, including 4,212 species and 23,086 reactions, is substantially larger and less manageable than the MK model adopted in this work and was used for verification purposes only. While the MK model from Vernuccio et al. (2019) takes a few minutes to integrate, the extended MK model takes almost 6 hours to integrate. However, Figure 2 shows that the selectivity for butene and heavier olefins is comparable between the 2 MK models. Therefore, to reduce computation costs, we consider the MK model from Vernuccio et al. (2019) as the “ground truth” model in this work.

In this work, we use data generated by simulating the MK model (Vernuccio et al., 2019) in a PBR at high conversion range to reparametrize the ROK models. We divide data from 73 MK simulations into two groups: 1) 64 low pressure simulations at pressure 2.18 atm, temperatures 483–522 K, and space velocities 62 to 740 m o l − f e e d m o l − H + − 1 m i n − 1 previously published in Ganesh et al. (2020) and 2) nine high pressure simulations at pressures of 10–30 atm, temperatures of 523–573 K, and space velocity of 1123 m o l − f e e d m o l − H + − 1 m i n − 1 .

We postulate a library of six candidate ROK models. Model M0 is a simple kinetic model with Arrhenius-form rate expressions. Model M1 is an ROK model adopted from literature (Oliveira et al., 2010). Model M2 updates the adsorption enthalpy chain length dependence and adsorption isotherm model from M1. Model M3 adds temperature dependence to the frequency factor from model M2. Models M4 and M5 add activation energy chain length dependence to models M2 and M3, respectively. The sequential structural changes introduced in the ROK models from M0 to M5 are intended to ultimately match the MK model structure presented by Vernuccio et al. (2019). In this paper, we use the phrase “functional form” to refer to the dependence of kinetic parameters on reaction conditions, including temperature and adsorption phenomena, and reactant characteristics, including chain length and chain symmetry.

To calibrate parameters in each candidate model, we minimize a linear combination of the sum of squared errors (SSE) and sum of squared percent errors (SSPE) between olefins concentrations predicted by the MK and ROK models in a PBR: min θ S S E + w S S P E s . t . r a t e l a w f r o m T a b l e (1) m a t e r i a l b a l a n c e f o r P B R f r o m E q . ( 20 ) b o u n d s (21) where w is the weight assigned to SSPE. θ is the vector estimated kinetic model parameters, which are defined as follows for models M0 through M5: θ ₀ = [α _1,j , α _2,j , E _1,j , E _2,j ] for all j ∈ {4, …, N}, θ ₁ = [ α 1 ′ , α 2 ′ , γ ₁, γ ₂, ϕ ₂, λ, ΔH′, E 1 ′ , E 2 ′ ], θ ₂ = [ α 1 ′ , α 2 ′ , γ ₁, γ ₂, ϕ ₂, α _ads , β _ads , E 1 ′ , E 2 ′ ], θ ₃ = [ α 1 ′ ′ ′ , α 2 ′ ′ ′ , γ ₁, γ ₂, ϕ ₂, α _ads , β _ads , E 1 ′ , E 2 ′ ], θ ₄ = [ α 1 ′ ′ , α 2 ′ ′ , γ′, δ, E ₀, κ ₁, κ ₂, α _ads , β _ads ], and θ ₅ = [ α 1 ′ ′ ′ ′ , α 2 ′ ′ ′ ′ , γ′, δ, E ₀, κ ₁, κ ₂, α _ads , β _ads ]. Additionally, we analyze the effects of model scaling and objective function SSE formulation on model fit quality to choose the best modeling strategy. The formulation and details of this analysis are reported in Supplementary Material.

The models were implemented in Pyomo (Hart et al., 2017) and discretized using backward finite-difference scheme via Pyomo.DAE (Nicholson et al., 2018). The resulting regression problems were formulated using Pyomo.parmest (Klise et al., 2019) and contain 22,192 (M0), 23,214 (M1, M2, and M3), and 21,097 (M4 and M5) variables and 22,168 (M0), 23,205 (M1, M2, and M3), and 21,088 (M4 and M5) constraints. The problems were solved using with IPOPT (Wächter and Biegler, 2006) and the HSL linear solver MA27 (HSL).

Kinetic models are often only partially identifiable (Vajda et al., 1989). A model is said to be identifiable if there exists a unique set of model parameter values that produce a unique model response. Conversely, a model lacks identifiability if multiple distinct sets of model parameter values produce the same model response (Seber and Wild, 1989). In nonlinear regression problems, the model being reparameterized is not identifiable when the least-squares regression objective is flat in one or more directions (related to fitted parameters) and the Fisher Information Matrix (FIM), which is proportional to the reduced Hessian for the regression problem, is (near) singular. The FIM, M , is also approximately the inverse of the parameter covariance matrix, V . In this work, we ignore prior information (Franceschini and Macchietto, 2008) when approximating M as: M θ ̂ , ϕ ′ ≈ V θ ̂ , ϕ ′ − 1 (22a) M θ ̂ , ϕ ′ ≈ Q θ ̂ , ϕ ′ ⋅ Σ p − 1 ⋅ Q T θ ̂ , ϕ ′ (22b) M u , u ′ θ ̂ , ϕ ′ ≈ ∑ j ′ = 1 N − 1 ∑ j ″ = 1 N − 1 ∑ m c = 1 M c ∑ m c ′ = 1 M c q j ′ , m c , u σ ̃ j ′ , m c , j ″ , m c ′ q j ″ , m c ′ , u ′ (22c)

In Eq. 22, the FIM M is evaluated using nominal parameter valuesparameter values θ ̂ at operating conditions ϕ′ . Recall, the parameter covariance matrix V represents the uncertainty in parameter estimates obtained from the least-squares regression problem. Classical nonlinear regression theory assumes the measurements are corrupted with random independently and identically distributed (i.i.d.) zero-mean observation uncertainty (Bates and Watts (1988)). Let the measurement covariance matrix Σ _p correspond to such multivariate Gaussian uncertainty. Furthermore, the sensitivity matrix Q contains the partial derivatives of each model prediction with respect to regressed model parameter. Thus, Eq. 22b is a first order approximation to propagate the measurement covariance Σ _p through the regressed model to obtain the parameter covariance V . Likewise, when the model residuals are small, Eq. 22b approximates the inverse of the reduced Hessian of the regression problem (Bard, 1974). Eq. 22c is the equivalent calculation expressed with summations. Here the pairs (j′, mc) and (j ^″, mc’) index the partial pressures predicted by the model. Specifically, j′ and j ^″ indicate the olefin species with carbon numbers j′ + 1 and j ^″ + 1 (for example, j′ = 1 corresponds to ethylene). Likewise, mc and mc’ indicate discretization points along the reactor, in terms of scaled catalyst mass per Eq. 20, at which model responses are recorded. The regressed parameters are indexed by u and u′. σ ̃ ( j ′ , m c ) , ( j ′ ′ , m c ′ ) is the element of the inverse of Σ _p that corresponds to the covariance between measurement (j′, mc) and (j ^″, mc’). Likewise, q _(j′,mc),u corresponds to the element of Q for the sensitivity of model response p _j’ measured at point mc (or mc’) with respect to model parameter θ _u : q j ′ , m c , u = ∂ p j ′ , m c ∂ θ u θ ̂ (23)

These calculations are performed using the Pyomo.DOE package (Wang and Dowling, 2022). We approximate the inverse of the measurement uncertainties σ ̃ ( j ′ , m c ) , ( j ′ ′ , m c ′ ) using the sum of squared residuals SSE from the regression problem in Eq. 21 as follows: σ ̃ j ′ , m c , j ″ , m c ′ ≈ S S E N D − N θ − 1 j ′ = j ″ and m c = m c ′ 0 otherwise (24) where N _D is the total number of data points used to calibrate the model and calculate SSE. N _θ is the total number of model parameters and defined in Table 1 for each ROK model.

If the FIM M approximated using Eq. 22c is (near) singular, eigendecomposition of M , can be used to identify the direction of “sloppiness” (Chis et al., 2016), i.e., the parameter(s) which can assume multiple values without changing the model response and least-squares regression objective value.

In Eq. 25, FIM M is decomposed into eigenvectors U = [ u 1 , … , u N θ ] and corresponding eigenvalues λ 1 , … , λ N θ .¹ Λ is a diagonal matrix with the eigenvalues on the diagonal. By definition, the eigenvectors u k ∈ R 1 × N θ are unit vectors in the parameter space. Recall, the FIM M is proportional to the reduced Hessian of the least-squares optimization problem. A (near-)zero eigenvalue indicates that the regression objective function has (near-)zero curvature in the direction of the corresponding eigenvector, which is a consequence of the model predictions being insensitive to perturbations in said direction. When interpreting the eigendecomposition, if λ k ≈ 0 and the corresponding u k predominantly points in the direction of a single parameter, said parameter is sloppy.

In this work, we use a sequential approach to identify ROK model sloppy parameters. In this approach, once a sloppy parameter is identified using Eq. 25, the parameter is removed as a decision variable from the regression problem and is fixed to a previously-obtained value that preserves model feasibility. Then the regression problem is re-solved, the FIM is recomputed for eigenvalue analysis, and the process is repeated until the FIM is no longer singular and the model is identifiable. This algorithm is further explained in Section 5.3.

Although an ill-posed model with sloppy parameters can be used for parameter estimation, the uncertainty in resulting parameter estimates is extremely large. Adding measurements (data) or removing (i.e., fixing) sloppy parameters as decision variables often makes the model identifiable and reduces the uncertainty of parameter estimates. This is a local analysis because M is approximated at specific values of θ via Eq. 22; repeating the analysis at different values of θ may reveal different sloppy parameters.

Using the calibrated ROK models, we compute the optimal temperature profiles for a staged PBR. In place of one PBR with volume V, catalyst mass M _cat , and temperature T, we consider N _Z isothermal PBRs, each with catalyst loading M _cat /N _Z , in series with inter-stage cooling in the staged PBR design, as shown in Figure 3.

Following from the PBR governing equation, Eq. 20, the design equation for an olefin with j carbon atoms at each stage can be written as: N Z d p j d m cat T z − P p inert S V n cat p inert 0 ∑ n = 2 j − 2 ∑ m = 2 j − n − r n , m | T z = 0 , ∀ z ∈ 1 , … , N z (26)

Under these assumptions, we propose a staged-PBR temperature optimization problem to maximize the formation of high-value olefins, with carbon numbers j _min to N: max T z ∑ j min ≤ j ≤ N p j z = N Z s . t . d e s i g n e q u a t i o n f r o m E q . ( 26 ) m a t e r i a l b a l a n c e b o u n d s (27)

In Eq. 27, the system operates at constant pressure and therefore, through ideal gas law, maximizing the partial pressure of select olefins maximizes the molar flow, i.e., formation of those olefins. The temperatures T _z for stages z ∈ {1, 2, …, N _z } are the design degrees of freedom. The number of reactors in series is varied while keeping the total available reactor volume constant. The problem is initialized using the 1-stage or previous best solution. For example, the 2-stage solution is used to initialize the 4-stage problem. Likewise, 1-, 3-, and 5-stage solutions are used to initialize the 5-, 6-, and 10-stage problems, respectively. The differential-algebraic equation (DAE) system was discretized along the reactor length using backward finite difference via Pyomo.DAE (Nicholson et al., 2018) and the resulting nonlinear optimization problem consists of 1,674 variables and N _z degrees of freedom. The problem is defined in Pyomo (Hart et al., 2017) and solved using IPOPT (Wächter and Biegler, 2006) with the HSL linear solver MA27 (HSL). Each optimization problem solved in under 14 seconds, including initialization.

The uncertainty in calibrated ROK model parameters induces uncertainty in the objective function of the staged reactor problem. We use the Uncertainty Propagation Toolbox in the IDAES-PSE framework (Lee et al., 2021) to propagate this parametric uncertainty through Eq. 28 using the following first-order error propagation formula (Rooney and Biegler, 2001): Var f x ∗ , θ ̂ ≈ ∂ f ∂ θ + ∂ f ∂ x ∂ x ∂ θ ⋅ V ⋅ ∂ f ∂ θ + ∂ f ∂ x ∂ x ∂ θ T (28)

In Eq. 28, V is the parameter covariance matrix calculated in Section 4.4, x^∗ is the optimal solution and f ( x ∗ , θ ̂ ) is the corresponding objective function value for the staged reactor problem obtained using the parameters θ ̂ . The nonlinear programming sensitivity code k_aug (Thierry, 2019) is used to calculate the gradients in Eq. 28.

To systematically train, evaluate, and compare ROK models from the library (Table 1), we performed weighted nonlinear regression (Eq. 21) using the MK simulation data described in Section 4.1. Table 2 reports quality of fits metrics R ², mean squared error (MSE), mean absolute error (MAE), and mean squared logarithmic error (MSLE) for the six ROK model alternatives.

For ROK model training with MK simulation data, MSLE is the most informative metric. MK simulation data used to train the ROK models consists of olefin partial pressure values ranging from O ( 10 − 8 ) to O ( 1 0 7 ) Pa. As Table 2 shows, R ², MSE, and MAE values for models M0 to M5 are almost identical and vary by less than 0.2%, 0.1%, and 0.05%, respectively, for the six candidate ROK models. The R ², MSE, and MSE metrics, therefore, are most influenced by olefins with higher partial pressures (propylene, butene, pentene, and hexene); olefins with lower partial pressures (ethylene, heptene, octene, and nonene) do not significantly affect these metrics. Moreover, the geometrical interpretation of R ² does not hold for nonlinear regression (Miaou et al., 1996; Spiess and Neumeyer, 2010). The species-wise MSE and MAE values for the remaining species are two and four orders of magnitude smaller, respectively, than the total values and, therefore, a log-scaled error metric such as MSLE is needed to assess model fit quality for olefins with low and high partial pressures. The MSLE values have the same order of magnitude for all olefin species, regardless of partial pressure. MSLE shows up to 33.5% difference in fit quality between the ROK models M5 (lowest MSLE) and M1 (highest MSLE). The MSLE of models M0, M2, and M3 are within 3% of each other and within 25% of models M4 and M5. The literature model M1, with its distinct kinetic functional forms (compared to the MK model as discussed in Section 4), has the largest MSLE value of all six ROK models which indicates it has the worst fit to MK simulation data. The fitted parameter values for the six ROK models are reported in Supplementary Tables S1-S6 in Supplementary Material. In ROK modeling literature, it is most common to qualitatively assess the quality of fit (Oliveira et al., 2010; Ying et al., 2015; Wang et al., 2021; Zhang and Sarathy, 2021) or report relative error metrics using experimental data (Radmanesh et al., 2006; Ebrahimi et al., 2018; Sani et al., 2018; Sun et al., 2018; John et al., 2019). In some cases, a metric such as R ² (Yan et al., 2019) is reported, but it is challenging to meaningfully compare such fit metrics across different models and datasets.

Figures 4, 5 confirm that M1 fails to emulate MK model behavior and highlight the importance of model functional form and model assumptions (see Table 1). By virtue of having unique pre-exponential factor and activation energy values fitted for each olefin species and reaction type, M0 emulates MK simulations with similar accuracy as models M2 to M5. Although model M1 was extended to consider transformations up to nonene, instead of octene, the activation energy chain length dependence (Eqs. (9)-(10)), adsorption enthalpy (fixed), and adsorption isotherm model assumption (Eq. 14) are different compared to the MK model (Vernuccio et al., 2019); consequently, M1 struggles to emulate MK simulations. The introduction of adsorption enthalpy chain length dependence and change in adsorption isotherm model assumptions in models M2 to M5 (and additionally updated activation energy functional form of models M4 and M5), discussed in detail in Sections 4.2.6 and 4.2.7, respectively, and elaborated in Table 1, enable these models to capture MK model conversion trends and magnitudes better than model M1. As a result, Model M1 is not considered a viable model for MK model emulation and excluded from further analysis. There is also a significant gap between the conversions predicted by MK and ROK models M0 and M2 to M5. We hypothesize this is because the ROK models are missing some of the pertinent information (e.g., mathematical expressions in the rate law) and purpose model refinement as future work.

Next, the calibrated ROK models were simulated with a mixed feed composition mirroring the Bakken shale basin (Ridha et al., 2018) and a reactor size estimated to create a space velocity of 25 m o l − f e e d m o l − H + − 1 m i n − 1 (which is reasonable for typical oligomerization reactors) over process-relevant temperatures of 473–1073 K and pressures of 2–40 atm.

Sensitivity analysis shown in Figure 6 indicates that models M3 to M5 are most reliable at process conditions. Since oligomerization is exothermic and the reverse cracking is endothermic, a higher conversion to heavier olefins is expected at lower temperatures and a lower conversion to heavier olefins is expected behavior at higher temperatures. Despite having a good in-sample fit, model M0 is unable to emulate this thermodynamic behavior; it demonstrates the opposite behavior by predicting higher conversions with increasing temperatures. Additionally, IPOPT struggles to converge material balances for the olefins at high temperatures and pressures (top-right corner of Figure 6A). The apparent good fit quality of M0 to MK simulation data, discussed previously, is, therefore, a result of over-parameterization of the ROK model. Referring back to Table 1, M0 is defined using pairwise unique activation energies and pre-exponential factors for each olefin species and reaction type (Eq. 2) and consists of 24 degrees of freedom (number of fitted parameters). This improves the fit quality of M0 to training data (MK simulations) but, leads to poor scaling with temperature and space velocity due to the lack of operating condition dependence of the parameters. Models M2 to M5, on the other hand, each consist of only nine fitted parameters but, capture oligomerization thermodynamics quite well (through the parametric dependencies on operating conditions) across the range of temperatures and pressures and predict maximum propylene conversion at temperatures below 600 K. However, IPOPT once again struggles to converge the olefin material balances in M2 and M3 at higher temperature and pressure conditions (top-right corner of Figure 6B). For models M0, M2, and M3, IPOPT also fails to converge the olefin material balances at higher temperatures and pressures. These convergence issues persist even with alternate model scaling strategies, which are discussed in Section 5.5 and Supplementary Section S1 of the Supporting Material. IPOPT converges successfully across the entire range of operating conditions analyzed for models M4 and M5.

Model identifiability analysis is important to properly characterize the uncertainty in estimated parameters. For models M2 to M5, one or more parameters are sloppy. We use the following sequential approach to identify and remove (fix) the sloppy model parameters:

1) Compute the FIM using Eq. 22

Using this procedure, we analyze model M5. We will use Eqs 29, 30, which are the expanded-form of the rate law from Eq. 8 and Table 1, to provide physical interpretations for the identifiability analysis results. log k 1 , n , m ′ = log α 1 ⁗ + log T T ref + n β ads + α ads − E 0 + κ 1 δ − Δ j γ ′ R T (29) log k 2 , n , m ′ = log α 2 ⁗ + log T T ref + n β ads + α ads − E 0 − 1 − κ 2 δ + Δ j γ ′ R T (30)

Identifiability Iteration 1 for M5: Without fixing any parameter, the condition number of the FIM is 3.49, ×, 10²¹. Supplementary Table S16 in Supplementary Material reports the eigen decomposition of the FIM for M5 with nine fitted parameters; these results show that parameter γ′ is sloppy. γ′ cannot be identified uniquely since Δj in Eqs 29, 30 denotes the change in number of C-atoms between reactants and products and, therefore, takes a value of zero (law of conservation).

Identifiability Iteration 2 for M5: With γ′ fixed, the condition number of the FIM decreases to 3.16 × 10¹⁸. With eight fitted parameters, the eigen decomposition reported in Supplememtary Table S17 in Supplementary Material shows δ is the next sloppy parameter. Referring to Eqs 29, 30, δ forms a bilinear combination with parameters κ ₁ and κ ₂. At most only one parameter in a bilinear combination can be uniquely estimated.

Identifiability Iteration 3 for M5: With γ′ and δ fixed (and seven remaining parameters), the condition number of the FIM decreases to 7.21 × 10¹⁷. The eigen decomposition reported in Supplementary Table S18 shows that E ₀ is the next sloppy parameter. In Eqs 29, 30, E ₀ forms a linear combination with two other fitted parameters—α _ads and κ _i from the bilinear term κ _i δ—out of which only one parameter can be estimated uniquely. E ₀ is, therefore, fixed in this iteration.

Identifiability Iteration 4 for M5: With γ′, δ, and E ₀ fixed, the condition number of the FIM decreased to 7.95 × 10¹⁶. Supplementary Table S19 reveals that κ ₁ is the next sloppy parameter. Inspecting Eq. 29 shows that despite fixing E ₀, κ ₁ in the bilinear term κ ₁ δ forms a linear combination with α _ads (fitted parameter) and, therefore, could not be identified uniquely.

Identifiability Iteration 5 for M5: With γ′, δ, E ₀, and κ ₁ fixed, the condition number of the FIM decreased to 1.56 × 10⁶. Supplementary Table S20 in Supplementary Material shows the eigen decomposition. This is smaller than the threshold of O ( 1 0 8 ) and the identifiability analysis procedure is terminated.

Per Table 1, model M4 is structurally similar to M5. As expected, repeating the identifiability analysis also reveals γ′, δ, E ₀, and κ ₁ as the sloppy parameters for model M4. Fixing these parameters improved the condition number of the FIM from 2.55 × 10²¹ to 2.40, ×, 10⁶. For completeness, Supplementary Tables S11—S15 in Supplementary Material tabulate the eigenvalue analysis for model M4 in sequential order.

We performed similar analysis for M2 and M3. Fixing parameters α 1 ′ (M2) and α 1 ′ ′ ′ (M3) reduced the FIM condition numbers from 3.87 × 10⁹ to 1.47 × 10⁷ and 7.70, ×, 10⁹ to 1.30, ×, 10⁷, respectively. Supplementary Tables S7–S10 tabulate the eigenvalue analysis for models M2 and M3, respectively.

Recalibration of models M2 to M5 with sloppy parameters fixed does not affect model fit quality. This is because the presence of sloppy parameters reduces the confidence of parameter estimates without affecting the overall model fit quality. In this work, all results for models M2 to M5 were generated after fixing the sloppy parameters. As we will see in Section 6, the fixing of sloppy parameters to make models identifiable not only reduces parameter estimate uncertainties but, as a consequence, improves the confidence in results obtained using the updated models.

For identifiability analysis described above, we assume the measurement covariance matrices Σ _p are identity matrices times the following scalars: 9.47 × 10⁻⁴ atm² (M2), 9.46 × 10⁻⁴ atm² (M3), 9.93 × 10⁻⁴ atm² (M4), and 9.92 × 10⁻⁴ atm² (M5).

Figure 7 shows that none of the ROK models, including M4 and M5, accurately recapitulate product distribution at both high and low pressures. The ROK models under-predict propylene conversion at low pressures (Figure 7A) and models M2 and M3 over-predict propylene conversion at high pressure (Figure 7B). Models M4 and M5, which have a low in-sample residual, match propylene conversions well only at high pressure but, over-predict the formation of nonene. The overproduction of C₉ species suggests that the flux of the MK model might be artificially truncated by the termination criterion adopted for the automated network generation. The exclusion of C_>9 from the reaction network results in an accumulation of C₉ species which are not allowed to further oligomerize. As a result, C₉ intermediates in the MK model are forced to deprotonate or generate smaller products via β-scission. It is worth noting that the number of species and reactions that are automatically generated for an oligomerization network increase exponentially as a function of the termination criterion (Vernuccio and Broadbelt, 2019). Thus, the inclusion of heavier species in the reaction mechanism would result in unmanageable models and impracticable computational costs.

The differences observed in Figure 7 for low and high pressure operation can be also ascribed to the different MK model assumptions that were used to generate the low and high-pressure data. In order to investigate this aspect, we used the extended MK model described in 3.1 to simulate the oligomerization process at 523 K and 10–30 atm. Figure 8 shows a parity plot for the product selectivity calculated with the extended model and the MK model used in this work. At 10 atm, the rate of hydride transfer which leads to the formation of paraffins is negligible and the omission of the corresponding reactions from the MK model results in minor deviations from the model performance. However, at 20 and 30 atm, the presence of these additional reactions alters the rates of oligomerization and the product selectivity. Thus, the omission of hydride transfer, cyclization as well as aromatization steps, not included in the MK model for numerical tractability, has an impact on the overall model performance. Due to this implicit issue with the data source, the overall in-sample fit quality of the ROK models is a trade-off across the range of pressure conditions over which the models were calibrated. In addition to simplified chemistry, this trade-off makes it very challenging to train a “one-size-fits-all” ROK model.

Finally, we studied the effect of weight w in the objective of the optimization problem in Eq. 21. Supplementary Figure S2 in Supplementary Material shows the trade-off in SSE and SSPE with varying values of w. A weight w of 0.1 balances the trade-off and produces calibrated models that best emulate the entire product distribution (Supplementary Figure S3). Supplementary Figure S1 in Supplementary Material shows that ROK model predictions are robust to the choice of model scaling, i.e., partial pressures as shown in Eq. 1, or mole fractions as shown in Supplementary Equation S2. Supplementary Figure S2 in Supplementary Material shows that model fit for the ROK models is also robust to SSE objective function scaling options shown in Supplementary Equation S3,S4.