A Unified Framework for Kinetic Parameter Estimation Based on

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Process Systems Engineering

A Unified Framework for Kinetic Parameter Estimation Based on Spectroscopic Data w/ or w/o Unwanted Contributions Weifeng Chen, Lorenz T. Biegler, and Salvador Garcia-Munoz Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b05273 • Publication Date (Web): 13 Feb 2019 Downloaded from http://pubs.acs.org on February 14, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A Unified Framework for Kinetic Parameter Estimation Based on Spectroscopic Data with or without Unwanted Contributions Weifeng Chen a, Lorenz T. Biegler b* and Salvador García Muñoz c a

b

College of Information Engineering, Zhejiang University of Technology, Hangzhou, 310023, P.R. China Center for Advanced Process Decision-making, Department of Chemical Engineering, Carnegie Mellon University, 5000 Forbes Ave. Pittsburgh, PA 15213, USA. c Small Molecule Design and Development, Lilly Research Labs, Indianapolis, IN 46285, USA

Kinetic estimation of chemical reactions is considered to be a crucial step prior to design a robust, controllable and safe production process. Spectroscopic measurements are widely used for kinetic analysis. However, there may be unwanted contributions in the measured spectra, which may come from the instrumental variations (such as the baseline shift or distortion) or from the presence of inert absorbing interferences with no kinetic behavior. Here, the kinetic estimation problem with time invariant contributions is studied in depth, and we derive conditions where the estimation accuracy is not affected by the time invariant contributions. Moreover, kinetic parameter estimation and separation of time invariant contributions can be performed simultaneously under proper conditions. Also, an approach for kinetic parameter estimation based on spectra with time variant contributions is proposed. Finally, a novel unified framework is developed for kinetic parameter estimation when there is no prior information for unwanted contributions.

1. Introduction Spectroscopic measurements are widely used for chemical reaction monitoring1-4. Enormous amounts of spectroscopic data could be gathered for the subsequent analysis. This obtained information describes the time evolution of different chemical components present in the reaction system simultaneously. The traditional strategy for kinetic parameter estimation is to recover concentration profiles from spectroscopic data by means of soft-modeling approaches5-8, such as multivariate curve resolution (MCR), followed by fitting to kinetic model in a separate step. This usually does not work well because of the problem of rotational ambiguity associated with the obtained profiles9-10. Although rotational ambiguity can be considerably reduced by means of natural constraints such as non-negativity, unimodality or balance-closure equations for concentration and absorption spectra profiles, the accuracy of estimated kinetic parameters is still severely affected. Another widely applied approach is Gauss-Newton-Levenberg-Marquardt (GNLM)11-12 based hard-modeling method. Several studies13-18 claim that hard-modeling type methods are more appropriate than soft-modeling approaches in chemical kinetics. In the course of this approach, the kinetic parameters are refined to minimize

*Correspondence

to: L. T. Biegler, E-mail: [email protected]

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the sum of squares of the residuals between the measured and modeled spectroscopic data. Once an initial guess for kinetic parameters is given, concentration profiles can be obtained by integrating the system of ordinary dynamic equations which represents the reaction mechanism. The pure component spectra can be regarded as linear parameters and calculated as the pseudo-inverse of concentration matrix. GNLM then updates the kinetic parameters (nonlinear parameters) after evaluating the fitting error of the spectroscopic data, and corresponding derivative information with respect to kinetic parameters. The latter are obtained by finite difference or sensitivity analysis. The drawback of hard-modeling methods is that they are poorly suited to unstable and ill-conditioned dynamic systems, which can lead to unbounded state profiles and convergence difficulties19. Instead of using the kinetic model as a hard constraint, a hybrid approach was developed20, which combines model-free multivariate curve resolution with model-based kinetic regularization. A multi-resolution approach for accelerating the convergence of multivariate curve resolution methods was subsequently proposed in Sawall et al.21, based on the idea that good initial values result in a fast solution at the next refined level. One disadvantage of the proposed approach is that singular value decomposition (SVD) or principal component analysis (PCA) of the rotation matrix is used to approximate the measured spectroscopic data. With this approximation, the residuals may not be equal to that of MCR with natural constraints. Moreover, direct rotation of PCA solutions might not give solutions that satisfy required constraints, like non-negativity22. In addition, it is difficult to determine the appropriate weight factors corresponding to regularization terms in the objective function; this impacts the accuracy of kinetic parameter estimation. Vaia and Sahinidis23 developed a global mixed-integer nonlinear programming algorithm for addressing the problem of simultaneous model structure determination, covariance estimation and parameter estimation in infrared spectroscopy. Nevertheless, this approach can be quite expensive even for medium-scale spectroscopic data sets. Finally, we proposed an approach based on maximum likelihood principles for simultaneous estimation of reaction kinetics and curve resolution24. However, variances of the noise in the system variables and spectral measurements must be estimated beforehand and the iterative optimization-based variance estimation approach may be time-consuming for large-scale cases. All of these proposed approaches for kinetic parameter estimation are based on the assumption that there are no unwanted contributions in the instrumental response or the measured spectroscopic data are well preprocessed. Many preprocessing strategies including multiplicative scatter correction, extended standard normal variate, detrending, Norris-Williams derivation and Savitzky-Golay derivation, and their combinations25 are commonly used. However, selection of the preprocessing method depends on the experience of the user, and it is possible to use the wrong preprocessing method for a given situation. Another approach proposed by de Juan et al.26 can handle the unwanted contributions in the instrumental response without preprocessing. This approach combines aspects of hard-modeling and soft-modeling methods, by using a multivariate curve resolution-alternating least squares method (MCR-ALS) along with the kinetic model treated as a hard constraint to fit the concentration data. Rotational ambiguity can be significantly reduced by the introduction of the kinetic model with reasonable

2


Page 2 of 27

Page 3 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


kinetic estimation also provided. This combined hard-model MCR-ALS was successfully applied to many processes including enzymatic reactions of different mixtures of hypoxanthine, xanthine and uric acid with xanthine oxidase, photochemical degradation process of decabromodiphenyl ether, and time-dependent kinetic processes in colloidal systems27-29. However, choosing suitable soft-modelling constraints for the concentration and pure spectra matrix still depends on experience. In addition, without derivative information and explicit search directions towards minimizing the fitting error, the combined hard-model MCR-ALS converges very slowly and the precision of estimated kinetic parameters may be insufficient. Finally, a dynamic optimization problem of fitting the concentration data needs to be solved at each iteration, which is time-consuming. In this work, we focus on the measured spectroscopic data with unwanted time invariant contributions (e.g. baseline shift, distortion and the presence of inert absorbing interferences with no kinetic behavior) and unwanted time variant contributions (e.g. drift). A novel simultaneous approach is proposed for kinetic parameter estimation based on measured spectroscopic data with unwanted contributions, and a unified framework for handing spectroscopic data with or without unwanted contributions is established.

