Model-based experimental screening for DOC parameter estimation

Computers and Chemical Engineering 74 (2015) 144–157

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Model-based experimental screening for DOC parameter estimation Björn Lundberg a,∗ , Jonas Sjöblom b , Åsa Johansson c , Björn Westerberg d , Derek Creaser a a

Department of Chemical and Biological Engineering, Chemical Reaction Engineering, Chalmers University of Technology, Göteborg SE-412 96, Sweden Department of Applied Mechanics, Division of Combustion, Chalmers University of Technology, Göteborg SE-412 96, Sweden Johnson Matthey AB, Västra Frölunda SE-421 31, Sweden d Scania AB, Södertälje SE-151 87, Sweden b c

a r t i c l e

a b s t r a c t

i n f o

Article history: Received 4 August 2014 Received in revised form 4 December 2014 Accepted 7 January 2015 Available online 19 January 2015 Keywords: Parameter estimation D-optimal design Diesel Oxidation Catalyst Multivariate Data Analysis Engine rig experiments

In the current study a parameter estimation method based on data screening by sensitivity analysis is presented. The method applied Multivariate Data Analysis (MVDA) on a large transient data set to select different subsets on which parameters estimation was performed. The subset was continuously updated as the parameter values developed using Principal Component Analysis (PCA) and D-optimal onion design. The measurement data was taken from a Diesel Oxidation Catalyst (DOC) connected to a full scale engine rig and both kinetic and mass transport parameters were estimated. The methodology was compared to a conventional parameter estimation method and it was concluded that the proposed method achieved a 32% lower residual sum of squares but also that it displayed less tendencies to converge to a local minima. The computational time was however significantly longer for the evaluated method. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Common kinetic mechanisms involve a large number of parameters. These mechanisms and the complex interaction between different reactions often result in a high correlation between the kinetic parameters (Bernaerts et al., 2000; Franceschini and Macchietto, 2008b). This severely complicates the parameter estimation process and increases the importance of Design of Experiments (DoE) and the parameter estimation algorithm. DoE in reaction kinetics has a long tradition reaching back to the pioneering work by Box and Lucas (1959), Box and Hunter (1965), Draper and Hunter (1966, 1967a,b), Box (1968), Hill et al. (1968), and Hunter et al. (1969). The focus has however been generally more directed toward steady-state systems rather than transient ones (Franceschini and Macchietto, 2008a) even though the need and benefits of transient experiments have been demonstrated (Goodwin, 1977; Berger et al., 2008). The system studied in the current paper is a reactor model of a full scale Diesel Oxidation Catalyst (DOC) connected to a Heavy Duty Diesel engine rig. The system is transient and the model parameters highly correlated; DoE is therefore of highest importance. A wide range of catalyst models with varying level of detail can be found in literature. Generally, the kinetics are either described

∗ Corresponding author. Tel.: +46 317722942. E-mail address: [email protected] (B. Lundberg). http://dx.doi.org/10.1016/j.compchemeng.2015.01.004 0098-1354/© 2015 Elsevier Ltd. All rights reserved.

