A joint model-based experimental design approach for the identiﬁcation of kinetic models in continuous ﬂow laboratory reactors

Continuous ﬂow laboratory reactors are typically used for the development of kinetic models for catalytic reactions. Sequential model-based design of experiments (MBDoE) procedures have been proposed in literature where experiments are optimally designed for discriminating amongst candidate models or for improving the estimation of kinetic parameters. However, the effectiveness of these procedures is strongly affected by the initial model uncertainty, leading to suboptimal design solutions and higher number of experiments to be executed. A joint model-based design of experiments (j-MBDoE) technique, based on multi-objective optimization, is proposed in this paper for the simultaneous solution of the dual problem of discriminating among competitive kinetic models and improving the estimation of the model parameters. The effectiveness of the proposed design methodology is tested and discussed through a simulated case study for the identiﬁcation of kinetic models of methanol oxidation over silver catalyst.


Introduction
Model-based design of experiments (MBDoE) techniques represent a consolidated tool for the rapid assessment and identification of fundamental kinetic models by optimally designing a set of experiments yielding the most informative data to be used for model identification (Franceschini and Macchietto, 2008). As per conventional model building procedures (Asprey and Macchietto, 2002;Blau et al., 2008), experiments are optimally designed with the following purposes: i) discriminating between structurally identifiable candidate models, in order to identify the most suitable model structure representing a system (Hunter and Reiner, 1965;Forzatti, 1983, 1984;Schwaab et al., 2006); ii) improving the precision of parameter estimates, once a suitable model structure is determined (Galvanin et al., 2007;Bandara et al., 2009). Whilst the first objective is achieved based on the maximisation of the discriminating power (i.e. a function for quantitatively evaluating the deviation between model predictions), the second is based on the maximisation of the expected information, given as a measurement function of the Fisher Information Matrix (FIM), allowing to increase the confidence on parameter estimation. The sequential iteration of steps i) (MBDoE for model discrimination) and ii) (MBDoE for improving parameter precision) leads to the detection of the best model structure and to a sta-tistically reliable estimation of the model parameters, minimising the experimental trials required by the model identification task. This optimal design approach has been recently applied also to the design of steady-state (Reizman and Jensen, 2012) as well as transient experiments (Schaber et al., 2014) for the development of kinetic models in microfluidic devices, underlining the potential of MBDoE techniques in the identification of reaction kinetics. Nevertheless, the conventional sequential MBDoE approach used for model building is affected by several limitations due to the intrinsic nature of the optimal design problem. In fact, at the beginning of the MBDoE procedure, when both the model structure and the set of model parameters are unknown, the design for model discrimination could be highly ineffective for discriminating amongst candidate kinetic models when the optimally designed experimental conditions are applied to the actual system. Furthermore, due to model uncertainty, the planned discriminating experiments could provide a very low level of information for the estimation of the kinetic parameters, and this fact could severely affect the reliability of model predictions. Finally, the need of sequentially performing the design for model discrimination and the design for improving parameter precision procedures leads to the execution of a large number of experiments for obtaining reliable kinetics, prolonging time and effort required by the entire modelling activity.
In order to overcome these issues, Hill et al., 1968 introduced the concept of joint experimental design, i.e. a design for both establishing the form of an adequate model representing a system (i.e. the model structure) and to obtain a precise estimation of its set of parameters. A multi-objective design criterion was proposed and F. Galvanin et al. / Computers and Chemical Engineering xxx (2016) Petrov et al. (1991) and Akiti et al. (1997) applied a joint design approach for the investigation of reaction kinetics. However, in both studies the design objective functions were evaluated on a grid of experimental conditions and no direct multi-objective optimisation algorithm was applied for the design of the optimal experimental conditions. In this paper, a joint model-based design of experiments (j-MBDoE) procedure is proposed for the development of kinetic models for simultaneously discriminating amongst candidate kinetic models and improving the estimation of kinetic parameters. Preliminary data from the reactors is used for model discrimination and to screen the most informative regions of the design space. According to j-MBDoE, trade-off solutions between metrics of Fisher information and discriminating power are computed using a multi-objective optimisation algorithm and are used to design a sequence of steady-state experiments. The effectiveness of the proposed design strategy is tested and critically discussed through a case study for the identification of kinetic models of methanol oxidation over silver catalyst, where the effectiveness of different optimal design configurations is compared and quantitatively assessed. The rest of the paper is organised as follows.
In Section 2, mathematical formulations suitable for the optimal design of experiments are introduced. MBDoE techniques are introduced for improving parameter estimation (Section 2.1), for model discrimination (Section 2.2) and a new formulation (j-MBDoE) (Section 2.4) is presented for the joint improvement of parameter estimation and model discrimination. In Section 2.3 techniques for information analysis and ranking of experiments from available preliminary data are introduced. In Section 3 the optimal design procedure for MBDoE and j-MBDoE is described. In Section 4 the oxidative dehydrogenation of methanol to formaldehyde over silver catalyst case study is presented and results from the application of the proposed design techniques are discussed. Concluding remarks are presented in Section 5.