2. Background Theory 2.1 Beer-Lambert’s Law Beer-Lambert's law predicts the absorbance of a solution d(t, λ) at a particular wavelength λ and given time t is the sum over all contributions of dissolved absorbing species, i.e.,

d (t ,  )  c1 (t ) s1 ( )  c2 (t ) s2 ( )  ...  cnc (t ) snc ( )

(1)

where ck(t) is the concentration for k-th species at time t, sk(λ) is the absorbance for k-th species with one-unit concentration at wavelength λ, and nc is the number of absorbing species. Considering the measurement error, Beer-Lambert's law can be expressed with the following model,

D =CS T  E

(2)

where D is the ntp by nwv spectra matrix, ntp is the number of sampling points in time and nwv is the number of variables for wavelength dimension. C is the ntp by nc concentration matrix, the row is sampling time and the column is the absorbing species. S is the nwv by nc absorbance matrix, the row is wavelength and the column is the absorbing species. E is the ntp by nwv measurement error matrix and each element of E is assumed to be mutually independent and follow the Gaussian distribution N(0, δ2), where δ2 is the variance. 2.2 Simultaneous Approach for Kinetic Parameter Estimation The reaction process is modeled as system of ordinary differential equations (ODE), 3



Page 4 of 27

dz (t )  f ( z (t ), θ ) dt

(3)

where z(t) (nz×1) is the differential state variable vector (nz is the number of differential state variables), the kinetic parameter vector θ (ntheta×1) is defined as [θ1, θ2, ..., θntheta]T and f is a mapping from Rnz+ntheta to Rnz. We discretize the ODE (3) with Radau collocation, which is algebraically stable [30-31], and choose K+1 interpolation points in finite element j to represent the state in this element [19] as K

K

z (t )   m ( ) z j ,m , t  tp j 1  h j ,  [0,1], m ( )  m0



m'  0,  m

(   m' )

( m   m' )

(4)

where tpj-1 and tpj represent the two endpoints of finite element j. zj,m is the value of z(t) at the j-th finite element [tpj-1 , tpj] and the m-th collocation point  m . K is the number of collocation points in a finite element. This leads to the following collocation equations, K



 m ( a ) z j ,m  h j  f ( z j ,a , θ )  0, j  1...N , a  1..K

(5)

m0

where N is the number of finite elements. We assume that each component i,l of E is mutually independent and has the same variance,

 i ,l  N (0,  2 ) , i=1…ntp, l=1…nwv

(6)

and then we have ntp nwv

p ( D  CS )   p ( i ,l ) T

(7)

i 1 l 1

where p(.) is the probability calculation operator. Under the assumption that all the absorbing species are known and the corresponding concentration profiles are independent, we apply the maximum likelihood estimation approach to Eq. (7), and obtain the following optimization problem, nc   min   di ,l   ck (ti ) sk (l )  i 1 l 1  k 1  ntp nwv

s.t.

K

2



 m ( a ) z j ,m  h j  f ( z j ,a , θ )  0, j  1..N , a  1..K

m 0

z j 1, K  z j ,0 , j  2..N K

ck (ti )   m ( i ) zn' , j ,m , i  (ti  tp j 1 ) (tp j  tp j 1 ) , n '  v(k ) m 0

i  1,.., ntp, l  1,.., nwv, k  1,.., nc

4


(8)

Page 5 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where di,l is the element of D in the i-th row and l-th column. zn(ti) is the n-th element (n=1,..,nz) of z(ti) and zn,j,m is the n-th element of zj,m. v(k) is a known mapping from the k-th absorbing species to the n’-th state variable in the ODE (3). The decision variables in this formulation are ck(ti), sk(λl), zn,j,m, and θ.

3. Kinetic Parameter Estimation with Unwanted Contributions in Spectroscopy The chemical reaction mechanism of isochoric process can be uniquely defined by the reactant coefficient matrix Re (nr×nz) and the product coefficient matrix Pr (nr×nz), where nr is the number of reactions and nz is the total number of reactants and products. The stoichiometric coefficient matrix St (nr×nz) can be calculated as Pr-Re. As in Ref. 32, we define a time invariant pseudo-equivalence matrix Φ as

 St       z0T  ZT   in  where z0 (nz×1) is the initial concentration vector, Zin (nz×nf) is the dosing concentration matrix with nf dosing times (t1, t2, …, tnf). As shown in Ref. 32, Φ has the same kernel as the (calculated) concentration matrix Z (ntp×nz). Hence, we can determine the rank deficiency of Z and linear dependent relationships among the columns of Z by analyzing the matrix Φ. According to the mapping relationship in problem (8), C is composed of the v(1)-th, v(2)-th, …, v(nc)-th columns of Z. Let Φsub, Stsub, z0,subT, and Zin,subT be composed of the v(1)-th, v(2)-th, …, v(nc)-th columns of Φ, St, z0T, and ZinT. Based on the theorem in References 32 and 33, we can deduce the rank deficiency of C and linearly dependent relationships among the columns of C can be determined by analyzing the matrix Φsub. We define the rank of kernel of Φsub as rkp. 3.1 Independent Absorbing Species In the following, we consider the case where rkp = 0, i.e. the concentration profiles of all the absorbing species with kinetic behavior are linearly independent. The unwanted contributions can be divided into the time invariant and the time variant instrumental variations. Beer-Lambert's law is modified as,

D =CS T  G  E

(9)

a) Time Invariant Instrumental Variations The time invariant instrumental variations include baseline shift, distortion and presence of inert absorbing interferences without kinetic behavior. In this work, the estimation problem with unwanted contributions is studied in depth. The time invariant instrumental variations can be uniformly described through the addition of a rank-one matrix G, given by the outer product:

5



G = eg T , e = [1 1 ... 1]T ( ntp1) , g =  g1 g 2 ... g nwv  , G  R ntpnwv T

Page 6 of 27

(10)

where the vector g (nwv×1) represents the unwanted contribution at each sampling time. If problem (8) is directly applied to handle spectroscopic data described as Eqs. (9)-(10), the correctness of the results depends on the structure of the reaction model, which can be based on the stoichiometric coefficient matrix Stsub and Zin,sub . According to the above definition, we have T d [c1 (t ), c2 (t ),..., cnc (t )]T d [ zv (1) (t ), zv (2) (t ),..., zv ( nc ) (t )]   St subT  [r1 (t ), r2 (t ),..., rnr (t )]T (11) dt dt

where ri (t) (i=1, …, nr) are the reaction rates. We define Ωsub as

 St sub   sub   T   Z in , sub  and the rank of kernel of Ωsub is rko. The kernel basis is denoted by {bi, i=1,…,rko} and the dimension of bi is nc×1. Since the following equations are satisfied,

bi T

d [c1 (t ), c2 (t ),..., cnc (t )]T  bi T St T  [r1 (t ), r2 (t ),..., rnr (t )]T  0 dt

(12)

and also biTZin,sub = 0, we can conclude that

[c1 (t ), c2 (t ),..., cnc (t )]bi  z0T bi , i  1,..., rko

(13)

If rko > 0, the unwanted contributions G can be decomposed as g G =CS gT，S g  [ s1g , s2g ,..., snwv ]  R ncnwv

(14)

The l-th column (slg) of Sg can be calculated as

slg 

gl bi , i  1,..., rko, l  1,..., nwv z0T bi

(15)

z0Tbi is not equal to zero because of the assumption that concentration profiles of all the absorbing species with kinetic behavior are linearly independent. Then, Eq. (10) can be rewritten as

6


Page 7 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


D =C ( S + S g )T  E

(16)

In this case, the estimated pure spectra matrix is (S + Sg) with problem (8). The unwanted contributions are absorbed completely in the estimation of (S + Sg) and there is no impact on the kinetic parameter estimation. Moreover, additional information is needed to separate S and Sg. If rko = 0, the objective of problem (8) is modified as follows, nc   d ck (ti ) sk (l )  g (l )     i , l  i 1 l 1  k 1  ntp nwv

2

(17)

where the time-invariant unwanted contributions are considered as a function of wavelength ( g(λ) ) and gl=g(λl). In order to decrease the number of decision variables, g(λ) is approximated by an interpolation polynomial with respect to λ. Assuming the wavelength span is equally divided into Ng finite elements leads to the following equation, K

g (l )   m ( lg ) g jg ,m , lg  (l   p jg 1 ) ( p jg   p jg 1 ) , l  1,.., nwv m0

(18)

g jg 1, K  g jg ,0 , jg  2,.., N g where λpjg-1 and λpjg are the endpoints of the jg-th finite element. g jg ,m are the variables that need to be estimated at the collocation points. Then, problem (8) is changed to nc   min   di ,l   ck (ti ) sk (l )  g (l )  i 1 l 1  k 1  s.t. Eq. (18) ntp nwv

K





m 0

m

2

( a ) z j ,m  h j  f ( z j ,a , θ )  0, j  1..N , a  1..K

(19)

z j 1, K  z j ,0 , j  2..N K

ck (ti )   m ( i ) zn' , j ,m , i  (ti  tp j 1 ) (tp j  tp j 1 ) , n '  v(k ) m 0

i  1,.., ntp, l  1,.., nwv, k  1,.., nc It is concluded that when the rank of kernel of Ωsub is equal to zero, the influence of the time invariant unwanted contributions for kinetic parameters estimation can be greatly decreased based on formulation (19).. The decision variables in this formulation are ck(ti), sk(λl), zn,j,m, θ, and gjg,m. b) Time Variant Instrumental Variations 7