by a microkinetic model (Dumesic et al., 1993; Olsson et al., 2002; Crocoll et al., 2005; Salomons et al., 2006), where every reaction is described by elementary steps, or by more empirical global kinetic models (Voltz et al., 1973; Ansell et al., 1996; Lafossas et al., 2011; Watling et al., 2012). However, the catalyst model is not only defined by the kinetics but the influence of heat and mass transport must also be taken into account. The general trend for monolith converters is still to describe heat and mass transfer with simple one-dimensional models (Güthenke et al., 2007), especially for the purpose of parameter estimation, but different levels of complexity ˇ epánek et al., 2012) up to three-dimensional CFD-based models (Stˇ exist. The preference for simple, and thereby fast, catalyst models in parameter estimation tends to make the kinetic parameters case specific, resulting in kinetic models only giving reliable predictions for a catalyst similar to the one used in the experiments from which the model was derived (Wang et al., 2008). This means that efficient methods to tune model parameters from existing models to measurement data from new catalysts are highly desirable. To tune the catalyst model to measurement data a method of parameter optimization is needed. The methods can broadly be divided into either global or local optimization algorithms where the latter historically have been far more favored in automotive catalysis parameter estimation. The local optimization methods are generally based on an evaluation of the parameter sensitivity to identify the direction of steepest reduction of the residual. The advantage of these algorithms is that they are fast and accurate but the drawback is that they have a tendency to converge to a

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

local minimum close to the initial guess and thereby failing to converge to the global minima (Blockley, 2010). The global methods generally evaluate a larger parameter space which means that they are more likely to find the global optima. The main drawback with the global methods is very high computational cost which means that they often are not realistic to use for complex systems with many parameters. To be successful the global algorithms also are more dependent on robust models since a wide parameter span will be needed to be evaluated. Examples from literature where global optimization algorithms have been applied in catalysis modeling include Glielmo and Santini (1999) who used the Genetic Algorithm (an evolutionary algorithm) and Ramanathan and Sharma (2011) who showed that an improved fit could be achieved if a global algorithm was followed by a local algorithm to increase the accuracy. A method of reducing the tendencies of the local methods to converge to a local minimum close to the initial guess may therefore be the more desirable way of improving parameter optimization for the current application. In our previous study (Lundberg et al., 2014), an experimental plan that included several different DOC catalyst configurations, was applied to a full scale engine rig system and used for parameter tuning. The best model fit was achieved with a model that included internal transport resistance with the effective diffusivities of the species being tuned in addition to the kinetic parameters. In the current study, the data and catalyst model formulation will be re-used but with a new approach to parameter tuning will be presented and evaluated. The method is based on a study by Sjoblom and Creaser (2008) who developed methods for model-based DoE including transient data. The methodology presented by Sjöblom and Creaser includes Multivariate Data Analysis (MVDA) to reduce the dimensionality of the largely correlated Jacobian and subsequent selection of a D-optimal sub-set of parameters to estimate based on the parameter sensitivity of simulation data. The relatively simple global kinetic model in the current study makes the number of parameters to estimate less of a concern. Instead the method was revised to select a set of time points from the large transient data set with high sensitivity for parameter changes and at the same time retaining a good representation of the entire data set. To achieve these properties the D-optimal onion design, introduced by Olsson et al. (2004), was selected as an algorithm for data point selection. The algorithms combine the D-optimal design, which by itself tends to only select the most extreme points (Pinto et al., 1990; Zullo, 1991), with a space filling approach making the selection more evenly distributed over the entire data set. Parameter estimation was performed based only in the selected time points with the aim of making use of improved statistical properties to improve the overall fit to measurement data. Since the parameter sensitivity, and thereby the data selection, will depend on the parameter values, the parameter estimation was performed in an iterative process where the data point selection was updated as new parameter values were found. This change of data points used for parameter tuning may also help to avoid local optima. As an evaluation, the efficiency and resulting fit to measurement data achieved by the presented method will be compared to the results of a traditional method of parameter estimation.

2. Experimental The most important characteristics of the reactor model and adjustable parameters are summarized in the current section and for completeness the defining equations are given in Appendix. A thorough description of the experimental set-up, measurements and engine operation design is given in the preceding study (Lundberg et al., 2014).

145

2.1. Reactor model A full scale catalyst, connected to an engine with varying inlet properties displays a highly dynamic behavior. This means that the catalyst outlet conditions will not only be influenced by the current inlet conditions but also those at previous time points. To describe this behavior a transient catalytic reactor model is needed. In this work a uniform radial flow and concentration distribution over the catalyst cross section was assumed which makes it sufficient to model only one channel. The single channel model, closely based on the model presented by Ericson et al. (2008), was discretized as tanks in series where the catalyst washcoat was discretized both radially and axially while the gas phase was only discretized axially. This 1D/2D (gas phase/washcoat) structure was chosen since it was considered a good compromise between accuracy and computational speed (Mladenov et al., 2010). A film theory model was used to model the external heat and mass transport between gas bulk and washcoat surface. Axial diffusion and radial temperature gradients in the washcoat were neglected. The global kinetics of the DOC has been described in literature with several different levels of detail. The simplest kinetic models only describe the oxidation of CO, HC and NO (Voltz et al., 1973; Oh and Cavendish, 1982; Wang et al., 2008) (where HC is represented as one molecular species) but more detailed models with a more complex description of HC (Stamatelos et al., 1999; Lafossas et al., 2011) and other additional reactions (Ansell et al., 1996; Salomons et al., 2006; Pandya et al., 2009) have also been developed. The kinetic model used in this study is of Langmuir–Hinshelwood type and was originally suggested in the classical work by Voltz et al. (1973) and later modified by Oh and Cavendish (1982). The model, which has been widely and frequently used in DOC modeling over the years, only includes three reactions of which one is equilibriumlimited: CO + 0.5O2 → CO2 C3 H6 + 4.5O2 → 3CO2 + 3H2 O NO + 0.5O2  2NO2 The reaction rates were calculated according to Eqs. (1)–(5). r1 = r2 = r3 = K =

k1 yCO yO2

(1)

G(yi , Ts ) k2 yC3 H6 yO2 G(yi , Ts ) k3 yNO yO2 G(yi , Ts )

(2)

 1−

K Kp

 (3)

yNO2

(4)

1/2

yNO yO

2

2 G(yi , Ts ) = Ts (1 + K4 yCO + K5 yC3 H6 )2 (1 + K6 yCO yC2

3 H6

0.7 )(1 + K7 yNO )

(5) where Kj is the reaction rate coefficient for the inhibition terms in the denominator G and Kp is the equilibrium constant for NO oxidation. At thermodynamic equilibrium, Kp will be equal to K and reaction rate r3 will be equal to zero. Both the rate coefficients kj and the inhibition terms Kj were described by Arrhenius expressions: kj = Aj e−EA,j /RTs

(6)

The start values for estimation of kinetic parameters were taken from the start values in Wang et al. (2008) where results from several studies (Chen et al., 1988; Koltsakis et al., 1997; Dubien

146

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

et al., 1998; Koltsakis and Stamatelos, 1999; Kandylas and Koltsakis, 2002; Tsinoglou and Koltsakis, 2003; Guojiang and Song, 2005) were compiled. The initial values for kinetic parameter estimation used in this study are shown in Section 3. The kinetic parameters in Eqs. (1)–(5) are highly correlated and since the initial parameter values were taken from different studies, the fit of the model to experimental data was expected to be poor before any parameter tuning was performed. However, the parameters were successfully used as a starting point for parameter tuning of a DOC model against engine rig data in Wang et al. (2008) which was also the intended application in the present study. 2.2. Adjustable parameters To tune the model to the measurement data a number of parameters in scaled and centered form (Sjoblom and Creaser, 2008) were adjusted. These parameters can be divided into kinetic parameters, mass transfer parameters and heat transfer parameters. In the current section only a brief introduction will be given to the adjustable parameters, a more detailed description can, however, be found in Section A.2 in Appendix. Three types of kinetic parameters are adjusted; pre-exponential factors, activation energies and activity scaling factors. The activity scaling factor (actq in Eq. (A2)) is simply a scale factor for all reaction rates on a certain site on a certain catalyst (parameters were tuned to several different catalysts simultaneously) with the purpose of accounting for different properties of catalyst samples, such as active noble-metal dispersion. Previously (Lundberg et al., 2014), the effective diffusivities were evaluated as a method of tuning the transport resistance with the aim of reducing the correlation between mass transport and kinetic parameters. The species were divided into two groups, where the first group contained O2 , NO, NO2 , and CO and the second group contained HC. In the first group all species are well defined with similar diffusivities (Massman, 1998) and could be expected to have similar mass transport properties in the washcoat with presumably the same bias from their true values. To reduce the number of parameters to tune the same scale factor, fDscl ,S in Eq. (A21), was used for all species in this group. The second group contained HC which was represented as C3 H6 but in reality it is a wide range of hydrocarbons with different mass transport properties. The scale factor used for HC was denoted fDscl ,L . The heat transfer parameters including thermal mass (sheat ), environmental temperature (T∞ ), and a lumped heat transfer coefficient (˛tot ) were earlier tuned to a satisfactory fit in (Lundberg et al., 2014) for the system at hand and will therefore not be retuned in the current study. 2.3. Experimental set-up A 13 L heavy duty diesel engine was used as the exhaust source and Swedish MK1 diesel, a commercial low-sulfur (approximately 5 ppm S) diesel, was used as fuel. The engine was equipped with a dynamometer control system enabling independent control of load and speed. A change in operating point for the engine could be quickly implemented, although the dynamics of the engine rig (as well as the monolith) itself were much slower due to the large thermal mass. To further widen the experimental range, four different catalyst configurations with different noble metal loading, lengths, and washcoat thicknesses were used. All catalysts, shown in Table 1, were model catalysts provided by Johnson Matthey. All configurations have a 30.5 cm diameter (12 inch), 62 channels/cm2 (400 cpsi) cell density and a wall thickness of 0.152 mm. The substrate material used was cordierite. The

Fig. 1. Temperature engine map and selected operating points.

catalysts were thermally conditioned for one hour at 600 ◦ C prior to use. 2.4. Engine operation design As mentioned in Section 1 the available exhaust composition, flow and temperature is limited by the operating points of the engine. It also takes several minutes for the catalyst inlet conditions to stabilize when switching between operating points, which means that experimental time will be important when deciding the number of different operating points when the full transient behavior is of interest. Since the number of catalyst configurations was large (Table 1) and some replicates also were necessary only 8 different operating points were selected. The operating points were selected to make as large steps as possible in the different variables including concentrations of NO, NO2 , HC, CO and O2 as well as temperature and flow rate with the purpose of making the experimental space as large (and orthogonal) as possible. Fig. 1 shows how the eight operating points span over the engine map for the case of temperature. Every engine operating point was run for 15 min before switching to another, thereby generating both transient data as well as data for operation close to steady-state. Further details about the engine operation design were reported earlier (Lundberg et al., 2014). 3. Parameter estimation method With the engine operating points selected above, both transient and stationary data are available over a wide range of temperatures, flow rates and concentrations which will aid the parameter tuning. There is however a risk that some parts of the data set will dominate the parameter tuning simply because there still is a large number of data points at similar conditions. This similarity in the data may result in the parameter tuning being less efficient since some parameters may have a very low sensitivity for a large part of the data points. Data reduction by Multivariate Data Analysis (MVDA) is a means by which to make the influence from data points at different conditions more equal. In the current study a method of data reduction is evaluated that can be divided into three steps. In the first step the sensitivity for all parameters, time points and residuals is calculated by finite differences. The resulting matrix is commonly denoted the Jacobian matrix. In the second step a principal component analysis (PCA) is performed on the Jacobian matrix (J) resulting in a linear model consisting of a loadings matrix (P) and a scores matrix (T) according to Eq. (7). J ≈T ×P

(7)

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

147

Table 1 Catalyst configuration used for parameter estimation (a–d). Catalyst configuration c consists of two catalysts in series. The Pt-loading is shown in both mass% of washcoat and g/ft3 total catalyst volume. Configuration

Pt-loading [mass% washcoat]

Pt-loading [g/ft3 monolith]

a b c d

0.30 0.59 0.30 0.59

15 30 15 15

The scores matrix gives a simplified description of the parameter sensitivity of all data points given in the Jacobian. The scores matrix has fewer columns than the Jacobian but the same number of rows, see Appendix for details. In the third step a D-optimal onion design is performed on the resulting scores matrix giving a selection of data points that will have a more equal influence from the different parameters. The reduced number of columns in the scores matrix compared to the Jacobian will reduce the computational time of finding a D-optimal design which would be very high if the Jacobian would be used directly. The method is described in detail in Appendix and is schematically summarized in Fig. 2.

Length [cm]

Washcoat thickness [mm]

10.2 10.2 2 × 10.2 10.2

0.110 0.110 0.110 0.055

account for differences in the average concentration the residuals were also weighted against the inverse of the average outlet mole fractions for the respective component. In the current study the data points used for parameter estimation will change and a residual calculation that is more adaptive to a change data is needed. To achieve a more balanced influence from the different components, the residual was calculated as the difference in conversion instead of the difference in mole fraction, as demonstrated in Eq. (8).



yi,in −  yi,in − yi,out yi,out − yi,in yi,in



yi,out − yi,out

3.1. Residual calculation

res =

The definition of the residuals and their weighting are of utmost importance to obtain a balanced description of the parameter sensitivity at every time point of the data set and thereby also for the effectiveness of the parameter search algorithm itself. In the previous study and parameter estimation on the presented data (Lundberg et al., 2014), the difference in outlet mole fraction between measured and simulated data was used as the residual. To

yi,out is the where wi is a weighting factor (further described below),  simulated outlet mole fraction, yi,out is the measured outlet mole fraction and yi,in is the measured inlet mole fraction. To compensate for the fact that the NO conversion is limited by thermodynamic constraints the residual for NO was instead calculated relative to the

×w =−

yi,in

× wi (8)

Original parameters

Parameter sensivity calculaons (at all me points)

Jacobian matrix

PCA

Residual sum of squares calculaon (at all me points)

Scores matrix

D-opmal onion design Selected me points

Parameter tuning (at selected me points)

No

Improved parameter values?

Yes

Fig. 2. Summary of the method used for parameter estimation. Improved parameter values means that parameters were found that reduced the residual sum of squares at the selected time points which will result in a right hand loop.

148

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

equilibrium conversion possible at the current condition as shown in Eq. (9).



res =

(yNO,in − yNO,out,lim ) − ( (yNO,in − yNO,out,lim ) − (yNO,out − yNO,out,lim ) yNO,out − yNO,out,lim ) − yNO,in − yNO,out,lim yNO,in − yNO,out,lim

where yout,lim is the lowest molar fraction of NO possible via NO oxidation at the inlet conditions due to thermodynamic constraints. Since the model did not account for NOx reduction, and since the NOx reduction was below 3% for more than 99% of the data set, the fit for NOx is considered to be described well enough by only the NO residual. In addition to NO, the residuals of both HC and CO were used for parameter estimation. It is not desirable to fit parameters to data points with close to 100% conversion, since the kinetic parameter sensitivity for these points will be very low. Therefore it was necessary to set the residuals to zero for certain components at such points in the data set. For the current data set the following limits for outlet mole fractions, based on measurement accuracy, were selected to generate a zero residual for the three different components: yHC,out < 2 ppm yCO,out < 5 ppm yNO,out − yNO,out,lim < 5 ppm When starting simulations with the catalyst model the only known conditions were those of the inlet and outlet flow. This means that starting properties such as concentrations and temperature inside the catalyst needed to be estimated which resulted in the first seconds of every simulation being unreliable. To avoid parameters being tuned against unreliable simulation results the residuals from the first 120 s of every experimental data set was set to zero. The weighting factor wi is calculated according to Eq. (10): 1 wi = √ ni

(10)

where ni is the number of data points used for residual calculation for species i after points with close to 100% conversion and points from the first 120 s of the experimental data sets have been removed. 3.2. Parameter tuning After data selection using D-optimal onion design a parameter estimation with the gradient search method is undertaken for only the selected time points. The gradient search method is very efficient for linear systems but can also be applied for nonlinear systems such as catalyst models. For a nonlinear system the residual function is first linearized for all parameters and then a step in the parameter space is made in the direction of the steepest descent (lowest residual). This process is usually repeated until the change in residual is below a certain tolerance. The gradient search method of choice in this work is the trust-region-reflective method (Coleman and Li, 1996) which is the standard method in Matlab, the software used for parameter estimation in this project, for over-determined nonlinear least square problems (function call lsqnonlin). In theory the proposed methodology with PCA and D-optimal design should make a gradient search more efficient compared to if all data points were used, since the correct direction for parameter adjustment should be better defined by a well-balanced subset of more sensitive data points. Since the properties of the sensitivity matrix (Jacobian in all data points) are dependent on the actual parameter values, the best selection of time points for parameter tuning will change as the parameters are tuned. To make the parameter tuning efficient it

 × wNO = −

yNO,out − yNO,out wNO yNO,in − yNO,out,lim

(9)

was therefore necessary to select new time points for parameter tuning as new parameter values developed. In the current work, the parameter tuning was limited to only a few steps in the parameter space which decreased the risk of obtaining a too large influence from the selected time points on the overall fit. To investigate the influence of the number of steps on the final fit, two parameter tuning methods were evaluated. The first method (referred to as 1 step) selected new time points after every step in parameter space while the second method (referred to as 5 steps) took five steps in parameter space before new time points were selected. The 1 step and 5 steps methods will together be referred to as the MVDA methods. The motivation for the “1 step” approach was that the points selected with D-optimal onion design were selected by evaluating the sensitivity for all data points (Jacobian matrix) for a defined set of parameter values. When the parameter values are changed, as during a gradient search, the sensitivity also changes and the point selection should be updated to ensure that the parameter estimation is performed based on the statistically best data points. By limiting the number of steps to one, new points will be selected every time a new set of parameter values is found. The drawback with the method is that the sensitivity for all data points needed to be frequently calculated. This is computationally expensive since all time points are used for these calculations and more computational time was spent selecting new data points than performing actual parameter estimation with the selected data points. In the second approach (“5 steps” case), the number of steps taken with the gradient search method (with the points selected with D-optimal onion design) was instead limited to 5. Since the steps with the gradient search method are still limited, the selected time points may still be relevant, despite the changes in parameter values. As a result more computational time is be spent on parameter tuning and less on selecting data points compared to the 1 step method. As a reference, a conventional parameter tuning with a gradient search method was also performed based on the entire data set using the same data points as for the sensitivity analysis described in Section A.3 in Appendix. The parameter tuning in the preceding study was performed with this method and data points but the changed formulation of the residual calculation between the two studies makes a new tuning necessary for comparison with the MVDA-cases. 3.3. Stop criteria Two stop criteria were applied in the study. The first stop criterion was applied to the reduction of the residual sum of squares for the full data set. 1. The parameter estimation was stopped if a new set of parameters that reduced the residual sum of squares for the full data set by more than 1 percent could not be found after six consecutive steps in parameter space (six right hand loops in Fig. 2). The second criterion was applied to the parameter estimation performed with the data points selected with D-optimal design. If no reduction of the residual sum of squares was achieved (referred

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

to as improved parameter values in Fig. 2) with the gradient search method in the data points selected with D-optimal design, a new Doptimal design will be performed. Even though the scores matrix is unchanged (since the parameters were unchanged) the D-optimal design is likely to select different data points since the selection of a design matrix with D-optimal design uses a random seed in the selection process. A new parameter estimation with the gradient search method was then performed on the new data set. This process was repeated until parameter values were found that reduce the residual sum of squares of the points selected with D-optimal design or until the second stop criteria was fulfilled. 2. The parameter estimation was stopped if no improved parameter values were found after more than six selections of data points for the same parameter values (corresponding to more than six left hand loops in Fig. 2) If none of the criteria were fulfilled a new sensitivity calculation on the whole data set (Jacobian matrix calculation) was performed and the cycle returned to “Parameter Sensitivity Calculations” in Fig. 2. 4. Results and discussion 4.1. Data point selection with D-optimal onion design The data point selection with D-optimal onion design was updated as the parameters were estimated according to the method described previously in Section 3. Fig. 3 shows one example of a selection from which it is clear that certain areas were frequently selected from the onion layers with high sensitivity, whereas other areas were mainly selected from the onion layers with low sensitivity. The points selected were similarly for the different catalyst configurations even though points from catalyst configuration c (data point 28,818–43,240) were less frequently selected for the depicted case. The points selected were taken from both transient and stationary sections of the data. As expected the data points selected in layers with high sensitivity were mainly taken from transient areas while data points selected in layers with low sensitivity were mainly taken from more stationary areas. 4.2. Evolution of residuals In Fig. 4 the evolution of the residual sum of squares of the gas phase components (NO, HC, and CO) for the whole data set is shown. The figure also shows the final residual sum of squares for the reference method as a benchmark for the two MVDA methods. The residual sum of squares for the whole data set is only calculated after a completed parameter estimation in the time points selected by D-optimal onion design. Neither of the methods (1 step or 5 steps) fulfilled either of the two stop criteria in Section 3.3 and were instead stopped after about 3 × 104 core hours of simulation. For the reference case, the first stop criterion was fulfilled after about 5 × 103 core hours of simulation with a residual sum of squares of 3580. The reference parameter tuning method was then approaching parameter values where the catalyst model was unable to find a stable solution and thereby converging to a local minimum. After 5 × 103 core hours of simulation the 5 steps method had a residual sum of squares of 3256 but the 1 step method on the other hand had a residual sum of squares of 8919, it is therefore not possible to conclude that the MVDA-methods enabled a faster residual reduction. Instead the advantage with the methods appears to be a reduced risk of converging to a local minimum. Even after the comprehensive parameter tuning already performed,

149

Table 2 Residual sum square of each component and complete data set together with the summation of residual sum square (rightmost column) for the different methods.

Start 1 step 5 steps Ref

NO

HC

CO

Sum

18,986 515 558 1238

9736 1161 1193 1173

16,954 730 742 1169

45,675 2405 2494 3580

the residuals in Fig. 4 still show a declining trend even when both methods have reached residual sum of squares below 2500. Although the rate of decline appears to be decreasing, no clear optima appear to have been reached. 4.3. Comparison to measurement data The parameter values that gave the lowest residual sum of squares, shown in Fig. 4, shall be regarded as the final tuning results of the methods. A compilation of the residual sum of sqares for the final results and the starting paramter values are shown for every gas phase component in Table 2. From these results it appears that HC was the most difficult component to fit to measurement data for the MVDA methods which may have several different reasons. Firstly HC and CO concentrations are highly correlated in the exhaust of a HDD engine which means that it can be difficult to generate experimental data, let alone subsets of data, where the parameters that influence each component can be distinguished. Since both components have high concentrations at similar time points, the larger residual for CO at the starting parameters may have influenced the parameter estimation to initially focus more on reducing the CO residual at the cost of decreased sensitivity for HC concentration. HC also differs from NO and CO since it is modeled as an average of a wide range of hydrocarbons. The physical properties are more difficult to model and the HC composition at the catalyst inlet will even change between different engine operating points. For the reference case, the HC residual was at about the same level as for the MVDA methods but both NO and CO residuals are noticeably higher. The results in Table 2 also indicate that both MVDA methods achivieved a very similar final fit, with HC being the largest residual, but that the reference method differs because the total residual is more equally distibuted over the three components. Two examples of the final results of the tuning of the kinetic and mass transport parameters for the different parameter estimation methods are shown in Figs. 5 and 6. Fig. 5 shows a very similar fit for the two MVDA methods for all components. The reference method shows a similar fit for CO and HC but for NO the transient was poorly captured. The NO production before the change in engine operating point is likely a result of NO2 reduction by CO or HC. These reactions are not included in the kinetic model and can therefore not be captured by any of the models. Fig. 6 also shows a very similar fit for the two MVDA methods, with the most noticeable difference compared to the reference method, again mainly being the better fit for NO. However the NO concentrations after the transient did not agree as well with measurement data as for the case in Fig. 5. The fit of HC and CO is somewhat better for the MVDA methods, but there were still notable differences compared to measurement data. The above figures show only two of 16 transients used for parameter estimation (4 catalyst configurations and 4 operating point changes for every configuration). The transients were, as previously mentioned, selected to give a wide temperature, flow and concentration window. The combinations of torque and speed that

150

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

Fig. 3. Measured outlet mole fraction of NO for an example of a selected data set colored by relative sensitivity. The order of the engine operating points is 1, 7, 2, 8, 3, 4, 5, 6 which were run sequentially for catalyst configurations a, b, c and d (separated by blue vertical lines). The relative sensitivity is an indication of how sensitive the data point is to changes in parameter value relative to the rest of the data set. A value of 0.9 for example, means that the data point has a higher sensitivity than 90% of the data set. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Residual sum of squares for every set of parameters for the two evaluated methods. The dashed line is the residual sum of squares for the final set of parameters from the reference method. (a) 1 step method (b) 5 steps method (c) 1 step method, magnified (d) 5 steps method, magnified.

Table 3 Residual sum of squares for the experiment displayed in Figs. 5 and 6 for every component. NO

HC

CO

Figure

1 step 5 steps Ref

40.9 42.3 85.5

40.4 38.7 43.9

94.9 79.4 68.2

5 5 5

1 step 5 steps Ref

66.2 76.6 212.8

112.6 111.9 112.2

80.2 101.9 157.2

6 6 6

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

b)

