# State estimation of a solid-state polymerization reactor for PET based on improved SR-UKF.

код для вставкиСкачатьASIA-PACIFIC JOURNAL OF CHEMICAL ENGINEERING Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 Published online 19 August 2009 in Wiley InterScience (www.interscience.wiley.com) DOI:10.1002/apj.306 Research Article State estimation of a solid-state polymerization reactor for PET based on improved SR-UKF Ji Liu and Xing Sheng Gu* Research Institute of Automation, East China University of Science and Technology, Shanghai 200237, China Received 7 September 2008; Revised 18 February 2009; Accepted 23 February 2009 ABSTRACT: A state estimator for the continuous solid-state polymerization (SSP) reactor of polyethylene terephthalate (PET) is designed in this study. Because of its invalidity in the application to some of the practical examples such as SSP processes, the square-root unscented Kalman filter (SR-UKF) algorithm is improved for the state estimation of arbitrary nonlinear systems with linear measurements. Discussions are given on how to avoid the filter invalidation and accumulating additional error. Orthogonal collocation method has been used to spatially discretize the reactor model described by nonlinear partial differential equations. The reactant concentrations on chosen collocation points are reconstructed from the outlet measurements corrupted with a large noise. Furthermore, the error performance of the developed ISR-UKF is investigated under the influence of various initial parameters, inaccurate measurement noise parameters and model mismatch. Simulation results show that this technique can produce fast convergence and good approximations for the state estimation of SSP reactor. ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. KEYWORDS: state estimation; solid-state polymerization reactor; improved square-root unscented Kalman filter; error performance INTRODUCTION The widespread application of polyethylene terephthalate (PET) in the production of automotive tire cords and soft drink bottles has led to a multibillion dollar market for the resin.[1] The PET resin obtained by the traditional synthetic routes (melt process) generally has an average molecular weight of 15 000?20 000 kg/kmol and an intrinsic viscosity (IV) of 0.55?0.65 dl/g, which is unsuitable for bottle production (bottle grade IV 0.72?0.85 dl/g) or high-strength fibers (tire-cord resin IV 0.95?1.05 dl/g). The continuous solid-state polymerization (SSP) process is the most feasible way to increase the molecular weight of PET and other stepgrowth polymers, such as poly(butylene terephthalate) (PBT), poly(ethylene naphthalate) (PEN) and nylon 6.6.[2] The PET melt is extruded to form about 2-mmlong solid pellets. As shown in Fig. 1, these pellets are fed continuously from the top of the reactor, which is essentially an empty cylinder, and an inert gas enters *Correspondence to: Xing Sheng Gu, Research Institute of Automation, East China University of Science and Technology, Box. 303, No. 130, Meilong Road, Shanghai 200237, China. E-mail: xsgu@ecust.edu.cn ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. from the opposite end. This inert gas removes the reaction byproducts so that high molecular weights can be obtained from the reaction. Meanwhile, it also acts to partly regulate the reactor temperature. It is assumed that both the solid flow and the gas flow are close to ideal plug-flow conditions.[3] The principal polycondensation chain-building reaction within the solid PET pellets is of the form Pn + Pm ? Pm+n + C (1) where C is a byproduct, such as ethylene glycol (EG) or water. Pn , Pm and Pn+m represent the polymer with various polymerization degrees. The mechanism of the SSP can be summarized in the following three steps[3] : 1. The reaction step involving the reaction of the corresponding end groups; 2. The diffusion of the byproducts (water and EG) through the polymer matrix (internal diffusion step) to the surface; 3. The diffusion of these byproducts from the surface to the inert gas (interphase diffusion step). The reaction process involves both radial mass transfer in the solid particles and axial mass transfer in Asia-Pacific Journal of Chemical Engineering PET pellet N2 and byproducts pure N2 PET pellet Figure 1. Diagram of an SSP reactor. the reactor. This is a multidimensional and multiphase process. The underlying continuous dynamics of the reactor is described by nonlinear partial differential equations (PDEs) taking into account the reactant concentrations, comprising the reactions of the corresponding end groups, reactant convection, and internal and interphase diffusion of byproducts. Hence, the overall dynamics of the SSP reactor can be described by a nonlinear distributed parameter system. The quality indices of polymer, for example intrinsic viscosity, average molecular weight and polymerization degree, are highly dependent on the operating conditions such as the reaction residence time and the reaction temperature. The knowledge of the time-varying spatial concentration profiles in the reactor is mandatory for the production quality control, monitoring and diagnosis of the SSP process for PET. However, it is technically difficult to measure the time-varying concentration profiles in the whole reaction tube. Instead, most of PET resin manufacture will use the above indices for quality control because they can be analyzed by on-line instrumentation. The concentration profiles of the reactants have to be estimated indirectly from the few outlet measurements. While the current status of the linear estimation theory is a mature topic, challenges exist for the estimation of more complex systems (nonlinear, distributed, hybrid, large-scale physical models), in particular goaloriented modeling issues (estimation for control, diagnosis, detection, monitoring, etc).[4] The SSP reaction process is such a complex system with the characteristic of nonlinear distributed parameter, which is hybrid and time varying. The extended Kalman filter (EKF), which used to be the first choice for state estimation of nonlinear models,[5,6] uses a local slope linearization approach for calculating the mean and covariance of the random variables. Such approximations are valid ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. STATE ESTIMATION OF A SSP REACTOR FOR PET only for ?small? perturbations around the mean. Currently, unscented Kalman filter (UKF) and squareroot UKF (SR-UKF) are the widely used nonlinear filtering strategies, in which the state distribution is approximated by a small set of carefully chosen support points called sigma points.[7] It is accurate to the second (extendable to third)order Taylor series expansion for arbitrary nonlinearity, while the EKF can only achieve first-order accuracy. Furthermore, it is not necessary to compute the Jacobian matrices in the UKF. As a result, the UKF can avoid substantial computational demand and may become a powerful candidate for real-time implementation. Applications of UKF and SR-UKF can be found in open literature, mainly for navigation purposes, e.g. the tracking of a vehicle entering the atmosphere,[8,9] navigation of an unmanned helicopter[10] and the estimation of spacecraft attitude.[11] However UKF has been barely applied to solve chemical engineering problems especially for nonlinear distributed parameter systems.[12] One of the successful applications has been the state estimation of a lumped parameter CSTR system by Romanenko and Castro.[13] Another one is for the stratified domestic hot water storage tank described by nonlinear PDEs.[12] The estimation used both a distributed parameter observer and the standard UKF posterior to the model discretization. A major drawback of the UKF is its instability due to the numerical calculation error.[4] The covariance matrices may lose their positive semidefiniteness during the filtering procedure. To improve the numerical properties, the SR-UKF was proposed by Wan and van der Merwe.[14] In their implementation, the Cholesky factors of the covariance matrices will be propagated directly, similar to the standard square-root Kalman filter (SRKF), avoiding the requirement of re-factorizing at each time step. The SR-UKF algorithm, making use of powerful linear algebra techniques, orthogonal matrix triangularization (QR decomposition) and Cholesky factor update, improves numerical stability and can usually guarantee positive semidefiniteness of covariance matrices. Although more robust than the UKF, it is found that the SRUKF cannot be applied for the state estimation of SSP processes and some other nonlinear stochastic systems[15,16] because the updated matrix of Cholesky factor is still negative definite, which is caused by two problems. One is the negative zeroth weight of sigma points, probably causing Cholesky factor update of the predicted state and measurement covariance matrices failing.[17] The other is the numerical computing error probably causing the Cholesky factor update of the updated covariance matrix failing. This paper contributes to the improvement of the existing SR-UKF Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj 379 380 J. LIU AND X. GU Asia-Pacific Journal of Chemical Engineering algorithm so that it can be applied to arbitrary nonlinear systems with linear measurements such as the SSP process. The rest of the paper is organized as follows: First of all, a simplified one-dimensional dynamic model of the SSP reactor, described by nonlinear PDEs, is introduced as well as the orthogonal collocation method (OCM)[18] used for spatial discretization. The improved SR-UKF (ISR-UKF) algorithm is then developed with theoretical explanation on how to avoid filter instability and bringing in additional error. Next, the state estimator of SSP reactor for PET is implemented on the basis of the ISR-UKF. The time-varying concentrations within the reactor are reconstructed from the outlet measurements. Finally, conclusions and future work will be addressed. MODELING Dynamic model of SSP reactor Although there are extensive publications on modeling of melt reactors for PET, few have been found for SSP reactors. Most of them are particle models,[1,19,20] describing the dynamics within micro particles. Modeling on SSP macro reactor can be found in a few publications.[3,21] It was recognized that these SSP models are very complex with many uncertainties, consequently leading to calculation and solution problems. The model selected is for the estimation of concentration profiles used as quality soft-sensor for quality control. This type of model only deals with the reactant concentrations. Operating conditions, such as flow velocity and reaction temperature, will be chosen as input variables. A simplified macro reactor model based on the mass balance, proposed earlier by the authors,[22] will be adopted. This distributed parameter model consists of a set of nonlinear PDEs with initial and boundary conditions. Boundary conditions: g(t, 0) = g0 , e(t, 0) = e0 ?e(t, z ) ?g(t, z ) |z =1 = |z =1 = 0 ?z ?z e(0, z ) = g(0, z ) = 0 (3) ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. (5) where g, e and Z are, respectively, the average concentrations (in the unit of mol/ ) of EG, hydroxyl end groups and diester end groups; t is the dimensionless time with respect to reaction residence time ? ; and z is the dimensionless axial coordinate with respect to reactor length L. Dz is the average diffusion coefficient of EG. ? is the effective factor of the interphase diffusion of EG. ? and K are, respectively, the reaction rate constant and equilibrium constant. It is evident that model parameters Dz , ? and ? are related to operating conditions. The model assumes constant feed concentrations g0 and e0 : that is, the reactant concentrations in the PET pellet from the melt process are assumed to be constants. For a complete study, the model must be able to simulate the start-up procedure when the reactor is gradually filled with solid particles. Hence, an empty reactor is presumed as the initial condition. The pellets are continuously fed into the reactor with a flow velocity v which is a control variable determining the reaction residence time, defined as ? = vL. In an industrial application, the reaction residence time ? is a discrete variable controlling the polymerization degree. The larger the ? , the higher the molecular weights of the PET obtained. In the state estimation of the SSP reactor, ? is chosen as a control variable, the concentrations g(t, z ) and e(t, z ) are the state variables and the outlet concentration measurements g(t, 1) and e(t, 1) are the observation data. Model discretization The state estimator is subsequently designed on the basis of the discretized model. For chemical processes, the OCM has been widely adopted. One of the advan- Dz ? ? 2 g(t, z ) 4g(t, z )Z (t, z ) ?g(t, z ) ?g(t, z ) 2 + = 2 + ??? e (t, z ) ? ?t ?z K L ?z 2 ?e(t, z ) ?e(t, z ) 4g(t, z )Z (t, z ) + = ?2?? e 2 (t, z ) ? ?t ?z K e(t, z ) Z (t, z ) = 1 ? 2 Initial conditions: (4) (2) tages of this method is that the collocation points are chosen automatically, which tend to converge much faster than the method of line.[23] The less number of Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering STATE ESTIMATION OF A SSP REACTOR FOR PET collocation points required reduces the computational effort in terms of calculation time and data storage. In the OCM, concentration is approximated by a nonsymmetrical polynomial of the form C (t, z ) = C0 (t) + zC1 (t) + z (1 ? z ) N ai (t)Pi ?1 (z ) i =1 (6) where C (t, z ) represents g(t, z ) or e(t, z ); C0 (t) = C (t, 0), the feed concentration; C (t, 1) is the outlet concentration. {Pn (z )} are a set of polynomials such that Pn (z ) is orthogonal to Pm (z ) when m = n. The collocation points are the roots of Pn (z ) : z1 = 0, z2 , z3 , . . . , zN +1 , zN +2 = 1, of which z1 and zN +2 are the boundary points. Here, the transformation of Eqn (2) to ordinary differential equations (ODEs) can be expressed as N +2 N +2 Dz ? dgi (t) + Aij gj (t) = 2 Bij gj (t) dt L j =1 j =1 4gi (t)(1 ? ei (t)/2) 2 + ??? ei (t) ? K (7) N +2 dei (t) Aij ej (t) + dt j =1 4gi (t)(1 ? ei (t)/2) = ?2?? ei2 (t) ? (8) K Aij and Bij are obtained from the roots of the orthogonal polynomials; gi and ei are the timevarying concentration at the i th collocation point zi . i goes from 2 to N + 1. For the boundaries, we have: g1 (t) = g0 , e1 (t) = e0 N +2 j =1 A(N +2)j gj (t) = N +2 As discussed in last section, the measurement equation is determined by y(t) = H x(t) + v(t) where H is a constant 2 О (2N + 2) matrix. y(t) is the two-dimensional observation vector, which is equal to the outlet concentration measurements of EG and the hydroxyl end groups, contaminated by an uncorrelated zero-mean Gaussian white noise v(t). THE IMPROVED SQUARE-ROOT UNSCENTED KALMAN FILTER Since its proposal by Julier and Ulhman,[7] the UKF is mainly used for an arbitrary nonlinear system. However, numerical instability often causes the covariance matrix P losing its positive definiteness during the filtering procedure. ? Consequently, the sigma points ?i ,k ?1 = x?k ?1 ▒ ( (n + ?)P)i cannot be correctly calculated. Furthermore, the most computationally expensive operation in the UKF is the calculation of the square root of the covariance matrices at each time step for the new set of sigma points. Although proposed by Wan and van der Merwe as a countermeasure,[14] the SR-UKF algorithm seems to encounter the same problem as the UKF. It will be discussed further in subsequent sections. The SR-UKF cannot always be successfully implemented because the updated matrices are probably negative definite while updating the Cholesky factors (available as ?cholupdate? in Matlab) such as (c) Sk? = cholupdate{Sk? , ?0,k |k ?1 ? x?? k , W0 } (14) (c) Syk = cholupdate{Syk , Y0,k |k ?1 ? y?? k , W0 } (15) (9) or A(N +2)j ej (t) = 0 (10) j =1 So the distributed parameter model described by Eqns (2)?(5) is discretized to the 2N + 2 ODEs, which can be rewritten as in state space as x?(t) = F (x(t), u(t)) (11) where F is the known nonlinear functions defined by Eqns (7),(8) and (10). X is the unknown time-varying state vector defined as x = [x1 , x2 , . . . x2N +2 ] = [g2 , g3 , . . . , gN +2 , e2 , e3 , . . . , eN +2 ]T Sk = cholupdate{Sk? , Kk Syk , ?1} (12) u represent the control variables. ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. (16) where Sk? , Syk and Sk , respectively, are the Colesky factors of the predicted state and measurement covariance matrices and the updated covariance matrix. Kk is the filter gain matrix. ?0,k |k ?1 and Y0,k |k ?1 represent, respectively, the time update of the zeroth sigma point ? of the state vector and its observation. x?? k and y?k are, respectively, the mean estimate of state vector and the mean observation. The zeroth weight of sigma points, W0(c) , is defined as W0(c) = T (13) ? + 1 ? ?2 + ? n +? (17) where n is the state dimension and ?, ? and ? are the adjustable parameters of the filter, of which ? and ? are very small. Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj 381 382 J. LIU AND X. GU Asia-Pacific Journal of Chemical Engineering It has been found from our investigation that if designed on the basis of the standard UKF or SR-UKF, the state estimator for the SSP reactor would encounter the above problems. The algorithm broke down at second iteration. This work intends to improve the SR-UKF performance. The state-space model after spatial discretization for a SSP reactor is x(k + 1) = F (x(k ), u(k )) + w(k ) y(k ) = H (k )x(k ) + v(k ) (18) where x(k ) is the n-dimensional state vector of the system at time step k , u(k ) is the input vector, y(k ) is the measurement vector, w(k ) and v(k ), respectively, represent the process noise and measurement noise, F Pk? = 2n and the updated matrix from Eqn (15) is ? P ▒ v uu T = Syk Syk T + sgn(W0(c) ) |W0(c) |(Y0,k |k ?1 ? T ? y?? k )(Y0,k |k ?1 ? y?k ) (21) If the zeroth weight W0(c) is negative to ensure a fine tuning of the weights so that the overall prediction error can be reduced,[17] these two matrices are most possibly negative definite. According to Julier,[7] a useful heuristic is to select n + ? = 3 in Eqn (17) when the state is assumed Gaussian. Therefore, when the algorithm approximates higher order systems or probability density distributions, for a fine tuning of the weights, ? is probably much smaller than zero and the above updated matrices will become negative definite. Improvement 1: A modification of the covariance matrix used by Julier and Uhlmann[17] is as follows: Wi(c) [?i ,k |k ?1 ? ?0,k |k ?1 ][?i ,k |k ?1 ? ?0,k |k ?1 ]T + Qk i =1 Pyk = 2n Wi(c) [Yi ,k |k ?1 ? Y0,k |k ?1 ][Yi ,k |k ?1 ? Y0,k |k ?1 ]T + Rk (22) i =1 is the nonlinear state update function and H is the measurement matrix. Moreover, w(k ) and v(k ) are the uncorrelated zero-mean Gaussian white noise sequences and their covariance matrices are, respectively, Qk and Rk . Calculation of Colesky factors of predicted state and measurement covariance matrices Problem 1: According to linear algebra theory, if S is the original Cholesky factor of P , that is P = SS T , rank1 update matrix of the Cholesky factor S is available in MATLAB as S ? = cholupdate{S , u, ▒v } (19) The process is actually the Colesky factorization of ? the matrix P ▒ v ? uu T ,[14] which requires that the updated matrix P ▒ v uu T must be positive semidefinite. But in the SR-UKF, the updated matrix inferred from Eqn (14) is ? (c) ? ?T T P ▒ v uu = Sk Sk + sgn(W0 ) |W0(c) |(?0,k |k ?1 ? T ? x?? k )(?0,k |k ?1 ? x?k ) (20) ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. ?Covariance? is evaluated on ?0,k |k ?1 or Y0,k |k ?1 . 1 , i = 1, Since all of the weights Wi(c) = 2(n + ?) и и и , 2n are positive, the covariance matrices must be positive semidefinite. So we modify the Cholesky factorization of the above covariance matrices to ? Sk = qr [ W1(c) (?1,k |k ?1 ? ?0,k |k ?1 ) и и и (c) W2n (?2n,k |k ?1 ? ?0,k |k ?1 ) T Qk ] Syk = qr [ W1(c) (Y1,k |k ?1 ? Y0,k |k ?1 ) иии T (c) W2n (Y2n,k |k ?1 ? Y0,k |k ?1 ) Rk ] (23) This improvement avoids the negative definite covariance matrices caused by negative zeroth weight and substitutes the Cholesky factor rank-1 update (denoted as cholupdate) by the QR decomposition. Calculation of Colesky factor of the updated covariance matrix Problem 2: It has been mentioned that the updated matrix must be positive semidefinite while updating the Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering STATE ESTIMATION OF A SSP REACTOR FOR PET Cholesky factor. In the SR-UKF technique, the updated matrix inferred from Eqn (16) is as follows: ? (24) P ▒ v uu T = Sk? Sk?T ? Kk Syk (Kk Syk )T The subtraction of two positive semidefinite matrices easily results in a negative definite matrix. It is easy for the updated covariance matrix in Eqn (24) to lose its positive semidefiniteness owing to the accumulated computing error. It is well known that UKF-based algorithms are sensitive to the computing error accumulation,[4] similar to the standard Kalman filter. Improvement 2: As the measurement equation for a nonlinear system is linear as described in Eqn (18), the updated covariance matrix can be derived with the same equations as in classical Kalman filtering ?1 T T ) Pxy Pk = Pk? ? Kk Pyy KTk = Pk? ? Pxy (Pyy (25) Meanwhile, from the UKF, it can be known that the matrices Pyy and Pk? are symmetrical. So the above equation can be rewritten as Pk = Pk? ? Pk? HkT (Hk Pk? HkT + Rk )?1 Hk Pk? (26) In order to improve the numerical stability, the standard SRKF algorithm[24] modifies Eqn (26) by using the Cholesky factors Sk? and Sk instead of Pk? and Pk , which are propagated directly at each iteration. The modification comes from the square-root covariance filtering (SRCF).[25] Sk = Sk? ? Sk? Fk (UkT )?1 (Uk + Rk )?1 FkT (27) where Fk and Uk are subject to Fk = and (Sk? )T HkT 1 Uk = [FkT Fk + Rk ] 2 Step 1: Initialization x?0 = E [x0 ], S0 = chol{E [(x0 ? x?0 )(x0 ? x?0 )T ]} Step 2: Sigma point calculation and time updating ?k ?1 = [x?k ?1 x?k ?1 + (n + ?)Sk x?k ?1 ? (n + ?)Sk ]T ?k |k ?1 = F [?k ?1 , uk ?1 ] x?? k = (29) Hence, originated from the SRKF, Eqns (27)?(29) are used to calculate the Cholesky factor Sk to improve the numerical stability of the SR-UKF. It is easy to prove that the nonzero matrix FkT Fk is positive definite. As a result, the square root of the matrix FkT Fk + Rk in Eqn (29) can be solved by Cholesky factorization. Of course, the matrix FkT Fk + Rk may also lose positive definiteness during the filtering process. In this case the QR decomposition is again imposed. Rk ]T ) (30) Uk = qr([FkT Therefore, the ISR-UKF algorithm for state estimation of arbitrary nonlinear system with linear measurements is summarized as follows: ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Wi(m) ?i ,k |k ?1 i =0 Sk? = qr{[ W1(c) (?1,k |k ?1 ? ?0,k |k ?1 ) T (c) W2n (?2n,k |k ?1 ? ?0,k |k ?1 ) Qk ] } иии Yk |k ?1 = Hk ?k |k ?1 y?? k = 2n Wi(m) Yi ,k |k ?1 i =0 Step 3: Measurement updating Syk = qr{[ W1(c) (Y1,k |k ?1 ? Y0,k |k ?1 ) (c) W2n (Y2n,k |k ?1 ? Y0,k |k ?1 ) Rk ]T } иии Pxk yk = (28) 2n 2n ? T Wi(c) [?i ,k |k ?1 ? x?? k ][Yi ,k |k ?1 ? y?k ] i =0 Kk = (Pxk yk /SyTk )/Syk ? x?k = x?? k + Kk yk ? y?k Fk = (Sk? )T HkT 1 Uk = [FkT Fk + Rk ] 2 Sk = Sk? [I ? Fk (UkT )?1 (Uk + Rk )?1 FkT ] Error of the improved filter In Improvement 1, the modified covariance matrix, taking Pyy for example, is equivalent to calculating the covariance using the original form plus a term [y?? k ? T [17] Y0,k |k ?1 ][y?? k ? Y0,k |k ?1 ] . According to Julier et al ., the calculated covariance is Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj 383 384 J. LIU AND X. GU Pyk = ?fx Pk? ?fx T +E Asia-Pacific Journal of Chemical Engineering 3 2 D 2 f (Dx D 3 f (Dx f )T Dx f (Dx f )T f )T + x + x 3! 2 О 2! 3! where ?fx is the Jacobian of F [?] with respect to x(k ), and the operator Dx f evaluates the total differential of F [?] when x = x ? x. The modified covariance ensures positive semidefiniteness and does not truncate any order. Thus the new filter is still correct up to the second order with error in the fourth or higher order. In Improvement 2, the Cholesky factor Sk of covariance matrix Pk is strictly derived from the SRCF. The improvement should not bring in further error with the improved numerical properties similar to those of the SRKF. Therefore, it can be concluded that the modification on the SR-UKF algorithm introduces no additional error. + ... (31) Thus, the state-space equations will consist of Eqn (11) with N = 5. According to Eqn (12), the state vector is 12 dimensional and defined as x(t) = [g(t, z2 ), и и и g(t, z7 ), e(t, z2 ) и и и , e(t, z7 )]T (34) The feed concentrations[22] on z1 = 0 are g(t, z1 ) = 10?5 mol/ and e(t, z1 ) = 0.0187 mol/ . For integration, the sampling interval T = 0.002 is selected. So the model can be approximated in discrete time. (35) xk +1 = xk + TF (xk , uk ) and restricted by the measurement equation yk = H xk + vk (36) 0 0 0 0 0 1 0 0 0 0 0 0 with H = 0 0 0 0 0 0 0 0 0 0 0 1 where vk is the white noise sequence and its covariance matrix is assumed as Rk = diag[10?10 , 10?4 ]. For simulation, the initial concentration profile in the SSP reactor is assumed to be unknown. The initial state vector and its covariance matrix for the filter are arbitrarily chosen as follows: ESTIMATOR IMPLEMENTATION FOR THE SSP REACTOR In this section an earlier lumping state estimation method combining OCM with ISR-UKF is applied to the SSP reactor. The SR-UKF is explained to encounter the second problem above, which causes the updated covariance matrix unable to be updated. Next an ISR-UKF-based estimator is implemented and its performance will be evaluated. In the following, we will present the detailed design, especially on how to choose the number of collocation points, initial profiles and the filter parameters, etc. Design and simulation The state estimation of the concentration profiles in SSP reactor is based on the dynamic model described by the Eqns (2)?(5). First of all, the model is transformed to a set of ODEs with the OCM. Here the time?spatial concentration profile C (t, z ) is approximated by the time-varying concentrations Ci (t, zi ) on five internal collocation points and two boundary points, that is N = 5. The internal points are generated from orthogonal polynomials? roots, which are z2 = 0.0469, z3 = 0.231, z5 = 0.769, z6 = 0.953 z4 = 0.5, (32) and the boundary points are z1 = 0, z7 = 1. (33) ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. x?0 = [10?4 10?4 10?4 10?4 10?4 10?4 10?4 10?4 10?4 10?4 P0 = block ? diag{10?10 I6 , 10?4 I6 } 10?4 10?4 ]T (37) Using the standard SR-UKF algorithm to the state estimate model, the program ceases while carrying out the Cholesky factor update in Eqn (16) at the second time step. The error information is that the updated matrix is not positive definite. It is evident that the pause is caused by the second problem described in the Section Calculation of Colesky Factors of the Updated Covariance Matrix. Besides, if ? = ?9, in order to ensure n + ? = 3, the zeroth weight W0(c) becomes negative. So the first problem described in the Section Calculation of Colesky Factors of Predicted State and Measurement Covariance Matrices also probably happens. However, the ISR-UKF avoids the filtering invalidation while applied to the SSP process. The ISR-UKF-based state estimator of SSP reactor achieves success. The start-up process states are first estimated for the reactor from initial empty to gradually filled with solid particles. The observation data corrupted by the noise are shown in Fig. 2, representing the concentration measurements of EG and hydroxyl at the outlet. The states are estimated from the observation serials based on the ISR-UKF. Figure 3 depicts the state and estimation Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering STATE ESTIMATION OF A SSP REACTOR FOR PET x 10-5(mol/O) (a) (b) Observation of hydroxyl concentration Observation of EG concentration 7 6 5 4 3 2 1 0 -1 -2 -3 0 0.5 1 1.5 2 2.5 3 3.5 4 (mol/O) 0.05 0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 0 0.5 1 1.5 t 2 2.5 3 3.5 4 t Figure 2. Plots of the available (simulated) observations of (a) EG and (b) hydroxyl in start-up process. of hydroxyl concentrations on all internal collocation points during this start-up process. It shows the fact that solid particles gradually fill the reactor from the inlet. Although large at beginning due to the arbitrary initial parameter settings of the filter, the estimation error rapidly converges. Even with the observation data highly corrupted by noise (shown in Fig. 2), the filter can overcome the disturbance and produce very good results. Next, states of the dynamic process will be estimated when the operating conditions are changed. In the SSP process, the reaction residence time is one of the operating variables determined by the flow velocity of the pellets. Figure 4 displays the observation data from on-line measurement when the residence time decreases from 30 to 10 h. The estimation results at certain collocation point are shown in Fig. 5. The filter parameters are kept the same. The same conclusion can be drawn on the filter performance as the start-up simulation. The selection of initial state of the filter will be discussed later. x 10-3 (mol/O) z1 Concentration of hydroxyl 18 z2 16 z3 14 z4 z5 12 10 z6 8 z7 6 4 2 0 -2 0 0.5 1 1.5 2 t 2.5 3 3.5 4 The state (solid) and its estimation (dash) of hydroxyl on collocation points during the start-up process. This figure is available in colour online at www.apjChemEng.com. Figure 3. x 10-5 (mol/O) (a) (b) Observation of hydroxyl concentration Observation of EG concentration 3.2 3 2.8 2.6 2.4 2.2 2 0 0.5 1 1.5 2 2.5 3 t (mol/O) 0.018 0.017 0.016 0.015 0.014 0.013 0.012 0.011 0.01 0.009 0 0.5 1 1.5 2 2.5 3 t Figure 4. Observations of (a) EG and (b) hydroxyl (b) when ? decreases to 10 h from 30 h. ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj 385 J. LIU AND X. GU Asia-Pacific Journal of Chemical Engineering (a) (b) x 10-5 (mol/O) 3.5 Concentration of EG on z5 Concentration of hydroxyl on z5 state estimate 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 x 10-3 (mol/O) 16 3 0 0.5 1 1.5 2 2.5 14 10 8 6 4 2 0 -2 3 state estimate 12 0 0.5 1 1.5 t (c) (d) Outlet concentration of hydroxyl 3.1 3 2.9 state estimate 2.8 2.7 2.6 2.5 2.4 2.3 2.2 0 0.5 1 1.5 2 2.5 3 t x 10-5 (mol/O) 3.2 Outlet concentration of EG 2 2.5 3 (mol/O) 0.015 0.0145 simulate estimate 0.014 0.0135 0.013 0.0125 0.012 0.0115 0.011 0.0105 0 0.5 1 1.5 2 2.5 3 t t Figure 5. The states and estimations of EG and hydroxyl concentrations on the fifth collocation points (a) and (b), and the right boundary in (c) and (d) when ? decreases to 10 h from 30 h. Performance evaluation x 10-5 (mol/O) 3 It is well known that the performance of the Kalman filter is largely determined by such factors as initial parameters of the filter, process and measurement noise parameters and uncertain model parameters. Here the performance of the designed estimator is evaluated under the impact of various initial parameters, inaccurate measurement noise parameters and model mismatch. Initial parameters The initial parameters of a filter largely determine the estimation error at the initial stage. The influence will become asymptotically negligible as more and more data are processed. These initial parameters mainly include initial state vector and its covariance matrix. In reality, the initial state should be a zero-vector, not as defined in Eqn (37). x?0 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]T (38) Its covariance matrix decreases to P0 = block ? diag {10?10 I6 , 10?6 I6 }. The estimation results at the third ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. 2.5 Concentration of EG on z4 386 state estimate 2 1.5 1 0.5 0 -0.5 0 0.5 1 1.5 2 t 2.5 3 3.5 4 Figure 6. The EG concentration and its estimate when X?0 is accurately set. internal points are shown in Fig. 6, from which we can see that it is much better estimated. When the operating variables are changed, the unknown initial profile of the state estimator is approximated by incorporating the known inlet concentrations Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering STATE ESTIMATION OF A SSP REACTOR FOR PET x 10-5 (mol/O) (a) Concentration of EG on z5 Concentration of hydroxyl on z5 16 3 simulate estimate 2.5 2 1.5 1 0.5 0 -0.5 -1 x 10-3 (mol/O) (b) 3.5 0 0.5 1 1.5 2 2.5 3 14 state estimate 12 10 8 6 4 2 0 -2 0 0.5 1 1.5 t 2 2.5 3 t Figure 7. The state and its estimate of EG (a) and hydroxyl (b) when X?0 is properly chosen. with the available measurement information. This can be done by a linear correlation of inlet and outlet concentrations. So the initial state vector of the filter in the dynamic simulation is chosen as follows: ?5 ?5 ?5 x?0 = [1.09 О 10 , 1.46 О 10 , 2 О 10 , 2.54 О 10 , 2.91 О 10?5 , 3 О 10?5 , 0.0184, 0.0172, 0.0154, 0.0135, 0.0123, 0.012]T (39) For the same conditions (the residence time decreases from 30 to 10 h), the simulation results are shown in Fig. 7. Compared to Fig. 5 (a) and (b), the justification in the initial parameter setting of the filter can remarkably improve the estimation performance. Measurement noise The performance of a filter largely relies on the process and measurement noise parameters of the system. The performance improves with the increasing of signal-tonoise ratio. In the previous section, we have demonstrated that ISR-UKF can handle the signal with heavy disturbance. However, the use of Kalman filters requires accurate information on the noise distribution. This is rather difficult since disturbance is unpredictable and uncertain. In the EKF estimation, error can be easily accumulated and even diverge owing to the inaccurate noise parameters presumed by the user. In order to evaluate the performance of the ISRUKF under the impact of inaccurate measurement noise parameters, we carried out another trial for Rk = diag[10?12 , 10?6 ] different from the actual value Rk = diag[10?10 , 10?4 ]. Simulation results are depicted in Fig. 8, where ?estimate 1? represents the estimation when the approximate value of noise covariance is taken as Rk = diag[10?10 , 10?4 ] and ?estimate 2? represents the estimation when Rk = diag[10?12 , 10?6 ]. It can be seen that, in the case of a wrong noise parameter, although the estimation error increases at the beginning, it converges faster and tracks the states well. It ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. state estimate1 estimate2 6 Concentration of EG on z5 ?5 x 10-5 (mol/O) 8 4 2 0 -2 -4 0 0.5 1 1.5 2 2.5 3 t Figure 8. The state and its estimates of EG concentration under different noise covariances Rk = diag[10?10 , 10?4 ] (estimate1) and Rk = diag[10?12 , 10?6 ](estimate2). proves that the ISR-UKF works well under inaccurate measurement noise estimation. Model mismatch The above simulations are based on an ideal situation where the exact system model is available to the estimator. This is often unrealistic in industrial practice. So, investigation has been conducted on the performance of the ISR-UKF when the filter model is mismatched to the system. It can happen in reality when uncertainty comes from diffusion coefficient Dz and effective factor ?. First, the effective factor ? is reduced 5%. From Fig. 9 it can be seen that the mismatched effective factor would significantly lower the capability of the estimator. Besides, the estimation error increases along the axial Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj 387 J. LIU AND X. GU Asia-Pacific Journal of Chemical Engineering reaction output and interphase diffusion. The sensitivity of Dz to the observation is smaller than that of ?. Many researchers have carried out sensitivity analysis.[26] This simulation result is in agreement with the fact that uncertain model parameters with various sensitivities affect the Kalman filtering performance differently.[27] Hence, the factors for the model mismatch play different roles on the estimation performance. The effective factor ? lowers significantly the capability of the estimator, whereas the diffusion coefficient Dz has only slight impact. A future work could develop an adaptive filter algorithm or a joint estimation algorithm for the state and parameter ? based on the ISRUKF. x 10-5 (mol/O) Concentrations of EG 3 2.5 z7 2 z4 1.5 z2 1 0 0.5 1 1.5 t 2 2.5 3 CONCLUSIONS Figure 9. The state (solid) and its estimate (dash) of EG concentration when ? reduces 5% at the second, fourth and right boundary points. This figure is available in colour online at www.apjChemEng.com. x 10-5 (mol/O) 3.1 state estimate1 estimate2 3 Outlet concentration of EG 388 2.9 2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1 0 0.5 1 1.5 t 2 2.5 3 Figure 10. The state and its estimate of EG concentration when ? reduces 5% (estimate1) and 10% (estimate2). coordinate of the reactor because the influence of the reduction of ? on the reaction and interphase diffusion gradually increases when the reaction proceeds. If ? reduced further, the estimation error will become larger, as shown in Fig. 10 where ? is reduced by 10%. However, the estimation of hydroxyl is only slightly affected since ? does not directly affect its mass balance. The change of diffusion coefficient shows that the estimator can work well with the mismatch of Dz to the extent of 10 times the actual value because the internal diffusion of EG is much smaller than the ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. This paper presents the design of a novel state estimator for a typical distributed chemical reactor, the continuous SSP reactor for PET. An improved ISR-UKF is developed on the basis of the standard SR-UKF algorithm. Simulation results in both the start-up and dynamic process show that the ISR-UKF used for the state estimation of SSP reactor converges fast and reconstructs the concentration profiles well even with highly contaminated observation data and inaccurate filter parameters. It can be further concluded from the study that the shape of initial profile and accuracy of measurement noise covariance do not affect the long-term convergence behavior of the ISR-UKF, although they have some influence at the initial stage. It is also found that the factors for the model mismatch play different roles in the estimation performance. The effective factor lowers significantly the capability of the estimator, while the diffusion coefficient affects slightly. The successful estimation of the reactant concentration in SSP reactor will be useful as a quality index soft-sensor for production quality control. Acknowledgements This work was supported by National Natural Science Foundation of China (Grant No.60774078), Shanghai Commission of Science and Technology (Grant No.08JC1408200) and Shanghai Leading Academic Discipline Project (Project No. B504). REFERENCES [1] F.K. Mallon, W.H. Ray. J. Appl. Polym. Sci., 1998; 69, 1233?1250. Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering [2] S.N. Vouyiouka, E.K. Karakatsani, C.D. Papaspyrides. Prog. Polym. Sci., 2005; 30, 10?37. [3] C. Algeri, M. Rovaglio. Ind. Eng. Chem. Res., 2004; 43, 4253?4266. [4] G.E. Hovland, T.P. von Hoff, E.A. Gallesteyc, M. Antoined, D. Farruggiod, A.D.B. Paice. Control Eng. Pract., 2005; 13, 1341?1355. [5] M. Boutayeb, D. Aubry. IEEE Trans. Automat. Contr., 1999; 44, 1550?1556. [6] K. Judd. Physica D, 2003; 183, 273?281. [7] S.J. Julier, J.K. Uhlmann, H.F. Durrant-Whyte. A new approach for filtering nonlinear system, in Proceedings of the American Control Conference, Washington, DC 1995; 1628?1632. [8] S.J. Julier, J.K. Uhlmann. Proc. IEEE, 2004; 92, p. 1958. [9] S.J. Julier, J.K. Uhlmann. Proc. IEEE, 2004; 92, 401?422. [10] R. van der Merwe, E.A. Wan. Sigma-point Kalman filters for integrated navigation, in Proceedings of the 60th Annual Meeting of the Institute of Navigation, Dayton, OH 2004. [11] M.C. VanDyke, J.L. Schwartz, C.D. Hall. Unscented Kalman filtering for spacecraft attitude state and parameter estimation, in 2004 AAS/AIAA Space Flight Mechanics Meeting, Maui, Hawaii, February, 2004. [12] T. Kreuziner, M. Bitzer, W. Marquardt. Control Eng. Pract., 2008; 16, 308?320. [13] A. Romanenko, J.A.A.M. Castro. Comput. Chem. Eng., 2004; 28, 347?355. ? 2009 Curtin University of Technology and John Wiley & Sons, Ltd. STATE ESTIMATION OF A SSP REACTOR FOR PET [14] E.A. Wan, R. van der Merwe. Proc. Int. Conf. Acoust. Speech Signal Process., 2001; 6, 3461?3464. [15] K. Xiong, H.Y. Zhang, C.W. Chan. Automatica, 2006; 42, 261?270. [16] M. Athans, R.P. Wishner, A.B. Bertolini. IEEE Trans. Automat. Contr., 1968; 3, 504?514. [17] S.J. Julier, J.K. Uhlmann, H.F. Durrant-Whyte. IEEE Trans. Automat. Contr., 2000; 45, 477?482. [18] B.A. Finlayson. The Method of Weighted Residuals and Variational Principles, Academic Press: New York, 1972. [19] G. Qiu, N.X. Huang, Z.L. Tang, L. Gerking. Chem. Eng. Sci., 1997; 52, 371?376. [20] C.K. Kang. J. Appl. Polym. Sci., 1998; 68, 837?846. [21] F.K. Mallon, W.H. Ray. J. Appl. Polym. Sci., 1998; 69, 1775?1788. [22] J. Liu, X.S. Gu, S.Z. Zhang. J. Chem. Ind. Eng. (China), 2008; 59, 1727?1731. [23] L.T. Fan, G.K.C. Chen, L.E. Erickson. Chem. Eng. Sci., 1971; 26, 379?387. [24] P.G. Kaminski, A.E. Bryson. IEEE Trans. Automat. Contr., 1971; AC-16(6), 727?735. [25] R. Battin. Astronautical Guidance, McGraw-Hill Book Company: New York, 1964. [26] R.E. Griffin, A.P. Sage. AIAA J., 1969; 7, 1890?1897. [27] K.T. Koussoulas, C.T. Leondes. Int. J. Syst. Sci., 1986; 17, 937?941. Asia-Pac. J. Chem. Eng. 2010; 5: 378?389 DOI: 10.1002/apj 389

1/--страниц