Page 8 of 27

The main time variant instrumental variation is drift in the spectroscopic data. Again, G in Eq. (9) can be expressed as a rank-1 matrix:

G = qg T , q =  q1 q2 ... qntp  , g =  g1 T

g 2 ... g nwv 

T

(20)

and the objective of problem (8) is modified as follows, nc   d    i ,l  ck (ti ) sk (l )  q (ti ) g (l )  i 1 l 1  k 1  ntp nwv

2

(21)

where the time-variant unwanted contributions are considered as a function of time and wavelength ( q(t)g(λ) ) and qi=q(ti), gl=g(λl). In order to decrease the number decision variables, g(λ) is approximated by Eq. (18) and q(t) is approximated by an interpolating polynomial with respect to t. Assume the time span is equally divided into Nq elements. Hence, we have the following equation, K

q (ti )   m ( iq )q jq ,m , iq  (ti  tp jq 1 ) (tp jq  tp jq 1 ) , i  1,.., ntp m0

(22)

q jq 1, K  q jq ,0 , jq  2,.., N q where q jq , m are the variables that need to be estimated at the collocation points. In addition, since there are no constraints except bounds to restrict ambiguity q (ti ) g (l ) 

1



q(ti ) and g (l ) , the issue of non-uniqueness occurs; e.g. there is a scale

q (ti ) g (l ) , where α is a nonzero scalar. This leads to nonunique values of qr and

gl and convergence difficulty in solving optimization problem with objective (21). In order to resolve the convergence problem, q(ti) is redefined as qr (ti ) :

q (ti )  qr (tntp )  1 , under the assumption that q(tntp) is q (tntp )

not equal to zero. Then Eq. (22) is changed to, K

qr (ti )   m ( iq )q rjq ,m , iq  (ti  tp jq 1 ) (tp jq  tp jq 1 ) , i  1,.., ntp m 0

q

r jq 1, K

q

r jq ,0

, jq  2,.., N q , q

r Nq , K

=1

Then, problem (8) is changed to

8


(23)

Page 9 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nc   min   di ,l   ck (ti ) sk (l )  qr (ti ) g (l )  i 1 l 1  k 1  s.t. Eq. (18), Eq. (22), Eq. (23) ntp nwv

K





m 0

m

2

( a ) z j ,m  h j  f ( z j ,a , θ )  0, j  1..N , a  1..K

(24)

z j 1, K  z j ,0 , j  2..N K

ck (ti )   m ( i ) zn' , j ,m , i  (ti  tp j 1 ) (tp j  tp j 1 ) , n '  v(k ) m0

i  1,.., ntp, l  1,.., nwv, k  1,.., nc It is concluded that the influence of unwanted time variant contributions for kinetic parameter estimation can be greatly decreased based on formulation (24).. The decision variables in this formulation are ck(ti), sk(λl), zn,j,m, θ, gjg,m, and qrjq,m. In addition, we observe that problem (19) is a special case of problem (24). 3.2 Dependent Absorbing Species In this case, rkp is greater than zero, and kernel basis of matrixΦsub is calculated. The strategy of defining rkp absorbing species as uncolored [32, 34] can be applied. We denote ncr (= nc-rkp) as the corresponding rank of the independent concentration profiles, and we assume the first concentration profiles c1(t), c2(t), ..., cncr(t) are independent. Then, the dependent profiles cncr+1(ti), cncr+2(ti), ..., cnc(ti) can be expressed as follows, [cncr+1(ti), cncr+2(ti), ...,cnc(ti)]T=H . [c1(ti), c2(ti), ...,cncr(ti)]T , i=1,...,ntp, HRrkp×ncr where H matrix is defined as Φsub_uT (Φsub_cT)-1, the subscript sub_c refers to the ncr colored species and the subscript sub_u refers to the rkp uncolored species. Defining ci = [ci(t1), ci(t2), ..., ci(tntp)]T, and the spectra as si=[si(λ1), si(λ2), …, si(λnwv)], i=1,...,nc, we have C.ST = [c1, c2, ..., cncr, cncr+1, ..., cnc] . [s1, s2, ..., sncr, sncr+1, ..., snc]T = [c1, c2, ..., cncr] . ([s1, s2, ..., sncr]+[sncr+1, ..., snc]H)T.