a) 1600

1200 [ppm]

1000

60 NO in NO out measured NO out sim, 1 Step NO out sim, 5 Steps NO out sim, Ref

C H in 3 6 50

[ppm]

1400

151

800 600

C3H6 out measured

40

C3H6 out sim, 1 Step

30

C3H6 out sim, Ref

C H out sim, 5 Steps 3 6

20

400 10

200 0 0

200

400

600

0 0

800 1000 1200 1400 1600 1800 Time [s]

200

400

600

c)

d)

250

500 CO in CO out measured CO out sim, 1 Step CO out sim, 5 Steps CO out sim, Ref

200 150

450 400 350 [°C]

[ppm]

800 1000 1200 1400 1600 1800 Time [s]

100

300 250 Tin

200

50

T

150 0 0

200

400

600

100 0

800 1000 1200 1400 1600 1800 Time [s]

out

measured

Tout simulated 200

400

600

800 1000 1200 1400 1600 1800 Time [s]

Fig. 5. Measured and simulated outlet concentrations of NO (a), C3 H6 (b), CO (c) and temperature (d) for a change in operating point from 1 to 7 for catalyst configuration c.

result in high HC and CO are few and therefore only two of the eight engine operating points used in the study generate significant HC and CO concentrations at the catalyst outlet. To demonstrate the performance of the different methods, transients with engine operating points with high HC and CO have been chosen in both Figs. 5 and 6. In addition, the high temperature at the start of the experiments provide an interesting transient change in both CO and HC conversion for the first 900 s and the change in engine operating point provides transients in NO and temperature. To put the fit of

