A Method of Online Safety Assessment for Industrial Process

ARTICLE pubs.acs.org/IECR

A Method of Online Safety Assessment for Industrial Process Operations Based on Hopf Bifurcation Analysis Lubin Ye, Zhengshun Fei, and Jun Liang* State Key Laboratory of Industrial Control Technology, Institute of Industrial Control, Zhejiang University, Hangzhou 310027, People’s Republic of China ABSTRACT: This paper proposes an approach of online safety assessment for industrial process operations. The core idea is to apply the Hopf bifurcation analysis of nonlinear systems in safety assessment. Herein, the Hopf bifurcation equations are used to characterize the critical stability boundary or manifold of the system, which is taken as the safety limit. Then, a safety index (SI) is constructed to denote the integrated exponential distances between each parameter and their bifurcation points. In the online implementation, the parameters are estimated by the extended Kalman filter (EKF), and the Hopf bifurcation points are generated by solving the nonlinear bifurcation equations. Afterward, the value of the SI can be online calculated. The introduced approach was then applied to a simulated gas phase polyethylene reactor process, in which the efficiency of the proposed method was verified in indicating the distance to the potential unsafe oscillation and in the early identification of potential threats.

1. INTRODUCTION Industrial processes are often operated under high temperatures or pressures with complex physical changes or chemical reactions. Improper operations in industrial processes may cause economic loss, harm the environment, or even threaten human health. Thus, safety is crucial for industrial processes and has attracted tremendous attentions from both academia and industry. To obtain the online quantitative characterization of process safety, some safety/security indices were proposed.1,2 The indexbased approaches are actually space mappings, which relate the output safety/security index with the input operation pattern. However, the probabilistic safety index1 for multimode processes was under the premise that critical safety limits have already been obtained. And since the safety assessment was obtained by a datadriven method, the consequent safety level generation was restricted in historical operation regions. The security index2 used the reactor temperature as the indication of security and then derived the critical points for the operation parameters, which categorize the operation space into several security zones. It was later pointed out3 that this approach requires extensive, time-consuming dynamic simulations and is impractical for largerscale processes. Some studies focused on the safety regions of manipulated variables when analyzing operation safety.4-9 They used the concept of dynamic security region (DSR), the area in the variable space where the system is still transient stable when it suffers from a given fault or disturbance. In these approaches, the characterization of DSR boundary6-9 is crucial. To achieve this, it is required that all the system’s unstable equilibriums are found, which is extremely difficult for high-order systems. Therefore, the applications of DSR methods were mainly in the power systems, and its use in process system engineering was scarcely reported. However, the idea of using system stability property in analyzing operation safety can be utilized for industrial processes. As nonlinear dynamic systems, the industrial processes, especially r 2011 American Chemical Society

chemical reactors, may have complex stability characteristics. A real accident in an ammonia synthesis reactor with a rapid temperature oscillation was reported, and it was concluded that the system instability was caused by the decrease of the reactor pressure.10,11 The calculation results showed that a pair of complex conjugate poles crossed the imaginary axis (a Hopf bifurcation occurred) when the pressure decreased to a low value. In the study of operating conditions in a polyethylene reactor, it was found that instability happened from a Hopf bifurcation point, which resulted in an oscillation of the reactor temperature and runaway toward unacceptable values.12 Literature surveys revealed that there exist bifurcation phenomena in many chemical or biochemical reactors.13-17 Some researchers applied bifurcation analysis when considering process safety,18-20 in which they conveyed that the Hopf bifurcation might result in unstable oscillations of the process and thus it was a safety threat. However in industrial processes, safety issues of the process operation are much more comprehensive than the Hopf bifurcations and oscillations. Potential risks of leakage, toxicity, fire, explosion, and so on should be considered in safety assessment, which is also related to loss prevention, hazard identification, risk evaluation, and reliability analysis.21,22 Practical process accidents may be caused by unexpected disturbances or operation changes that result in the runaways and safety violations of key variables like pressure and temperature. Hopf bifurcation is one type in all possible causes of unsafe conditions. It may cause unstable oscillation and lead to disastrous accidents sometimes, but at other times might be tolerable. Whatever, the oscillation behavior is an unexpected condition, and the Hopf bifurcation can be considered as one of the potential threats to operation safety. In this paper, the concerns of process operation safety are focused Received: September 11, 2010 Accepted: January 20, 2011 Revised: January 20, 2011 Published: February 18, 2011 3403

dx.doi.org/10.1021/ie1018854 | Ind. Eng. Chem. Res. 2011, 50, 3403–3414

Industrial & Engineering Chemistry Research on Hopf bifurcation phenomena, and the safety assessment is based on the Hopf bifurcation analyses, as well. The usual method of practical numeric bifurcation analysis involves an iterative procedure known as the “continuation” algorithm,23,24 which can be implemented by some software packages, like AUTO25,26 and MATCONT.27,28 The “continuation” is a time-consuming algorithm and thus the above safety analysis approaches18-20 based on bifurcations only generated the equilibrium trajectory with respect to one parameter. If more parameters need to be considered simultaneously, the work will become intractable due to the extensive computation burden. Therefore, the “continuation” algorithm is inappropriate for online use. The critical value for one parameter in the above bifurcation analyses can be seen as the stable boundary or safety limit of the parameter. To characterize stable boundary for more parameters, a group of Hopf bifurcation equations were used.29-34 These equations showed that there might be infinite bifurcation points, locating in a high-dimension curved surface (critical manifold). This manifold limits the parameter stable region and was utilized to ensure system stability, robustness, and flexibility in process designing and optimization.30-32,34 Compared with the onedimension bifurcation analyses, 12-20 the method based on bifurcation equations is able to characterize the bifurcation manifold of multiple parameters. The bifurcation values can be obtained by solving these algebraic equations, which is obviously faster than the continuation algorithm. Therefore, the method based on the bifurcation equations is more appropriate for online calculation. In this paper, an online process safety assessment method is developed based on Hopf bifurcation analysis. Herein, a safety index (SI) is constructed to measure the distances between the current operation point and the critical manifold of Hopf bifurcations (the stable boundary). The region inside the manifold is stable, but the SI quantity is not uniformly distributed. Those points nearer to the manifold are regarded to have less value of SI, and vice versa. This index is an integrated function of the exponential parametrical distances, which is designed to be continuous, smooth, and decreasing monotonically with any distance. Although the constructed SI does not include all safety issues concerned in practical operation, it quantifies the risk of Hopf bifurcation and the subsequent unexpected oscillations. Since the industrial process parameters are drifting during practical running and may affect the location of the bifurcation manifold and system stability properties, they are estimated by an extended Kalman Filter (EKF), which was described detailedly in the literature.35-39 An online calculation of the bifurcation points is implemented by solving the algebraic bifurcation equations, which is actually a nonlinear programming (NLP) problem. The gas phase polyethylene (PE) reactor12,13,17 is a typical industrial process that has bifurcation phenomena, which requires safety considerations, so it is taken as the simulation example. The remainder of the paper is organized as follows. The preliminaries of EKF and equations for critical manifold of Hopf bifurcation are introduced in section 2. The proposed approach, including solving for critical Hopf bifurcation points, the safety index (SI) construction, SI sensitivity analysis, and the online safety assessment implementation procedures based on SI are presented in section 3. The efficiency of the approach is testified in a simulated gas phase polyethylene reactor example in section 4 along with result discussions. Lastly, a conclusion is given.