We also redefine [rs1, rs2, ..., rsncr] = [s1, s2, ..., sncr]+[sncr+1, ..., snc]H so that C.ST still satisfies the Beer-Lambert law. Also, problems (8) and (24) can still be used by changing nc to ncr, and linear dependence of concentration profiles will not affect the kinetic parameter estimation. So far, kinetic parameter estimation based on spectroscopic data with unwanted contributions can be properly performed under the assumption that the type of unwanted contributions (time-invariant or time-variant) is known beforehand by choosing problem (8), (19) or (24). In practice, although we may not know whether there are unwanted contributions in the measured spectroscopic data and which type it is, we can determine which

9



problem formulation should be used by checking the number of components with independent concentration profiles from the measured spectroscopic data matrix (of rank ncms); this can be done with singular value decomposition (SVD) [35] under the assumption that the pure spectra for all the absorbing species are independent. If ncr is equal to ncms, there are no unwanted contributions in the measured spectroscopic data or the unwanted contributions can be decomposed into C . SgT; this is similar to the case in section 3.1, where problem (8) is used. If ncr = ncms-1 there exist unwanted contributions in the measured spectroscopic data and then problem (24) is used; this covers the scenarios with unwanted time-invariant or time-variant contributions. In this work we assume the main time variant instrumental variation is due to drift. Other types of time variant disturbance may lead to ncr < ncms-1 and formulation (24) should be modified correspondingly. The approach for kinetic parameter estimation based on spectroscopic data with unwanted contributions (or not) can be summarized as follows, Step 1 Construct the time invariant pseudo-equivalence matrix Φsub, and calculate the corresponding rank of null space rkp. Step 2 If rkp > 0, choose only ncr absorbing species to link to the kinetic model, i.e. change nc to ncr in the objective function of problem (8) and (24). Step 3 Obtain the number of components existing in the measured spectroscopic data (ncms) by SVD. Step 4 If ncr = ncms, choose problem (8) to estimate kinetic parameters; otherwise, use problem (24) to perform kinetic parameter estimation. Step 5 Calculate standard deviations for the estimated kinetic parameters based on the approach proposed in Chen et al. [32].

4. Numerical Studies To demonstrate how to determine kinetic parameters from spectra with or without unwanted contributions, we consider four case studies. The first two provide a detailed presentation of the approach by using two simulated data sets obtained from a two-step consecutive reaction and an aspirin production process with independent absorbing species. The third case study considers a reaction with dependent absorbing species. The fourth case study applies our approach to an actual experiment (Michael reaction). K, Ng, Nq in Eq. (18) and (22) are all set to 3. All the involved problems are solved with IPOPT 3.8.1 (an interior point filter line search strategy based nonlinear programming solver36) on a Lenovo S230u Twist running Windows 8 with Intel® Core™ i7-3537U CPU, 2.00 GHz, and 8.0 GB RAM.

10


Page 10 of 27

Page 11 of 27

4.1 Reaction (AB and BC) [37] We consider a two-step consecutive reaction (AB and BC) with an initial concentration of species A of 1.0×10-3mol/L. The data matrices were obtained by simulating the reaction mechanism with rates k1=2.0s-1 and k2=0.2s-1. Three IR spectra with considerable overlaps were used as pure species spectra of the three assumed components. Measurement error was introduced by adding normally distributed random numbers. The dimension of the data matrix D1 is 300×100. A, B and C are absorbing species. The time invariant pseudo-equivalence matrix Φsub for this case is

 sub

 1   0 1.0 103

1 0 1 1  0 0 

The rank of Φsub is 3 and rkp = 0, nc=ncr=3. Hence, the concentration profiles for A, B and C are independent and problem (8) is solved. The estimated values for k1 and k2 are 1.97266 and 0.20122. The standard deviations for k1 and k2 are 1.2×10-2 and 7.4×10-4. The estimated pure absorbance profiles S1 are shown as Figure 1. 160 A

B

C

140

Absorbance

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


120

100

80

60 2600

2500

2400

2300

2200

2100

2000 -1

1900

1800

1700

1600

Wavenumber (cm ) Figure 1. Estimated pure absorbance for A, B and C

Consider the spectra D2 with time-invariant unwanted contributions as D2(i, l) = D1(i, l) + (0.01 + 2.5×10-4×λl2.15), l=1, 2, … , 100. The matrix Ωsub for this case is

 1 1 0   sub     0 1 1  11



rko = 1 and the kernel basis of Ωsub is b1=[1, 1, 1]T. If we know the unwanted contributions are time-invariant, based on Eq. (15), the absorbance matrix S2calc = (S + Sg) can be calculated as 2.15  0.01  2.5 104  12.15 0.01  2.5 104  22.15  0.01  2.5 104  100 S2calc  S   b1， b1 ,..., b1  0.001, 0, 0 b1 0.001, 0, 0 b1 0.001, 0, 0 b1  

however, we do not know the exact S which is used to generate D1. S1 is obtained by solving Problem (8) with D1, which is an estimation for S. S1 is used to replace S for calculating S2calc. According to the proposed approach, problem (8) is used for handling D2 and also yields S2 for comparison in Figure 2. The estimated values for k1 and k2 are 1.97266 and 0.20122 and relative errors are 1.37% and 0.61%. The standard deviations for k1 and k2 obtained by problem (8) are 1.2×10-2 and 7.4×10-4. They are the same as those of D1. The profiles of estimated absorbance matrix S2 and S2calc (using (15)) are shown as Figure 2. 160

140

Absorbance

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 27

120

100

80

60 2600

Estimated A Estimated B Estimated C Calculated A Calculated B Calculated C 2500

2400

2300

2200

2100

2000 -1

1900

1800

1700

1600

Wavenumber (cm ) Figure 2. Estimated and calculated pure absorbance for A, B and C

From Figure 2, we can see that the estimated and calculated pure absorbance profiles completely overlap for all the absorbing species. Eqs. (15)-(16) are verified. Moreover, if time-invariant property of the unwanted contributions is not known beforehand for this case, problem (8) can also be solved correctly by checking the equivalence between ncms and ncr using SVD. In the following, consider the spectra D3 with time variant drift as D3(i, l) = D1(i, l) + (0.2×ti1.2) × (0.01 + 2.5×104×λ 2.15), l

i=1, 2, …, 300, l=1, 2, …, 100. ncms of D3 is determined as 4 by SVD. ncms - 1 = ncr is satisfied.

According to the proposed approach, problem (24) is used for handling D3. Successful termination without regularization of the Hessian matrix in IPOPT indicates that second order sufficient conditions are satisfied and

12


Page 13 of 27

that all the parameters are estimated uniquely. The estimated values for k1 and k2 are 1.96938 and 0.19978. The relative errors are 1.53% and 0.11%, respectively. The standard deviations for k1 and k2 are calculated in this case and they are 4.1×10-3 and 9.1×10-5, respectively. The concentration profiles for species A, B, and C with the estimated parameter values and the profile bands corresponding to the standard deviations are shown as Figure 3. However, if time variant drifted is ignored and problem (8) is used, the (poorly) estimated values for k1 and k2 are 1.67077 and 0.0001 (the lower bound for k2). -3