Optimal design of experiments in flow reactors: mathematical formulation
As a first approximation, an isothermal plug flow reactor can be modelled as a dynamic plug flow reactor (PFR) in the form: where c i is the species concentration, r j and ij are the reaction rate and the stoichiometric coefficient of the i-th species in the j-th reaction respectively, z is the axial coordinate, v z is the speed of fluid flow in the z-direction and t is the integration time. Eq. (1), together with the reaction rate expressions (of which mathematical structure needs to be identified), represent a system of differential and algebraic equations (DAEs) which can be written in a general form as: with the set of initial conditions x(0,0) = x 0 , where x(z,t) is the N x -dimensional vector of state variables (i.e. the concentrations c i in (1)), u(z,t) is the N u -dimensional vector of manipulated input variables (i.e. flow rate, temperature, pressure of the system) and is the N Â -dimensional set of unknown kinetic model parameters to be estimated. The symbolîs used to indicate the estimate of a variable (or of a set of variables): thus, y(z,t) is the N y -dimensional vector of measured values of the outputs, whileŷ(z,t) is the N ydimensional vector of the corresponding values estimated by the model. Note that, in a general case, state variables and input variables can be characterised in both time and space.
The optimisation of experimental conditions is carried out by computing the n ϕ -dimensional experiment design vector = [y 0 , u, t sp , ] T , which includes the N y -dimensional set of initial conditions y 0 of the measured variables (which is a subset of x 0 , since x 0 contains all the state variables defining the system, not only the measured ones), the set of manipulated inputs u, the set of time instants at which the output variables are sampled (t sp ) and (possibly) the experiment duration . These conditions can be optimised for improving parameter precision or for discriminating among candidate kinetic models depending on the MBDoE objective function used in the optimisation.

Optimal design of experiments for improving kinetic parameter estimation
Conventional MBDoE techniques for improving parameter estimation aim at decreasing the parameter uncertainty region predicted by model through the solution of the optimisation problem: Eq. (4) is subject to Eqs.
(2) and (3) (the model equations) and to a n -dimensional set of constraints on design variables, defining the design space D (i.e. the operating range of experimental decision variables). According to Eq. (4), the experiment is designed so as to minimise a measurement function of V , representing the design criterion. The most common design criteria are the A-, D-, E-optimal criteria, minimising the trace, the determinant and the maximum eigenvalue of V Â respectively (Pukelsheim, 1993), or its singular values (Galvanin et al., 2007). In (4) V and H are the variancecovariance matrix of model parameters and the Fisher Information Matrix (FIM), respectively. H is defined by where n sp is the number of sampling points of measured variables, s ij is the ij-th element of the N y × N y inverse matrix of measurements error and H 0 Â is the prior dynamic information matrix, taking into account preliminary statistical information about the parametric system before each trial is carried out.

Optimal design of experiments for kinetic model discrimination
MBDoE techniques for model discrimination are based on the maximisation of the model discrimination function, i.e. a function MD evaluating the deviations between predictions of candidate models (discriminating power). Multi-response MBDoE criteria for model discrimination are usually expressed in the following form (Hunter and Reiner, 1967): Several modifications of Eq. (6) have been proposed in literature, in order to take into account the measurements uncertainty Forzatti, 1983, 1984) and/or the a-priori statistical reliability of each model (Box and Hill, 1966;Schwaab et al., 2002). In particular, the following expression was proposed by Schwaab et al. (2002): In Eq. (7) the discriminating power MD is a function ofŷ M ,i andŷ N ,i (the i-th predicted responses of model M and N), of 2 y,i (the variance of measurement errors for the i-th response) and P i , which is the relative probability of the i-th model to be the "true" (best) model computed from a-posteriori statistics obtained after parameter estimation. In particular, P i [%] can be computed for the i-th model from the expression where 2 i is the chi-square value of the i-th model evaluated after parameter estimation (Bard et al., 1977). According to Eq. (8), which follows the formulation for relative probability suggested by Schwaab et al. (2006), the model with the lowest chi-square value (i.e. the model providing the best fitting performance) is the most reliable model for describing the system. A series of experiments provide a good discrimination between candidate models when P i is significantly higher (ideally close to 100%) for at least one model.