the depicted experiments in Figs. 5 and 6 in relation to the fit of the remaining 14 experiments the residual sum of squares for the presented experiments are shown in Table 3. The trends in Table 3 for two of the experiments agree with the results for all experiments shown in Table 2. The fit for NO is better for the MVDA methods in both experiments while the fit for HC is about the same for all methods. Even though the fit for CO is better for the reference case in one of the experiments, the fit of CO over both experiments in the Table 3 is still better for the MVDA methods.

b)

a) 1000 800

[ppm]

700 600

NO in NO out measured NO out sim, 1 Step NO out sim, 5 Steps NO out sim, Ref

C3H6 out measured C3H6 out sim, 1 Step

50

500 400 300

C H out sim, 5 Steps

40

3 6

C3H6 out sim, Ref

30 20

200

10

100 0 0

C H in 3 6

60

[ppm]

900

70

200

400

600

0 0

800 1000 1200 1400 1600 1800 Time [s]

200

400

600

c)

d)

250

450 CO in CO out measured CO out sim, 1 Step CO out sim, 5 Steps CO out sim, Ref

150

400 350 [°C]

[ppm]

200

100

300 250 Tin

200 50 0 0

800 1000 1200 1400 1600 1800 Time [s]

T

150 200

400

600

800 1000 1200 1400 1600 1800 Time [s]