1

Concentration (mol/L)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


x 10

-4

x 10

A B C

7.72

0.8

7.7

-4

1.22

0.6

1.24

6.48

x 10

6.46 0.4

-4

x 10 3.2

0.2

3.18 0.56 0 0

5.73 5.74 5.75

1

2

3

0.58 4

5

Time (s)

6

7

8

9

10

Figure 3. Profiles for A, B, C and profile bands corresponding to the standard deviations

4.2 Aspirin Synthesis with Dissolution, Reaction and Crystallization Processes A large parameter estimation problem with complex kinetics is considered. The reaction mechanism of aspirin production (a non-isochoric process) is described in Ref. 17 with the following reactions, k1 SA  AA   ASA  HA k2 ASA  AA   ASAA  HA k3 ASAA  H 2O   ASA  HA k4 AA  H 2O   2 HA kd SA( s )   SA kc ASA   ASA( s )

The stoichiometric coefficients matrix St is

13


(25)


 0 1 1  0 0 1  0 0 0 St    0 0 1  1 1 0  0 0 0

1 0 1 1 2 0 0

0 0 0 0 1

Page 14 of 27

0 1 1 0  1 1 1  0 0 1 0 0 0  1 0 0  1

0

(26)

each column of (26) corresponding to SA(solid), SA, AA, HA, ASA(solid), ASA, ASAA and H2O. The reaction, dissolution and crystallization rate laws are shown as follows,

r1  k1cSA (t )c AA (t ) r2  k2 c ASA (t )c AA (t ) r3  k3c ASAA (t )cH 2O (t ) r4  k4 c AA (t )cH 2O (t ) k  c (T )  c (t )  , if m (t )  0 SA SA rd   d 0, if mSA (t )  0

(27)

d

sat SA



sat rg  kc max  c ASA (t )  c ASA (T ), 0 



c

The concentration and mass profiles for the reaction system can be described by the following system of ODEs, 





I) m SA   M SAVrd , II) c SA

V  rd  r1  cSA V





V  r1  r2  r4  c AA V



V  r1  r2  r3  2r4  cHA V

III) c AA



IV) c HA 



V) m ASA  M ASAVrg , VI) c ASA



V  r1  r2  r3  rg  c ASA V





VII) c ASAA

V  r2  r3  c ASAA V 

f V VIII) c H 2O  r3  r4  cHin2O  cH 2O V V ns 4    f IX) V  V i    i , j rj   i ,d rd   i ,c rc   i cHin2O  V i 1  j 1  

14


(28)

Page 15 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Deionized water was delivered by the syringe pump in the process. This is represented by (VIII) in Eq. (28). The rate f(t) was 1mL/min, and we introduce rp = f(t)cinH2O/V(t) as a pseudo reaction with no reactant and water as the product. Hence the St matrix can be expanded to

 0 1 1   0 0 1 0 0 0  St ex   0 0 1  1 1 0  0 0 0 0 0 0 

1 0 1 1 2 0 0

0 0 0 0 1

0 0

0 1 1 0  1 1 1  0 0 1 0 0 0  1 0 0  0 0 1  1

0

(29)

We generate simulated UV/Vis absorbance data for known parameters k1=0.036031, k2=0.15961, k3=6.8032, k4=1.8029, kc=0.75669, kd=7.1109 and cSAsat=2.0627 and SA, ASA and ASAA as the absorbing species. The pure component spectra S for species SA, ASA and ASAA are the same as in [17] and measurement error was introduced with normally distributed numbers. The dimension of the simulated spectra matrix D1 is 471×111. Here ncms = 3 for D1 as determined by SVD and ncms = ncr is satisfied. According to the proposed approach, formulated (8) is used for handling D1. The estimated values are k1=0.036025, k2=0.159578, k3=6.97449, k4=1.87378, kc=0.755951, kd=7.10329 and cSAsat=2.06281. Consider the spectra D2 with unwanted contributions as D2(i, l) = D1(i, l) + (0.001 +10-2 log(λl)), i=1,2,…,471, l=1, 2, … , 111. If we know the unwanted contributions are time-invariant, problem (19) is chosen based on the analysis of the kernel basis of Ωsub, which is null (rko = 0). The estimated values are k1=0.0360256, k2=0.15959, k3=6.96646, k4=1.87101, kc=0.756297, kd=7.10288 and cSAsat=2.06281. The comparisons between the true values and the estimated values are shown as Table 1. Table 1. The statistics of the estimated values obtained by problem (19) Parameter k1 k2 k3 k4 kc kd cSAsat

True Value 0.036031 0.15961 6.8032 1.8029 0.75669 7.1109 2.0627

Estimated Value 0.0360256 0.15959 6.96646 1.87101 0.756297 7.10288 2.06281

Absolute Error 5.4×10-6 2.0×10-5 -0.16 -6.8×10-2 3.9×10-4 8.0×10-3 1.1×10-4

Relative Error 1.5×10-4 1.3×10-4 2.4×10-2 3.8×10-2 5.2×10-4 1.1×10-3 5.3×10-5

Standard Deviation 4.3×10-6 1.1×10-4 9.7×10-2 3.9×10-2 1.9×10-3 9.2×10-3 3.5×10-4

From Table 1 we can see that relative errors of all the estimated parameters are small and the corresponding standard deviations are also small. If the unwanted contributions are not known to be time invariant for this case, the proposed approach uses problem (24) after checking that ncms-1 = ncr=3 (using SVD for D2). The solution process terminates successfully and IPOPT validates that second order sufficient conditions are satisfied, so that

15



all the parameters are estimated uniquely. The estimated values are k1=0.0360256, k2=0.15959, k3=6.96646, k4=1.87101, kc=0.756297, kd=7.10288 and cSAsat=2.06281. The comparisons between the true values and the estimated values obtained by problem (24) are shown as Table 2. The profiles bands corresponding to the standard deviation for absorbing species SA, ASA and ASAA are shown as Figure 4. Table 2. The statistics of the estimated values obtained by problem (24) Parameter k1 k2 k3 k4 kc kd cSAsat

True Value 0.036031 0.15961 6.8032 1.8029 0.75669 7.1109 2.0627

Estimated Value 0.0360228 0.159765 6.96336 1.86905 0.756992 7.10366 2.06277

Absolute Error 8.2×10-6 -1.5×10-4 -0.16 -6.6×10-2 -3.0×10-4 7.2×10-3 -7.0×10-5

Relative Error 2.3×10-4 9.7×10-4 2.4×10-2 3.7×10-2 4.0×10-4 1.0×10-3 3.4×10-5

Standard Deviation 4.3×10-6 1.1×10-4 9.7×10-2 3.9×10-2 1.9×10-3 9.2×10-3 3.5×10-4

3 SA


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 27

ASA

ASAA

2.5 1.958 2

1.957

9.72 9.74 9.76