Ranking of experiments and information analysis
In some cases, prior to the optimal design of experiments, information from preliminary (historical) experiments might be available. Examples of preliminary experiments include: a) Unplanned experiments: experiments where the conditions are fixed by the experimenter without following a precise rationale; b) Screening experiments: experiments based on a black-box design of experiments (DoE) approach where response surface methodologies (Montgomery, 2013) are used, once the relevant experimental factors are identified together with the system (measured) responses; c) Exploratory experiments: experiments focused on specific points of the design space aiming at confirming an observed or expected behaviour (dictated by experience or knowledge on the system).
The b) approach should be preferred, because of the high level of randomisation realised by adopting DoE techniques, and their effectiveness on reducing the expected variance of observations (Pukelsheim, 1993).
The availability of preliminary data is extremely valuable for two main reasons: 1. it allows obtaining a preliminary parameter estimates for each candidate model, characterised by an estimated value 0 , and a-posteriori statistics expressed in terms of variance-covariance matrix V ; 2. given 0 , it allows for the analysis of the information acquired from each experiment for the identification of each model parameter according to the FIM (5).
In order to evaluate the relative amount of information which can be obtained for the estimation of the j-th model parameters from i-th experiment a new index, called the Relative Fisher Information (RFI) index, is introduced here: In Eq. (9) H ij is the FIM for the i-th experiment for the j-th candidate model evaluated from Eq. (5) and H j is the global information obtained from the N exp experiments for the identification of the j-th model according to a norm || . || (trace, determinant or maximum eigenvalue). In this study, the trace of FIM has been used as a suitable matrix norm. The utility of Eq. (9) is that RFI allows for a ranking of the available experiments, underlining the most informative regions of the design space D to be exploited for parameter identification for each candidate model. An illustrative example of a ranking of experiments procedure is given in Fig. 1a screening of the design space according to a conventional D-optimal design (Fig. 1a) based on response surface models (Montgomery, 2013) is carried out for evaluating the most informative regions of the experimental design space (molar fraction of reactants in this case) by using RFI for two candidate kinetic models, Model 1 (M1) and Model 2 (M2) (Fig. 1b). Note that further details regarding the experi- mental design space and the candidate kinetic models are given in Section 4.
In the figure, experiment 20 (EXP20, indicated in blue) is characterised by the highest value of RFI for both models. Hence, the experimental region with high concentration of reactants allows for a significant improvement in the estimation of M1 and M2 model parameters. Conversely, experiments 1 and 16 (EXP1, EXP16) represent poorly informative experiments for the estimation of Model 2 (red bars in the plot) and Model 1 (white bars in the plot) parameters (respectively). Even if repeated, these experiments will not improve the estimation of the parameters in an appreciable way, because they are carried out under experimental conditions for which a very low sensitivity to the model parameters is realised, according to Eq. (5). The value of the design vector providing the highest RFI for each model is used in the initialisation of optimal design algorithms in order to avoid non-informative regions of the design space and regions where the FIM is nearly singular, potentially hindering the parametric identifiability of each candidate model (Galvanin et al., 2013).
Please cite this article in press as: Galvanin . Points indicated by "MD" and "PE" represent, respectively, the optimal solution for model discrimination and the optimal solution for improving parameter estimation.