ARTICLE

2. PRELIMINARIES 2.1. Extended Kalman Filter. For a nonlinear system, the mathematical model can be represented by

_x ¼ f ðx, u, θÞ y ¼ gðx, u, θÞ

ð1Þ

where x ∈ Rn is the vector of state variables, u ∈ Rm is the vector of manipulated variables, θ ∈ Rl is the vector of process parameters and y ∈ Rc is the vector of measurable output variables; f( 3 ) and g( 3 ) are vector functions with appropriate dimensions. The extended Kalman Filter (EKF) is a recursive algorithm based on local linearization of the system model (1) to provide the minimum-variance state estimates.38 If the drift of the parameters is very slow compared with the estimation frequency, the parameter can be seen as being constant at every estimation step. The unknown parameters are then regarded as the extended state variables; the process model is then modified as35,36 " " # f ðx, u, θÞ _x _z ¼ _ ¼ 0 θ y ¼ gðx, u, θÞ

ð2Þ

Here, we denote z = [xT,θT]T as the expanded vector of state variables. To make sure all unmeasured parameters are observable without estimation bias, it is required that the number of estimated parameters does not exceed the number of measured output variables.36 After linearization and discretization, the EKF formulations of the system (2) are as follows.35,36 Model prediction. 2 3 ^xðk þ 1jkÞ 5 ^zðk þ 1jkÞ ¼ 4 ^ θðk þ 1jkÞ 2

^xðkjkÞ þ ¼4 ^ θðkjkÞ

R ðk þ 1ÞTs kTs

^ f ðxðtÞ, uðkÞ, θðkjkÞÞ dt

R

3 5 ð3Þ

where Ts is the sampling interval of EKF and the term kTs(kþ1)Tsf(x(t), u(k), θ̂(k|k)) dt is obtained by numerical integration. Measurement Correction. 2 3 2 3 ^xðk þ 1jk þ 1Þ ^xðk þ 1jkÞ 5¼4 5 ^zðk þ 1jk þ 1Þ ¼ 4 ^ ^ þ 1jkÞ θðk þ 1jk þ 1Þ θðk ^ þ 1jkÞÞ þ Kðk þ 1Þ½yðk þ 1Þ - gð^xðk þ 1jkÞ, uðkÞ, θðk

ð4Þ

where K(k þ 1) is the Kalman gain matrix as Kðk þ 1Þ ¼ R z ðk þ 1jkÞΨk þ 1 T ½Ψk þ 1 R z ðk þ 1jkÞΨk þ 1 T þ R 2 -1

ð5Þ

Here, R2 is the covariance of measurement noise and Rz is the estimated covariance matrix of state variables, which are calculated iteratively by R z ðk þ 1jkÞ ¼ Φk R z ðkjkÞΦk T þ R w R z ðk þ 1jk þ 1Þ ¼ ½I - Kðk þ 1ÞΨk þ 1 R z ðk þ 1jkÞ 3404

ð6Þ ð7Þ

dx.doi.org/10.1021/ie1018854 |Ind. Eng. Chem. Res. 2011, 50, 3403–3414

Industrial & Engineering Chemistry Research

ARTICLE

where Rw is the covariance matrix of state disturbance. Here, " # Mk Dk Φk ¼ ð8Þ 0 I Ψ k ¼ ½ Hk

Fk

A 1 3 ðω 1 þ jω 2 Þ ¼ jλðω1 þ jω 2 Þ

ð9Þ

R where Mk = exp(AkTs) and Dk= kTs(kþ1)Tsexp[Ak 3 ((k þ 1)Ts - τ)] dτ Bk. The matrices Ak, Bk, Hk, and Fk are obtained by linearization of the system equations (2) at the sampling instant, generated by38 ∂f Ak ¼ ∂x ^ ^ x ¼ xðkjkÞ, θ ¼ θðkjkÞ , u ¼ uðkÞ ∂f Bk ¼ ∂θ ^ ^ x ¼ xðkjkÞ, θ ¼ θðkjkÞ , u ¼ uðkÞ ∂g Hk ¼ ∂x ^ ^ x ¼ xðkjk - 1Þ, θ ¼ θðkjk - 1Þ, u ¼ uðkÞ ∂g Fk ¼ ð10Þ ∂θ ^ ^ x ¼ xðkjk - 1Þ, θ ¼ θðkjk - 1Þ, u ¼ uðkÞ The EKF procedures of eqs 3-10 provide online estimates for both states and parameters of the system in eq 1. 2.2. Critical Manifold of Hopf Bifurcation. Bifurcation is a particular feature of nonlinear dynamic systems, differing from linear systems. When parameters vary, the properties of system equilibrium may have a sudden change, causing phenomena like multiple steady states, hysteresis behaviors, limit cycles, or chaos.15,30 To be brief, the appearance of a topologically nonequivalent phase portrait under variation of a parameter is called a bifurcation.15,23,24 There are many types of bifurcations, including fold, Hopf, cusp, isola.15,30 Hopf and fold are the most common bifurcations present in nonlinear systems.15 Hopf bifurcations are responsible for instable oscillation phenomena and fold bifurcations are usually the causes of multiple steady states and hysteresis behaviors. For chemical reactors, Hopf bifurcations were more frequently reported and have more severe hazards. Therefore, the current work focuses on Hopf bifurcations in process safety assessment. Now consider the nonlinear system of eq 1. Since there is no difference between u and θ in mathematics, they are both taken as the system parameters. Thus, an extended parameter vector is adopted that β = [uT θT]T, and then the state equation of system (1) can be rewritten by x_ ¼ f 1 ðx, βÞ ð11Þ A Hopf bifurcation point (x*, β*) satisfies the following equations:30,32-34 f 1 ðx, βÞ ¼ 0 ð12Þ A 1 ω 1 þ λω 2 ¼ 0

ð13Þ

A 1 ω 2 - λω 1 ¼ 0

ð14Þ

ω1 T ω1 þ ω2 T ω2 - 1 ¼ 0

ð15Þ

ω1 ω2 ¼ 0

ð16Þ

T

corresponding eigenvalue (the real part is zero). Equation 12 assures that any Hopf bifurcation point is the equilibrium of the system, eqs 13 and 14 are derived from the eigenvector relation:

where A1 = ∂f1/∂x|x=x*,β=β* is the Jacobian matrix of f1. The vectors ω1 and ω2 are the real and imaginary part of eigenvector of A1 when x = x*, β = β*, and λ is the imaginary part of the

and eqs 15 and 16 normalize the eigenvector. It can be easily obtained that x* ∈ Rn, β* ∈ Rmþl, ω1 ∈ Rn, ω2 ∈ n R , and λ ∈ R1. So, there are (3n þ m þ l þ 1) variables. Equations 12-16 totally have (3n þ 2) equations. Therefore, the degree of freedom (DOF) is (m þl - 1), which implies that the Hopf bifurcation points are located in a (m þ l - 1)-dimensional manifold.