0.727

1.5

0.726 0.725 148.7

1.78

1

1.775

0.5

0.96 0.98 0 0

148.75

50

100

1

1.02

Time (min)

150

200

250

Figure 4. The profiles bands corresponding to the standard deviation for absorbing species SA, ASA and ASAA

The results in Table 2 are nearly the same as Table 1. The fact that Problem (19) can be regarded as a special case of problem (24) is verified by this case. The estimated unwanted contributions g obtained by solving problems (19) and (24) are shown as Figure 5.

16


Page 17 of 27

0.0275

Unwanted Contribution

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Exact Obtained by formulation (19) Obtained by formulation (24)

0.027

0.0265

0.026

0.0255 0.0255 0.0254 0.0254 0.0253 270

275

0.0267

280

0.0266 0.0255

0.025 260

0.0265

280

300

320

Wavelength (nm)

340

360 362 364 366 368 360

380

Figure 5. The exact and estimated unwanted contributions obtained by solving problem (19) and (24)

From Figure 5, we can see that both the estimated unwanted contributions obtained by problems (19) and (24) are properly estimated compared with the exact values. 4.3 Rank Deficient Simulation Case The reaction mechanism is shown as, k1 A  B  C k2 C  D  E

(30)

k3 D  E  F

The initial concentrations for A, B, C, D, E and F are 0.5 mol/L, 0.5 mol/L, 0, 0.5 mol/L, 0, 0. The absorbing species are A, C, D and F, and four IR spectra were simulated for the four Gaussian peaks. Simulated IR absorbance data for known k1=0.3, k2=0.1, and k3=0.05 were generated. Normally distributed white noise was added to the data matrix. The dimension of the simulated data matrix D1 is 301×101 and is shown as Figure 6

17



Page 18 of 27

0.3 0.25 0.2 0.15 0.1 0.05 0 -0.05 2500

2450

2400

2350

2300

2250

-1

2200

2150

2100

Wavenumber (cm ) Figure 6. The simulated spectra with white noise

The rank of the matrix Φsub is 3. Hence the rank of the concentration profiles of A, C, D and F is 3. The concentration profiles of A, C, and D are independent by checking the rank of the corresponding submatrix of Φsub. In this work, A, C, and D are defined as the colored species, i.e. only these species are linked to the kinetic model, while F is defined as the uncolored species. Here ncms = 3 for D1 as determined by SVD and ncms = ncr is satisfied. According to the proposed approach, formulated (8) is used for handling D1. The estimated values are k1=0.296515, k2=0.100143, k3=0.0498703 and the corresponding standard deviations are 2.3×10-3, 5.2×10-4, and 1.9×10-3. Now consider the time variant spectra D2 with unwanted contributions as D2(i, l) = D1(i, l) + (5×105×t 1.6)×(0.01×log(λ )0.6), l i

i=1,2,…,301, l=1, 2, … , 101. Here ncms= 4 for D2 as determined by SVD and ncms -

1 = ncr is satisfied. According to the proposed approach, problem (24) is used for handling D2. IPOPT validates that sufficient second order conditions hold, so all the parameters are estimated uniquely. The estimated values and relative errors are k1=0.306701 (2.23%), k2=0.0966451 (3.35%), k3=0.0511816 (2.36%) and the corresponding standard deviations are 1.5×10-3, 2.8×10-4, and 1.7×10-4. These estimated values are acceptable compared with true values. The profiles bands corresponding to the standard deviation for absorbing species A, C, D and F are shown as Figure 7. Note that if time variant drift is ignored and Problem (8) is used, the (poorly) estimated values are k1=0.267331, k2=0.168616, k3=0.0297447.

18


Page 19 of 27

0.5 A


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.4

C

0.385

F

0.1406

0.3845

0.1404 13

0.3

D

13.5

193.1 193.2 193.3 0.132

0.114 0.113 0.112

0.2

0.1318 22

288.6 288.7 288.8 288.9

22.5

0.1

0 0

50

100

150

200

Time (s)

250

300

Figure 7. The profiles bands corresponding to the standard deviation for absorbing species A, C, D and F

4.4 Michael Reaction As the last case study, we consider the Michael Reaction38 with actual industrial data. The experimental steps are shown as follows, 

Add solid A and solvent 1 and agitate until clear solution is observed



Increase the reaction temperature to desired value



Transfer base/THF solution to the reactor



Transfer the solution of C/THF over 3.5 h at constant flow rate.



After addition is complete, add 3.5 ml of THF to rinse the syringe



Age the solution

The reaction initial conditions are shown as Table 3, Table 3 Initial conditions for Michael reaction Species AH B C Solvent 1 Solvent 2 Amount

4.83 g

0.378 g

4.18

59.2 ml

12.3 ml

Solvent 1 is used to dissolve A and solvent 2 is used to dissolve C. Base (B) is dissolved in 1.3 ml of solvent 2. C is added via controlled addition over 3.5 h. We consider FT-IR spectra obtained from the following reaction mechanism.

19



Page 20 of 27

k1 AH  B   A   BH  k2 A   C   AC k 2 AC  A  C

(31)

3 AC  AH   P  A

k

k4 AC  BH   P  B

where AH (Michael Donor) and C (Michael Acceptor) are starting materials, B is a base, BH+, A-, and AC- are reaction intermediates, and P is the reaction product. The rate laws are shown as follows, r1  k1c AH cB r2  k2 c A cC r2  k2 c AC 

(32)

r3  k3c AC  c AH r4  k4 c AC  cBH 

The change of concentration for each liquid of the reaction system can be described by, 

dV const flowrate, t  3.5h dc AH V  ,  r1  r3  c AH dt 0, t  3.5h dt V 



dc  dcB V V  r1  r4  cB , A  r1  r2  r2  r3  c A dt V dt V dcBH  dt dc AC  dt



(33)



mC _ add / V / 3.5, t  3.5h dc V V  r1  r4  cBH  , C  r2  r2  cC   V dt V 0, t  3.5h 



dc V V  r2  r2  r3  r4  c AC  , P  r3  r4  cP V dt V

In this case study, AH, B, C and P are absorbing. The rank of the matrix Φsub is 4. Hence, ncr is 4. ncms = 5 for the measured spectroscopic data as determined by SVD and ncms - 1 = ncr is satisfied. According to the proposed approach, problem (24) is used for estimating kinetic parameters, IPOPT validates that sufficient second order conditions are satisfied and all the parameters are estimated uniquely. The estimated values are shown as Table 4. Table 4. Estimated values for kinetic parameters obtained by problem (24) and (8) Kinetic Parameter

k1

k2

k-2

k3

k4

Lower Bound

1.0×10-5

1.0×10-5

1.0×10-5

1.0×10-5

1.0×10-5

Upper Bound

50

Infinity

150

Infinity

Infinity

Estimated by Problem (24)

1.42258×10-2

3.69329×103

150

1.82117×101

6.30813×102

20


Page 21 of 27