Joint model-based design of experiments (j-MBDoE)
A new joint Model-Based Design of Experiment (j-MBDoE) formulation is proposed here with the purpose of simultaneously discriminating between N M candidate models and ensuring at the same time the best possible estimation of kinetic parameters.
The multi-objective j-MBDoE can be formulated as: In (10) JD is the optimal design vector for the joint design problem, w MD M,N is the MN-th element of a N M × N M -dimensional selection matrix W MD , a matrix constituted by binary values [0,1] specified by the user to tailor the discrimination to a specific subset of candidate models. In Eq. (11) « PE is the design optimality function andw PE j is the j-th element of the N M -dimensional selection vector W PE a matrix constituted by binary values [0,1] specified by the user in order to minimise the parameter uncertainty of selected models. For example, if w MD M,N = 1 this implies that the user is interested in distriminating between the M-th and N-th model, while w MD M,N = 0 implies that the user is not interested in discriminating between the M-th and N-th model. Similarly, w PE N = 1 implies that the user is interested in improving the parameter estimation for the N-th model. If w PE N = 0 the same model won't be considered in the optimal design. The optimisation (10-11) is carried out for a number of values of , ranging from the minimum value ( = MIN , where MIN is the design optimality PE computed from MBDoE for model discrimination) to the maximum value ( = MAX , where MAX is the design optimality PE computed from the A-optimal MBDoE for improving parameter estimation) using an epsilon constraint methodology for multi-objective optimisation (Mavrotas, 2009).
An illustrative example on the application of j-MBDoE for the design of a sequence of steady-state experiments where a temperature profile is optimised is shown in Fig. 2. The optimisation Eq. (10) subject to Eqs (11) and Eqs. (2) and (3) provides a set of Pareto solutions as illustrated in Fig. 2a. Each solution of the Pareto front is represented by a different temperature and a different capability of discriminating among rival models (evaluated from MD ) or improving the estimation of parameters (evaluated from « PE ). Starting from point 1 (PE experiment for improving parameter estimation, characterised by the maximum value of « PE and the minimum value of MD ) to point 4 (MD solution, experiment for discriminating among candidate models characterised by the maximum value of MD and the minimum value of PE ) intermediate Pareto solutions (trade-off between the two main objectives) are used to design a sequence of N exp = 4 steady state experiments (Fig. 2b). Note that, after the first experiment (EXP1), if the system needs to be stabilised before proceeding with the second experiment (EXP2) a transient time ( 1 ) has to be taken into account. If the transient time between consecutive experiments is sufficiently quick (i.e. i → 0), the sequence of experiments EXP1-4 can be designed according to a single dynamic experiment (Asprey and Macchietto, 2000) where the switching times t sw i of design variables are optimised.
The analysis of the shape of the Pareto front and the distribution of optimal solutions can provide useful insights on the potential efficiency of j-MBDoE. Possible shapes and distributions of Pareto points are illustrated in Fig. 3. If intermediate solution points between PE and MD solutions follow a "C" type scenario, a good trade-off between the two objectives can be realised. Solutions are evenly distributed and there is nearly the equivalent probability of reducing the uncertainty of model parameters as well as of discriminating among competing models. If solutions follow the "A" type scenario, the situation is even more beneficial, because both discriminating power (« MD ) and design A-optimality (« PE ) can significantly be improved as compared to Scenario C by exploiting intermediate design solutions. Conversely, the "B" type scenario is associated with systems where a poor trade-off between objectives can be realised. If a high discrimination needs to be realised, Please cite this article in press as: Galvanin (11) and (2,3). Square points indicate intermediate optimal design solutions; extreme solution points indicated by "MD" and "PE" diamonds represent, respectively, the optimal solution for model discrimination and the optimal solution for improving parameter estimation. experiments tend to provide a low amount of information for the estimation of kinetic parameters (i.e. under these conditions models become resistant to an improvement in parameter precision).