100 0

out

measured

Tout simulated 200

400

600

800 1000 1200 1400 1600 1800 Time [s]

Fig. 6. Measured and simulated outlet concentrations of NO (a), C3 H6 (b), CO (c) and temperature (d) for a change in operating point from 2 to 8 for catalyst configuration b.

1.00 12.3 10.9 4.85 7.71 1.00 0.706 0.706 0.793 1.00 1.00 0.740 0.730 0.828 0.701

actd a [%] actb a [%]

3.98 50.6 3.96 1.35 1.93 2.08 × 10 3.88 × 105 2.59 × 105 5.33 × 104 8.25 × 104

1.00 1.00 1.00 1.00 0.766 3.10 × 104 2.28 × 104 2.68 × 104 2.87 × 104 3.05 × 104 −9.65 × 104 −9.67 × 104 −9.81 × 104 −9.74 × 104 −9.69 × 104

acta , actc a [%] EA ,7 [J/mol] EA,6 [J/mol]

65.5 39.0 0.151 77.7 95.8 4.50 × 10 3.18 × 1013 1.77 × 1014 2.34 × 1017 2.71 × 1014 4.00 × 10 4.70 × 1017 1.06 × 1018 1.67 × 1022 1.25 × 1019 1.00 × 10 2.17 × 1015 2.59 × 1016 2.78 × 1020 4.89 × 1016

−3.00 × 103 4.54 × 103 2.21 × 103 −3.30 × 103 −1.64 × 103 −7.99 × 103 −5.54 × 103 −1.02 × 104 −8.50 × 103 −8.06 × 103 7.00 × 104 7.82 × 104 8.18 × 104 7.29 × 104 7.46 × 104 Start 1 step 5 steps Ref Preceding study

EA,5 [J/mol] EA,4 [J/mol] EA,3 [J/mol]

Start (Wang et al., 2008) 1 step 5 steps Ref Preceding study (Lundberg et al., 2014)

a acta is activity scaling factor for catalyst configuration a (0.30 wt% Pt), actb is activity scaling factor for catalyst configuration b (0.59 wt% Pt), and actd is activity scaling factor for catalyst configuration d (0.59 wt% Pt with thin washcoat). b Effective diffusivity scaling for small (S) and large (L) molecules.

[–] 1.00 1.49 1.34 1.62 2.11

b

fDscl ,L [–] b

fDscl ,S

8.00 × 10 7.33 × 104 7.85 × 104 7.90 × 104 7.67 × 104 4.79 × 10 4.85 × 105 2.99 × 106 4.84 × 1010 4.71 × 107

EA,2 [J/mol] EA,1 [J/mol]

5

A7 [–] A6 [–] A4 [–]

A5 [–]

3

14

17

20

A2 [mol, K/m2 , s]

A3 [mol, K/m2 , s] Table 4 Results of kinetic parameter tuning.

Table 4 shows the final parameter values for both approaches and the reference case together with the starting parameters values from the study by Wang et al. (2008) and the parameters of the best fit from the preceding study by the authors (Lundberg et al., 2014). Since catalyst configuration c (see Table 1) consisted of the catalyst in configuration a in series with another identical catalyst, the same activity scaling factor was assumed for both configurations a and c. The correlation between activity scaling factors and preexponential factors in the numerator (A1–3 ) were linear. An increase by a certain factor of the activity scaling factor gives the same increase in reaction rate as an increase in the pre-exponential factor by the same factor. This was true since the kinetic model only uses one type of reaction site. To make the kinetic parameters more comparable, the activity scaling factors for every set of parameters (1 step, 5 steps, Ref and preceding study) were scaled with a number making the highest activity scaling factor of that set of parameters equal to 1% at the same time as the pre-exponential factors in the numerator were divided by the same number. This did not change the behavior of the kinetic model; but merely shifted the factor from one parameter to another. The activity scaling factor was in other words only a handle to tune differences between the catalysts and should only be interpreted as a relative difference and not an absolute value for activity scaling. In general it can be expected that a high noble metal loading should correspond to a low dispersion and thereby also a low activity scaling factor. It is therefore not unexpected that the lowest loading, i.e. 0.30 wt% Pt loading (acta and actc ), also has highest activity scaling for all cases of the current study. The complexity of the kinetic model and the catalyst model complicates a thorough analysis of every single parameter value but some general remarks can be made. The original parameters gave a too high conversion for all components and a very poor fit in general. The tuned parameters have therefore been changed in a direction where the reaction rates in general are slower. This is not always obvious when examining the parameter values since a large pre-exponential factor in the numerator (A1–3 ) may be compensated by a large pre-exponential factor in the denominator (A4–7 ) and differences in activation energies (EA,j ) will make different preexponential factors significant at different temperatures. As indicated by the residual analysis in the previous, section the parameter values are very similar for the two MVDA methods, while the reference case stands out by having large pre-exponentials in the numerator (A1–3 ). The reference case however also has a value for A7 that is about five orders of magnitude larger than for the other models. The fact that A7 is in the denominator expression for NO inhibition means that it will always have an effect since the NO levels are always above zero due to thermodynamic limitations. The difference between the reference case and the MVDA methods is thus a stronger dependence on NO concentration for the reference case more than an overall large difference in reaction rates due to large pre-exponentials in the numerator (A1–3 ) for the reference case. In the preceding study it was concluded that the sign of EA,7 made the NO inhibition term mimic internal transport resistance since an increase in temperature increases the term and thereby affects all reaction rates negatively. The sign of EA,7 is still positive for all evaluated methods but for both MVDA methods EA,5 has changed from a negative starting value to a positive final parameter value. This means that the MVDA model has two terms mimicking internal transport resistance, one as a function of NO concentration and the other a function of HC concentration, compared to one for the reference case. These values together with the decreased mass transport resistance due to increased diffusivity scaling parameters mean that the actual transport resistance in the model is reduced and replaced by a fictitious one connected to

4

4.4. Analysis of final parameter values

10.0 × 104 9.51 × 104 9.52 × 104 9.86 × 104 9.80 × 104

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

A1 [mol, K/m2 , s]

152

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

the kinetics. The result of the parameter tuning with MVDA shows the importance of having a well formulated model and also displays how an efficient parameter tuning method can tune the parameters to less realistic values, and still achieve a good fit, even if the model formulation is deficient. In general the parameter values from the preceding study are closer to the values of the MVDA methods than the reference case. The parameters from the preceding study also gave an overall better fit with a residual sum of squares of 3094 compared to 3580 for the reference case. There are two main differences between the parameter estimation method of the preceding study and the reference case. Firstly the residuals were calculated differently, see Section 3.1, and secondly a ten times higher parameter weighting for kinetic parameters (w in Eqs. (A15) and (A17)) was used in the current study to make the Jacobian less sensitive to numerical noise. It appears that this increased parameter weighting has resulted in estimated parameter values further away from the original ones with an overall poorer fit. This could be an effect of some parameters, such as A1 –A3 , having a large influence on the fit at an early stage of the parameter estimation which makes them change fast and thereby reach values where the sensitivity for the other parameters are low. In other words the risk of finding a local minimum may be increased by the increased parameter weighting. One of the aims of the MVDA method was to make the influence of parameters more equal and it appears that this has prevented the parameters from reaching extreme values for both the 1 step and the 5 steps method even though the parameter weighting was high. 4.5. Final remarks By continuously changing data points on which the parameter estimation is performed the risk of finding a local minimum appears to have been reduced at the cost of longer computational time. As mentioned in the introduction, these characteristics are typically associated with global parameter estimation methods such as simulated annealing. The difference with the method at hand is, however, that since a gradient search method is applied, the parameters will evolve without evaluating a large number of parameter sets (as is characteristic of a global search method) which will put less demand on the model stability. This will also reduce computational cost which already is high for such a complex system of equations and large data set. The method of parameter estimation should be useful for complex dynamic models where model stability limits what parameter values are possible and where the number of data points is large. The model is not likely to be suitable if the diversity of the data points is low, for example if the model should be tuned to only a few stationary data points. With low diversity in the data it is not likely to achieve a notable improvement in parameter sensitivity by changing the data point selection. The computational effort is then better spent on tuning the model to the few available data points directly. For models where the stability is highly dependent on starting conditions, such as detailed kinetic models, the large number of short simulation sections where starting conditions are needed may make the method less suitable. For example it may be very difficult to find stable starting conditions if the first data point is in the middle of a sharp transient such as a catalyst ignition point. 5. Conclusions