Estimated by Problem (8)

4.2283×101

1.6163×10-1

1.0×10-5

1.0×10-5

9.74021×10-3

The concentration profiles for absorbing species AH, B, C and P with the estimated values obtained by problem (24) and (8) are shown in Figure 8 Formulation (24)

Formulation (8)

0.4

0.4 AH

0.35

Mol Fraction

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


B

C

P

0.35

0.3

0.3

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0 0

200

400

600 800 Time (min)

1000 1200 1400

0 0

AH

200

400

B

600 800 Time (min)

C

P

1000 1200 1400

Figure 8. The concentration profiles for AH, B, C and P with the estimated values obtained by problem (24) and (8)

From Table 4 and Figure 8, we see that kinetic parameters k3 and k4 estimated by problem (8) are very small, which is not consistent with the product yield. In addition, we can find that k-2 is active at its upper bound obtained by problem (24). The optimal objective value for this case is 3.2323066×10-1. The multiplier for k-2 is 1.85×10-7, which is very small compared with the optimal objective. Moreover, if the upper bound for k-2 increases to 300, the optimal objective value is 3.2321665×10-1 and the estimated values are 1.46101×10-2, 7.11941×103, 300, 1.79975×101, and 6.29248×102. There is a little improvement for the optimal objective value with large kinetic parameters variations. The solution for k-2 still hits the upper bound and the multiplier for k-2 is 4.71×10-8. From these data, we can see that the objective function near the optimal solution is flat, which means the estimation for Michael reaction is nonunique without some suitable constraint on the kinetic parameters. Hence, it is reasonable to fix some parameters, in order to improve the estimability. Here, a similar approach to the one in Ref. 32 is used. First, we sort the five free kinetic parameters using a Gram-Schmidt orthogonalization procedure39 with the estimated kinetic parameter values obtained by problem (24) shown in Table 4. We find that the sensitivity vector corresponding to k2 has the smallest projected norm. Hence, the upper bound for k-2 is changed to 300 and k2 is fixed to 3.69329×103. Problem (24) is solved again and IPOPT confirms a unique solution. The standard deviations for the estimated kinetic parameters are calculated with fixed k2. The statistical results are shown in Table 5. Table 5. Variances and standard deviations for the estimated kinetic parameters with fixed k2 Kinetic Parameter

k1

k2

k-2

k3

k4

Upper Bound

50

Infinity

300

Infinity

Infinity

21



Estimated Value (EV)

1.44836×10-2

3.69329×103

1.5302×102

1.79586×101

6.22744×102

Variance

2.2×10-10

------

1.2

6.9×10-3

3.5×101

Standard Deviation (SD)

1.5×10-5

------

1.1

8.3×10-2

5.9

Ratio of SD to EV

0.10%

------

0.72%

0.46%

0.95%

From Table 5, we can see that the ratios of standard deviation to estimated values are less than 1%, which indicates that estimability is satisfied and the problem is well-posed. The concentration profiles for species AH, B, C and P with the estimated values for all the parameters in Table 5 and the bands corresponding to the standard deviations are shown in Figure 9. 0.4 0.35

Mol Fraction

0.3

0.2688

AH

0.2686

B

C

P

0.2684 71.8 72 72.2 72.4

0.25 0.283

0.0835

0.2

0.0248

0.2828

0.083

0.15

695 695.5 696

0.0247

202.8203203.2

0.1

1107.4

1107.6

0.05 0 0

200

400

600

800

1000

Time (min)

1200

1400

Figure 9. Profiles for AH, B, C and P with the estimated values obtained by problem (24) and profile bands corresponding to the standard deviations

The obtained absorbance profiles for AH, B, C and P are shown as Figure 10. 15 AH

Absorbance

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 27

B

C

P

10

5

0 2600

2400

2200

2000

1800

1600

1400 -1

Wavenumber (cm ) 22


1200

1000

800

600

Page 23 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10. Absorbance profiles for AH, B, C and P obtained by solving problem (24) with fixed k2

5. CONCLUSIONS Unwanted contributions commonly appear in the measured spectroscopic data because of the baseline shift or distortion, the presence of inert absorbing interferences with no kinetic behavior, and drifts. These unwanted contributions may have significant impact on the accuracy of the kinetic parameters estimation. In this work, estimation problem with time invariant contributions is studied in depth. We derive conditions under which the time invariant contributions have no influence on the estimation accuracy and provide rigorous equations (14)(15) for calculating the estimated S matrix. These were demonstrated and verified by the two-step consecutive reaction case.. The estimated S matrix is perfectly consistent with the theoretical calculation. When time invariant contributions affect estimation accuracy, problem (19) is proposed to estimate kinetic parameters and separate time invariant contribution from measured spectra simultaneously. In addition, problem (19) is considered as a special case of problem (24), which is designed for handing spectra with time variant contributions. The results for solving the problem of aspirin production process with problems (19) and (24) are nearly the same. The estimated kinetic parameters and the separated time invariant contributions are very close to the true values. The relative errors and the standard deviations show that the accuracy is acceptable. The above conclusions are under the assumption that all absorbing species are independent. With the assistance of a colored strategy, kinetic parameter estimation for reaction with dependent absorbing species can also be performed similarly. Simulation spectra with drift of a reaction with dependent absorbing species demonstrated good performance of the proposed approach. Sometimes, it is not easy to know whether the unwanted contribution is time invariant, even if there exist unwanted contributions in the measured spectra. For this reason, a unified framework for kinetic parameter estimated is established in this paper. All of the numerical results including the experimental case (Michael reaction) support the advantages of the proposed framework. Finally, we note that our handling of drift in this unified formulation rests on the linear subspace properties for Beer’s Law. In future work we plan to consider extensions to these properties, in order to handle nonlinear features such as peak drift40 and nonlinear deviations from Beer’s Law.

ACKNOWLEDGEMENT This research was supported by Lilly Research Award Program. We thank Carla Luciani, Michael Frederick and Zhenqi (Pete) Shi (Eli Lilly) for their assistance collecting experimental data.

23