3. ONLINE PROCESS SAFETY ASSESSMENT UTILIZING HOPF BIFURCATION ANALYSIS In this section, a safety index (SI) is proposed to measure the distance of the current operation point to the Hopf bifurcation manifold, which evaluates the risk of potential bifurcation phenomena. An online implementation of safety assessment is then formulated, consisting of online parameter estimation, bifurcation points solving, and SI calculation. The requirements of the application for the proposed approach are introduced as follows: (1) The process model is correct without mismatch. (2) The process sensors of output variables are functioning well without faults and the measurements are unbiased. (3) The functions in eqs 12-16 are continuous and differentiable, and their derivatives are also continuous at most interested states. (4) The parameters are drifting much more slowly than the frequency of EKF estimation and SI calculation. So, the parameters are regarded as constants in one period of SI implementation. The first and the second assumptions are required in the EKF estimation and bifurcation analysis. Further, the current work is based on the locations of bifurcation points, which are obtained through the model equations. Thus, these two assumptions are necessary for the applicability of the model-based method. Assumption (3) assures that eqs 12-16 are quite smooth, which makes it easy to find the bifurcation root for these equations. Assumption (4) implies that at each online stage of SI calculation, the parameters are approximate constants and the constructed SI is not involved with dynamic problems, and this condition is also necessary in the EKF for parameter estimation. 3.1. Online Solving for Critical Hopf Bifurcation Values of Parameters. To obtain the online safety assessment, the para-

meters of process system in eq 11 have to be estimated by the EKF in section 2.1. Some parameters may be the manipulated variables or measurable disturbance, the values of which are directly provided without estimation. At every online implementation instance, in order to get the distances from the current value of each parameter to the critical manifold, the Hopf bifurcation points of each parameter have to be obtained. In traditional bifurcation analyses,23-28 the continuation algorithm is usually adopted, in which the computation cost is very large and the computation burden increases exponentially with the dimension of considered parameters. However in the current work, the Hopf bifurcation values are obtained by solving bifurcation equations 12-16, which is a problem of nonlinear least-squares. 3405



ARTICLE

To find the root, DOF of the equations should be 0, and (m þ l - 1) free variables have to be configured. Since the parameter drifts are very slow, the difference between the values of two adjacent instances is small. Thus, when searching for the Hopf bifurcation point of the ith parameter (denoted by β*) i at an instance, the other (m þ l - 1) parameters are assigned to be at their currently estimated values. Then, eqs 12-16 can be rewritten by f 1 ðx, ½β1 , :::, βi - 1 , βi , βi þ 1 , :::, βl þ m T Þ ¼ 0 A 1 ω 1 þ λω 2 ¼ 0 A 1 ω 2 - λω1 ¼ 0 ω1 T ω1 þ ω2 T ω2 - 1 ¼ 0 ω1 T ω2 ¼ 0 ^ ðtÞ β1 ¼ β 1 ::: ^ βi - 1 ¼ β i - 1 ðtÞ ^ βi þ 1 ¼ β i þ 1 ðtÞ ::: ^ βl þ m ¼ β l þ m ðtÞ

ð17Þ

^1(t), ..., β î-1(t), β îþ1(t), ..., β ^lþm(t) are the estimated where β values of parameters at time t, or the directly measured values of input variables and disturbances. In eqs 17, DOF = 0 and the critical Hopf bifurcation value β*i can be determined by the nonlinear least-squares problem, which is actually an application of nonlinear programming (NLP). Herein, the Gauss-Newton method is used here to search for the solution. The solving procedures can be easily executed in some computation software like MATLAB. To reduce the searching time, the root generated in the previous step is taken as the initial point for the Gauss-Newton algorithm of the current step. The Hopf bifurcation values of other (m þ l - 1) parameters can be generated in the same way. Under the assumptions, the equations in eqs 17 are continuous and differentiable, and the derivatives are also continuous at most interested states in space. So the objective function of the nonlinear least-squares problem is quite smooth, which is very beneficial to the numerical solving method of NLP. Since the parameter drift is very slow in every step, the NLP’s optimal point is not far away from the previous optimality result, which is taken as the initial point in the current step. So, this implementation avoids the local minimum problem to a great extent. Therefore, the searching algorithm can easily find the optimal result in every step of NLP. In the studied experiments, it was found that the Gauss-Newton algorithm always gained the optimal points. This approach by solving algebraic equations is much faster than the continuation algorithm and is more appropriate for online use when considering multiple-parameter drift since the impact of the parameter dimension effect on the NLP computation burden is not so obvious. 3.2. Safety Index (SI) Construction. In the existing indexbased approach of operation safety assessment, the safety limit or critical boundary is usually adopted.1,2 Most safety limits were determined directly from the process’s physical or chemical constraints, for example, tolerated safe ranges of temperatures, pressures, or liquid levels. In the current work, dynamic features of the process system is involved, and the concept of “safety” is specifically focused on the issue of Hopf bifurcation, which is a type of unsafe threat to the process operation. The above-

discussed Hopf bifurcation manifold is then utilized here as the safety limit. The following presented safety index (SI) is a normalized measure of the spatial distance to the Hopf bifurcation manifold. The original spatial distance between the current value and the Hopf bifurcation point of a parameter is inappropriate for assessment since those parameters denote different physical or chemical quantities and have diverse dimensions. To get a uniform measure of distance, all parameters have to be scaled to be dimensionless and the distance measures are also expected to be normalized to the range of [0, 1]. For the ith parameter, the î(t) and β*i at time t is formulated as distance between β ! " # ^ ðtÞ - β Þ2 vi 2 ðβ i i ¼ 1 - exp ð18Þ si ¼ 1 - exp hi σ i 2 hi σ i 2 î(t) and βi*, and where vi is the original spatial distance between β 2 σi is the estimation variance of βi; hi is a positive scaling factor, which determines the curve shape of the exponential function. Remark 1. Herein, the negative exponential function is formulated to characterize the distance because it can restrict the value range to [0, 1]. Further, it is continuous, smooth, and monotonic, and these properties remain at any order of derivative. It can be concluded that the value of si in eq 18 approaches 1 î(t) is far from its critical point βi*, and the value when β î(t) is close to βi*. approaches 0 when β The exponential expression in eq 18 only provides the distance measure for one parameter. In order to give a general assessment on the distance of process operating condition to the Hopf bifurcation, a safety index (SI) has to be constructed by integrating all si of the (m þ l) parameters. Herein, the SI is formulated as the weighted geometric mean of the (m þ l) exponential distances, as mY þl si wi ð19Þ SI ¼ i¼1