Optimal experimental design procedure for the identification of kinetic models
A sequential procedure for the optimal design of experiments (i.e. application of MBDoE or of the newly proposed j-MBDoE) is illustrated in Fig. 4. The procedure starts once a number of candidate kinetic models is proposed, given preliminary information in terms of historical data (i.e. previous data from the reacting systems from unplanned experiments) and design space D (i.e. the space of variability of the elements of design vector ).
Five main consecutive steps are carried out: i) a preliminary model discrimination based on a-posteriori statistics obtained after parameter estimation is carried out; this step allows the preliminary evaluation of kinetic parametersÂ 0 and probabilities P 0 i from 2 values according to Eq. (8); ii) a ranking of preliminary experiments is realised in order to detect the experimental regions yielding the maximum amount of Fisher information (through RFI, Eq. (9)); the most informative experimental conditions are used for the initialisation of optimal design algorithms through the "initial guess" vector 0 ; iii) a MBDoE or j-MBDoE optimisation is carried out to obtain a set of solutions opt of the problem Eq. (4) or Eq. (7) subject to Eqs.
(2) and (3) (MBDoE) or Eq. (10) subject to Eq. (11) and Eqs. (2) and (3) (j-MBDoE); iv) the optimally determined experimental conditions opt are implemented in the experimental system, providing a new set of experimental data y to be used for the identification of kinetic models; v) a new parameter estimation is carried out where new estimated valuesÂ are determined for each model; a-posteriori statistics allow the model discrimination in terms of P i .
The iterative procedure involving steps iii) to v) leads to the detection of the best model structure representing the system, elucidating the most plausible kinetic mechanism for a given catalyst, ensuring at the same time a precise estimation of kinetic parameters. The procedure stops when both a single model is found to be adequate for representing the system and a statistically satisfactory estimation of its kinetic parameters is achieved.
The model adequacy has been assessed for each candidate kinetic model by using a chi-square ( 2 ) test (Stewart et al., 1998). For each model the chi-square has been computed and compared with 2 Ref , a reference value from a 2 distribution with (N − N Â ) degrees of freedom (N is the total number of experimental points, N Â is the number of model parameters). In (12) y ij is the j-th observation of the i-th measured response,ŷ ij is the relative model prediction, while 2 y i is the expected variance for the i-th measured response. The best model in terms of fitting performance is the model with the lowest value of 2 and, if 2 < 2 Ref , the chi-square test is passed and the model provides an adequate representation of experimental data.
The precision in the estimation of kinetic parameters after the execution of the designed experiments has been evaluated in terms of t-test. For a statistically precise estimation the t-value of the i-th kinetic parameter can be evaluated from In Eq. (13) t i is the t-value related to the i-th model parameter, Â i and v Â i are, respectively, the estimated value and the estimated variance of the i-th kinetic parameter obtained from maximum likelihood parameter estimation (Bard et al., 1977), t(·) is the t-value distribution with a [1/2 + (1 − ␣)/2]% confidence level and (N − N Â ) degrees of freedom, where ␣ is the statistical level of significance and N the total number of experimental points (in the system under investigation for 95% significance, ␣ = 0.05 and t(·) is approximately 2). For a statistically precise parameter estimation of the i-th model parameter, t i should be higher than t ref , a reference t-value given by a Student t-distribution with (N − N Â ) degrees of freedom.
gPROMS ModelBuilder (Process Systems Enterprise, 2012) was used as simulation software for the estimation of kinetic parameters and for the statistical assessment of model adequacy. The SQP optimiser of gPROMS has been used for the solution of the NLP optimisation problem for both MBDoE and j-MBDoE design configurations. SQP initialisation was carried out by using different starting points randomly distributed around the initial guess design vector 00 , determined in step 2 of the proposed procedure, in order to avoid the possibility of incurring into local minima. However, given the reduced number of design variables involved in each single steady-state experiment, no numerical issues were observed.

Case study: identification of kinetic models of methanol oxidation over silver catalyst
The effectiveness of the procedure proposed in Section 3 is tested on a simulated case study concerning the identification of kinetic models for the oxidative dehydrogenation of methanol to formaldehyde on silver catalyst (Galvanin et al., 2015). Data from preliminary experiments, carried out on a specifically designed silicon-glass microreactor where the silver catalyst is deposited on the channel wall by sputtering (Cao and Gavriilidis, 2005), were used for a preliminary estimation of model parameters, as well as for the initialisation of optimal design algorithms. The concentration of reactants, temperature, pressure and volumetric flow rate can be controlled and monitored in the system; concentrations are analysed online at the inlet and outlet of the microreactor using gas chromatography.
Please cite this article in press as: Galvanin  In this case study the elements of the design vector which can be optimised by experimental design are: 1. Composition of reactants in terms of molar fractions: methanol (0.07-0.14), oxygen (0.03-0.10) and water (0.02-0.22) modelled as initial conditions y 0 ; 2. Temperature T (720K < T < 830K) modelled as input u; 3. Pressure P (1-2atm) modelled as input u; 4. Flow rate F (10-100 mL/min) modelled as input u.
The ranges of operability shown in parenthesis in the above represent the design space D. Concentration measurements are available as molar fractions of CH 3 OH, O 2 , CH 2 O, H 2 , H 2 O and CO 2 at the inlet (z = 0) and outlet (z = l, where l is the reaction channel length [mm]) of the reactor and they are assumed to be corrupted by Gaussian noise with zero mean and a standard deviation of 3% on the readings. The microreactor was modelled as a plug flow reactor in the Eq. (1) form, where the reaction channel length containing the catalyst was 12.5 mm long, 0.12 mm high and 0.6 mm wide.

Formulation of candidate kinetic models
In the identification study the simplified model proposed by Andreasen et al. (2005) was used as a reference model (Model 1). According to Model 1, the following reactions constitute the base (global) mechanism: Two additional simplified kinetic models of increasing level of complexity were considered: 1. Model 2: including the following total oxidation reactions for both CH 3 OH and CH 2 O: 2. Model 3: including total oxidation reactions Eqs. (16) and (17) Table 1 Preliminary model discrimination results: chi-square statistics ( 2 ), number of parameters (N ) Akaike information criterion (AIC) and probability to be the best model (P 0 i ).