Reference (1) Roggo Y, Chalus P, Maurer L, Lema-Martinez C, Edmond A, Jent N. A review of near infrared spectroscopy and chemometrics in pharmaceutical technologies. J Pharm Biomed Anal. 2007; 44(3): 683-700. (2) Quinn AC, Gemperline PJ, Baker B, Zhu M, Walker DS. Fiber-optic UV/visible composition monitoring for process control of batch reactions. Chemom Intell Lab Syst. 1999; 45(1-2): 199-214. (3) Ng DJW, Assirelli M. Mixing study in batch stirred vessels using a fibre optic UV-VIS monitoring technique: a novel method. Chem Eng Res Des. 2007; 85(10): 1348-1354. (4) Gurden SP, Westerhuis JA, Smilde AK. Monitoring of batch processes using spectroscopy. AIChE J. 2002; 48(10): 2283-2297. (5) Lawton WH, Sylvestre EA. Self modelling curve resolution. Technometrics.1971; 13(3): 617-633. (6) Lawton WH, Sylvestre EA, Maggio MS. Self modelling nonlinear regression. Technometrics. 1972; 14(3): 513-532. (7) Jiang J, Liang Y, Ozaki Y. Principles and methodologies in self-modeling curve resolution. Chemom Intell Lab Syst. 2004; 71(1): 1-12. (8) De Juan A, Tauler R. Multivariate curve resolution (MCR) from 2000: progress in concepts and applications. Crit Rev Anal Chem. 2006; 36(3-4): 163-176. (9) Abdollahi H, Tauler R. Uniqueness and rotation ambiguities in multivariate curve resolution methods. Chemom Intell Lab Syst. 2011; 108(2): 100-111. (10)Skvortsov AN. Estimation of rotation ambiguity in multivariate curve resolution with charged particle swarm optimization (cPSO-MCR). J Chemometr. 2014; 28(10): 727-739. (11)Levenberg K. A method for the solution of certain problems in least squares. Q J Appl Math. 1944; 2(2): 164-168. (12)Marquardt DW. An algorithm for least-squares estimation of nonlinear parameters. J Soc Ind Appl Math. 1963; 11(2): 431-441. (13)Maeder M, Zuberbuhler AD. Nonlinear least-square fitting of multivariate absorption data. Anal Chem. 1990; 62(20): 2220-2224. (14)Bijlsma S, Boelens HFM, Smilde AK. Determination of rate constants in second-order kinetics using UV-visible spectroscopy. Appl Spectrosc. 2001; 55(1): 77-83. (15)Puxty G, Maeder M, Hungerbuhler K. Tutorial on the fitting of kinetics models to multivariate spectroscopic measurements with non-linear least squares regression. Chemom Intell Lab Syst. 2006; 81(2): 149-164. (16)Kompany-Zareh M, Khoshkam M. Applying constraints on model-base methods: Estimation of rate constants in a second order consecutive reaction. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy. 2013; 102: 319-326. (17)Joiner DE, Billeter J, McNally MEP, Hoffman RM, Gemperline PJ. Comprehensive kinetic model for the dissolution, reaction, and crystallization process involved in the synthesis of aspirin. J Chemometr. 2014; 28(5): 420-428.

24


Page 24 of 27

Page 25 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18)Pla FFP, Baeza JJB, Llopis E, Baeza MP, Fernandez L. Comparative study of hard- and soft-modeling algorithms for kinetic data processing. International Journal of Chemical Kinetics. 2016; 48(8): 449463. (19)Biegler LT. Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes. Philadelphia: Society for Industrial and Applied Mathematics; 2010. (20)Sawall M, Borner A, Kubis C, Selent D, Ludwig R, Neymeyr K. Model-free multivariate curve resolution combined with model-based kinetics: algorithm and applications. J Chemometr. 2012; 26(10): 538-548. (21)Sawall M, Kubis C, Borner A, Selent D, Neymeyr K. A multiresolution approach for the convergence acceleration of multivariate curve resolution methods. Analytica Chimica Acta. 2015; 891: 101-112. (22)Tauler R. Application of non-linear optimization methods to the estimation of multivariate curve resolution solutions and of their feasible band boundaries in the investigation of two chemical and environmental simulated data sets. Analytica Chimica Acta. 2007; 595: 289-298. (23)Vaia A, Sahinidis NV. Simultaneous parameter estimation and model structure determination in FTIR spectroscopy by global MINLP optimization. Comput Chem Eng. 2003; 27: 763-779. (24)Chen W, Biegler LT, Garcia-Muñoz S. An approach for simultaneous estimation of reaction kinetics and curve resolution from process and spectral data. J Chemometr. 2016; 30(9): 506-522. (25)Rinnan A, Berg FVD, Engelsen SB. Review of the most common pre-processing techniques for nearinfrared spectra. Trends in Analytical Chemistry. 2009; 28(10): 1201-1222. (26)De Juan A, Maeder M, Martinez M, Tauler R. Combining hard- and soft-modelling to solve kinetic problems. Chemom Intell Lab Syst. 2000; 54(2): 123-141. (27)Amigo JM, de Juan A, Coello J, Maspoch S. A mixed hard- and soft-modelling approach to study and monitor enzymatic systems in biological fluids. Analytica Chimica Acta. 2006; 567: 245-254. (28)Mas S, de Juan A, Lacorte S, Tauler R. Photodegradation study of decabromodiphenyl ether by UV spectrophotometry and a hybrid hard- and soft-modelling approach. Analytica Chimica Acta. 2008; 1828. (29)Mas S, Bendoula R, Agoda-Tandjawa G, de Juan A, Roger JM. Study of time-dependent structural changes of laponite colloidal system by means of near-infrared spectroscopy and hybrid hard- and softmodelling multivariate curve resolution-alternating least squares. Chemom Intell Lab Syst. 2015; 142: 285-292. (30)Biegler LT, Cervantes AM, Wächter A. Advances in simultaneous strategies for dynamic process optimization. Chemical Engineering Science. 2002; 57(4): 575-593. (31)Kameswaram S, Biegler LT. Convergence rates for direct transcription of optimal control problems using collocation at Radau points. Computational Optimization and Applications. 2008; 41(1): 81-126. (32)Chen W, Biegler LT, Garcia-Muñoz S. Kinetic parameter estimation based on spectroscopic data with unknown absorbing species. AIChE J. 2018; 64(10): 3595-3613. (33)Rodrigues D, Srinivasan S, Billeter J, Bonvin D. Variant and invariant states for chemical reaction systems. Comput Chem Eng. 2015; 73: 23-33.

25



(34)Maeder M, Neuhold YM. Practical data analysis in chemistry (Data Handing in Science and Technology, 26). Elsevier Science, Amsterdam, 2007. (35)Golub GH, Reinsch C. Singular value decomposition & least squares solutions. Num Math. 1970; 14(5): 403-420. (36) Wächter A, Biegler LT. On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Mathematical Programming. 2006; 106(1): 2557. (37)Jaumot J, Juan A, Tauler R. MCR-ALS GUI 2.0: new features and applications. Chemometrics and Intelligent Laboratory Systems. 2015; 140(15): 1-12. (38)Mather BD, Viswanathan K, Miller KM, Long TE. Michael addition reactions in macromolecular design for emerging technologies. Prog Polym Sci. 2006; 31: 487-531. (39)Kravaris C, Hahn J, Chu Y. Advances and selected recent developments in state and parameter estimation. Comput Chem Eng. 2013; 51: 111-123. (40)Tulpan D, Léger S, Belliveau L, Culf A, Čuperlović-Culf M. MetaboHunter: an automatic approach for identification of metabolites from 1H-NMR spectra of complex mixtures. BMC Bioinformatics. 2011; 12:400.

26


Page 26 of 27

Page 27 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


27


A Unified Framework for Kinetic Parameter Estimation Based on

Recommend Documents