where wi is the weighting factor for si, which denotes P the importance of the ith parameter in safety and satisfies i = mþl wi = 1, "i, 0 e wi e 1. If all (m þ l) parameters have the 1 1/(mþl) . same weight, wi = 1/(m þ l) and SI = (Πmþl i = 1si) Remark 2. The SI constructed in eq 19 has the following properties: (1) The value of SI is in the range of [0, 1]. î(t) = β*, (2) If any si = 0 or β i which means βi reaches its Hopf bifurcation point, the value of SI = 0, and if "i, vi f ¥, SI = 1. (3) The value of SI increases monotonically when any vi increases. (4) The value of first order partial derivative of SI with respect to any squared distance vi2, that is, ∂SI/∂(vi2), increases monotonically when vi decreases. Properties (1 and 2) are obvious. Property (1) assures that the value of SI is in the normalized range of [0, 1], which is necessary for giving a general assessment. If any parameter βi reaches its bifurcation point (Hopf), the process system will have an unstable oscillation, which is considered as an unsafe condition in practical operations. Contrarily, when all parameters are far from their critical bifurcation points, the process operation is regarded to be “quite safe”. From property (2), it can be seen that the worst and best safety conditions correspond with the case of SI=0 and SI=1, respectively. 3406



ARTICLE

Properties (3-4), which are proved in the Appendix, show the monotonic features of the constructed SI. Property (3) reveals that an increment in any spatial distance vi will help increasing the SI value, which is expected in intuitive cognitions. Property (4) demonstrates that the marginal decrement in SI with respect to Δ(v2i ) becomes more obvious when the vi value is closer to 0. This property implies that the SI curve decreases more and more steeply when the operating point approaches the critical manifold. In most plant operators’ minds, the safety degree (denoted by the SI in this article) should be more sensitive to the same perturbation in the parametric distance when the process is closer to the critical manifold. Contrarily, when the process operating point is far from the critical manifold, the safety degree should be less affected by spatial distance variations. So, property (3 and 4) of the SI well emulates the human’s cognitions in safety characterizations. In eq 18, hi is the scaling factor that determines the curve shape of si. If hi value is too large, the curve variation of eq 18 is too mild and si has insufficient sensitivities to the variations of vi. If hi value is too small, the curve variation is too severe and si is oversensitive to vi. It is expected the value of hi is properly chosen that the value of si has an appropriate variation range when vi drifts in a reasonable scope. Practically, the state that process parameters are at their nominal values is considered to be an ideal safe condition. So under the nominal condition, the value of si is expected to be at a predefined criterion, which should be a large value close to 1. Herein, the criterion is chosen to be 0.99, which means that the nominal parameter conditions are deemed to have all si = 0.99. It can be demonstrated by the curve of eq 18 that si covers a proper range for efficient SI assessment. î(t) with its nominal value β0i (t) in Substituting the parameter β eq 18, one can get " 1 - exp

ðβ0 ðtÞ - βi , 0 Þ2 - i hi σ i 2

l X ∂log10 ðSIÞ ∂log10 ðSIÞ ∂si ∂log10 ðSIÞ ∂sj ¼ ¼ ^ ^ ^ ∂s ∂sj i ∂βj ∂ βj ∂β i¼1 j l ^ - β Þ X ð1 - sj Þ 2ðβ ∂log10 ðSIÞ ∂si 1 j j þ ¼ 3 wj 3 3 hσ2 ^ logð10Þ ∂s s i j j j ∂ βj i ¼ 1, i6¼ j

þ

^ , :::, β ^ ^ ^ Fðx , βi , ω1 , ω 2 , λ, β 1 i - 1 , βi þ 1 , :::, βl þ m Þ 3 2 T ^ , :::, β ^ ^ ^ f 1 ðx, ½β 1 i - 1 , βi , βi þ 1 , :::, βl þ m Þ 7 6 7 6 A 1 ω 1 þ λω 2 7 6 7 6 A ω λω ¼6 1 2 1 7 7 6 ω Tω þ ω Tω - 1 2 2 5 4 1 1 T ω1 ω2

ð20Þ

ðβ0i ðtÞ - βi , 0 Þ2 σi 2 3 lnð0:01Þ

ð23Þ

^ , :::, β ^ ^ ^ Fðx , βi , ω 1 , ω 2 , λ, β 1 i - 1 , βi þ 1 , :::, βl þ m Þ ¼ 0 ð24Þ ^j can be formulated as The partial derivative of F with respect to β ∂F ∂x ∂F ∂βi ∂F ∂ω1 ∂F ∂ω 2 þ þ þ ^ ^ ^ ^ ∂x ∂β ∂ω ∂ω 1 ∂ βj 2 ∂β ∂βi ∂βj j j

∂F ∂λ þ ^ ∂λ ∂β j

l X k ¼ 1, k6¼ i

^ ∂F ∂β k ¼0 ^ ∂β ^ ∂β k

ð25Þ

j

^k, k 6¼ j is not affected by with β ^j, Since the estimated parameter β ^ ^ ∂βk/∂βj = 0, k 6¼ j. Thus the above equation can be written as

∂F ∂x ∂F ∂βi ∂F ∂ω 1 ∂F ∂ω 2 ∂F ∂λ ∂F þ þ þ þ þ ¼0 ^ ^ ∂β ^ ^ ^ ^ ∂x ∂β ∂ω ∂ω ∂λ ∂β 1 ∂βj 2 ∂β ∂βi ∂βj j j j j

ð26Þ

ð21Þ or

As presented above, the constructed SI of eq 19 is a normalized measure of closeness to Hopf bifurcation phenomena, which will result in unstable oscillations in most cases and is one of the potential safety threats of process operation. There might be other causes or disturbances that bring unsafe conditions or accidents. However, it is impractical to integrate all concerned safety issues in one formulation, and the current work based on the SI only focuses on the risk of Hopf bifurcation. Since the SI in eq 19 is related to the estimated and bifurcation value of the process parameters, the calculation result is dependent on the parameter values and might be affected by some uncertain factors, like biased estimation or measurements. A sensitivity analysis is introduced in order to give the SI reliability with respect to those changes. For the SI, multiplicative changes are concerned because it is a measure of exponential distance. The partial derivative of log10(SI)

ð22Þ

then F( 3 ) ∈ R3nþ2 is a vector of functions, and eqs 17 can be rewritten by

þ ¼ 0:99

^ Þ∂β 1 ð1 - si Þ 2ðβi - β i i wi 3 3 3 2 ^ logð10Þ s h σ i i i ∂ β i ¼ 1, i6¼ j j l X

^j is determined by relations in eqs 17. Let Here the term ∂βi*/∂β

#

where β*i ,0 is the bifurcation value of βi0(t) under the nominal condition. Subsequently, it is hi ¼ -

^j is generated as respect to β

3 ∂x 7 6 ∂β 6 ^j 7 6 7 6 ∂β 7 6 i 7 6 ^ 7 7 6 " #6 ∂β j 7 6 ∂F ∂F ∂F ∂F ∂F 6 ∂ω 1 7 7 ∂F , , , , ¼0 6 ^ 7þ ^ 6 ∂x ∂βi ∂ω 1 ∂ω2 ∂λ ∂β 7 6 j 7 ∂ βj 6 ∂ω 2 7 7 6 6 ^ 7 6 ∂β j 7 7 6 6 ∂λ 7 5 4 ^ ∂β 2

