- Email: [email protected]
A convex optimization framework for the identification of homogeneous reaction systems
Automatica 114 (2020) 108823
Contents lists available at ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
A convex optimization framework for the identification of homogeneous reaction systems✩ ∗
Ali Al-Matouq a , , Tyrone Vincent b a b
Department of Engineering Management, Prince Sultan University, Riyadh 11586, Saudi Arabia Department of Electrical Engineering and Computer Science, Colorado School of Mines 1600 Illinois St., Golden, CO 80401, USA
article
info
Article history: Received 7 May 2018 Received in revised form 4 October 2019 Accepted 30 December 2019 Available online xxxx Keywords: System identification Identification of kinetic models Non-negative Lasso Homogeneous reactions
a b s t r a c t A new convex optimization framework for the simultaneous identification of reaction stoichiometries and kinetics for homogeneous reaction systems is developed. The identification problem is posed as a non-negative Lasso program that incorporates both a set of hypothesized reaction stoichiometries and a set of hypothesized reaction rate laws. Estimation error bounds are established and an iterative technique for finding a proper value of the Lasso tuning parameter for model selection and parameter identification is also developed. The technique was tested via simulation on a continuous stirred tank reactor process demonstrating recovery of stoichiometry, reaction rate laws and near recovery of model parameters using noisy reactor species concentration measurements. This was demonstrated even for a case that was previously described to be structurally unidentifiable. The same was also demonstrated for the case when a subset of reactor species concentration measurements is only available. © 2020 Elsevier Ltd. All rights reserved.
1. Introduction Detailed first principle models of reaction processes are desired in many applications including process analysis, design, control and optimization. These models require a structure that is based on first principles and model parameters identified from experimental data. Such models are useful for providing chemical insight into the process dynamics in comparison with general (black box) models. The chemist provides a list of postulated reactions and reaction rate laws based on a set of chemical and thermodynamic principles. Experimental measurements are then used to validate different hypothesis using parameter identification techniques. However, identifying chemical reaction models from experimental data can become a very tedious task even for simple reaction processes. The difficulty arises because the identification of chemical reaction processes is both a model structure and model parameter identification problem combined that also deals with model structures that are nonlinear in the parameters (Maria, 2004). Statistical estimation theory and direct deterministic estimation ✩ This project was supported by Prince Sultan University, Saudi Arabia under Grant Number IBRP-ENG-2014-12-24. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Juan C. Aguero under the direction of Editor Torsten Söderström. ∗ Corresponding author. E-mail addresses: [email protected] (A. Al-Matouq), [email protected] (T. Vincent). https://doi.org/10.1016/j.automatica.2020.108823 0005-1098/© 2020 Elsevier Ltd. All rights reserved.
approaches are commonly used to solve these identification problems, see for example the review in Maria (2004). Parameter estimation is usually performed for each model candidate and a model is then chosen using established model selection techniques (Maria, 2004). These conventional methods will require, in most cases, large amounts of experimental data to identify the best model. Moreover, the solution is likely to be non-unique and may require solving complex minimization problems that can be very difficult to solve. Efforts to simplify the identification problem was made by Bonvin and his group (see for example Bhatt, Kerimoglu, Amrhein, Marquardt, & Bonvin, 2012; Bonvin & Rippin, 1990; Brendel, Bonvin, & Marquardt, 2006; Rodrigues, Billeter, & Bonvin, 2015) where an incremental (step by step) approach for the identification of reaction stoichiometries and kinetics is followed. Postulated reaction stochiometries are first analyzed using target factor analysis (Bonvin & Rippin, 1990) and then the most probable reaction rate laws are isolated by estimating the values of the reaction fluxes and the reaction rates before identifying the reaction parameters. These multi-stage techniques, although can reduce the computations involved, but are likely to introduce errors in each stage. Methods to overcome problems in multi-stage techniques are addressed in Willis and von Stosch (2017) using mixed integer programming which can become very difficult to solve. Particularly, non-convexity problems become computationally intractable when very large number of experiments and candidate reaction models are to be investigated. Hence, efficient simultaneous identification techniques are needed for the identification of chemical reaction models. A
2
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
few recent studies have exploited the use of convex optimization in identifying chemical reaction networks. The study in Papachristodoulou and Recht (2007) used a convex formulation with the ℓ1 norm penalty for identifying simultaneously the chemical reaction network (reaction stoichiometry) and kinetic parameters from concentration measurements. The same study, however, demonstrated through an example that multiple reaction networks can be consistent with the same measurements using this formulation and hence, there is no unique solution. This phenomenon was theoretically analyzed in Santosa and Weitz (2011) where a quantitative assessment of distinguishability of chemical reaction networks from concentration measurements was made and showed that the ℓ1 norm inverse problem is not unique and multiple reaction networks can exist that are consistent with the same given measurements. In fact this study concluded that additional a priori information is needed to address non-uniqueness problems. Another attempt using the ℓ1 norm was made in Kugler, Gaubitzer, and Muller (2009) with additional a priori information in the form of plausible reaction stoichiometry. However, the reaction rate laws were predetermined to be fixed and the resulting problem was non-convex and therefore can be subject to non-uniqueness and convergence problems as with previous techniques. In this study, both a set of hypothesized reaction stoichiometries and a set of hypothesized reaction rate laws will be incorporated in the identification problem as compared to previous studies. Also, the identification problem is formulated as a convex optimization problem that can identify the reaction stochiometry, the reaction rate laws and kinetic parameters simultaneously. We also establish the associated estimation error bounds and develop an iterative algorithm for tuning the non-negative Lasso operator that is used to search for the most plausible reaction model. The main feature of the newly developed identification framework compared to previous techniques is that the problem formulation is convex and can provide unique solutions when an appropriate tuning parameter is selected. This is particularly advantageous for problems involving vast amounts of species and postulated reaction mechanisms/reaction rate laws. In addition, the model structure and model parameters are determined simultaneously which helps to reduce uncertainties as compared with multi-stage techniques (Brendel et al., 2006). The objective of this study is to introduce both the algorithmic development and the theoretical treatment of the new identification framework and to demonstrate the potential for using it in practice through simulation experiments. The following is the outline of this study. Section 2 will introduce the linear time varying state space representation of the species mass balance equations for a continuous stirred tank reactor process. Section 3 will introduce the identification problem formulation that incorporates the species mass balance equations, the hypothesized reaction stochiometries and the hypothesized reaction rate laws. Section 4 will provide the associated parameter estimation error bounds and the assumptions required. An algorithm for searching for the most plausible reaction model is developed also in this section in view of theoretical results. Finally, Section 5 will introduce three different simulation examples demonstrating the potential of the new identification framework. The following notation is used in this study: R represents the set of real numbers; A ∈ Rm×p is an m × p matrix with real values; ∥z ∥l1 is the ℓ1 norm of vector z and |z | represents the number of non-zero elements in the vector z. row s(A) and col(A) represent the set of row vectors and column vectors of matrix A, respectively. The support of z ∈ Rp , denoted as S(z) is the set of indices associated with the non-zero elements of z while the complement support (indexing zero elements of z) is denoted by Sc (z). The vector composed of the elements of z ∈ Rp indexed by
Fig. 1. Process flow diagram of CSTR.
S(z) is denoted by zS . Similarly, the matrix that is composed of the column vectors indexed by the index set S(z) is denoted by AS ∈ Rm×|z | . Estimated variables are denoted with a hat (e.g. xˆ ), while noisy measurement vectors or matrices that contain noisy ˜ measurements use tilde symbol (e.g. y˜ and A). 2. Reactor model equations Consider reactions occurring in a continuous stirred tank reactor process (CSTR) with a single inlet and a single outlet stream as shown in the simple process flow diagram in Fig. 1. The following equations describe the species mole balance equations assuming well mixed homogeneous operation of the reactor (Rawlings & Ekerdt, 2002): d
(v (t)ci (t)) = qin (t)cin,i (t) − qout (t)ci (t) + v (t)fi (t), dt ci (0) = ci,0 , i = 1, . . . , ns
(1)
where, v (t), (L) denotes the reactor volume at time t in liters, qin (t), qout (t), (L/s) denote the inlet and outlet volumetric flow rates respectively; c1 (t), . . . , cns (t), (mol/L) represent the ns reactor species concentrations and cin,1 (t), . . . , cin,ns , (mol/L) represent the ns reactor inlet species concentrations. Here, f1 (t), . . . , fns (t), (mol/L min) represent the production (or consumption) rate (reaction flux) for each species which is the net molar amounts produced (or consumed) per unit time per unit volume by all reactions. The above species mole balance equations can also be used for reactions occurring in a batch or semi-batch reactor as explained in Rawlings and Ekerdt (2002). To simplify presentation, we will assume in this study that the reactants and products are incompressible (i.e. density is assumed constant). A different model can be used for the case when density variations are present (Bhatt et al., 2012). As a result, from the mass balance equation we may write the following: dv (t) dt
= qin (t) − qout (t),
v (0) = v0
(2)
After substituting (2) into (1), we obtain: dci (t) dt
=
) qin (t) ( cin,i (t) − ci (t) + fi (t) v (t)
(3)
We now define the reaction flux vector f(t) = [f1 (t), . . . , fns (t)]T as follows: f(t) = NT r(t)
(4) nr ×ns
where, N ∈ R is the reaction stoichiometric matrix with elements defining the stoichiometric coefficients of nr reactions involving ns species and is represented as follows:
ν11 ⎢ .. N=⎣ . νnr 1
⎡
··· .. . ···
⎤ ν1ns .. ⎥ . ⎦ νnr ns
(5)
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
where, νji gives the stoichiometric coefficient for the ith species in the jth reaction (νji ≤ 0 for reactants and νji ≥ 0 for products) for j = 1, . . . , nr and i = 1, . . . , ns . The following gives an example of how the stoichiometric matrix can be formed. Example 2.1. Consider the following mechanism for thermal decomposition of acetone as proposed in Rice, Johnston, and Evering (1932): k1
reaction#1: CH3 COCH3 −→ CH3 + CH3 CO k2
3
Assumption 3.1. The availability of a set of m sampled noisy measurements, with a sample period of Ts , for the reactor inlet ˜ 0 , . . . , u˜ m−1 , u˜ k ∈ Rns , and species concentrations, denoted as u a subset of reactor species concentrations denoted as y˜ 1 , . . . , y˜ m , where, y˜ k ∈ Rn˜ s and n˜s ≤ ns represents the number of species measured. To simplify presentation, we will assume that n˜ s = ns . The case when a subset of reactor species concentrations is measured is covered in the simulation section of this study.
reaction#2: CH3 CO −→ CH3 + CO k3
reaction#3: CH3 + CH3 COCH3 −→ CH4 + CH3 COCH2 k4
reaction#4: CH3 COCH2 −→ CH2 CO + CH3 k5
reaction#5: CH3 + CH3 COCH2 −→ CH3 COC2 H5 We order the species as follows: {CH3 COCH3 , CH3 , CH3 CO, CO, CH3 COCH2 , CH2 CO, CH4 , CH3 COC2 H5 }. Consequently, the corresponding stoichiometric matrix using the same order of species and reactions given above is as follows:
⎡
−1
⎢ 0 ⎢ N = ⎢ −1 ⎣ 0 0
1 1 −1 1 −1
1 −1 0 0 0
0 1 0 0 0
0 0 1 −1 −1
0 0 0 1 0
0 0 1 0 0
0 0 0 0 1
⎤
ν˜ 11 . ˜N = ⎢ ⎣ .. ν˜ n˜ r 1
(6)
where, kj > 0 represents the reaction rate constant for reaction j and gj (c(t), t) represents a continuous function that depends on the reactor species concentration vector c(t) defined by c(t) := [c1 (t), . . . , cns (t)]T . Temperature effects can be taken into account by including temperature as an additional function variable and using, for example, the Arrhenius law with a constant pre-exponential factor (Rawlings & Ekerdt, 2002). We can now represent our model in state space form as follows: x˙ (t) = A(t)x(t) + B(t)u(t) + f(t),
x(0) = x0
y(t) = C x(t)
x(t) = c(t), u(t) = cin (t) qin (t)
v (t)
··· .. . ···
⎤ ν˜ 1ns .. ⎥ . ⎦ ν˜ n˜ r ns
(8)
where ν˜ ℓ,i gives the stoichiometric coefficient for the ith species in the ℓth hypothesized reaction and n˜ r ≥ nr . Furthermore, we ˜ assume that row s(N) ⊆ row s(N). The above assumption implies that the correct hypothesis (or most plausible hypothesis) for the reaction stochiometry (represented by N) is part of the set of all hypothesized stochiometries ˜ (represented by N). Assumption 3.4. The availability of a set of sl hypothesized reaction rate law functions for each hypothesized reaction l given by: Rl := {g˜l,1 (y(t), t), . . . , g˜l,sl (y(t), t)}, l = 1, . . . , n˜ r
(9)
where, g˜l,q (y(t), t) represents the qth hypothesized reaction rate law function for the lth hypothesized reaction. Furthermore, we assume that the true reaction rate law corresponding to the reaction is within this set.
(7)
where,
A(t) = −
Assumption 3.3. The availability of a conjectured reaction ˜ ∈ Rn˜ r ×ns that reflects our knowledge stochiometry matrix N about the feasible reactions that may occur inside the reactor. The postulated stochiometry matrix is represented by:
⎡
⎥ ⎥ ⎥ ⎦
Going back to Eq. (4), the vector signal r(t) is defined as r(t) := [r1 (t), . . . , rnr (t)]T , where rj (t), (mol/L min) is the reaction rate for the jth reaction at time t. To simplify presentation, we will also assume that the reaction rate is governed by the following general nonlinear relationship: rj (t) = kj gj (c(t), t)
Assumption 3.2. The availability of a set of m sampled noisy measurements, with a sample period of Ts , for the reactor volume, given by v˜ 1 , . . . , v˜ m and for the inlet volumetric flow rate, given by q˜ in,0 , . . . , q˜ in,m−1 .
Ins , B(t) = −A(t), C = Cm
cin (t) := [cin,1 (t), . . . , cin,ns (t)]T Here, Ins is the identity matrix of size ns × ns . The matrix Cm ∈ Rnm ×ns is defined according to the available reactor measurements. 3. Problem formulation The problem under study is to determine the kinetic model for the CSTR reactor process which consists of finding the stoichiometry of the nr occurring reactions (given by matrix N); the correct model structure of the associated reaction rate laws (given by gj (c(t), t)) and the corresponding kinetic parameters kj for j = 1, . . . , nr from limited species concentration measurements of a CSTR process. The solution space will be limited by postulated reactions and reaction rate laws for the chemical reaction under study. We first introduce the assumptions required.
Note, that Assumption 3.4 requires that the postulated reaction rate law functions, g˜l,q (y(t), t), depend on measured reactor species concentrations y(t) only. Feasibility of postulated reaction stochiometries and reaction rate laws are based on chemical and thermodynamic principles, including, for example, consistency with element balances and transition state theory (Rawlings & Ekerdt, 2002). ˜ as We denote the jth row of N as nTj and the lth row of N
˜ Tl . The vector sequence of fluxes f˜k derived from noisy measuren ments is defined as follows: ˜ T R˜ k α∗ f˜k = N
(10)
where, R˜ 1,k ⎢ 0 ⎢
0 R˜ 2,k
⎣ 0
0 0
⎡ ˜k = ⎢ R
0 R˜ l,k =
[
{ kl :=
g˜l,1 (y˜ k , kTs ) kj eqj , 0,
0 0
··· ··· .. .
0
0
R˜ n˜ r ,k
⎤
k1 ⎢ k2 ⎥
⎡
⎥ ⎥ ⎥ , α∗ = ⎢ ⎣ ⎦
g˜l,2 (y˜ k , kTs )
⎤
.. .
⎥ ⎦
kn˜ r
···
g˜l,sl (y˜ k , kTs )
˜ Tl ∈ rows(N) and gj = g˜l,qj if n otherwise
]
4
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
for j ∈ {1, . . . , nr } and for l = 1, . . . , n˜ r . Here, eqj is the element column vector of dimension sl with 1 positioned at index qj ; i.e. the index of the correct reaction rate law in the set Rl and 0 is the zero column vector of dimension sl . Consequently, the nonzero elements of α∗ are pointing to both the true process reactions and associated reaction rate law functions describing the process. Furthermore, the values of the nonzero elements of α∗ are nothing but the kinetic parameters kj for j = 1, . . . , nr . We recognize here that α∗ can be sparse, particularly when the number of hypothesized reactions and associated reaction rate laws is significantly larger than the true number of reactions nr . Hence, identifying the vector elements of α∗ will determine both the model structure and parameters among all hypothesized model structures for the reaction process (i.e. among all hypothesized reaction stoichiometries and associated reaction rate laws). At this stage, we would like to formulate an optimization problem for identifying the most plausible reaction model, among all possible models hypothesized, that can incorporate all experimental noisy measurements while exploiting a sparse and least squares estimation framework for recovering the model structure and parameters in view of the assumptions mentioned earlier. Towards this goal, we first approximate (7) using noisy measurements y˜ k as a discrete time state space system as follows: xk = Ak−1 xk−1 + Bk−1 uk−1 + Dk−1 fk−1 y˜ k = Ck xk + ϵk
k = 1, . . . , m − 1,
x∗0 = x(0)
(11)
where, xk = x(kTs ) is the discrete state vector. The time varying matrices Ak , Bk , Dk and Ck can be found using the complete discretization method, assuming zero order hold approximation (Tóth, 2010): Ak =eA(kTs )Ts , Bk = A−1 (kTs )(Ak − I)B(kTs ) Dk =A−1 (kTs )(Ak − I), Ck = C
(12)
The above discretization method gives small discretization errors compared to other techniques as studied in Tóth (2010). Model (11) can be expanded and written as a matrix equation as follows: Y˜m = Xβ∗ + Tm Um + ζ Y where, X = [Ψ Θm ], β ∗ =
[
α∗ x∗0
]
, Ψ α∗ = Hm Fm
(13)
The above matrices and vectors are defined as follows:
⎡
⎤
⎡
C1 A 0 C2 A1 A0
⎤
⎡ ⎤ u0 ⎢ ⎥ ⎢ . ⎥ . ⎥ ⎥ , Um := ⎢ Y˜m := ⎣ .. ⎦ , Θm := ⎢ .. ⎣ .. ⎦ , ⎣ ⎦ . ∏m y˜ m um−1 Cm k=0 Am−k ⎡ ⎤ C1 B0 0 ··· 0 ⎢ ⎥ .. ⎢ ⎥ . C2 A1 B0 C2 B1 0 ⎢ ⎥, Tm := ⎢ ⎥ .. .. ⎣ ⎦ . . 0 ∏m−1 Cm k=0 Am−k B0 · · · · · · Cm Bm−1 ⎡ ⎤ C1 D 0 0 ··· 0 ⎢ ⎥ .. ⎢ ⎥ . C2 A1 D0 C2 D 1 0 ⎥, Hm := ⎢ ⎢ ⎥ .. .. ⎣ ⎦ . . 0 ∏m−1 Cm k=0 Am−k D0 · · · · · · Cm Dm−1 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ˜ T R0 f0 N ϵ0 ⎢ ⎥ ⎢ . ⎥ ⎢ . ⎥ .. Fm := ⎣ .. ⎦ , Ψ = Hm ⎣ ⎦ , ζ Y := ⎣ .. ⎦ . ˜ T Rm−1 ϵm fm−1 N y˜ 1
∑n˜
r We note here that Ψ ∈ Rm×stot , where, stot = l=1 sl represents the total number of reaction rate laws proposed for all hypothesized reactions. The above data equation can be viewed as a linear equation with a noisy observation vector Y˜m ; a regression matrix X; an unknown sparse + dense coefficient vector β∗ ; a bias term Tm Um and a measurement noise vector ζ Y . However, in our problem neither Tm Um nor X are known for certainty. Instead, a noisy version of Tm Um , denoted by T˜m U˜ m and a noisy version ˜ = [Ψ˜ Θ ˜ m ] can be formed using the noisy of X denoted by X ˜ k , y˜ k , v˜ k and q˜ in,k for k = 0, . . . , m − 1. Hence, the observations u following alternative data equation can be written:
˜ m = X˜ β∗ + ζ m Y where,
˜ m = Y˜m − T˜m U˜ m Y ˜ = X + Z˜ m X ζ m = (ζ Y + ζ U − Z˜ m β∗ )
(14)
Here, ζ U is a vector that represents the uncertainty in T˜m U˜ m and Z˜ m is an m × (stot + ns ) matrix that represents the uncertainty ˜ The estimation problem can now be phrased as follows: in X. ˜ m , the noisy regression given the vector of noisy measurements Y ˜ as defined above, find an estimate of the sparse + dense matrix X coefficient vector β∗ . The above problem can be classified as a sparse recovery estimation problem under regression matrix uncertainty in which estimators have been proposed in Loh and Wainwright (2011), Rosenbaum, Tsybakov, et al. (2010) and Zhu, Leus, and Giannakis (2011) that take into account perturbations in the regression ma˜ and the observation vector Y˜ m . However, these techniques trix X require the solution of non-convex optimization problems, not mentioning the additional a priori information required in these methods that can be difficult to obtain in practice. The estimator proposed in Rosenbaum et al. (2010) uses a convex program, however, the method requires precise knowledge about the noise present. Alternatively, it was found in Zhu et al. (2011) that using a Lasso operator to find β∗ can give relatively good results even for the case when the number of regression parameters is much higher than the number of observations. The same was also demonstrated in Rosenbaum et al. (2010) when Lasso with thresholding is used. Hence, in this study, we propose to make use of the following Lasso formulation:
ˆ σˆ 2 = β, arg min
β⪰0,σ 2
log(σ ) +
β = [αT , xT0 ]T
1 2mσ 2
∥α∥1 ∥Y˜ m − X˜ β∥22 + λ σ (15)
The above Lasso formulation was proposed in Städler, Bühlmann, and Geer (2010) using the negative log-likelihood estimation framework to estimate both the coefficient vector β and the variance σ 2 of ζ m assuming noise independence. The division of the ℓ1 norm penalty by σ was made to make the estimates equivariant under scaling of the measurements. This also helps to find a suitable value for λ approximately independent of the noise parameter σ as will be shown later. To transform the above non-convex problem to a convex one, the change of variables ρ = σ1 and φ = ρβ can be used (Städler et al., 2010). With the additional inclusion of positivity constraints, the final form of the estimator to be used and analyzed
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
˜ m = X˜ β∗ +ζ m and φ ∗ = ρβ Theorem 4.2. Suppose that Y ˆ ∗ , where:
in this study is the following:
φˆ βˆ = [αˆ T , xˆ T0 ]T = , ρˆ
where,
ˆ ρˆ = arg min −log(ρ ) + φ,
1
5
ρˆ = ∥ρ Y˜ m − X˜ φ∥22 + λwT φ
˜ Tm X˜ φˆ + Y
√
˜ Tm X˜ φˆ )2 + 4∥Y˜ m ∥2 m (Y 2 ˜ m ∥2 2∥Y 2
(18)
(16)
Suppose also that 4.1 is satisfied for some κ ∗ > 0 and that λ ≥ 1 ˜T ∥X ζ m ∥∞ , then the following error bound applies: m
Some additional remarks relating to (16) are worth mentioning:
(19)
φ⪰0,ρ
2m
w = [1Tstot 0T ]T
(1) The positivity constraints on φ were imposed since the kinetic parameters and the initial reactor concentration measurements x0 must be positive. (2) The third penalty is equivalent to the ℓ1 norm penalty in (15) since positivity constraints on φ are imposed. This penalty exploits solution sparsity in α and is used here for searching for the most plausible model structure among all hypothesized reactions and reaction rate laws. A method for selecting a suitable value for the tuning parameter λ will be introduced later. (3) Additional constraints can be imposed including, for example, relationships between kinetic parameters and upper and lower bounds if needed.
√ nr + ns ∥βˆ − β∗ ∥2 ≤ 3λ √ ρκ ˆ ∗
Proof. See Appendix A. Theorem 4.2 is applicable for both low dimensional settings, when m × ns ≫ (stot + ns ) and high dimensional settings; i.e. when (stot + ns ) ≫ m × ns . The result demonstrates that the error can be made small, for example, by (1) increasing the number of measurements and (2) increasing the value of κ ∗ . This can be done, for example, using experiment design that can reduce ˜ As a special case, we the correlations between the columns of X. present the following corollary. Corollary 1. Suppose that Assumption 4.1 holds and φ ∗ = ρβ ˆ ∗. Suppose also that the elements of ζ m can be approximated as zero mean i.i.d. noise with v ar(ζ m ) ≤ σ 2 , if:
√
4. Estimation performance and model selection
λ=2
stot + ns m
Cmax ρσ ˆ
(20)
then with probability 1/m the error bound becomes:
In the following, we are interested in analyzing the performance of the estimator in (16) and developing a model selection criterion. We will analyze the important assumptions needed for bounding the error and then analyze the factors that influence estimation performance with the objective of designing a suitable model selection criterion.
ρˆ (stot + ns )(nr + ns ) Cmax σ (21) mκ ∗ ˜ j ∥2 for j ∈ {1, . . . , stot + ns } and X˜ j represents where, Cmax = maxj ∥X ˜ the jth column of X.
4.1. Error bounds
Proof. See Appendix A.
Our approach for performance evaluation will be bounding ˆ − β∗ . Small error will indicate the rethe ℓ2 norm error of β covery of reactions associated with large kinetic parameters yet with possibly false inclusion or exclusion of reactions associated with small kinetic parameters. Our approach will make use of methods developed in Hastie, Tibshirani, and Wainwright (2015) and Wainwright (2009, 2019) for finding ℓ2 norm error bounds for unconstrained Lasso using the restricted eigenvalue assumption. We first present this key assumption. Assumption 4.1 (Restricted Eigenvalue Condition Hastie et al., ˜ satisfies the following: 2015). The noisy regression matrix X 1 ∥X˜ v∥22 ≥ κ ∗ ∥v∥22 , κ ∗ > 0, ∀v ∈ C , m where, C := {∥vSc∗ ∥1 ≤ 3∥vS ∗ ∥1 }
(17)
Here, vS ∗ represents the subvector of v with support S , which is the set of non-zero elements of the true coefficient vector φ ∗ . Similarly, vSc∗ represents the subvector of v with support index Sc∗ corresponding to the zero elements in φ ∗ . The above condition is called the restricted eigenvalue condition since without the constraint on v , κ will represent the minimum eigenvalue of the ˜ T X. ˜ The above assumption can be satisfied symmetric matrix X ˜ is rank deficient for some special cases as discussed even when X in Wainwright (2019). We now present the theorem for bounding ˆ. the error on β ∗
√
∥βˆ − β ∥2 ≤ 6 ∗
4.2. Model selection and parameter identification As evident from Theorem 4.2, the selection of λ can influence the estimation error and consequently influence both the resulting model structure and estimated kinetic parameters. In light of Theorem 4.2, we desire to find a technique that repeatedly solves (16) for different values of λ in order to find the most plausible reaction model structure and parameters. Using first ˆ ≻ 0 when: order optimality conditions, we can easily show that φ
λ < λmax = max √ j
˜ T Y˜ m X j ˜ m ∥2 m∥Y
,
j = 1, . . . , stot
(22)
The technique developed for finding λ is shown as a flow chart in Fig. 2 and is explained next. The maximum number of iterations imax is first specified and a sequence of values λ1 = 0 < λ2 < · · · , < λimax = λmax is generated on a logarithmic scale to have more iterations near zero as follows:
λi = log10 10(iλmax )/imax , i = 2, 3, . . . , imax − 1
(23)
Consequently, at each iteration i, the minimization problem (16) is solved using λi and the ith active support index is found ˆ j ≥ ϵmin }, where ϵmin ≪ 1 represents a user as S i = {j|φ
ˆ j . Next, defined threshold for considering non-zero elements of φ Assumption 4.1 is verified by finding a predicted value for κ ∗ assuming the candidate active support index S i by solving the following optimization problem: κˆ i := min v∈Ci
1 ∥X v∥22 m ∥v∥22
6
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
more input excitation are required. Generally, more measurements can help to satisfy Assumption 4.1 and reduce the estimation errors as implied in (19), while more input excitation can help to increase the value of κ ∗ . We note here that imax should be selected large enough to properly cover the candidate model structures within the solution path. In this study, trial and error was used to select imax . 5. Simulation study In the following, we demonstrate our results using the same simulation model used in Brendel et al. (2006) and in Bhatt et al. (2012) for comparison. The example involved a homogeneous reaction system for the acetoacetylation of pyrrole with diketene. The reaction stoichiometry is given by: K
P + D −→ PAA K D + D −→ DHA K D −→ OL K PAA + D −→ G where, (P) represents pyrrole, (D) diketene, (PAA) 2-acetoacetyl pyrrole, (DHA) dehydroacetic acid, (OL) oligomers and (G) is a byproduct. The corresponding true stoichiometric matrix using the same order is given by:
⎡ −1 ⎢−2 N=⎣ −1 −1
−1
+1
0 0 0
0 0 −1
0 +1 0 0
0 0 +1 0
⎤
0 0 ⎥ 0 ⎦ +1
and the true reaction rates are described by: Fig. 2. Algorithm for finding a suitable value for λ and final model identification.
where, Ci := {∥vS i ∥1 ≤ 3∥vS i ∥1 }
(24)
c
The above minimization problem can be solved using, for example, sequential quadratic programming with several iterations until convergence. This will provide an estimate for κ for the sake of verifying Assumption 4.1 at iteration i and for calculating a threshold value. The threshold value kmin is found based on the error bound given in (19) as follows:
√
n˜ r + ns kmin = 3λi √
(25)
ρˆ κˆ i i
This threshold value is then used to find the active support Sˆ ∗ ˆ j ≥ kmin }. Based on this active support, the defined as Sˆ ∗ = {j|β solution of the following minimization problem is found: ∗ 1 φˆ Sˆ ∗ , ρˆ ∗ = arg min −log(ρ ) + ∥ρ Y˜ m − X˜ Sˆ ∗ φSˆ ∗ ∥22 φ⪰0,ρ
(26)
2m
∗
ˆ Sˆ ∗ is augmented with zeros according to Sˆc∗ to form The solution φ ∗ φˆ . This solution is then tested for feasibility by testing whether ˆ ∗ , . . . kˆ ∗ ]T = the identified kinetic parameter vector, given by [k 1 n˜ r
ˆ ρˆ ∗ , contains at most one kinetic parameter per reaction wT φ/ (i.e. one reaction rate law per reaction). The first model in the solution path that passes the above feasibility test is flagged as the most plausible model based on the criteria to select the model with the largest number of feasible reactions. The solution path can be further investigated to determine if the reaction process can be modeled with fewer number of reactions, if desired. If on the other hand, no model was found that satisfies the above selection criterion, than more experimental data and/or
ra (t) =ka cP (t)cD (t)cK (t), rb (t) = kb cD (t)2 cK (t) rc (t) =kc cD (t),
rd (t) = kd cPAA (t)cD (t)cK (t)
with kinetic parameters given by ka = 0.0523, kb = 0.1279, kc = 0.0281 and kd = 0.1. As in Brendel et al. (2006), the catalyst is assumed to deactivate according to the formula: cK (t) =
v0 cK v (t) 0
where cK0 is the initial concentration of catalyst in the reactor. We now present three different identification scenarios to demonstrate the potential of the technique. 5.1. Scenario I: kd = 0 For scenario 1 of the present simulation study the input parameters used in this experiment are identical to the simulation scenario given in Bhatt et al. (2012) and Brendel et al. (2006) for comparison, presented again here. A set of 32 reaction rate law candidates are assumed for the four reactions described earlier, corresponding to feasible power-law kinetic functions g˜i,j (c(t), t) shown in Table 1. The rate constant of the fourth reaction was set to zero; i.e. kd = 0 which makes the fourth reaction not occurring to test whether the identification method can also determine the correct reaction stoichiometry. Simulation of the true reactor model with TS ∗ = 10 s, m = 361 and tf = 60min was first done with the following initial conditions: v (0) = 0.5 L, cD (0) = 0.14 mol/L, cP (0) = 0.3 mol/L, cPAA (0) = 0.08 mol/L, cDHA (0) = 0.01 mol/L, cOL (0) = 0.01 mol/L, cG (0) = 0.01 mol/L and cK 0 = 2 mol/L. Furthermore, the inlet feed was set to qin (kTS ∗ ) = 5 ml/min and cDin (kTS ∗ ) = 6 mol/L, for k = 1, . . . , m while all other species concentrations set to zero. The true concentration measurement vector was formed as: yk = [cD (tk ), cP (tk ), cPAA (tk ), cDHA (tk ), cOL (tk ), cG (tk )]T
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823 Table 1 Hypothesized reaction rate laws. Reaction a: K
Reaction b:
P + D −→ PAA
K
D + D −→ DHA
Reaction c:
Reaction d:
D −→ OL
PAA + D −→ G
K
g˜a,1 = 1
g˜b,1 = 1
g˜c ,1 = 1
g˜d,1 = 1
g˜a,2 = cD
g˜b,2 = cD
g˜c ,2 = cD
g˜d,1 = cD
g˜a,3 = cP
g˜b,3 = cD2
g˜c ,3 = cD2
g˜d,3 = cPAA
g˜a,4 = cK
g˜b,4 = cD cK
g˜c ,4 = cD cK
g˜d,4 = cK
g˜a,5 = cD cP
g˜b,5 = cD2 cK
g˜c ,5 = cD2 cK
g˜d,5 = cD cPAA
g˜a,6 = cP cK
g˜b,6 = cK
g˜c ,6 = cK
g˜a,7 = cD cK
g˜d,6 = cPAA cK g˜d,7 = cD cK
g˜a,8 = cD cP cK
g˜d,7 = cD cPAA cK
g˜a,9 = cD cP2
2 g˜d,9 = cD cPAA
g˜a,10 = cD2 cP
g˜d,10 = cD2 cPAA
where, tk = kTS ∗ . Consequently, noisy concentration measurements were then formed as: y˜ k = yk + ϵk
ϵk = σ · ηk ,
σ = 0.01 × [0.38, 0.45, 0.45, 0.63, 0.52, 0.05]T
and ηk ∈ R6 is random vector sequence generated using Matlab’s normally distributed pseudorandom function (randn). Following our proposed algorithm in Fig. 2, the dictionary Ψ˜ ∈ R361×32 was first formed using simulated noisy concentration measurements for all measured species y˜ k and using the candidate reaction rate laws shown in Table 1. Afterwards, the ˜ m and matrices Θ ˜ m and X˜ were formed according to their vector Y definitions in (14). The maximum tuning parameter λmax was calculated using (22) and was found to be 42. A set of imax = 20 tuning parameters was then generated using (23). The thresholds ϵmin = 1 × 10−6 was used. Consequently, the algorithm in Fig. 2 was implemented using Matlab CVX (Grant & Boyd, 2012) over the solution path starting with λ1 = 1 × 10−6 . The solution time for each iteration in the solution path was found to be approximately 1.5 s using a Macbook pro with 2.7 GHz Intel Core i5 and 8 GB RAM. In every iteration i we found an estimate for κi by solving (24) using Matlab’s ‘‘fmincon’’ function with the option of sequential quadratic programming solver and maximum of 500 iterations. The results for the identification step for each iteration of the algorithm in 2 for i = 1, 2, 3, 4, 5 are shown in Fig. 3. The values of kˆ ∗i for each iteration are plotted as a stem plot. As shown in the figure, after 4 iterations, it was found that κˆ 4 = 7 × 10−5 and associated threshold was kmin = 0.021. At this iteration the solution first became feasible. The associated estimation ℓ2 norm ˆ S ∗ − βˆ ∗∗ ∥2 = 1 × 10−4 , shown in the error was found to be ∥β S last plot in the same figure, demonstrating near recovery of the true parameters. Fig. 4 shows the true, measured and estimated reactor species concentrations for this experiment. In comparison with the method in Brendel et al. (2006), the method gave a corresponding error of 7 × 10−4 . 5.2. Scenario II: kd ̸ = 0 Reference to scenario I explained above, it was mentioned in Bhatt et al. (2012) that the fourth reaction was excluded in the analysis due to lack of structural identifiability. We simulated the exact case described above when the fourth reaction was included by setting kd = 0.1 to test this scenario using the proposed technique in this study. The other simulation details are exactly similar to the details discussed in scenario I above except that the inlet concentration for species D was set randomly as cDin (kTS ∗ ) = 6+ηk mol/L, where ηk ∈ R is a zero mean normally
7
distributed pseudorandom number with variance of 0.01 generated in Matlab in order to introduce random input excitation. We have also introduced zero mean random noise with variance of 0.001 in the reactor volume flow rate measurements q˜ in,k to make the problem more difficult to estimate. The algorithm in Fig. 2 was then implemented as before using the same settings mentioned in scenario I above. After 5 iterations, a feasible model was found, see Fig. 5. The associated estimation errors were found and shown in the bottom right plot demonstrating near exact recovery of the true parameters. The example demonstrates that the new technique is capable of estimating a previously described unidentifiable problem by introducing more reactor inlet excitation. 5.3. Scenario III: Subset of reactor species measured In this simulation experiment, we would like to examine the case when a subset of reactor concentration measurements is available. Using the same homogeneous reaction system example explained earlier in scenario 1, we tested the case when only concentration measurements for species D, P , PAA and OL are available (i.e. cDHA and cG are not measured). In this case, the dictionary matrix Ψ˜ will be similar to scenario 1 since the hypothesized reaction rate laws depend on the first three measured species only. Consequently, matrix C in (7) was redefined as follows:
⎡
1 ⎢0 C =⎣ 0 0
0 1 0 0
0 0 1 0
0 0 0 0
0 0 0 1
⎤
0 0⎥ 0⎦ 0
To make this simulation experiment even more difficult to identify, we have introduced random white noise in the volume measurements v˜ k and the volume flow rate measurements q˜ in,k using zero mean white noise with variance of 0.001. The other simulation experiment conditions (model parameters, inlet conditions, initial condition, time step and simulation time) are identical to the conditions described earlier for scenario 1 except that λ1 was set to 1 × 10−5 . The plot in Fig. 6 shows the result of implementing the algorithm in Fig. 2 for the fourth iteration. The figure demonstrates again near exact recovery of reaction stoichiometry, the correct reaction rate laws and near exact recovery of the kinetic parameters despite not measuring two of the reactor species and despite having noisy reactor volume and flow rate measurements. 6. Conclusion A new identification framework for the identification of reaction stoichiometries and kinetic parameters for homogeneous reaction systems is presented. The identification problem makes use of hypothesized reaction stochiometries and hypothesized reaction rate laws. The new identification framework is demonstrated for a CSTR process. A dictionary of plausible reaction rate trajectories is formed from noisy concentration measurements. Subsequently, a special non-negative Lasso formulation that is equi-variant under scaling was used to determine the correct model structure among all possible model candidates and the associated kinetic parameters simultaneously. Estimation error bounds were found and are applicable for both the low and high dimensional settings. A simple solution path searching technique to find the most plausible reaction model structure and parameters in view of the theoretical treatment was also introduced. The methodology was illustrated on a simple example used in the literature demonstrating exact recovery of reaction stoichiometry, model structure and near recovery of model parameters using noisy measurements under three different scenarios. In comparison with previously studied identification frameworks, we mention the following:
8
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
Fig. 3. Identified vs. true kinetic parameters for simulation scenario 1 (when kd = 0). Each plot shows one iteration of the algorithm in Fig. 2.
Fig. 4. True, measured and estimated reactor species concentrations for simulation scenario #1.
(1) The new technique is capable of recovering the reaction
true reaction stoichiometry from a postulated stoichiomet-
stoichiometry, the associated reaction rate laws and near
ric matrix that may include extraneous reactions. Also, the
exact recovery of the kinetic parameters simultaneously
new technique does not require the estimation of reac-
using one convex optimization problem. Hence, techniques
tion fluxes, extents of reactions nor reaction rates prior to
for elucidating the correct reaction stoichiometry are not
identification. Direct identification is achieved from exper-
needed because the method can concurrently recover the
imental data which minimizes errors.
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
9
Fig. 5. Identified vs. true kinetic parameters for simulation scenario 2 (when kd = 0.1). Each plot shows one iteration of the algorithm in Fig. 2.
Extensions to the technique are needed for other types of reactors, including, plug flow reactors (Bermúdez, Esteban, Ferrín, Rodríguez-Calo, & Sillero-Denamiel, 2017) and distributed reaction systems (Rodrigues, Billeter, & Bonvin, 2017). Acknowledgments The authors would like to thank the anonymous reviewers for their valuable comments that helped to improve this work. Appendix A. Proof of Theorem 4.2 From optimality of the solution of (16) we may write:
− log(ρˆ ) + Fig. 6. Kinetic parameters identified vs. true kinetic parameters for simulation scenario 3.
(2) Since the new technique requires the solution of one convex optimization problem, the method has the potential for identifying large kinetic models that involve large number of reactor species; experimental data and hypothesized reaction rate laws and stochiometries compared to previous techniques. (3) The solution path, determined by the iterative technique shown in Fig. 2, can help to determine if the reaction process can be modeled with fewer number of reactions, if desired. The method can also be extended to other cases, for example when density are and/or temperature variations are present. These cases require additional measurements and some modifications to the mole balance and reaction rate equations used.
≤ − log(ρˆ ) +
1 2m 1 2m
ˆ 22 + λwT φˆ ∥ρˆ Y˜ m − X˜ φ∥ ∥ρˆ Y˜ m − X˜ φ ∗ ∥22 + λwT φ ∗
Defining vˆ := φˆ − φ ∗ and expanding the above inequality we obtain the following: 1 { 2m
} ∥X˜ vˆ ∥22 − 2ρˆ Y˜ Tm X˜ vˆ + 2(φ ∗ )T X˜ T X˜ vˆ + λwT vˆ ≤ 0
˜ m = X˜ β ∗ + ζ m and the substitution Using the data equation Y φ ∗ = ρ ∗ β∗ , then after rearranging we obtain: 1 2m
∥X˜ vˆ ∥22 ≤
ρˆ
(
m
ζm +
ρˆ − ρ ∗ ˜ ∗ Xφ mρ ∗
)T
˜ vˆ − λwT vˆ X
Using Holder’s inequality, we can write: 1 2m
∥X˜ vˆ ∥22 ≤
1 m
∥X˜ T ζ φ ∥∞ ∥ˆv ∥1 − λwT vˆ
where,
ζ φ = ρζ ˆ m+
ρˆ − ρ ∗ ˜ ∗ Xφ ρ∗
(A.1)
10
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
Hence, if ρ ∗ ≈ ρˆ , then ζ φ = ρζ ˆ m . Now considering that ∥ˆv ∥1 = ∥ˆvS ∗ ∥1 + ∥ˆvSc∗ ∥1 and assuming:
λ
2 m the right hand side of (A.1) can be written as:
(A.2)
≤
3λ
(∥ˆvS ∗ ∥1 + ∥ˆvSc∗ ∥1 ) + λ(∥ˆvS ∗ ∥1 − ∥ˆvSc∗ ∥1 )
2
∥ˆvS ∗ ∥1 ≤
3λ √ 2
nr + ns ∥ˆv ∥2
∥ˆvSc∗ ∥1 ≤ 3∥ˆvS ∗ ∥1 and hence vˆ ∈ C . Consequently, following Assumption 4.1 we obtain:
√
nr + ns
κ∗
Finally, assuming that ρ ∗ ≈ ρˆ , we obtain the error bound in (19). For the expression for ρˆ , we first note that the Lagrangian of (16) can be written as (Boyd & Vandenberghe, 2004): 1 ∥ρ Y˜ m − X˜ φ∥22 + (λwT − µT )φ 2m where, µ is a Lagrange multiplier corresponding to the positivity constraint on φ . Now, ρ = ρˆ become a minimizer of (16) if the following KKT condition is satisfied (Boyd & Vandenberghe, 2004): L(φ, ρ, µ) = − log(ρ ) +
∇ρ L = −
1
ρˆ
+
ρˆ ˜ 2 1 ∥Ym ∥2 − Y˜ Tm X˜ φˆ = 0
m
m
(A.3)
Accordingly, (A.3) is a quadratic equation in ρˆ which has a positive solution given by (18). Appendix B. Proof of Corollary 1 We desire to select λ so that (A.2) is satisfied with high probability. Since ρ ∗ ≈ ρˆ is assumed then ζ φ = ρζ ˆ m . First, we define z as follows:
ρˆ
˜ T ζm X m Consider an arbitrary j ∈ {1, . . . , stot + ns }, then the jth element of z (denoted by zj ) can be written as: z :=
zj =
ρˆ ˜ T Xj ζ m
m
˜ T denotes the jth row of X˜ T . Assuming that ζ m can be where X j approximated as zero mean, white noise with variance σ 2 , we can say that: var(zj ) ≤
ρˆ 2 σ 2 ˜ 2 ∥Xj ∥2 2
m Consequently, using Chebyshev’s inequality we have:
P{|zj | ≥ t } ≤
ρˆ 2 σ 2 ˜ 2 ∥X j ∥2 2 2
m t Since ∥z ∥∞ corresponds to the maximum over stot + ns variables, we have:
P{∥z ∥∞ ≥ t } ≤ (stot + ns )
Cmax ρσ ˆ
and λ = 2t we have the following tail bound:
λ
}≤
1
References
where the first inequality is true because vˆ Sc∗ = φˆ Sc∗ ≥ 0. Now 1 since 2m ∥X˜ vˆ ∥22 ≥ 0, from the first inequality above we have:
∥φˆ − φ∗ ∥2 ≤ 3λ
m
2 m and the result follows.
m
2
stot + ns
P{∥z ∥∞ ≥
ρˆ ˜ T ∥X ζ m ∥∞ ∥ˆv ∥1 − λwT vˆ λ
√ t=
ρˆ ≥ ∥X˜ T ζ m ∥∞
≤
˜ j ∥2 for j ∈ {1, . . . , stot + ns }. Setting t as: where, Cmax = maxj ∥X
ρˆ 2 σ 2 m2 t 2
2 Cmax
Bermúdez, Alfredo., Esteban, N., Ferrín, J. L., Rodríguez-Calo, J. F., & SilleroDenamiel, M. R. (2017). Identification problem in plug-flow chemical reactors using the adjoint method. Computers & Chemical Engineering, 98, 80–88. Bhatt, Nirav, Kerimoglu, Nimet, Amrhein, Michael, Marquardt, Wolfgang, & Bonvin, Dominique (2012). Incremental identification of reaction systems: A comparison between rate-based and extent based approaches. Chemical Engineering Science, 83, 24–38. Bonvin, D., & Rippin, D. W. T. (1990). Target factor analysis for the identification of stoichiometric models. Chemical Engineering Science, 45(12), 3417–3426. Boyd, Stephen Poythress, & Vandenberghe, Lieven (2004). Convex optimization. Cambridge University Press. Brendel, Marc, Bonvin, Dominique, & Marquardt, Wolfgang (2006). Incremental identification of kinetic models for homogeneous reaction systems. Chemical Engineering Science, 61(16), 5404–5420. Grant, Michael, & Boyd, Stephen (2012). CVX: Matlab software for disciplined convex programming, version 2.0 beta. http://cvxr.com/cvx. Hastie, Trevor, Tibshirani, Robert, & Wainwright, Martin (2015). Statistical learning with sparsity: The lasso and generalizations. CRC Press. Kugler, Philipp, Gaubitzer, Erwin, & Muller, Stefan (2009). Parameter identification for chemical reaction systems using sparsity enforcing regularization: A case study for the chlorite iodide reaction. The Journal of Physical Chemistry A, 113(12), 2775–2785. Loh, Po-Ling, & Wainwright, Martin J. (2011). High-dimensional regression with noisy and missing data: Provable guarantees with non-convexity. In Advances in neural information processing systems (pp. 2726–2734). Maria, G. (2004). A review of algorithms and trends in kinetic model identification for chemical and biochemical systems. Chemical and Biochemical Engineering Quarterly, 18(3), 195–222. Papachristodoulou, Antonis, & Recht, Ben (2007). Determining interconnections in chemical reaction networks. In Proc. American control conference (pp. 4872–4877). Rawlings, James Blake, & Ekerdt, John G. (2002). Chemical reactor analysis and design fundamentals. Nob Hill Pub, Llc. Rice, F. O., Johnston, W. R., & Evering, B. L. (1932). The thermal decomposition of organic compounds from the standpoint of free radicals. ii. experimental evidence of the decomposition of organic compounds into free radicals. Journal of the American Chemical Society, 54(9), 3529–3543. Rodrigues, Diogo, Billeter, Julien, & Bonvin, Dominique (2015). Incremental model identification of distributed two-phase reaction systems. IFAC-PapersOnLine, 48(8), 266–271. Rodrigues, Diogo, Billeter, Julien, & Bonvin, Dominique (2017). Generalization of the concept of extents to distributed reaction systems. Chemical Engineering Science. Rosenbaum, Mathieu, Tsybakov, Alexandre B., et al. (2010). Sparse recovery under matrix uncertainty. The Annals of Statistics, 38(5), 2620–2651. Santosa, Fadil, & Weitz, Benjamin (2011). An inverse problem in reaction kinetics. Journal of Mathematical Chemistry, 49(8), 1507–1520. Städler, Nicolas, Bühlmann, Peter, & Geer, Sara Van De (2010). ℓ 1-penalization for mixture regression models. Test, 19(2), 209–256. Tóth, Roland (2010). Modeling and identification of linear parameter-varying systems, vol. 403. Springer. Wainwright, Martin J. (2009). Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1 -constrained quadratic programming (lasso). IEEE Transactions on Information Theory, 55(5), 2183–2202. Wainwright, Martin J. (2019). High-dimensional statistics: A non-asymptotic viewpoint, vol. 48. Cambridge University Press. Willis, Mark J., & von Stosch, Moritz (2017). Simultaneous parameter identification and discrimination of the nonparametric structure of hybrid semi-parametric models. Computers & Chemical Engineering, 104, 366–376. Zhu, Hao, Leus, Geert, & Giannakis, Georgios B. (2011). Sparsity-cognizant total least-squares for perturbed compressive sampling. IEEE Transactions on Signal Processing, 59(5), 2002–2016.
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823 Ali Al-Matouq received the B.Sc. and M.Sc. degrees with distinction in electrical engineering from King Fahd University of Petroleum and Minerals, Dhahran, in 2002 and 2007, respectively and an M.Sc. degree in Chemical Engineering and a Ph.D. degree in Electrical Engineering from Colorado School of Mines in 2012 and 2014, respectively. He also has 5 years of professional experiance as a process control and automation engineer with an ISA certification. He is currently an Assistant Professor in the Engineering Management department of Prince Sultan University, Riyadh. His research interests include process control and optimization with applications in chemical and biomedical engineering.
11
Tyrone Vincent received the B.S. degree in electrical engineering from the University of Arizona, Tucson, in 1992, and the M.S. and Ph.D. degrees in electrical engineering from the University of Michigan, Ann Arbor, in 1994 and 1997, respectively. He is currently a Professor in the Department of Electrical Engineering and Computer Science, at the Colorado School of Mines, Golden. His research interests include system identification, estimation, and fault detection with applications in materials processing, and energy systems.
Contents lists available at ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
A convex optimization framework for the identification of homogeneous reaction systems✩ ∗
Ali Al-Matouq a , , Tyrone Vincent b a b
Department of Engineering Management, Prince Sultan University, Riyadh 11586, Saudi Arabia Department of Electrical Engineering and Computer Science, Colorado School of Mines 1600 Illinois St., Golden, CO 80401, USA
article
info
Article history: Received 7 May 2018 Received in revised form 4 October 2019 Accepted 30 December 2019 Available online xxxx Keywords: System identification Identification of kinetic models Non-negative Lasso Homogeneous reactions
a b s t r a c t A new convex optimization framework for the simultaneous identification of reaction stoichiometries and kinetics for homogeneous reaction systems is developed. The identification problem is posed as a non-negative Lasso program that incorporates both a set of hypothesized reaction stoichiometries and a set of hypothesized reaction rate laws. Estimation error bounds are established and an iterative technique for finding a proper value of the Lasso tuning parameter for model selection and parameter identification is also developed. The technique was tested via simulation on a continuous stirred tank reactor process demonstrating recovery of stoichiometry, reaction rate laws and near recovery of model parameters using noisy reactor species concentration measurements. This was demonstrated even for a case that was previously described to be structurally unidentifiable. The same was also demonstrated for the case when a subset of reactor species concentration measurements is only available. © 2020 Elsevier Ltd. All rights reserved.
1. Introduction Detailed first principle models of reaction processes are desired in many applications including process analysis, design, control and optimization. These models require a structure that is based on first principles and model parameters identified from experimental data. Such models are useful for providing chemical insight into the process dynamics in comparison with general (black box) models. The chemist provides a list of postulated reactions and reaction rate laws based on a set of chemical and thermodynamic principles. Experimental measurements are then used to validate different hypothesis using parameter identification techniques. However, identifying chemical reaction models from experimental data can become a very tedious task even for simple reaction processes. The difficulty arises because the identification of chemical reaction processes is both a model structure and model parameter identification problem combined that also deals with model structures that are nonlinear in the parameters (Maria, 2004). Statistical estimation theory and direct deterministic estimation ✩ This project was supported by Prince Sultan University, Saudi Arabia under Grant Number IBRP-ENG-2014-12-24. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Juan C. Aguero under the direction of Editor Torsten Söderström. ∗ Corresponding author. E-mail addresses: [email protected] (A. Al-Matouq), [email protected] (T. Vincent). https://doi.org/10.1016/j.automatica.2020.108823 0005-1098/© 2020 Elsevier Ltd. All rights reserved.
approaches are commonly used to solve these identification problems, see for example the review in Maria (2004). Parameter estimation is usually performed for each model candidate and a model is then chosen using established model selection techniques (Maria, 2004). These conventional methods will require, in most cases, large amounts of experimental data to identify the best model. Moreover, the solution is likely to be non-unique and may require solving complex minimization problems that can be very difficult to solve. Efforts to simplify the identification problem was made by Bonvin and his group (see for example Bhatt, Kerimoglu, Amrhein, Marquardt, & Bonvin, 2012; Bonvin & Rippin, 1990; Brendel, Bonvin, & Marquardt, 2006; Rodrigues, Billeter, & Bonvin, 2015) where an incremental (step by step) approach for the identification of reaction stoichiometries and kinetics is followed. Postulated reaction stochiometries are first analyzed using target factor analysis (Bonvin & Rippin, 1990) and then the most probable reaction rate laws are isolated by estimating the values of the reaction fluxes and the reaction rates before identifying the reaction parameters. These multi-stage techniques, although can reduce the computations involved, but are likely to introduce errors in each stage. Methods to overcome problems in multi-stage techniques are addressed in Willis and von Stosch (2017) using mixed integer programming which can become very difficult to solve. Particularly, non-convexity problems become computationally intractable when very large number of experiments and candidate reaction models are to be investigated. Hence, efficient simultaneous identification techniques are needed for the identification of chemical reaction models. A
2
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
few recent studies have exploited the use of convex optimization in identifying chemical reaction networks. The study in Papachristodoulou and Recht (2007) used a convex formulation with the ℓ1 norm penalty for identifying simultaneously the chemical reaction network (reaction stoichiometry) and kinetic parameters from concentration measurements. The same study, however, demonstrated through an example that multiple reaction networks can be consistent with the same measurements using this formulation and hence, there is no unique solution. This phenomenon was theoretically analyzed in Santosa and Weitz (2011) where a quantitative assessment of distinguishability of chemical reaction networks from concentration measurements was made and showed that the ℓ1 norm inverse problem is not unique and multiple reaction networks can exist that are consistent with the same given measurements. In fact this study concluded that additional a priori information is needed to address non-uniqueness problems. Another attempt using the ℓ1 norm was made in Kugler, Gaubitzer, and Muller (2009) with additional a priori information in the form of plausible reaction stoichiometry. However, the reaction rate laws were predetermined to be fixed and the resulting problem was non-convex and therefore can be subject to non-uniqueness and convergence problems as with previous techniques. In this study, both a set of hypothesized reaction stoichiometries and a set of hypothesized reaction rate laws will be incorporated in the identification problem as compared to previous studies. Also, the identification problem is formulated as a convex optimization problem that can identify the reaction stochiometry, the reaction rate laws and kinetic parameters simultaneously. We also establish the associated estimation error bounds and develop an iterative algorithm for tuning the non-negative Lasso operator that is used to search for the most plausible reaction model. The main feature of the newly developed identification framework compared to previous techniques is that the problem formulation is convex and can provide unique solutions when an appropriate tuning parameter is selected. This is particularly advantageous for problems involving vast amounts of species and postulated reaction mechanisms/reaction rate laws. In addition, the model structure and model parameters are determined simultaneously which helps to reduce uncertainties as compared with multi-stage techniques (Brendel et al., 2006). The objective of this study is to introduce both the algorithmic development and the theoretical treatment of the new identification framework and to demonstrate the potential for using it in practice through simulation experiments. The following is the outline of this study. Section 2 will introduce the linear time varying state space representation of the species mass balance equations for a continuous stirred tank reactor process. Section 3 will introduce the identification problem formulation that incorporates the species mass balance equations, the hypothesized reaction stochiometries and the hypothesized reaction rate laws. Section 4 will provide the associated parameter estimation error bounds and the assumptions required. An algorithm for searching for the most plausible reaction model is developed also in this section in view of theoretical results. Finally, Section 5 will introduce three different simulation examples demonstrating the potential of the new identification framework. The following notation is used in this study: R represents the set of real numbers; A ∈ Rm×p is an m × p matrix with real values; ∥z ∥l1 is the ℓ1 norm of vector z and |z | represents the number of non-zero elements in the vector z. row s(A) and col(A) represent the set of row vectors and column vectors of matrix A, respectively. The support of z ∈ Rp , denoted as S(z) is the set of indices associated with the non-zero elements of z while the complement support (indexing zero elements of z) is denoted by Sc (z). The vector composed of the elements of z ∈ Rp indexed by
Fig. 1. Process flow diagram of CSTR.
S(z) is denoted by zS . Similarly, the matrix that is composed of the column vectors indexed by the index set S(z) is denoted by AS ∈ Rm×|z | . Estimated variables are denoted with a hat (e.g. xˆ ), while noisy measurement vectors or matrices that contain noisy ˜ measurements use tilde symbol (e.g. y˜ and A). 2. Reactor model equations Consider reactions occurring in a continuous stirred tank reactor process (CSTR) with a single inlet and a single outlet stream as shown in the simple process flow diagram in Fig. 1. The following equations describe the species mole balance equations assuming well mixed homogeneous operation of the reactor (Rawlings & Ekerdt, 2002): d
(v (t)ci (t)) = qin (t)cin,i (t) − qout (t)ci (t) + v (t)fi (t), dt ci (0) = ci,0 , i = 1, . . . , ns
(1)
where, v (t), (L) denotes the reactor volume at time t in liters, qin (t), qout (t), (L/s) denote the inlet and outlet volumetric flow rates respectively; c1 (t), . . . , cns (t), (mol/L) represent the ns reactor species concentrations and cin,1 (t), . . . , cin,ns , (mol/L) represent the ns reactor inlet species concentrations. Here, f1 (t), . . . , fns (t), (mol/L min) represent the production (or consumption) rate (reaction flux) for each species which is the net molar amounts produced (or consumed) per unit time per unit volume by all reactions. The above species mole balance equations can also be used for reactions occurring in a batch or semi-batch reactor as explained in Rawlings and Ekerdt (2002). To simplify presentation, we will assume in this study that the reactants and products are incompressible (i.e. density is assumed constant). A different model can be used for the case when density variations are present (Bhatt et al., 2012). As a result, from the mass balance equation we may write the following: dv (t) dt
= qin (t) − qout (t),
v (0) = v0
(2)
After substituting (2) into (1), we obtain: dci (t) dt
=
) qin (t) ( cin,i (t) − ci (t) + fi (t) v (t)
(3)
We now define the reaction flux vector f(t) = [f1 (t), . . . , fns (t)]T as follows: f(t) = NT r(t)
(4) nr ×ns
where, N ∈ R is the reaction stoichiometric matrix with elements defining the stoichiometric coefficients of nr reactions involving ns species and is represented as follows:
ν11 ⎢ .. N=⎣ . νnr 1
⎡
··· .. . ···
⎤ ν1ns .. ⎥ . ⎦ νnr ns
(5)
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
where, νji gives the stoichiometric coefficient for the ith species in the jth reaction (νji ≤ 0 for reactants and νji ≥ 0 for products) for j = 1, . . . , nr and i = 1, . . . , ns . The following gives an example of how the stoichiometric matrix can be formed. Example 2.1. Consider the following mechanism for thermal decomposition of acetone as proposed in Rice, Johnston, and Evering (1932): k1
reaction#1: CH3 COCH3 −→ CH3 + CH3 CO k2
3
Assumption 3.1. The availability of a set of m sampled noisy measurements, with a sample period of Ts , for the reactor inlet ˜ 0 , . . . , u˜ m−1 , u˜ k ∈ Rns , and species concentrations, denoted as u a subset of reactor species concentrations denoted as y˜ 1 , . . . , y˜ m , where, y˜ k ∈ Rn˜ s and n˜s ≤ ns represents the number of species measured. To simplify presentation, we will assume that n˜ s = ns . The case when a subset of reactor species concentrations is measured is covered in the simulation section of this study.
reaction#2: CH3 CO −→ CH3 + CO k3
reaction#3: CH3 + CH3 COCH3 −→ CH4 + CH3 COCH2 k4
reaction#4: CH3 COCH2 −→ CH2 CO + CH3 k5
reaction#5: CH3 + CH3 COCH2 −→ CH3 COC2 H5 We order the species as follows: {CH3 COCH3 , CH3 , CH3 CO, CO, CH3 COCH2 , CH2 CO, CH4 , CH3 COC2 H5 }. Consequently, the corresponding stoichiometric matrix using the same order of species and reactions given above is as follows:
⎡
−1
⎢ 0 ⎢ N = ⎢ −1 ⎣ 0 0
1 1 −1 1 −1
1 −1 0 0 0
0 1 0 0 0
0 0 1 −1 −1
0 0 0 1 0
0 0 1 0 0
0 0 0 0 1
⎤
ν˜ 11 . ˜N = ⎢ ⎣ .. ν˜ n˜ r 1
(6)
where, kj > 0 represents the reaction rate constant for reaction j and gj (c(t), t) represents a continuous function that depends on the reactor species concentration vector c(t) defined by c(t) := [c1 (t), . . . , cns (t)]T . Temperature effects can be taken into account by including temperature as an additional function variable and using, for example, the Arrhenius law with a constant pre-exponential factor (Rawlings & Ekerdt, 2002). We can now represent our model in state space form as follows: x˙ (t) = A(t)x(t) + B(t)u(t) + f(t),
x(0) = x0
y(t) = C x(t)
x(t) = c(t), u(t) = cin (t) qin (t)
v (t)
··· .. . ···
⎤ ν˜ 1ns .. ⎥ . ⎦ ν˜ n˜ r ns
(8)
where ν˜ ℓ,i gives the stoichiometric coefficient for the ith species in the ℓth hypothesized reaction and n˜ r ≥ nr . Furthermore, we ˜ assume that row s(N) ⊆ row s(N). The above assumption implies that the correct hypothesis (or most plausible hypothesis) for the reaction stochiometry (represented by N) is part of the set of all hypothesized stochiometries ˜ (represented by N). Assumption 3.4. The availability of a set of sl hypothesized reaction rate law functions for each hypothesized reaction l given by: Rl := {g˜l,1 (y(t), t), . . . , g˜l,sl (y(t), t)}, l = 1, . . . , n˜ r
(9)
where, g˜l,q (y(t), t) represents the qth hypothesized reaction rate law function for the lth hypothesized reaction. Furthermore, we assume that the true reaction rate law corresponding to the reaction is within this set.
(7)
where,
A(t) = −
Assumption 3.3. The availability of a conjectured reaction ˜ ∈ Rn˜ r ×ns that reflects our knowledge stochiometry matrix N about the feasible reactions that may occur inside the reactor. The postulated stochiometry matrix is represented by:
⎡
⎥ ⎥ ⎥ ⎦
Going back to Eq. (4), the vector signal r(t) is defined as r(t) := [r1 (t), . . . , rnr (t)]T , where rj (t), (mol/L min) is the reaction rate for the jth reaction at time t. To simplify presentation, we will also assume that the reaction rate is governed by the following general nonlinear relationship: rj (t) = kj gj (c(t), t)
Assumption 3.2. The availability of a set of m sampled noisy measurements, with a sample period of Ts , for the reactor volume, given by v˜ 1 , . . . , v˜ m and for the inlet volumetric flow rate, given by q˜ in,0 , . . . , q˜ in,m−1 .
Ins , B(t) = −A(t), C = Cm
cin (t) := [cin,1 (t), . . . , cin,ns (t)]T Here, Ins is the identity matrix of size ns × ns . The matrix Cm ∈ Rnm ×ns is defined according to the available reactor measurements. 3. Problem formulation The problem under study is to determine the kinetic model for the CSTR reactor process which consists of finding the stoichiometry of the nr occurring reactions (given by matrix N); the correct model structure of the associated reaction rate laws (given by gj (c(t), t)) and the corresponding kinetic parameters kj for j = 1, . . . , nr from limited species concentration measurements of a CSTR process. The solution space will be limited by postulated reactions and reaction rate laws for the chemical reaction under study. We first introduce the assumptions required.
Note, that Assumption 3.4 requires that the postulated reaction rate law functions, g˜l,q (y(t), t), depend on measured reactor species concentrations y(t) only. Feasibility of postulated reaction stochiometries and reaction rate laws are based on chemical and thermodynamic principles, including, for example, consistency with element balances and transition state theory (Rawlings & Ekerdt, 2002). ˜ as We denote the jth row of N as nTj and the lth row of N
˜ Tl . The vector sequence of fluxes f˜k derived from noisy measuren ments is defined as follows: ˜ T R˜ k α∗ f˜k = N
(10)
where, R˜ 1,k ⎢ 0 ⎢
0 R˜ 2,k
⎣ 0
0 0
⎡ ˜k = ⎢ R
0 R˜ l,k =
[
{ kl :=
g˜l,1 (y˜ k , kTs ) kj eqj , 0,
0 0
··· ··· .. .
0
0
R˜ n˜ r ,k
⎤
k1 ⎢ k2 ⎥
⎡
⎥ ⎥ ⎥ , α∗ = ⎢ ⎣ ⎦
g˜l,2 (y˜ k , kTs )
⎤
.. .
⎥ ⎦
kn˜ r
···
g˜l,sl (y˜ k , kTs )
˜ Tl ∈ rows(N) and gj = g˜l,qj if n otherwise
]
4
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
for j ∈ {1, . . . , nr } and for l = 1, . . . , n˜ r . Here, eqj is the element column vector of dimension sl with 1 positioned at index qj ; i.e. the index of the correct reaction rate law in the set Rl and 0 is the zero column vector of dimension sl . Consequently, the nonzero elements of α∗ are pointing to both the true process reactions and associated reaction rate law functions describing the process. Furthermore, the values of the nonzero elements of α∗ are nothing but the kinetic parameters kj for j = 1, . . . , nr . We recognize here that α∗ can be sparse, particularly when the number of hypothesized reactions and associated reaction rate laws is significantly larger than the true number of reactions nr . Hence, identifying the vector elements of α∗ will determine both the model structure and parameters among all hypothesized model structures for the reaction process (i.e. among all hypothesized reaction stoichiometries and associated reaction rate laws). At this stage, we would like to formulate an optimization problem for identifying the most plausible reaction model, among all possible models hypothesized, that can incorporate all experimental noisy measurements while exploiting a sparse and least squares estimation framework for recovering the model structure and parameters in view of the assumptions mentioned earlier. Towards this goal, we first approximate (7) using noisy measurements y˜ k as a discrete time state space system as follows: xk = Ak−1 xk−1 + Bk−1 uk−1 + Dk−1 fk−1 y˜ k = Ck xk + ϵk
k = 1, . . . , m − 1,
x∗0 = x(0)
(11)
where, xk = x(kTs ) is the discrete state vector. The time varying matrices Ak , Bk , Dk and Ck can be found using the complete discretization method, assuming zero order hold approximation (Tóth, 2010): Ak =eA(kTs )Ts , Bk = A−1 (kTs )(Ak − I)B(kTs ) Dk =A−1 (kTs )(Ak − I), Ck = C
(12)
The above discretization method gives small discretization errors compared to other techniques as studied in Tóth (2010). Model (11) can be expanded and written as a matrix equation as follows: Y˜m = Xβ∗ + Tm Um + ζ Y where, X = [Ψ Θm ], β ∗ =
[
α∗ x∗0
]
, Ψ α∗ = Hm Fm
(13)
The above matrices and vectors are defined as follows:
⎡
⎤
⎡
C1 A 0 C2 A1 A0
⎤
⎡ ⎤ u0 ⎢ ⎥ ⎢ . ⎥ . ⎥ ⎥ , Um := ⎢ Y˜m := ⎣ .. ⎦ , Θm := ⎢ .. ⎣ .. ⎦ , ⎣ ⎦ . ∏m y˜ m um−1 Cm k=0 Am−k ⎡ ⎤ C1 B0 0 ··· 0 ⎢ ⎥ .. ⎢ ⎥ . C2 A1 B0 C2 B1 0 ⎢ ⎥, Tm := ⎢ ⎥ .. .. ⎣ ⎦ . . 0 ∏m−1 Cm k=0 Am−k B0 · · · · · · Cm Bm−1 ⎡ ⎤ C1 D 0 0 ··· 0 ⎢ ⎥ .. ⎢ ⎥ . C2 A1 D0 C2 D 1 0 ⎥, Hm := ⎢ ⎢ ⎥ .. .. ⎣ ⎦ . . 0 ∏m−1 Cm k=0 Am−k D0 · · · · · · Cm Dm−1 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ˜ T R0 f0 N ϵ0 ⎢ ⎥ ⎢ . ⎥ ⎢ . ⎥ .. Fm := ⎣ .. ⎦ , Ψ = Hm ⎣ ⎦ , ζ Y := ⎣ .. ⎦ . ˜ T Rm−1 ϵm fm−1 N y˜ 1
∑n˜
r We note here that Ψ ∈ Rm×stot , where, stot = l=1 sl represents the total number of reaction rate laws proposed for all hypothesized reactions. The above data equation can be viewed as a linear equation with a noisy observation vector Y˜m ; a regression matrix X; an unknown sparse + dense coefficient vector β∗ ; a bias term Tm Um and a measurement noise vector ζ Y . However, in our problem neither Tm Um nor X are known for certainty. Instead, a noisy version of Tm Um , denoted by T˜m U˜ m and a noisy version ˜ = [Ψ˜ Θ ˜ m ] can be formed using the noisy of X denoted by X ˜ k , y˜ k , v˜ k and q˜ in,k for k = 0, . . . , m − 1. Hence, the observations u following alternative data equation can be written:
˜ m = X˜ β∗ + ζ m Y where,
˜ m = Y˜m − T˜m U˜ m Y ˜ = X + Z˜ m X ζ m = (ζ Y + ζ U − Z˜ m β∗ )
(14)
Here, ζ U is a vector that represents the uncertainty in T˜m U˜ m and Z˜ m is an m × (stot + ns ) matrix that represents the uncertainty ˜ The estimation problem can now be phrased as follows: in X. ˜ m , the noisy regression given the vector of noisy measurements Y ˜ as defined above, find an estimate of the sparse + dense matrix X coefficient vector β∗ . The above problem can be classified as a sparse recovery estimation problem under regression matrix uncertainty in which estimators have been proposed in Loh and Wainwright (2011), Rosenbaum, Tsybakov, et al. (2010) and Zhu, Leus, and Giannakis (2011) that take into account perturbations in the regression ma˜ and the observation vector Y˜ m . However, these techniques trix X require the solution of non-convex optimization problems, not mentioning the additional a priori information required in these methods that can be difficult to obtain in practice. The estimator proposed in Rosenbaum et al. (2010) uses a convex program, however, the method requires precise knowledge about the noise present. Alternatively, it was found in Zhu et al. (2011) that using a Lasso operator to find β∗ can give relatively good results even for the case when the number of regression parameters is much higher than the number of observations. The same was also demonstrated in Rosenbaum et al. (2010) when Lasso with thresholding is used. Hence, in this study, we propose to make use of the following Lasso formulation:
ˆ σˆ 2 = β, arg min
β⪰0,σ 2
log(σ ) +
β = [αT , xT0 ]T
1 2mσ 2
∥α∥1 ∥Y˜ m − X˜ β∥22 + λ σ (15)
The above Lasso formulation was proposed in Städler, Bühlmann, and Geer (2010) using the negative log-likelihood estimation framework to estimate both the coefficient vector β and the variance σ 2 of ζ m assuming noise independence. The division of the ℓ1 norm penalty by σ was made to make the estimates equivariant under scaling of the measurements. This also helps to find a suitable value for λ approximately independent of the noise parameter σ as will be shown later. To transform the above non-convex problem to a convex one, the change of variables ρ = σ1 and φ = ρβ can be used (Städler et al., 2010). With the additional inclusion of positivity constraints, the final form of the estimator to be used and analyzed
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
˜ m = X˜ β∗ +ζ m and φ ∗ = ρβ Theorem 4.2. Suppose that Y ˆ ∗ , where:
in this study is the following:
φˆ βˆ = [αˆ T , xˆ T0 ]T = , ρˆ
where,
ˆ ρˆ = arg min −log(ρ ) + φ,
1
5
ρˆ = ∥ρ Y˜ m − X˜ φ∥22 + λwT φ
˜ Tm X˜ φˆ + Y
√
˜ Tm X˜ φˆ )2 + 4∥Y˜ m ∥2 m (Y 2 ˜ m ∥2 2∥Y 2
(18)
(16)
Suppose also that 4.1 is satisfied for some κ ∗ > 0 and that λ ≥ 1 ˜T ∥X ζ m ∥∞ , then the following error bound applies: m
Some additional remarks relating to (16) are worth mentioning:
(19)
φ⪰0,ρ
2m
w = [1Tstot 0T ]T
(1) The positivity constraints on φ were imposed since the kinetic parameters and the initial reactor concentration measurements x0 must be positive. (2) The third penalty is equivalent to the ℓ1 norm penalty in (15) since positivity constraints on φ are imposed. This penalty exploits solution sparsity in α and is used here for searching for the most plausible model structure among all hypothesized reactions and reaction rate laws. A method for selecting a suitable value for the tuning parameter λ will be introduced later. (3) Additional constraints can be imposed including, for example, relationships between kinetic parameters and upper and lower bounds if needed.
√ nr + ns ∥βˆ − β∗ ∥2 ≤ 3λ √ ρκ ˆ ∗
Proof. See Appendix A. Theorem 4.2 is applicable for both low dimensional settings, when m × ns ≫ (stot + ns ) and high dimensional settings; i.e. when (stot + ns ) ≫ m × ns . The result demonstrates that the error can be made small, for example, by (1) increasing the number of measurements and (2) increasing the value of κ ∗ . This can be done, for example, using experiment design that can reduce ˜ As a special case, we the correlations between the columns of X. present the following corollary. Corollary 1. Suppose that Assumption 4.1 holds and φ ∗ = ρβ ˆ ∗. Suppose also that the elements of ζ m can be approximated as zero mean i.i.d. noise with v ar(ζ m ) ≤ σ 2 , if:
√
4. Estimation performance and model selection
λ=2
stot + ns m
Cmax ρσ ˆ
(20)
then with probability 1/m the error bound becomes:
In the following, we are interested in analyzing the performance of the estimator in (16) and developing a model selection criterion. We will analyze the important assumptions needed for bounding the error and then analyze the factors that influence estimation performance with the objective of designing a suitable model selection criterion.
ρˆ (stot + ns )(nr + ns ) Cmax σ (21) mκ ∗ ˜ j ∥2 for j ∈ {1, . . . , stot + ns } and X˜ j represents where, Cmax = maxj ∥X ˜ the jth column of X.
4.1. Error bounds
Proof. See Appendix A.
Our approach for performance evaluation will be bounding ˆ − β∗ . Small error will indicate the rethe ℓ2 norm error of β covery of reactions associated with large kinetic parameters yet with possibly false inclusion or exclusion of reactions associated with small kinetic parameters. Our approach will make use of methods developed in Hastie, Tibshirani, and Wainwright (2015) and Wainwright (2009, 2019) for finding ℓ2 norm error bounds for unconstrained Lasso using the restricted eigenvalue assumption. We first present this key assumption. Assumption 4.1 (Restricted Eigenvalue Condition Hastie et al., ˜ satisfies the following: 2015). The noisy regression matrix X 1 ∥X˜ v∥22 ≥ κ ∗ ∥v∥22 , κ ∗ > 0, ∀v ∈ C , m where, C := {∥vSc∗ ∥1 ≤ 3∥vS ∗ ∥1 }
(17)
Here, vS ∗ represents the subvector of v with support S , which is the set of non-zero elements of the true coefficient vector φ ∗ . Similarly, vSc∗ represents the subvector of v with support index Sc∗ corresponding to the zero elements in φ ∗ . The above condition is called the restricted eigenvalue condition since without the constraint on v , κ will represent the minimum eigenvalue of the ˜ T X. ˜ The above assumption can be satisfied symmetric matrix X ˜ is rank deficient for some special cases as discussed even when X in Wainwright (2019). We now present the theorem for bounding ˆ. the error on β ∗
√
∥βˆ − β ∥2 ≤ 6 ∗
4.2. Model selection and parameter identification As evident from Theorem 4.2, the selection of λ can influence the estimation error and consequently influence both the resulting model structure and estimated kinetic parameters. In light of Theorem 4.2, we desire to find a technique that repeatedly solves (16) for different values of λ in order to find the most plausible reaction model structure and parameters. Using first ˆ ≻ 0 when: order optimality conditions, we can easily show that φ
λ < λmax = max √ j
˜ T Y˜ m X j ˜ m ∥2 m∥Y
,
j = 1, . . . , stot
(22)
The technique developed for finding λ is shown as a flow chart in Fig. 2 and is explained next. The maximum number of iterations imax is first specified and a sequence of values λ1 = 0 < λ2 < · · · , < λimax = λmax is generated on a logarithmic scale to have more iterations near zero as follows:
λi = log10 10(iλmax )/imax , i = 2, 3, . . . , imax − 1
(23)
Consequently, at each iteration i, the minimization problem (16) is solved using λi and the ith active support index is found ˆ j ≥ ϵmin }, where ϵmin ≪ 1 represents a user as S i = {j|φ
ˆ j . Next, defined threshold for considering non-zero elements of φ Assumption 4.1 is verified by finding a predicted value for κ ∗ assuming the candidate active support index S i by solving the following optimization problem: κˆ i := min v∈Ci
1 ∥X v∥22 m ∥v∥22
6
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
more input excitation are required. Generally, more measurements can help to satisfy Assumption 4.1 and reduce the estimation errors as implied in (19), while more input excitation can help to increase the value of κ ∗ . We note here that imax should be selected large enough to properly cover the candidate model structures within the solution path. In this study, trial and error was used to select imax . 5. Simulation study In the following, we demonstrate our results using the same simulation model used in Brendel et al. (2006) and in Bhatt et al. (2012) for comparison. The example involved a homogeneous reaction system for the acetoacetylation of pyrrole with diketene. The reaction stoichiometry is given by: K
P + D −→ PAA K D + D −→ DHA K D −→ OL K PAA + D −→ G where, (P) represents pyrrole, (D) diketene, (PAA) 2-acetoacetyl pyrrole, (DHA) dehydroacetic acid, (OL) oligomers and (G) is a byproduct. The corresponding true stoichiometric matrix using the same order is given by:
⎡ −1 ⎢−2 N=⎣ −1 −1
−1
+1
0 0 0
0 0 −1
0 +1 0 0
0 0 +1 0
⎤
0 0 ⎥ 0 ⎦ +1
and the true reaction rates are described by: Fig. 2. Algorithm for finding a suitable value for λ and final model identification.
where, Ci := {∥vS i ∥1 ≤ 3∥vS i ∥1 }
(24)
c
The above minimization problem can be solved using, for example, sequential quadratic programming with several iterations until convergence. This will provide an estimate for κ for the sake of verifying Assumption 4.1 at iteration i and for calculating a threshold value. The threshold value kmin is found based on the error bound given in (19) as follows:
√
n˜ r + ns kmin = 3λi √
(25)
ρˆ κˆ i i
This threshold value is then used to find the active support Sˆ ∗ ˆ j ≥ kmin }. Based on this active support, the defined as Sˆ ∗ = {j|β solution of the following minimization problem is found: ∗ 1 φˆ Sˆ ∗ , ρˆ ∗ = arg min −log(ρ ) + ∥ρ Y˜ m − X˜ Sˆ ∗ φSˆ ∗ ∥22 φ⪰0,ρ
(26)
2m
∗
ˆ Sˆ ∗ is augmented with zeros according to Sˆc∗ to form The solution φ ∗ φˆ . This solution is then tested for feasibility by testing whether ˆ ∗ , . . . kˆ ∗ ]T = the identified kinetic parameter vector, given by [k 1 n˜ r
ˆ ρˆ ∗ , contains at most one kinetic parameter per reaction wT φ/ (i.e. one reaction rate law per reaction). The first model in the solution path that passes the above feasibility test is flagged as the most plausible model based on the criteria to select the model with the largest number of feasible reactions. The solution path can be further investigated to determine if the reaction process can be modeled with fewer number of reactions, if desired. If on the other hand, no model was found that satisfies the above selection criterion, than more experimental data and/or
ra (t) =ka cP (t)cD (t)cK (t), rb (t) = kb cD (t)2 cK (t) rc (t) =kc cD (t),
rd (t) = kd cPAA (t)cD (t)cK (t)
with kinetic parameters given by ka = 0.0523, kb = 0.1279, kc = 0.0281 and kd = 0.1. As in Brendel et al. (2006), the catalyst is assumed to deactivate according to the formula: cK (t) =
v0 cK v (t) 0
where cK0 is the initial concentration of catalyst in the reactor. We now present three different identification scenarios to demonstrate the potential of the technique. 5.1. Scenario I: kd = 0 For scenario 1 of the present simulation study the input parameters used in this experiment are identical to the simulation scenario given in Bhatt et al. (2012) and Brendel et al. (2006) for comparison, presented again here. A set of 32 reaction rate law candidates are assumed for the four reactions described earlier, corresponding to feasible power-law kinetic functions g˜i,j (c(t), t) shown in Table 1. The rate constant of the fourth reaction was set to zero; i.e. kd = 0 which makes the fourth reaction not occurring to test whether the identification method can also determine the correct reaction stoichiometry. Simulation of the true reactor model with TS ∗ = 10 s, m = 361 and tf = 60min was first done with the following initial conditions: v (0) = 0.5 L, cD (0) = 0.14 mol/L, cP (0) = 0.3 mol/L, cPAA (0) = 0.08 mol/L, cDHA (0) = 0.01 mol/L, cOL (0) = 0.01 mol/L, cG (0) = 0.01 mol/L and cK 0 = 2 mol/L. Furthermore, the inlet feed was set to qin (kTS ∗ ) = 5 ml/min and cDin (kTS ∗ ) = 6 mol/L, for k = 1, . . . , m while all other species concentrations set to zero. The true concentration measurement vector was formed as: yk = [cD (tk ), cP (tk ), cPAA (tk ), cDHA (tk ), cOL (tk ), cG (tk )]T
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823 Table 1 Hypothesized reaction rate laws. Reaction a: K
Reaction b:
P + D −→ PAA
K
D + D −→ DHA
Reaction c:
Reaction d:
D −→ OL
PAA + D −→ G
K
g˜a,1 = 1
g˜b,1 = 1
g˜c ,1 = 1
g˜d,1 = 1
g˜a,2 = cD
g˜b,2 = cD
g˜c ,2 = cD
g˜d,1 = cD
g˜a,3 = cP
g˜b,3 = cD2
g˜c ,3 = cD2
g˜d,3 = cPAA
g˜a,4 = cK
g˜b,4 = cD cK
g˜c ,4 = cD cK
g˜d,4 = cK
g˜a,5 = cD cP
g˜b,5 = cD2 cK
g˜c ,5 = cD2 cK
g˜d,5 = cD cPAA
g˜a,6 = cP cK
g˜b,6 = cK
g˜c ,6 = cK
g˜a,7 = cD cK
g˜d,6 = cPAA cK g˜d,7 = cD cK
g˜a,8 = cD cP cK
g˜d,7 = cD cPAA cK
g˜a,9 = cD cP2
2 g˜d,9 = cD cPAA
g˜a,10 = cD2 cP
g˜d,10 = cD2 cPAA
where, tk = kTS ∗ . Consequently, noisy concentration measurements were then formed as: y˜ k = yk + ϵk
ϵk = σ · ηk ,
σ = 0.01 × [0.38, 0.45, 0.45, 0.63, 0.52, 0.05]T
and ηk ∈ R6 is random vector sequence generated using Matlab’s normally distributed pseudorandom function (randn). Following our proposed algorithm in Fig. 2, the dictionary Ψ˜ ∈ R361×32 was first formed using simulated noisy concentration measurements for all measured species y˜ k and using the candidate reaction rate laws shown in Table 1. Afterwards, the ˜ m and matrices Θ ˜ m and X˜ were formed according to their vector Y definitions in (14). The maximum tuning parameter λmax was calculated using (22) and was found to be 42. A set of imax = 20 tuning parameters was then generated using (23). The thresholds ϵmin = 1 × 10−6 was used. Consequently, the algorithm in Fig. 2 was implemented using Matlab CVX (Grant & Boyd, 2012) over the solution path starting with λ1 = 1 × 10−6 . The solution time for each iteration in the solution path was found to be approximately 1.5 s using a Macbook pro with 2.7 GHz Intel Core i5 and 8 GB RAM. In every iteration i we found an estimate for κi by solving (24) using Matlab’s ‘‘fmincon’’ function with the option of sequential quadratic programming solver and maximum of 500 iterations. The results for the identification step for each iteration of the algorithm in 2 for i = 1, 2, 3, 4, 5 are shown in Fig. 3. The values of kˆ ∗i for each iteration are plotted as a stem plot. As shown in the figure, after 4 iterations, it was found that κˆ 4 = 7 × 10−5 and associated threshold was kmin = 0.021. At this iteration the solution first became feasible. The associated estimation ℓ2 norm ˆ S ∗ − βˆ ∗∗ ∥2 = 1 × 10−4 , shown in the error was found to be ∥β S last plot in the same figure, demonstrating near recovery of the true parameters. Fig. 4 shows the true, measured and estimated reactor species concentrations for this experiment. In comparison with the method in Brendel et al. (2006), the method gave a corresponding error of 7 × 10−4 . 5.2. Scenario II: kd ̸ = 0 Reference to scenario I explained above, it was mentioned in Bhatt et al. (2012) that the fourth reaction was excluded in the analysis due to lack of structural identifiability. We simulated the exact case described above when the fourth reaction was included by setting kd = 0.1 to test this scenario using the proposed technique in this study. The other simulation details are exactly similar to the details discussed in scenario I above except that the inlet concentration for species D was set randomly as cDin (kTS ∗ ) = 6+ηk mol/L, where ηk ∈ R is a zero mean normally
7
distributed pseudorandom number with variance of 0.01 generated in Matlab in order to introduce random input excitation. We have also introduced zero mean random noise with variance of 0.001 in the reactor volume flow rate measurements q˜ in,k to make the problem more difficult to estimate. The algorithm in Fig. 2 was then implemented as before using the same settings mentioned in scenario I above. After 5 iterations, a feasible model was found, see Fig. 5. The associated estimation errors were found and shown in the bottom right plot demonstrating near exact recovery of the true parameters. The example demonstrates that the new technique is capable of estimating a previously described unidentifiable problem by introducing more reactor inlet excitation. 5.3. Scenario III: Subset of reactor species measured In this simulation experiment, we would like to examine the case when a subset of reactor concentration measurements is available. Using the same homogeneous reaction system example explained earlier in scenario 1, we tested the case when only concentration measurements for species D, P , PAA and OL are available (i.e. cDHA and cG are not measured). In this case, the dictionary matrix Ψ˜ will be similar to scenario 1 since the hypothesized reaction rate laws depend on the first three measured species only. Consequently, matrix C in (7) was redefined as follows:
⎡
1 ⎢0 C =⎣ 0 0
0 1 0 0
0 0 1 0
0 0 0 0
0 0 0 1
⎤
0 0⎥ 0⎦ 0
To make this simulation experiment even more difficult to identify, we have introduced random white noise in the volume measurements v˜ k and the volume flow rate measurements q˜ in,k using zero mean white noise with variance of 0.001. The other simulation experiment conditions (model parameters, inlet conditions, initial condition, time step and simulation time) are identical to the conditions described earlier for scenario 1 except that λ1 was set to 1 × 10−5 . The plot in Fig. 6 shows the result of implementing the algorithm in Fig. 2 for the fourth iteration. The figure demonstrates again near exact recovery of reaction stoichiometry, the correct reaction rate laws and near exact recovery of the kinetic parameters despite not measuring two of the reactor species and despite having noisy reactor volume and flow rate measurements. 6. Conclusion A new identification framework for the identification of reaction stoichiometries and kinetic parameters for homogeneous reaction systems is presented. The identification problem makes use of hypothesized reaction stochiometries and hypothesized reaction rate laws. The new identification framework is demonstrated for a CSTR process. A dictionary of plausible reaction rate trajectories is formed from noisy concentration measurements. Subsequently, a special non-negative Lasso formulation that is equi-variant under scaling was used to determine the correct model structure among all possible model candidates and the associated kinetic parameters simultaneously. Estimation error bounds were found and are applicable for both the low and high dimensional settings. A simple solution path searching technique to find the most plausible reaction model structure and parameters in view of the theoretical treatment was also introduced. The methodology was illustrated on a simple example used in the literature demonstrating exact recovery of reaction stoichiometry, model structure and near recovery of model parameters using noisy measurements under three different scenarios. In comparison with previously studied identification frameworks, we mention the following:
8
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
Fig. 3. Identified vs. true kinetic parameters for simulation scenario 1 (when kd = 0). Each plot shows one iteration of the algorithm in Fig. 2.
Fig. 4. True, measured and estimated reactor species concentrations for simulation scenario #1.
(1) The new technique is capable of recovering the reaction
true reaction stoichiometry from a postulated stoichiomet-
stoichiometry, the associated reaction rate laws and near
ric matrix that may include extraneous reactions. Also, the
exact recovery of the kinetic parameters simultaneously
new technique does not require the estimation of reac-
using one convex optimization problem. Hence, techniques
tion fluxes, extents of reactions nor reaction rates prior to
for elucidating the correct reaction stoichiometry are not
identification. Direct identification is achieved from exper-
needed because the method can concurrently recover the
imental data which minimizes errors.
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
9
Fig. 5. Identified vs. true kinetic parameters for simulation scenario 2 (when kd = 0.1). Each plot shows one iteration of the algorithm in Fig. 2.
Extensions to the technique are needed for other types of reactors, including, plug flow reactors (Bermúdez, Esteban, Ferrín, Rodríguez-Calo, & Sillero-Denamiel, 2017) and distributed reaction systems (Rodrigues, Billeter, & Bonvin, 2017). Acknowledgments The authors would like to thank the anonymous reviewers for their valuable comments that helped to improve this work. Appendix A. Proof of Theorem 4.2 From optimality of the solution of (16) we may write:
− log(ρˆ ) + Fig. 6. Kinetic parameters identified vs. true kinetic parameters for simulation scenario 3.
(2) Since the new technique requires the solution of one convex optimization problem, the method has the potential for identifying large kinetic models that involve large number of reactor species; experimental data and hypothesized reaction rate laws and stochiometries compared to previous techniques. (3) The solution path, determined by the iterative technique shown in Fig. 2, can help to determine if the reaction process can be modeled with fewer number of reactions, if desired. The method can also be extended to other cases, for example when density are and/or temperature variations are present. These cases require additional measurements and some modifications to the mole balance and reaction rate equations used.
≤ − log(ρˆ ) +
1 2m 1 2m
ˆ 22 + λwT φˆ ∥ρˆ Y˜ m − X˜ φ∥ ∥ρˆ Y˜ m − X˜ φ ∗ ∥22 + λwT φ ∗
Defining vˆ := φˆ − φ ∗ and expanding the above inequality we obtain the following: 1 { 2m
} ∥X˜ vˆ ∥22 − 2ρˆ Y˜ Tm X˜ vˆ + 2(φ ∗ )T X˜ T X˜ vˆ + λwT vˆ ≤ 0
˜ m = X˜ β ∗ + ζ m and the substitution Using the data equation Y φ ∗ = ρ ∗ β∗ , then after rearranging we obtain: 1 2m
∥X˜ vˆ ∥22 ≤
ρˆ
(
m
ζm +
ρˆ − ρ ∗ ˜ ∗ Xφ mρ ∗
)T
˜ vˆ − λwT vˆ X
Using Holder’s inequality, we can write: 1 2m
∥X˜ vˆ ∥22 ≤
1 m
∥X˜ T ζ φ ∥∞ ∥ˆv ∥1 − λwT vˆ
where,
ζ φ = ρζ ˆ m+
ρˆ − ρ ∗ ˜ ∗ Xφ ρ∗
(A.1)
10
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823
Hence, if ρ ∗ ≈ ρˆ , then ζ φ = ρζ ˆ m . Now considering that ∥ˆv ∥1 = ∥ˆvS ∗ ∥1 + ∥ˆvSc∗ ∥1 and assuming:
λ
2 m the right hand side of (A.1) can be written as:
(A.2)
≤
3λ
(∥ˆvS ∗ ∥1 + ∥ˆvSc∗ ∥1 ) + λ(∥ˆvS ∗ ∥1 − ∥ˆvSc∗ ∥1 )
2
∥ˆvS ∗ ∥1 ≤
3λ √ 2
nr + ns ∥ˆv ∥2
∥ˆvSc∗ ∥1 ≤ 3∥ˆvS ∗ ∥1 and hence vˆ ∈ C . Consequently, following Assumption 4.1 we obtain:
√
nr + ns
κ∗
Finally, assuming that ρ ∗ ≈ ρˆ , we obtain the error bound in (19). For the expression for ρˆ , we first note that the Lagrangian of (16) can be written as (Boyd & Vandenberghe, 2004): 1 ∥ρ Y˜ m − X˜ φ∥22 + (λwT − µT )φ 2m where, µ is a Lagrange multiplier corresponding to the positivity constraint on φ . Now, ρ = ρˆ become a minimizer of (16) if the following KKT condition is satisfied (Boyd & Vandenberghe, 2004): L(φ, ρ, µ) = − log(ρ ) +
∇ρ L = −
1
ρˆ
+
ρˆ ˜ 2 1 ∥Ym ∥2 − Y˜ Tm X˜ φˆ = 0
m
m
(A.3)
Accordingly, (A.3) is a quadratic equation in ρˆ which has a positive solution given by (18). Appendix B. Proof of Corollary 1 We desire to select λ so that (A.2) is satisfied with high probability. Since ρ ∗ ≈ ρˆ is assumed then ζ φ = ρζ ˆ m . First, we define z as follows:
ρˆ
˜ T ζm X m Consider an arbitrary j ∈ {1, . . . , stot + ns }, then the jth element of z (denoted by zj ) can be written as: z :=
zj =
ρˆ ˜ T Xj ζ m
m
˜ T denotes the jth row of X˜ T . Assuming that ζ m can be where X j approximated as zero mean, white noise with variance σ 2 , we can say that: var(zj ) ≤
ρˆ 2 σ 2 ˜ 2 ∥Xj ∥2 2
m Consequently, using Chebyshev’s inequality we have:
P{|zj | ≥ t } ≤
ρˆ 2 σ 2 ˜ 2 ∥X j ∥2 2 2
m t Since ∥z ∥∞ corresponds to the maximum over stot + ns variables, we have:
P{∥z ∥∞ ≥ t } ≤ (stot + ns )
Cmax ρσ ˆ
and λ = 2t we have the following tail bound:
λ
}≤
1
References
where the first inequality is true because vˆ Sc∗ = φˆ Sc∗ ≥ 0. Now 1 since 2m ∥X˜ vˆ ∥22 ≥ 0, from the first inequality above we have:
∥φˆ − φ∗ ∥2 ≤ 3λ
m
2 m and the result follows.
m
2
stot + ns
P{∥z ∥∞ ≥
ρˆ ˜ T ∥X ζ m ∥∞ ∥ˆv ∥1 − λwT vˆ λ
√ t=
ρˆ ≥ ∥X˜ T ζ m ∥∞
≤
˜ j ∥2 for j ∈ {1, . . . , stot + ns }. Setting t as: where, Cmax = maxj ∥X
ρˆ 2 σ 2 m2 t 2
2 Cmax
Bermúdez, Alfredo., Esteban, N., Ferrín, J. L., Rodríguez-Calo, J. F., & SilleroDenamiel, M. R. (2017). Identification problem in plug-flow chemical reactors using the adjoint method. Computers & Chemical Engineering, 98, 80–88. Bhatt, Nirav, Kerimoglu, Nimet, Amrhein, Michael, Marquardt, Wolfgang, & Bonvin, Dominique (2012). Incremental identification of reaction systems: A comparison between rate-based and extent based approaches. Chemical Engineering Science, 83, 24–38. Bonvin, D., & Rippin, D. W. T. (1990). Target factor analysis for the identification of stoichiometric models. Chemical Engineering Science, 45(12), 3417–3426. Boyd, Stephen Poythress, & Vandenberghe, Lieven (2004). Convex optimization. Cambridge University Press. Brendel, Marc, Bonvin, Dominique, & Marquardt, Wolfgang (2006). Incremental identification of kinetic models for homogeneous reaction systems. Chemical Engineering Science, 61(16), 5404–5420. Grant, Michael, & Boyd, Stephen (2012). CVX: Matlab software for disciplined convex programming, version 2.0 beta. http://cvxr.com/cvx. Hastie, Trevor, Tibshirani, Robert, & Wainwright, Martin (2015). Statistical learning with sparsity: The lasso and generalizations. CRC Press. Kugler, Philipp, Gaubitzer, Erwin, & Muller, Stefan (2009). Parameter identification for chemical reaction systems using sparsity enforcing regularization: A case study for the chlorite iodide reaction. The Journal of Physical Chemistry A, 113(12), 2775–2785. Loh, Po-Ling, & Wainwright, Martin J. (2011). High-dimensional regression with noisy and missing data: Provable guarantees with non-convexity. In Advances in neural information processing systems (pp. 2726–2734). Maria, G. (2004). A review of algorithms and trends in kinetic model identification for chemical and biochemical systems. Chemical and Biochemical Engineering Quarterly, 18(3), 195–222. Papachristodoulou, Antonis, & Recht, Ben (2007). Determining interconnections in chemical reaction networks. In Proc. American control conference (pp. 4872–4877). Rawlings, James Blake, & Ekerdt, John G. (2002). Chemical reactor analysis and design fundamentals. Nob Hill Pub, Llc. Rice, F. O., Johnston, W. R., & Evering, B. L. (1932). The thermal decomposition of organic compounds from the standpoint of free radicals. ii. experimental evidence of the decomposition of organic compounds into free radicals. Journal of the American Chemical Society, 54(9), 3529–3543. Rodrigues, Diogo, Billeter, Julien, & Bonvin, Dominique (2015). Incremental model identification of distributed two-phase reaction systems. IFAC-PapersOnLine, 48(8), 266–271. Rodrigues, Diogo, Billeter, Julien, & Bonvin, Dominique (2017). Generalization of the concept of extents to distributed reaction systems. Chemical Engineering Science. Rosenbaum, Mathieu, Tsybakov, Alexandre B., et al. (2010). Sparse recovery under matrix uncertainty. The Annals of Statistics, 38(5), 2620–2651. Santosa, Fadil, & Weitz, Benjamin (2011). An inverse problem in reaction kinetics. Journal of Mathematical Chemistry, 49(8), 1507–1520. Städler, Nicolas, Bühlmann, Peter, & Geer, Sara Van De (2010). ℓ 1-penalization for mixture regression models. Test, 19(2), 209–256. Tóth, Roland (2010). Modeling and identification of linear parameter-varying systems, vol. 403. Springer. Wainwright, Martin J. (2009). Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1 -constrained quadratic programming (lasso). IEEE Transactions on Information Theory, 55(5), 2183–2202. Wainwright, Martin J. (2019). High-dimensional statistics: A non-asymptotic viewpoint, vol. 48. Cambridge University Press. Willis, Mark J., & von Stosch, Moritz (2017). Simultaneous parameter identification and discrimination of the nonparametric structure of hybrid semi-parametric models. Computers & Chemical Engineering, 104, 366–376. Zhu, Hao, Leus, Geert, & Giannakis, Georgios B. (2011). Sparsity-cognizant total least-squares for perturbed compressive sampling. IEEE Transactions on Signal Processing, 59(5), 2002–2016.
A. Al-Matouq and T. Vincent / Automatica 114 (2020) 108823 Ali Al-Matouq received the B.Sc. and M.Sc. degrees with distinction in electrical engineering from King Fahd University of Petroleum and Minerals, Dhahran, in 2002 and 2007, respectively and an M.Sc. degree in Chemical Engineering and a Ph.D. degree in Electrical Engineering from Colorado School of Mines in 2012 and 2014, respectively. He also has 5 years of professional experiance as a process control and automation engineer with an ISA certification. He is currently an Assistant Professor in the Engineering Management department of Prince Sultan University, Riyadh. His research interests include process control and optimization with applications in chemical and biomedical engineering.
11
Tyrone Vincent received the B.S. degree in electrical engineering from the University of Arizona, Tucson, in 1992, and the M.S. and Ph.D. degrees in electrical engineering from the University of Michigan, Ann Arbor, in 1994 and 1997, respectively. He is currently a Professor in the Department of Electrical Engineering and Computer Science, at the Colorado School of Mines, Golden. His research interests include system identification, estimation, and fault detection with applications in materials processing, and energy systems.