Model
Model 1  and a selective oxidation reaction The hydrogen oxidation reaction H 2 + 1/2O 2 → H 2 O (20), known to occur at higher temperatures (Schubert et al., 1998), has also been included in each reaction mechanism. The reaction rate expressions for Eqs. (14) amd (15) are the ones proposed by Andreasen et al. (2005), while power law kinetic expressions have been used for the other reactions in the proposed mechanisms. Details on reaction rate expressions are given in Appendix A.
In all these preliminary experiments He was used as an inert, the volumetric flow rate was kept at F = 26.5 mL/min and the pressure at P = 1.6 atm. A model discrimination was carried out from experimental data by assessing for each proposed model: i) the lack-of-fit in terms of 2 obtained after parameter estimation is carried out; ii) the Akaike information criterion (AIC) investigating the tradeoff between fitting capability and model complexity in terms of number of model parameters (N Â ); iii) the preliminary probability P i evaluated from Eq. (8). Results after model discrimination are shown in Table 1. Model 3 provides the most satisfactory results in terms of lack-of-fit, underlined by the lower 2 , and turns out to be the most probable "true" model (P 3 = 39%), even if with a significant uncertainty as the P i are quite similar and far from 100% (P 1 = 27%, P 2 = 34%). OFAT experiments E1-E21 do not ensure a clear discrimination among the different models and additional experiments need to be designed and carried out. Furthermore, note how Model 3 represents the most complex model in terms of number of model parameters, as clearly indicated by the relatively high AIC value, which would tend to promote the use of Model 1 (the simplest model).
From fitting results (Fig. 5b), it is clear that methanol profiles are always represented in a satisfactory way by all the proposed models, but Model 1 (solid line) is ineffective in representing oxygen, formaldehyde and CO 2 molar fractions. Interestingly, the representation of oxygen and methanol concentration as a function of the investigated temperature was significantly improved by including total oxidation reactions in the model formulation (Model 2 and 3) as compared to the original model formulation (Model 1). Furthermore, a better representation of both CO 2 and CH 2 O is realised if competitive dehydrogenation and selective oxidation pathways are included (Model 3). The model availability allows definition of the best experimental conditions in order to estimate the model parameters with the greatest precision. Fig. 5a shows the RFI evaluated from Eq. (6) for a ranking of experiments E1-E21. Each proposed model shows a different RFI response to a change in experimental conditions. In particular: i) an increment in temperature T would be beneficial for Model 2, but would be unhelpful for the estimation of Model 1 and Model 3 kinetic parameters; ii) an increment on oxygen concentration in the feed is beneficial for the identification of all the proposed models; iii) higher methanol concentration in the feed is beneficial for the identification of Model 2 and 3, while a maximum in the information level is realised for Model 1; iv) water concentration increases the information for Model 1 and 2, while does not particularly affect Model 3 identification. The experimental conditions giving the highest RFIs for the three candidate models were found to be: T = 808 K; F = 26.5 mL/min, P = 1.6 atm; maximum concentration of reactants y 0 = [y CH3OH y O2 y H2O ] T = [0.14 0.10 0.22] T . These values for the experimental conditions define 0 (initial guess design vector) and have been used in the initialization of optimal design of experiments algorithms. The optimal design solution opt provided by MBDoE or j-MBDoE is usually very sensitive to the initial guess on design variables. The availability of 0 , evaluated from preliminary experiments, allows for the robust initialisation of the optimisation algorithm, providing a sufficiently informative experiment in terms of design optimality PE .

Optimal design of experiments: proposed configurations
Two optimal design configurations are proposed for the identification of kinetic models from the execution of a maximum of N exp = 16 designed steady state experiments: 1. MBDoE: 8 experiments are used for model discrimination according to Eqs. (10) and (8) experiments are used for the improvement on parameter estimation according to a conventional A-optimal MBDoE following Eq. (4) subject to Eqs. (2) and (3); 1. j-MBDoE: 16 experiments are designed from the solution of (10) subject to (11) and to the model Eqs. (2) and (3).
Synthetic "experimental" data are obtained by simulation of Model 3 with = [A i E a,i ] T = [3.20E + 07 2.58E + 05 4.00E + 04 1.36E + 05 3.25E + 05 2.50E + 03 1.44E + 05 1.62E + 05 8.67E + 04 1.44E + 05 1.51E + 05 8.50E + 04] T as the "true" parametric set defining the reaction system and by adding normally distributed noise with a mean of zero and a standard deviation of 3% on the readings. In the optimal experimental design phase, the preliminary estimationÂ 0 reported in Appendix A has been used for each model.