ð27Þ

j

where [∂F/∂x*, ∂F/∂βi*, ∂F/∂ω1, ∂F/∂ω2, ∂F/∂λ] is actually the jacobian matrix of F at the root β*. i Denote Ji = [∂F/∂x*, ∂F/∂β*, i 3407


Industrial & Engineering Chemistry Research ∂F/∂ω1, ∂F/∂ω2, ∂F/∂λ]. It can be generated that 2 ∂f 1 0 0 6 A1 ∂β 6 i 6 6 ∂ðA 1 w 1 Þ ∂ðA 1 w 1 Þ 6 A1 λI 6 ∂x ∂βi 6 Ji ¼ 6 ∂ðA w Þ ∂ðA w Þ 1 2 1 1 6 -λI A 1 6 6 ∂x ∂β i 6 60 0 w1 T w2 T 4 0 0 w2 T w1 T

ARTICLE

3 0 7 7 7 7 w2 7 7 7 7 7 w1 7 7 7 0 7 5 0

If Ji is full-ranked (easily satisfied in practical applications), it can be generated according to the implicit function theorem (IFT): 20 1 0 1 0 1T 0 1T 0 1T 3T T T 6@∂x A @∂βi A @∂ω 1 A @∂ω 2 A @ ∂λ A 7 , , , , 4 5 ^ ^ ^ ^ ^ ∂β ∂β ∂β ∂β ∂β j

¼ - J-1 i

j

j

j

j

Figure 1. Schematic of an industrial gas phase polyethylene reactor.

∂F ^ ∂β

j

ð28Þ ^j is the (n þ 1) element of the obtained Then the value of ∂βi*/∂β ^j in eq 22 with the result of eq 28, vector. Substitute the ∂βi*/∂β ^j can be calculated. Further, and then the sensitivity ∂log10(SI)/∂β these partial derivatives are scaled to be dimensionless. To keep the value bias of the obtained SI in acceptable multiplicative bias range, it is required that ∂log ðSIÞ 10 ^ e Tol ð29Þ Δβ Δlog10 ðSIÞ j ^ ∂β j where Tol is the acceptable change in log10(SI). So, the tolerance ^j, denoting the reliability of SI with of perturbation (TOP) of β ^ respect to βj, can be generated as 2 3-1 ∂log ðSIÞ 10 5 Tol TOP ¼ 4 ð30Þ ^ ∂β j

Since the constructed SI is actually an exponential measure of distance, SI is insensitive when its value is close to 1 and becomes very sensitive when it is close to 0, as implied in property (4) of SI. Therefore, the TOP decreases steeply when SI value approaches 0, which is confirmed in the case study results. It is then suggested that the logarithm of TOP (LTOP) is used to denote the SI reliability in practical applications. 3.3. Online Implementation of Process Operation Safety Assessment Based on the SI. Before the online safety assessment implementation, the initial Hopf bifurcation points for each parameter have to be decided. Herein, this work can be done in the software package MATCONT27,28 in MATLAB. In real industrial processes, the process parameters may drift with time. But their variation rate is much slower than the process dynamics. Since the Hopf bifurcation points determined by eqs 17 are steady state points of process system in eq 11 and the SI formulation in eqs 18 and 19 is also based on the steady state values, the sampling frequency of solving for bifurcation points in eqs 17 and the consequent safety assessment should not be too fast. The sampling interval of safety assessment could be

empirically decided according to the practical dynamic features of the process. Generally, the proposed online implementation procedures of process operation safety assessment based on the constructed SI are as follows: (1) Obtain the nominal values of process system parameters (including the process parameters, disturbances and manipulated variables). (2) Under the nominal condition, find the initial Hopf bifurcation values of those parameters through the MATCONT software package in MATLAB. (3) Using the nominal values of parameters and their initial Hopf bifurcation points, determine the scaling factor hi for each exponential distance measure si by eq 21. (4) Begin the online implementations. Estimate the system parameters through EKF in eqs 3-10. (5) At every assessment sample, derive the Hopf bifurcation point for each parameter by solving eqs 17 through NLP algorithms. (6) Obtain the SI value through eqs 18 and 19, which provides a general safety assessment on the current operation condition.

4. CASE STUDIES The proposed online safety assessment method was applied to a gas phase polyethylene (PE) reactor process. Figure 1 shows the schematic of a gas phase PE reactor process.12,13,40,41 The feed to the reactor consists of ethylene, comonomer, hydrogen, inerts, and catalysts. A recycle stream of unreacted gases flows from the top of the reactor, and is then cooled in a heat exchanger in counter-current flow with cooling water. The cooling rates in the heat exchanger are adjusted by mixing cold and warm water streams while maintaining a constant total cooling water flow rate through the heat exchanger. The reactor process model was formulated based on mass and heat balances.12 It was pointed out the mass balances on hydrogen and comonomer were not considered as they have mild effects on the reactor dynamics.12 The heat exchange process was then modified to be characterized by a log-mean temperature difference (LMTD) model.13 Then, the mathematical model for 3408



ARTICLE

Table 1. State Variables of PE Reactor state variables

descriptions

nominal values

units

[In]

molar concentration of inerts in the gas phase

439.68

mol/m3

[M1]

molar concentration of ethylene in the gas phase

326.72

mol/m3 mol

Y1

moles of active site type 1

3.835

Y2

moles of active site type 2

3.835

mol

T

reactor temperature

356.21

K

Tw1

temperature of cooling water stream from exchanger

290.37

K

Tg1

temperature of recycle gas

294.36

K

the PE reactor has the following form:13,40,41 d½In ¼ dt d½M1 ¼ dt

FIn -

FM1

Table 2. Manipulated Input Variables of PE Reactor

½In bt ½M1 þ ½In Vg

input variables

½In bt - RM1 ½M1 þ ½In Vg

dY1 RM1 MW1 Y1 ¼ Fc ac - kd1 Y1 dt Bw dY2 RM1 MW1 Y2 ¼ Fc ac - kd2 Y2 dt Bw

ð31Þ

Hf þ Hg1 - Hg0 - Hr - Hpol dT ¼ dt Mr Cpr þ Bw Cppol dTw1 Fw UA ¼ ðTwi - Tw1 Þ ðTw1 - Tg1 Þ Mw Cpw dt Mw dTg1 Fg UA ¼ ðT - Tg1 Þ þ ðTw1 - Tg1 Þ Mg Cpg dt Mg where bt ¼ Vp Cv

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð½M1 þ ½InÞ 3 RR 3 T - Pv

RM1 ¼ ½M1 3 kp0 exp Cpg ¼

-Ea 1 1 3 ðY1 þ Y2 Þ R T Tf

½M1 ½In Cpm1 þ CpIn ½M1 þ ½In ½M1 þ ½In

Hf ¼ FM1 Cpm1 ðTf eed - Tf Þ þ FIn CpIn ðTf eed - Tf Þ

ð32Þ