153

method, referred to as the MVDA method, was applied to full scale catalytic exhaust gas data to estimate kinetic parameters. The MVDA method displays a better overall fit for all components resulting in a residual sum of squares 32% below what was achieved with a conventional method used as reference. The method displayed less tendencies to converge to local minima and to reach areas where the catalyst model is unstable; however, the method was more computationally expensive than the reference method. In general the method should be applicable on other complex dynamic models where model stability limits what parameter values are possible and where the number of data points is large. Both evaluated configurations of the MVDA method gave similar final fit to measurement data. Neither method could be identified as less computationally expensive than the other. Acknowledgements The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at C3SE. Financial support from Swedish Energy Agency (Grant Nos. 322821 and 32282-2) is gratefully acknowledged. Appendix. A.1. Reactor model Quasi-steady state was considered to prevail for gas phase species (Eqs. (A1) and (A2)) and gas temperature (Eq. (A6)) since the characteristic time constants are far shorter for these processes compared to the solid heat balance (Eq. (A7)) and variations in inlet conditions. The gas bulk mass balance is given by:



0=

Vk,0

dci,k,0 dt



= Ftot (yi,k−1,0 − yi,k,0 ) − i,k,0 (ci,k,0 − ci,k,1 ) (A1)

and the washcoat mass balance is given by:



0=

dci,k,n Vk,n ε dt



= i,k,n−1 (ci,k,n−1 − ci,k,n )

−i,k,n (ci,k,n − ci,k,n+1 ) +



vi,j actq rj,k,n Am,k,n

(A2)

j

for n≥1 where rj,k,n is the reaction rate in mol/m2 Pt/s and Am,k,n is the catalyst active surface area in channel segment k and wall layer n, Vk,n is the volume of channel segment k and layer n and ε is the washcoat porosity. Index k indicates the segment number (axial discretization) where k = 1 represents the first segment in the direction of the flow and k = K represents the last. The index n indicates the layer number (radial discretization) where n = 1 represents the first washcoat layer closest to the gas bulk, n = N represents the last layer and n = 0 represents the gas bulk. The index i indicates the species and index q represent catalyst configuration. The mass transfer coefficients  i,k,n are given by: i,k,0 =

Ak (1/kc,i,k ) + (0.5x1 /Deff,i,k )

(A3)

i,k,n =

Deff,i,k Ak 0.5xn + 0.5xn+1

(A4)

for n = 1,. . ., N − 1. And for n = N it is: The method of parameter estimation evaluated in the current paper is based on data screening by sensitivity analysis where MVDA has been applied to a large transient data set to select a subset on which parameter estimation is performed. The subset was continuously updated as the parameter values developed. The

i,k,N = 0

(A5)

where Ak is the radial mass and heat transfer area in channel segment k. For simplicity, Ak was assumed to be constant for all wall layers. Deff,i,k is the effective pore diffusion coefficient of component

154

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

i in channel segment k, kc,i,k is the film transfer coefficient of component i in channel segment k and xn is the thickness of washcoat layer n. The gas energy balance is given by:



0=

Vk,0 cp,g

dTg,k



dt

Deff,i,k = = Ftot cp,g (Tg,k−1 − Tg,k ) − hk Ak (Tg,k − Ts,k ) (A6)

where cp,g is the heat capacity for the gas, hk is the heat transfer coefficient in channel segment k, Tg,k and Ts,k are the temperatures in channel segment k in the gas bulk and of the catalyst respectively. Note that the solid temperature, Ts,k , was not discretized radially since the high solid conductivity and short conduction distance will result in very small radial temperature gradients. The main temperature transport resistance is in other words assumed to be in the film. The solid energy balance is given by: dTs,k dt



cp,s ms,k +

sheat K



 n

−

actq rj,k,n Am,k,n (−Hj )

(A7)

where ms,k is the mass of catalyst material and substrate in channel segment k, cp,s is the mass weighted average of washcoat and substrate heat capacity, As is the channel cross sectional area of the substrate and washcoat, Hj is the heat of reaction j. sheat is a heat sink originating from solid material not included in the catalyst (e.g. catalyst canning and insulation). The heat loss to the environment is modeled by the lumped heat transfer parameter ˛tot and the effective environmental temperature T∞ . K is the total number of channel segments. The solid axial heat flux was calculated as: Ts,k − Ts,k−1 0.5zk + 0.5zk−1

(A8)

for k = 2,. . ., K − 1, K qk = 0

(A9)

hk =

ShDi,k d

Nug d

Di,k = Dref,i

Ts,k Tref

dp 3



8RTs,k Mi

(A14)

where dp is the mean pore diameter. A.2. Adjustable parameters

(A15)

where w is a weight factor for every parameter and p is the parameter that is changed by the gradient search method. The superscript “o” in the equations above indicates that this is the original values taken from literature or previous successful parameter tuning. The reaction rate coefficients were centered around a reference temperature to decrease the parameter correlation (Box et al., 1973). o kj,ref = Aoj e

−(E o /RTref )

(A16)

A,j

The reaction rate coefficient kj,ref was scaled according to: o ln(kj,ref ) = ln(kj,ref ) + wa,j pa,j

(A17)

The tuned rate constant at a reference temperature was then used to calculate the pre-exponential factor (A18)

When Eq. (A18) is inserted into the expression for the reaction rate coefficient in Eq. (6), the following expression for the tuned reaction rate coefficient is obtained kj = kj,ref e−[(EA,j /R)(1/TS −1/Tref )]

(A19)

The mass transport was tuned by adjusting the effective diffusivities for the species taking part in the reactions.

(A10)

fDscl ,i = fDo

(A11)

This will change the expression for the effective diffusivity according to

where Di,k is the diffusion coefficient for component i in channel segment k, g is the gas heat conductivity and d is the channel dimension. Sh and Nu are the Sherwood and the Nusselt numbers, respectively. Only the asymptotic values for Sh and Nu were used and thus entrance effects were neglected. Asymptotic values were taken from Tronconi and Forzatti (1992). The gas diffusivities for each component i at reference temperature (Tref ) were calculated using the Fuller–Schettlet–Giddins equation (Murzin and Salmi, 2005). Their dependence on temperature was expressed as:



DKi,k =

Aj = kj,ref eEA,j /RTref

for k = 1, k = K + 1where s is the heat conductivity for the solid material (same conductivity assumed for the washcoat and the cordierite support) and zk is the length of channel segment k (the total heat transfer distance is half the distance of each segment). The mass and heat transfer is described with the film model and the mass and heat transfer coefficients are given by: kc,i,k =

(A13)

where fD is a factor that takes into consideration the porosity and the tortuosity of the porous material. DKi,k is the Knudsen diffusivity which was calculated as:

o EA,j = EA,j + WEa,j pEa,j

j

˛tot (Ts,k − T∞ ) K

fD (1/Di,k ) + (1/DKi,k )

Both the pre-exponential factors, Aj , and the activation energies, EA,j , can be tuned. The scaling of EA,j was straightforward and performed according to:

= hk Ak (Tg,k − Ts,k ) − As (qk+1 − qk ) +

qk = −s

The effective diffusivity for pore diffusion was initially estimated based on an additive resistance, also known as the Bosanquet formula (Froment et al., 2011):

1.75 (A12)

scl ,i

Deff,i,k =

+ wDscl ,i pDscl ,i

fD fD ,i (1/Di,k ) + (1/DK,i,k ) scl

The value for fDo

scl ,i

(A20)

(A21)