Model-based design of experiments (MBDoE)
Results from the application of a sequence of steady state experiments designed by MBDoE are reported in Table 2 and shown in Fig. 6 as elements of the design vector for each experiment: flow rate (F), pressure (P), temperature (T) (Fig. 6a) and concentration of reactants (y 0 ) (Fig. 6b). The simulated experimental response is given in Fig. 6c. The first eight experiments are used for discriminating among candidate models (MD series) while the remaining eight experiments are used for improving the estimation of kinetic parameters (PE series). The discrimination takes place between couples of models (M12, M23, M13) or between the full set of models (M123), while the improvement of parameter precision is realised for single models (M1, M2, M3), for couples of models (M12, M13, M23) and for the full set of candidate models (M123). In order to achieve the maximum discrimination between Model 1 and Model 2, the MBDoE optimisation pulls the system Please cite this article in press as: Galvanin Table 3 MBDoE: model discrimination results in terms of chi-square statistics ( 2 ), reference chi-square ( 2 Ref ) and probability of each model to be the best model (P i ). Models failing the 2 test are indicated in bold.  Fig. 6a). Conversely, low temperatures are used for discriminating between Model 1-2 and Model 3. Together with the A-optimal design criterion, also the D-and E-optimal design criteria have been tested. Results are not shown for the sake of conciseness, but only a relative difference on temperature profiles of 2% has been observed, and no difference in terms of all the other design variables. Results from model discrimination are given in Table 3. Note how MBDoE design configuration is very efficient on discriminating among candidate kinetic models. In fact, after the execution of the sequence of designed experiments and the subsequent estimation of kinetic parameters, Model 3 is uniquely recognised as the best model representing the reaction system (P 3 close to 98%). The model is also adequate on representing the system, as underlined by the chi-square test which is passed for Model 3 (i.e. 2 < 2 Ref ). Results obtained after parameter estimation are given in Table 4. The analysis of t-values and 95% confidence intervals shows a clear limitation of the conventional A-optimal MBDoE. After the execution of the optimally designed experiments, a particularly poor estimation of parameters A 4 , A 5 , A 6 , E a4 and E a5 is obtained. These parameters are estimated with a large uncertainty (particularly parameters A 5 and E a5 , describing the rate of reaction for the total oxidation of formaldehyde, a critical reaction limiting the selectivity to the desired product). Interestingly, a significant improvement in the precision of all the other parameters is achieved. In this case the experimental design is efficiently discriminates among candidate kinetic models, but does not provide enough information for a statistically satisfactory estimation of the full set of kinetic parameters.

Joint model-Based design of experiments (j-MBDoE)
According to j-MBDoE, the first four experiments are designed for the simultaneous discrimination and improvement of parameter precision of the three candidate models (M123) while the remaining experiments are applied to couples of models (M12, F. Galvanin et al. / Computers and Chemical Engineering xxx (2016) xxx-xxx    M13, M23). Fig. 7a and Fig. 7b show the results in terms of Pareto solutions of the j-MBDoE multi-objective design problem. Results are quite interesting. When j-MBDoE is applied to Model 1 and Model 2 only (M12), experiments providing high discriminating power can be executed ensuring at the same time a good level of information for the estimation of kinetic parameters (filled circles in the plot, experiments jMB5-jMB8). However, when Model 3 is involved in the design, experiments providing high discriminating power tend to provide a low amount of information. This is particularly noticeable when experiments are designed for discriminating between Model 1 and Model 3 (filled squares in the plot, experiments jMB9-jMB12) and Model 2 and Model 3, where a great variability on A-optimality is realised (filled triangles in the plot, experiments jMB13-jMB16). The optimal values of experimental design variables are given in Table 5 and illustrated in Fig. 7c and Fig. 7d. Fig. 7e shows the values of simulated experimental responses for the optimally designed experiments. Interestingly, flow rate and water concentration in the feed are manipulated in order to exploit the trade-off solutions when j-MBDoE is applied to Model 1 and Model 3, while methanol concentration and temperature are primarily manipulated when j-MBDoE is applied to Model 2 and Model 3. These excitation patterns generate a strong variability on water, formaldehyde and oxygen readings (experiments jMB9-jMB16, Fig. 7e).