Hg1 ¼ Fg ðTg1 - Tf ÞCpg Hg0 ¼ ðFg þ bt ÞðT - Tf ÞCpg Hr ¼ Hreac MW1 RM1 Hpol ¼ Cppol ðT - Tf ÞRM1 MW1 The reactor process system has seven states, that is, [In], [M1], Y1, Y2, T, Tw1, and Tg1. The five flow rate FIn, FM1, Fc, Fw, and Fg

descriptions

values

units

FIn

flow rate of inert feed

5

mol/s

FM1

flow rate of ethylene feed

190

mol/s

Fc

flow rate of catalyst feed

5.8/3600

kg/s

Fw

flow rate of cooling water

(3.11 105)

kg/s

Fg

flow rate of recycle gas

(18 10-3) 8500

mol/s

Vp

bleed stream valve position

0.5

and the bleed stream valve Vp can be manipulated by operators and are taken as the input variables. Further, the three temperatures and the overhead bleed flow rate bt in eqs 32 are regarded to be measurable and taken as the system’s output variables. Other variables are the process parameters or intermediate variables. The descriptions, nominal values, and units of the states, input variables, parameters, and intermediate variables40 are listed in Tables 1-4, respectively. Herein, the reactor model in eqs 31 and 32 consists of differential algebraic equations (DAEs) and it can be easily expressed by ordinary differential equations (ODEs) using the implicit function theorem (IFT).34 A simulation module of the PE reactor process was constructed in the SIMULINK platform of MATLAB. The values for some parameters, including Vg, R, RR, MW1, Cpm1, Cpw, CpIn, Cppol, Mw, Mg, MrCpr, Hreac, Tf, kp0, and Ea are determined by chemical/physical properties of the substances or the process equipments, and they are thus constants and not considered in the bifurcation analysis. Bw (mass of polymer in the fluidized bed) is dependent on other factors and is also excluded. Those parameters of Pv, Cv, UA, kd1, kd2, ac may drift with time when the operation conditions change or some faults occur. And Tfeed and Twi, that is, the temperature of feed flow and cooling water inlet flow, are regarded as the measurable disturbances, which may also change. Therefore, the eight variables Pv, Cv, UA, kd1, kd2, ac, Tfeed, Twi and the six input variables FIn, FM1, Fc, Fw, Fg, Vp were considered in the bifurcation analysis. To find the initial Hopf bifurcation points, the 14 variables were analyzed in the MATCONT software package. Each of them was set to change gradually from its nominal value while others remained constant. It was shown in MATCONT results that bifurcations (Hopf) were only found in the variations of ac, kd1, kd2, Fc, Fg, Twi. Thus, the subsequent SI calculation only considered these six parameters, and thus (m þ l) = 6 in the current study. In the six variables, Twi is a measurable disturbance and Fc and Fg are input variables. So, only the three parameters ac, kd1, kd2 had to be estimated by the EKF. In the calculation of SI, the scaling factor in eq 18 for each si was obtained in eq 21. To verify the effectiveness of the 3409



ARTICLE

Table 3. Parameters of PE Reactor parameters

description

values

units m3

Vg

volume of gas phase in the reactor

500

Pv

pressure downstream of bleed vent

17

atm

Bw

mass of polymer in the fluidized bed

7 104

kg

kp0

pre-exponential factor for polymer propagation rate

85 10-3

m3/(mol 3 s)

Ea

activation energy

(9000)(4.1868)

J/mol

Cpm1

specific heat capacity of ethylene

(11)(4.1868)

Cv

vent flow coefficient

7.5

J/(mol 3 K) atm-0.5 3 mol/s

Cpw CpIn

specific heat capacity of water specific heat capacity of inert gas

(103) (4.1868) (6.9) (4.1868)

Cppol

specific heat capacity of polymer

(0.85 103) (4.1868)

kd1

deactivation rate constant for catalyst site 1

0.0001

kd2

deactivation rate constant for catalyst site 2

0.0001

J/(kg 3 K) J/(mol 3 K) J/(kg 3 K) s-1 s-1

-3

MW1

molecular weight of monomer

28.05 10

Mw

mass holdup of cooling water in heat exchanger

3.31 104

Mg

mass holdup of gas stream in heat exchanger

6060.5

mol

MrCpr Hreac

product of mass and heat capacity of reactor walls heat of reaction

(1.4 107) (4.1868) (-894 103) (4.1868)

J/K J/kg J/(K 3 s)

kg/mol kg

UA

product of heat exchanger coefficient with area

(1.14 106) (4.1868)

Tf

reference temperature

360,

K

Tfeed

feed temperature

293

K

Twi

inlet cooling water temperature to heat exchanger

289.56

K

RR

ideal gas constant

8.20575 10-5

m3 3 atm/(mol 3 K)

R

ideal gas constant

8.314

J/(mol 3 K)

ac

active site concentration of catalyst

0.548

mol/kg

Table 4. Intermediate Variables of PE Reactor variables

descriptions

Hf

enthalpy of fresh stream

Hg0

enthalpy of total gas outflow stream from reactor

Hg1

enthalpy of cooled recycle gas stream to reactor

Hpol

enthalpy of polymer

Hr

heat liberated by polymerization reaction

constructed SI, two cases were studied, in which some parameters or disturbance changed gradually with time, simulating the drift processes that are usually encountered in practical industry. First, the parameter ac was set to decrease from its nominal value (0.548 mol/kg), that is, there was a reduction in the catalyst active concentration. The decrease rate was -0.0013(mol/kg) 3 h-1, which was quite slow. The EKF sampling time was set to 0.1 h and the SI calculation period was 2 h, and a total time length of 180 h was simulated. The case study was run on the computer of the following hardware and software platform: CPU, Intel Q8200 2.33 GHz; RAM, 3.25GB; Software platform, MATLAB R2009a. The average computational time of generating the SI result in one step was about 0.5888 s, which was very short compared with the online implementation interval (2 h). In the online implementation of the proposed method, the computational time was mainly cost in the solving of NLP, at about 0.5859 s, and the time used in algebraic calculation of SI through eq 19 and sensitivity analysis were only 5.1801 10-5 and 0.0028 s, respectively. The results of SI and the temperatures of the reactor, the cooling water from exchanger, and the recycle gas were plotted in

Figure 2. SI and temperature results when ac decreases.

Figure 2. The SI was at a high level at the beginning although ac had some drifts. Then, the SI began to decrease obviously after the sample of 88 h. In the sample of about 170 h, the SI value fell to zero, which implied a Hopf bifurcation point of the system was 3410



ARTICLE

Figure 3. Results for the exponential distances for each parameter when ac decreases.

Figure 4. The EKF estimation and Hopf bifurcation results for each parameter when ac decreases (blue, the estimated or measured values; green, the bifurcation values).

Figure 5. The reliability analysis: LTOP results when ac decreases.

reached and the reactor was going to have unstable oscillations or limit cycle. The prediction was soon proved by the temperature results in Figure 2. After about 1-2 h, the reactor temperature began to oscillate with the amplitude increasing exponentially. And the temperatures of cooling water from the heat exchanger and recycle gas also showed unstable oscillation phenomena. The severe oscillation may result in unsafe operation conditions of the process. Thus, the proposed SI can help early identification of those potential threats. To further analyze the causes in SI decrease, si values for each parameter were demonstrated in Figure 3. And the parameters