is typically around 1.

The heat loss to the environment was modeled as a part of the total heat balance in Eq. (A7) and the parameters sheat , ˛tot and T∞ could be adjusted. The parameter scaling and tuning was performed in the same way as the other parameters: o T∞ = T∞ + wT∞ pT∞

(A22)

˛tot = ˛otot + W˛tot p˛tot

(A23)

o sheat = sheat + Wsheat psheat

(A24)

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

A.3. PCA and parameter sensitivity Principal Component Analysis, PCA, is a mathematical method where large sets of observations of possibly correlated variables are orthogonally transformed into linearly uncorrelated variables. The number of uncorrelated variables, referred to as principal components, is usually chosen to be of a smaller magnitude compared to the number of variables in the untreated data. The PCA can be performed with several different purposes such as identification of classes of data and outliers, simplification, data reduction, modeling, variable selection, and prediction (Martens and Naes, 1989). A simple way of describing the method would be as a simplification of a data matrix (X) where the more similarity within the objects results in fewer terms needed for a satisfactory representation of the data. The resulting linear model can be described by Eq. (A25): X = TP + E

(A25)

where T is the scores matrix, P is the loading matrix, and E is the error. If the dimension of X is N × K then T has the dimension N × A and P has the dimension A × K, where A is the number of components in the PCA. The scores matrix can be used to interpret the relation between the N observations (rows in X) and the loading matrix can be used to interpret the relation between the K variables (columns in X). In the present study the Jacobian (J), or the parameter sensitivity matrix, was used as the X-matrix. The Jacobian is defined as: J(f, ) =

∂f ∂

155

data set and could even result in a sub-optimization if used for parameter estimation. In a D-optimal onion design the data set is divided into layers where every layer includes points in a specified parameter sensitivity range. The so-called onion D-optimal design is then applied to every onion layer and a number of data points will thereby be selected in every parameter sensitivity range, resulting in a more balanced data point selection. The D-optimal onion design method, as applied in the current study, can be divided into a number of steps. 1. The scores matrix from the PCA analysis (T) were column-wise centered and scaled. The new matrix is called Tcs . 2. The elements of Tcs were squared, summed row-wise and thereafter sorted from high to low. The resulting vector is denoted Tcs2,sort and has the data points with the highest parameter sensitivity first and the data points with the lowest sensitivity last. 3. The rows in Tcs were sorted in the same order as Tcs2,sort and is denoted Tcs,sort 4. The number of onion layers (nlayers ) was calculated from the number of data points (m) and number of components in the PCA analysis (n): nlayers = 0.01m/(n + 1) − 1, where nlayers was rounded up to the closest integer. 5. The number of data points in every layer was calculated according to npoints (j) = 100 npoints (j) = 10(n + 1)kj−2

j=1 j = [2, 3, . . ., nlayers , nlayers + 1]

(A26)

where f is some function value which, in the case of the current study, was the residuals for HC, CO and NO (as described in Eqs. (8) and (9)) and  was the parameter vector in the nonlinear model. The number of columns in the J matrix will then correspond to the number of residual types (different components in this study) times the number of parameters to be tuned. In the current study the Jacobian was estimated by finite differences and the number of rows was equal to the number of time points (observations). Before the PCA was performed on J, all the rows with all components set to zero were removed. This was either a result of the time points being within the first 120 s of an experiment or all conversions being close to 100% according to the limits given in the previous section. In the PCA model of the Jacobian matrix the resulting scores matrix contained information about the relation between the time points and the loadings matrix contained information about the relation between the residual sensitivities. By analyzing the scores matrix it was possible to identify the time points where the influence from parameter changes on the different residuals were most noticeable, that is the parameter sensitivity was high. These time points were thereby good candidates to use if the number of time points used for parameter tuning were to be reduced. A.4. D-optimal onion design For linear models traditional experimental designs such as full factorial designs, fractional factorial designs, and response surface designs are suitable when the factors are relatively unconstrained. For nonlinear models less traditional experimental designs such as D-optimal design are more favorable, especially when the data is defined (in a so called Candidate set) and also when the number of experimental runs is restricted (Martens and Naes, 1989). The advantage with applying a D-optimal onion design, instead of a regular D-optimal design, is that the selection of points from the data set will be more balanced. A regular D-optimal design will result in a selection where only the most extreme points are chosen for parameter estimation. These points will have the highest parameter sensitivity but may not be representative for the whole



k was calculated to fulfill npoints where  (1:nlayers ) ≤ m ≤ npoints (1:nlayers + 1) 6. A vector denoted int was created to divide the data set into intervals: int (j) = 0  npoints (1 : j − 1) int (j) =

j=1 j = [2, 3, . . ., nlayers , nlayers + 1]

7. The matrixes on which D-optimal design was applied then became Tonion,j = [Tcs,sort (int(j) + 1 : int(j + 1)), 2 (Tcs,sort (int(j) + 1 : int(j + 1))) ] Tonion,j = Tcs,sort (int(j) + 1 : int(j + 1))

j=1 j = [2, 3, . . ., nlayers ]

8. D-optimal design was applied to different layers of the onion design according to D-optimal (Tonion,j , 2n + 2) D-optimal (Tonion,j , n + 1)

j=1 j = [2, 3, . . ., nlayers ]

where the first input to D-optimal design was the rows in the selected interval and the second input was the number of points to be selected. The method results in a selection of 1% of the total number of data points on which PCA was applied. The exponentially increasing number of points in every onion layer resulted in the majority of the points being selected at high sensitivity since an equal number of points (number of components in the PCA +1) was selected from every layer except the first. An example of the sectioning of the data set into layers and the relative parameter sensitivity of the data points in the layers is shown in Fig. A1. In the first onion layer, where the parameter sensitivity is highest, the selected rows in the centered and scaled scores matrix was expanded by additional columns containing the same columns again but squared. The reason for this was to capture the nonlinear behavior in the parameter sensitivity.

156

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157

Fig. A1. D-optimal onion layer thickness colored by relative sensitivity. The relative sensitivity is an indication of how sensitive the data point is to changes in parameter value relative to the rest of the data set. The color for one bar is set to the average relative sensitivity of the data points in that onion layer.

The D-optimal design algorithm selects rows (time points) from the different onion layers that maximizes the parameter volume (e.g. the determinant of (Tonion,j )T Tonion,j ). The D-optimal algorithm used here was that contained in the Matlab function candexch (included in Matlab Statistics Toolbox). By this algorithm the time points with the most influence on the change in determinant (highest sensitivity) could be identified. Together with the scaling of the scores matrix (T), this ensured that the data points selected in every layer widely spread the sensitivity for the different parameters. In other words, the D-optimal design algorithm aided the selection of data to provide a balanced influence from parameters and the onion layer method provided a balanced influence from different parts of the data set. The fact that different parts of the data set had different parameter sensitivity becomes apparent when one set of selected data points is studied as was illustrated in Section 4.1.