Discussion
Analysing the results in terms of model discrimination (Table 6) it is apparent that the j-MBDoE design efficiency is comparable to the one realised by a conventional A-optimal MBDoE (Table 3). Results are similar both in terms of model discrimination (P i is still around 97% for Model 3, which is correctly detected as the "true" model) and model adequacy (i.e. 2 test is amply satisfied for Model 3) on representing the experimental points. However, results become more interesting by analyzing the parameter estimation results (Table 7) where the superiority of j-MBDoE becomes apparent. In this case, thanks to the exploitation of the trade-off design solutions, j-MBDoE allows the full set of kinetic parameters to be estimated in a statistically sound way. Note how the amount of Fisher information generated by a sequence of experiments designed by j-MBDoE is higher than the one realised by a conventional A-optimal MBDoE, and this was predicted also in the experimental design phase. This is clear from the design results shown in the last row of Table 2 (MBDoE) and 5 (j-MBDoE). In fact, the sum of design optimality functions PE predicted by MBDoE was 259.12 against 348.02 for j-MBDoE. At the same time, the sum of discriminating power MD obtained from the two configurations is very similar: 27.5 for MBDoE, 34.03 for j-MBDoE. As a result, both the configurations are very effective for discriminating among competitive kinetic models, as underlined by the very similar probabilities P i obtained in Table 3 and Table 6.
Please cite this article in press as: Galvanin

Conclusion
In this paper, a new joint model-based design of experiments (j-MBDoE) procedure has been proposed for the identification of kinetic models in continuous flow laboratory reactors. The procedure allows for a multi-objective optimisation of the information obtained from experiments for the purpose of discriminating among candidate kinetic models, improving at the same time the statistical reliability on the estimation of kinetic parameters. In the first stage of the proposed procedure, preliminary historical data are exploited for ranking the available experiments and for detecting the most informative experimental conditions to be used for the initialisation of optimal design algorithms. In a second stage, experiments are optimally designed by j-MBDoE ensuring the maximum discriminating power (i.e. the maximum difference between model responses) and the maximisation of the expected Fisher information, required for the reduction of the parametric uncertainty.
A case study for the discrimination of simplified kinetic models of methanol oxidation over silver catalyst clearly shows how the newly proposed j-MBDoE procedure outperforms the conventional A-optimal MBDoE. Whilst the two approaches show a very similar capability on discriminating among candidate kinetic models, only j-MBDoE allows for a statistically sound estimation of the full set of kinetic parameters. In fact, highly informative discriminating experiments are designed by j-MBDoE thanks to the exploitation of trade-off Pareto solutions, following a compromise between design optimality and discriminating power. The analysis of Pareto solutions gives precious insights on the identifiability of the parametric system, underlining the relative resistance of candidate models to an improvement on parameter estimation. The evaluation of Relative Fisher Information (RFI) obtained from experiments with different discriminating power, shows how the most informative regions can be exploited during the experimentation, opening new important perspectives for designing experiments for the identification of kinetics.

Appendix A. Preliminary parameter estimation
A maximum likelihood parameter estimation was carried out from one factor at a time (OFAT) experiments data (E1-E21, see the main text for details). Kinetic parameters, including preexponential factors A i (expressed in [s −1 (mol/m 3 ) a+b−1−c−d ] for the generic reaction rate r i = k i c a A c b B /c c C c d D ) and activation energies E ai (expressed in [J/mol]) were computed according to the Arrhenius equation evaluated for each reaction rate and each candidate model (see Table A1 for details). In the parameter estimation procedure Eq. (A1) was reformulated as for the purpose of minimising the existing structural correlation between pre-exponential factors and activation energies (Himmelblau, 1968;Buzzi-Ferraris and Manenti, 2010). Parameters Â 1,i and Â 2,i have been estimated accordingly, following a log-likelihood parameter estimation technique (Bard et al., 1977). Preliminary parameter estimation resultsÂ 0 are given in Table A2-A4 in terms of estimated values and a-posteriori statistics including 95% confidence intervals and t-values. As can be seen Table A1 Set of reactions involved in the proposed kinetic models and related kinetic expressions.

Reactions
Reaction Number of kinetic parameters 6 10 12 from the estimation results given in Table A2-A4, only a few parameters are estimated in a statistically sound way from OFAT design (E a1 and E a2 for Model 1; A 1 , E a1 and E a2 for Model 2; E a1 and E a6 for Model 3), whilst some parameters are affected by a large uncertainty, particularly in Model 2 (A 3 ,A 4 , E a4 , E a5 ) and Model 3 (A 4 ,E a2 ).