estimation and their critical bifurcation values were shown in Figure 4. It can be concluded from these two plots that the parameters ac, Fc, and Twi were approaching their critical bifurcation points and thus they were the main causes for the reduction in SI. If one wants to stop the further deterioration and improve the SI level, that is, keep the operation condition far away from the Hopf bifurcation, some preventative control methods should be implemented by shifting these parameters or input variables from their critical values. The reliability of the SI was analyzed in the LTOP results in Figure 5 with the Tol of SI at 0.0792 (0.0792 = log10(1.2)). It was considered that a multiplicative change of 1.2 is tolerable for the SI calculation. It can be seen that ac has a relatively large perturbation range when SI was high and close to 1. However, the tolerated bias range was getting less and less when SI approached 0, implying that SI was quite sensitive to the parametric disturbance when SI was low. Compared with existing safety analyses based on bifurcation,18-20 the proposed approach based on SI has the advantage of finding the bifurcation points online when multiple parameters/disturbances simultaneously change. Here, a designed scenario of two drifting parameters was simulated. The parameter ac was set to decrease at the rate of -0.0013(mol/kg) 3 h-1, similar to that of 3411


Industrial & Engineering Chemistry Research the previous case. Meanwhile, the disturbance kd1 was increasing slowly at a rate of 0.0125 (s 3 h)-1 (the deactivation speed of catalyst might change). The SI and temperature results were illustrated in Figure 6, and it can be seen that the decrease rate of SI is much faster when these two parameters drift simultaneously than in the case where only ac was decreasing. At about 70 h, the SI fell near to 0 and the subsequent unstable oscillation occurred after 2-3 h. The results of exponential distances si and curves for the estimated values of parameters and their Hopf bifurcation points were plotted in Figures 7 and 8, respectively. From Figures 7 and 8, it can be seen that the SI decreased faster because not only were the estimated values of ac and kd1

Figure 6. SI and temperature results when ac and kd1 drift.

ARTICLE

drifting, but also their bifurcation values were changing. It can be known that ac, kd1, Fc, and Twi were responsible for the approaching the Hopf bifurcation. The preventative control should be designed to adjust these variables away from their bifurcation points. Reliability analysis results were demonstrated in Figure 9, where a similar LTOP feature can be generated. In the example studies, it can be concluded that the proposed approach based on the SI can give an efficient online monitoring on the distance of the PE reactor operation to the Hopf bifurcation when parameters slowly drift. More importantly, it is able to provide early alerts of the potential threat of unexpected oscillations in the Hopf bifurcations of the reactor process. And the causes of deteriorations in SI can be found by analyzing each parameter, which will be quite helpful in designing further adjusting or preventative control strategies.

5. CONCLUSION In the current work, an approach is proposed to provide online safety assessment of industrial process operations based on bifurcation analyses. Hopf bifurcation may result in unstable oscillations, which are unexpected in practical operations. So it is considered as one of the potential threats to operation safety. Herein, the Hopf bifurcation manifold is characterized by a group of algebraic equations. Afterward, a safety index (SI) is constructed, which is an integration of exponential measures of distances between each parameter and their critical bifurcation points. In the online implementation of safety assessment, the parameters are estimated by the extended Kalman filter (EKF), and the Hopf bifurcations are obtained by solving the nonlinear equations by NLP algorithms. A gas phase polyethylene reactor simulation example was studied to test the efficiency of the proposed method. From the example results, it can be concluded that the constructed SI results accorded well with the Hopf bifurcation behaviors when multiple parameters drift simultaneously, and it can provide an early identification and alerts for the potential danger of severe oscillations in the reactor temperature. The detailed information in each parameter and their critical values can help find the causes of decreasing distance to the bifurcation and suggest the directions for further improving adjustment or preventative control.

Figure 7. Results for the exponential distances for each parameter when ac and kd1 drift. 3412



ARTICLE

Figure 8. The EKF estimation and Hopf bifurcation results for each parameter when ac and kd1 drift (blue, the estimated or measured values; green, the bifurcation values).

Since the terms wi, SI/si2, 1/(hi2σi4), exp(-2vi2/(hiσi2)), and exp(-vi2/(hiσi2)) are all positive and (wi - 1) is negative, ∂2SI/ ∂(vi2)2 < 0. Thus, the first order derivative, that is, ∂SI/∂(vi2), increases monotonically when vi2 decreases. Since the variation trend of vi is the same with that of vi2, ∂SI/∂(vi2) increases monotonically when vi decreases, which is property (4) of the SI.

’ AUTHOR INFORMATION Corresponding Author

*Tel.: (þ86-571)8795-1071-424. Fax: (þ86-571)8795-1068. E-mail: [email protected].

’ ACKNOWLEDGMENT Figure 9. The reliability analysis: LTOP results when ac and kd1 drift.

’ APPENDIX The properties (3-4) of the constructed SI in eq 19 are proved here. According to eqs 18 and 19, it is easily obtained that ! ∂SI SI 2vi vi 2 ¼ wi exp ðA1Þ ∂vi si hi σ i 2 hi σ i 2 In eq A1, ∂SI/∂vi > 0 because the terms wi, SI/si, 2vi/(hiσi2), and exp(-vi2/(hiσ2i )) are all positive. Therefore, SI increases monotonically when vi increases, which is just property (3) of the SI. The second order derivative of SI with respect to vi is quite complicated. Therefore, the derivative with respect to the squared spatial distance vi2 is analyzed as ∂2 SI SI 1 2vi 2 ¼ wi ðwi - 1Þ 2 2 4 exp 2 si hi σi hi σ i 2 ∂ðvi 2 Þ ! SI 1 vi 2 þ ð - 1Þwi exp si hi 2 σ i 4 hi σ i 2

!

ðA2Þ

This research was supported by the National Science Foundation of China under Grant No. 60574047, the Research Fund for the Doctoral Program of Higher Education in China under Grant No. 20050335018, and the National High Technology Research and Development Program of China under Grant No. 2007AA04Z168 and 2009AA04Z154.

’ REFERENCES (1) Ye, L.; Liu, Y.; Fei, Z.; Liang, J. Online Probabilistic Assessment of Operating Performance Based on Safety and Optimality Indices for Multimode Industrial Processes. Ind. Eng. Chem. Res. 2009, 48, 10912. (2) Lou, H. H.; Muthusamy, R.; Huang, Y. Process Security Assessment: Operational Space Classification and Process Security Index. Trans IChemE, Part B 2003, 81, 418. (3) Uygun, K.; Huang, Y.; Lou, H. H. Fast Process Security Assessment Theory. AIChE J. 2004, 50, 2187. (4) Kaye, R. J.; Wu, F. F. Dynamic Security Regions of Power Systems. IEEE Trans. Circuits Syst. 1982, 29, 612. (5) Wu, F. F.; TSAI, Y.-K. Probabalistic Dynamic Security Assessment of Power Systems: Part I-Basic Model. IEEE Trans. Circuits Syst. 1983, 30, 148. (6) Chiang, H.-D.; Hirsch, M. W.; Wu, F. F. Stability Regions of Nonlinear Autonomous Dynamical Systems. IEEE Trans. Autom. Control 1988, 33, 16. (7) Cheng, D.; Ma, J. In Calculation of Stability Region, Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii, USA, 2003, 5615. 3413