References Ansell GP, Bennett PS, Cox JP, Frost JC, Gray PG, Jones AM, et al. The development of a model capable of predicting diesel lean NOx catalyst performance under transient conditions. Appl Catal B: Environ 1996;10:183–201. Berger RJ, Kapteijn F, Moulijn JA, Marin GB, De Wilde J, Olea M, et al. Dynamic methods for catalytic kinetics. Appl Catal A: Gen 2008;342:3–28. Bernaerts K, Versyck KJ, Van Impe JF. On the design of optimal dynamic experiments for parameter estimation of a Ratkowsky-type growth kinetics at suboptimal temperatures. Int J Food Microbiol 2000;54:27–38. Blockley R. Encyclopedia of aerospace engineering. Chichester, West Sussex, UK/Hoboken, NJ: John Wiley & Sons; 2010. p. 5229–38. Box MJ. The occurrence of replications in optimal designs of experiments to estimate parameters in non-linear models. J R Stat Soc Ser B (Methodol) 1968;30:290–302. Box GEP, Hunter WG. The experimental study of physical mechanisms. Technometrics 1965;7:23–42. Box GEP, Lucas HL. Design of experiments in non-linear situations. Biometrika 1959;46:77–90. Box GEP, Hunter WG, MacGregor JF, Erjavec J. Some problems associated with the analysis of multiresponse data. Technometrics 1973;15:33–51. Chen DKS, Bissett EJ, Oh SH, Van Ostrom DL. A three-dimensional model for the analysis of transient thermal and conversion characteristics of monolithic catalytic converters. SAE, international technical paper number: 880282; 1988. Coleman TF, Li YY. An interior trust region approach for nonlinear minimization subject to bounds. Siam J Optim 1996;6:418–45. Crocoll M, Kureti S, Weisweiler W. Mean field modeling of NO oxidation over Pt/Al2 O3 catalyst under oxygen-rich conditions. J Catal 2005;229:480–9. Draper NR, Hunter WG. Design of experiments for parameter estimation in multiresponse situations. Biometrika 1966;53:525–33. Draper NR, Hunter WG. The use of prior distributions in the design of experiments for parameter estimation in non-linear situations. Biometrika 1967a;54:147–53.

Draper NR, Hunter WG. The use of prior distributions in the design of experiments for parameter estimation in non-linear situations: multiresponse case. Biometrika 1967b;54:662–5. Dubien C, Schweich D, Mabilon G, Martin B, Prigent M. Three-way catalytic converter modelling: fast- and slow-oxidizing hydrocarbons, inhibiting species, and steam-reforming reaction. Chem Eng Sci 1998;53:471–81. ˜ AA. The microkinetics of Dumesic JA, Rudd DF, Aparicio LM, Rekoske JE, Trevino heterogeneous catalysis. Washington, DC: American Chemical Society; 1993. Ericson C, Westerberg B, Odenbrand I. A state-space simplified SCR catalyst model for real time applications. SAE number: 2008-01-0616; 2008. Franceschini G, Macchietto S. Model-based design of experiments for parameter precision: state of the art. Chem Eng Sci 2008a;63:4846–72. Franceschini G, Macchietto S. Novel anticorrelation criteria for model-based experiment design: theory and formulations. AIChE J 2008b;54:1009–24. Froment GF, De Wilde J, Bischoff KB. Chemical reactor analysis and design. Hoboken, NJ: Wiley; 2011. Glielmo L, Santini S. A two-time-scale infinite-adsorption model of three way catalytic converters during the warm-up phase. J Dyn Syst Meas Control 1999;123:62–70. Goodwin GC. Dynamic system identification: experiment design and data analysis/Graham C. Goodwin and Robert L. Payne. New York: Academic Press; 1977. Guojiang W, Song T. CFD simulation of the effect of upstream flow distribution on the light-off performance of a catalytic converter. Energy Convers Manag 2005;46:2010–31. Güthenke A, Chatterjee D, Weibel M, Krutzsch B, Koˇcí P, Marek M, et al. Current status of modeling lean exhaust gas aftertreatment catalysts. In: Marin GB, editor. Advances in chemical engineering. Academic Press; 2007. p. 103–283. Hill WJ, Hunter WG, Wichern DW. A joint design criterion for the dual problem of model discrimination and parameter estimation. Technometrics 1968;10:145–60. Hunter WG, Hill WJ, Henson TL. Designing experiments for precise estimation of all or some of the constants in a mechanistic model. Can J Chem Eng 1969;47:76–80. Kandylas IP, Koltsakis GC. NO2 -assisted regeneration of diesel particulate filters: a modeling study. Ind Eng Chem Res 2002;41:2115–23. Koltsakis GC, Stamatelos AM. Modeling dynamic phenomena in 3-way catalytic converters. Chem Eng Sci 1999;54:4567–78. Koltsakis GC, Konstantinidis PA, Stamatelos AM. Development and application range of mathematical models for 3-way catalytic converters. Appl Catal B: Environ 1997;12:161–91. Lafossas F, Matsuda Y, Mohammadi A, Morishima A, Inoue M, Kalogirou M, et al. Calibration and validation of a diesel oxidation catalyst model: from synthetic gas testing to driving cycle applications; 2011. Lundberg B, Sjoblom J, Johansson Å, Westerberg B, Creaser D. Parameter estimation of a DOC from engine rig experiments with a discretized catalyst washcoat model. SAE Int J Engines 2014;7:1093–112. Martens H, Naes T. Multivariate calibration. John Wiley & Sons; 1989. Massman WJ. A review of the molecular diffusivities of H2 O, CO2 , CH4 , CO, O3 , SO2 , NH3 , N2 O, NO, and NO2 in air, O2 and N2 near STP. Atmos Environ 1998;32:1111–27. Mladenov N, Koop J, Tischer S, Deutschmann O. Modeling of transport and chemistry in channel flows of automotive catalytic converters. Chem Eng Sci 2010;65:812–26. Murzin D, Salmi T. Mass transfer and catalytic reactions. Catalytic kinetics. Amsterdam: Elsevier Science; 2005. p. 341–418.

B. Lundberg et al. / Computers and Chemical Engineering 74 (2015) 144–157 Oh SH, Cavendish JC. Transients of monolithic catalytic converters. Response to step changes in feedstream temperature as related to controlling automobile emissions. Ind Eng Chem Prod Res Dev 1982;21:29–37. Olsson L, Fridell E, Skoglundh M, Andersson B. Mean field modelling of NOx storage on Pt/BaO/Al2 O3 . Catal Today 2002;73:263–70. Olsson I-M, Gottfries J, Wold S. D-optimal onion designs in statistical molecular design. Chemom Intell Lab Syst 2004;73:37–46. Pandya A, Mmbaga J, Hayes RE, Hauptmann W, Votsmeier M. Global kinetic model and parameter optimization for a diesel oxidation catalyst. Top Catal 2009;52:1929–33. Pinto JC, Lobão MW, Monteiro JL. Sequential experimental design for parameter estimation: a different approach. Chem Eng Sci 1990;45:883–92. Ramanathan K, Sharma CS. Kinetic parameters estimation for three way catalyst modeling. Ind Eng Chem Res 2011;50:9960–79. Salomons S, Votsmeier M, Hayes RE, Drochner A, Vogel H, Gieshof J. CO and H2 oxidation on a platinum monolith diesel oxidation catalyst. Catal Today 2006;117:491–7. Sjoblom J, Creaser D. Latent variable projections of sensitivity data for experimental screening and kinetic modeling. Comput Chem Eng 2008;32:3121–9.

157

Stamatelos AM, Koltsakis GC, Kandylas IP, Pontikakis GN. Computer aided engineering in diesel exhaust aftertreatment systems design. Proc Inst Mech Engin Part D: J Automob Eng 1999;213:545–60. ˇ epánek J, Koˇcí P, Marek M, Kubíˇcek M. Catalyst simulations based on coupling of Stˇ 3D CFD tool with effective 1D channel models. Catal Today 2012;188:87–93. Tronconi E, Forzatti P. Adequacy of lumped parameter models for SCR reactors with monolith structure. AIChE J 1992;38:201–10. Tsinoglou DN, Koltsakis GC. Effect of perturbations in the exhaust gas composition on three-way catalyst light off. Chem Eng Sci 2003;58:179–92. Voltz SE, Morgan CR, Liederman D, Jacob SM. Kinetic study of carbon monoxide and propylene oxidation on platinum catalysts. Prod R&D 1973;12:294–301. Wang TJ, Baek SW, Lee JH. Kinetic parameter estimation of a diesel oxidation catalyst under actual vehicle operating conditions. Ind Eng Chem Res 2008;47:2528–37. Watling TC, Ahmadinejad M, T¸ut¸uianu M, Johansson Å, Paterson MAJ. Development and validation of a Pt–Pd diesel oxidation catalyst model. SAE technical paper number: 2012-01-1286; 2012. Zullo LC. (PhD thesis) Computer aided design of experiments: an engineering approach (PhD thesis). University of London; 1991.