Industrial & Engineering Chemistry Research (8) Xue, A.; Wu, F. F.; Lu, Q.; Mei, S. Power System Dynamic Security Region and Its Approximations. IEEE Trans. Circuits Syst. 2006, 53, 2849. (9) Xue, A.; Wu, F. F.; Ni, Y.; Mei, S.; Lu, Q. Power System Transient Stability Assessment Based on Quadratic Approximation of Stability Region. Electr. Power Syst. Res. 2006, 76, 709. (10) Morud, J. C.; Skogestad, S. Analysis of Instability in an Industrial Ammonia Reactor. AIChE J. 1998, 44, 888. (11) Mancusi, E.; Merola, G.; Crescitelli, S.; Maffettone, P. L. Multistability and Hysteresis in an Industrial Ammonia Reactor. AIChE J. 2000, 46, 824. (12) McAuley, K. B.; Macdonald, D. A.; McLellan, P. J. Effects of Operating Conditions on Stability of Gas-Phase Polyethylene Reactors. AIChE J. 1995, 41, 868. (13) Dadebo, S. A.; Bell, M. L.; McLellan, P. J.; McAuley, K. B. Temperature Control of Industrial Gas Phase Polyethylene Reactors. J. Process Control 1997, 7, 83. (14) Russo, L. P.; Bequette, B. W. Operability of Chemical Reactors: Multiplicity Behavior of a Jacketed Styrene Polymerization Reactor. Chem. Eng. Sci. 1997, 53, 27. (15) Zhang, Y.; Henson, M. A. Bifurcation Analysis of Continuous Biochemical Reactor Models. Biotechnol. Prog. 2001, 17, 647. (16) Zhang, Y.; Zamamiri, A. M.; Henson, M. A.; Hjortso, M. A. Cell Population Models for Bifurcation Analysis and Nonlinear Control of Continuous Yeast Bioreactors. J. Process Control 2002, 12, 721. (17) Salau, N. P. G.; Neumann, G. A.; Trierweiler, J. O.; Secchi, A. R. Dynamic Behavior and Control in an Industrial Fluidized-Bed Polymerization Reactor. Ind. Eng. Chem. Res. 2008, 47, 6058. (18) Molnar, A.; Markos, J.; Jelemensky, L. Safety Analysis of CSTR towards Changes in Operating Conditions. J. Loss Prev. Process Ind. 2003, 16, 373. (19) Molnar, A.; Markos, J.; Jelemensky, L. Some Considerations for Safety Analysis of Chemical Reactors. Trans. IChemE, Part A 2005, 83, 167. (20) Vandova, Z. S.; Jelemensky, L.; Markos, J.; Molnar, A. Steady States Analysis and Dynamic Simulation as a Complement in the HAZOP Study of Chemical Reactors. Trans IChemE, Part B 2005, 83, 463. (21) Crowl, D. A.; Louvar, J. F. Chemical Process Safety: Fundamentals with Applications, 2nd ed.; Prentice Hall: NJ, 2002. (22) Rahmana, M.; Heikkila, A.-M.; Hurme, M. Comparison of Inherent Safety Indices in Process Concept Evaluation. J. Loss Prev. Process Ind. 2005, 18, 327. (23) Kubecek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Spring- Verlag: New York, 1983. (24) Chen, Y. Bifurcation and Chaos Theory of Nonlinear Vibration System; Higher Education Press: Beijing, 1993. (25) Doedel, E. J. AUTO: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations; California Institute of Technology: Pasadena, CA, 1986. (26) Doedel, E. J.; Oldeman, B. E. AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations; Concordia University: Montreal, Canada, 2009. (27) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A.; MATCONT, A MATLAB Package for Numerical Bifurcation Analysis of ODEs. ACM Trans. Math. Software 2003, 29, 141. (28) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A.; Sautois, B. Matcont: A Matlab Package for Dynamical Systems with Applications to Neural Activity; Department of Applied Mathematics and Computer Science Ghent University, Mathematisch Instituut Universiteit Utrecht, Vakgroep Informatietechnologie KAHO Sint-Lieven: Belgium, 2006. (29) Dobson, I. Computing a Closest Bifurcation Instability in Multidimensional Parameter Space. J. Nonlinear Sci. 1993, 3, 307. (30) Monnigmann, M.; Marquardt, W. Normal Vectors on Manifolds of Critical Points for Parametric Robustness of Equilibrium Solutions of ODE Systems. J. Nonlinear Sci. 2002, 12, 85. (31) Monnigmann, M.; Marquardt, W.; Prozesstechnik, L. f. In Parametrically Robust Control-Integrated Design of Nonlinear Systems

ARTICLE

Proceedings of the American Control Conference, Anchorage, AK, 2002; p 4321. (32) Monnigmann, M.; Marquardt, W. Steady-State Process Optimization with Guaranteed Robust Stability and Feasibility. AIChE J. 2003, 49, 3110. (33) Marquardt, W.; Monnigmann, M. Constructive Nonlinear Dynamics in Process Systems Engineering. Comput. Chem. Eng. 2005, 29, 1265. (34) Monnigmann, M.; Marquardt, W. Steady-State Process Optimization with Guaranteed Robust Stability and Flexibility: Application to HDA Reaction Section. Ind. Eng. Chem. Res. 2005, 44, 2737. (35) Li, R.; Corripio, A. B.; Henson, M. A.; Kurtz, M. J. On-line State and Parameter Estimation of EPDM Polymerization Reactors Using a Hierarchical Extended Kalman Filter. J. Process Control 2004, 14, 837. (36) Ricker, N. L.; Lee, J. H. Nonlinear Modeling and State Estimation for the Tennessee Eastman Challenge Process. Comput. Chem. Eng. 1995, 19, 983. (37) Sirohi, A.; Choi, K. Y. On-Line Parameter Estimation in a Continuous Polymerization Process. Ind. Eng. Chem. Res. 1996, 35, 1332. (38) Lee, J. H.; Ricker, N. L. Extended Kalman Filter Based Nonlinear Model Predictive Control. Ind. Eng. Chem. Res. 1994, 33, 1530. (39) Jang, S.-S.; Joseph, B.; Muka, H. Comparison of Two Approaches to On-Line Parameter and State Estimation of Nonlinear Systems. Ind. Eng. Chem. Process Des. Dev. 1986, 25, 809. (40) Gani, A.; Mhaskar, P.; Christofides, P. D. Fault-Tolerant Control of a Polyethylene Reactor. J. Process Control 2007, 17, 439. (41) Ohran, B. J.; Pena, D. M. d. l.; Davis, J. F.; Christofides, P. D. Enhancing Data-Based Fault Isolation through Nonlinear Control. AIChE J. 2008, 54, 223.

3414


A Method of Online Safety Assessment for Industrial Process

Recommend Documents