Fault Diagnosis Using the Hybrid Method of Signed Digraph and

Ind. Eng. Chem. Res. 2006, 45, 9061-9074

9061

Fault Diagnosis Using the Hybrid Method of Signed Digraph and Partial Least Squares with Time Delay: The Pulp Mill Process Gibaek Lee* Department of Chemical and Biological Engineering, Chungju National UniVersity, Chungju 380-702, Korea

Thidarat Tosukhowong and Jay H. Lee School of Chemical and Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332

Chonghun Han* School of Chemical and Biological Engineering, Seoul National UniVersity, Seoul 151-742, Korea

Although the existence of a time delay worsens the diagnosis accuracy, most research efforts on fault diagnosis have not considered the time delay. The hybrid fault diagnosis method that is based on the signed digraph (SDG) and the partial least-squares (PLS) has the advantage of improving diagnosis resolution and accuracy, compared to those of previous qualitative methods, and is capable of enhancing the ability to diagnose multiplefault [Lee et al., Ind. Eng. Chem. Res. 2003, 42, 6145-6154]. The current research modified this method to include the information of the time delay between process variables, and the time delay can be obtained through the open-loop responses of set-point changes. The modified method is applied for the fault diagnosis of the pulp mill process, which produces pulp from wood chips. It is one of the biggest processes in the fault diagnosis area, and it will be a new benchmark process to evaluate the fault diagnosis method. Also, its units, such as the storage tank and the bleaching tower, have a very large time delay. Dynamic principal component analysis (DPCA)-based fault diagnosis was used to get a reference of the diagnosis performance comparison. Although fault detection of the proposed method is frequently slower than DPCA, it gives more accurate results. Also, it demonstrates much-better diagnosis capabilities, compared to the original hybrid method. Introduction Fault diagnosis system analyzes process data on-line and diagnoses faults to improve the safety of the chemical plant and their plant personnel. A variety of fault diagnosis methods are classified as model-driven techniques and data-driven techniques, and they include the expert system, state estimation such as observer and EKF (extended Kalman filter), signed digraph (SDG), fault tree, qualitative simulation, statistical methods, and neural network.1 The main difficulties in developing a diagnostic system of chemical processes are due to the exclusive features of chemical processes of nonlinearity, uncertainty, and complexity. For more than the past two decades, many research efforts have focused on nonlinear fault diagnosis and robust fault diagnosis for uncertain processes.2-4 However, although time delay exists widely in real engineering problems, the researches on fault diagnosis that consider time delay are few.5,6 This study combines the time delay information with our previous hybrid method, combining SDG and the partial least squares (PLS) approach.7 The hybrid method improves the diagnosis resolution and accuracy, compared to previous qualitative methods. Moreover, it enhances the reliability of the diagnosis for all predictable faults, including both novel faults and multiple-fault. As a testbed of the modified method, we chose a pulp mill process, which is a much larger benchmark process than the Tennessee Eastman process that has been used to evaluate fault diagnosis methods for more than a decade.8-10 Although fault * To whom correspondence should be addressed. Fax, +82-43-8415220; E-mail: [email protected] (for G.L.). Fax, +82-2-888-7295; E-mail: [email protected] (for C.H.).

diagnosis systems of pulp mills are in demand for an operational advantage in a very competitive global market, to our knowledge, currently no literature on the fault diagnosis of a pulp mill exists. A pulp mill process that was simulated by Castro and Doyle included 8200 states, 140 inputs, 114 outputs, and 6 operation modes, and it can be a realistic testbed for plantwide fault diagnosis methods.11,12 In addition, this process shows the common characteristics of the industrial processes such as long measurement delay, many components, several recycle streams, and large nonlinearity. Especially, its equipments, such as the storage tank and the bleaching tower, have very large time delays (>2 h), which are suitable to test the proposed method. Fault Diagnosis Method System Decomposition Based on SDG. The hybrid method that combines SDG and PLS starts by decomposing the system to focus on the measured variables in SDG. Each arc of the SDG represents the instantaneous effect originated from the source node to the target node. Only the source nodes that are connected to a particular target node by means of the arcs can affect the target node. Because unmeasured nodes are unable to demonstrate the direct effects from faults, they are removed to obtain the reduced digraph, which consists of the original SDG with the unmeasured nodes removed. Each decomposed subprocess includes a central measured variable (target variable) as well as measured variables (source variables) and faults that are connected to the target variable. For example, consider Figure 1a, which shows the SDG focusing on M1 in a process. The nodes M2, M3, M4, U1, and U2 have a direct impact on M1, whereas none of the faults shown in the figure connect directly

10.1021/ie060793j CCC: $33.50 © 2006 American Chemical Society Published on Web 11/11/2006

9062

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006

In this equation, BPLS is the coefficient of the PLS regression model, and the weight W is defined in the NIPALS algorithm:

T ) XW

Figure 1. (a) Partial signed digraph focusing on M1. (b) Reduced digraph of the decomposed subprocess for M1.

to M1. On the other hand, M2 is directly affected by two faults of F1 and F2. Also, U2 and U1 are directly affected by faults of F3, F4, and F5, respectively. However, U1 and U2 cannot reveal the direct effects made by F3, F4, and F5, because these nodes are not measured. After the unmeasured nodes of U1 and U2 are removed, the nodes that have a direct effect on M1 are the measured nodes of M2, M3, and M4 and the faults of F3, F4, and F5, as shown in Figure 1b. If these faults do not occur, only the source nodes of M2, M3 and M4 affect M1. Local fault diagnosis can be performed on each decomposed subprocess. Together with the system decomposition, this fault diagnosis method, when locally executed on each target variable, has the advantage of being able to diagnose all types of multiplefault, except for those that affect the same measured nodes. Fault Diagnosis Based on Dynamic PLS Models. A simple diagnosis method based on the decomposition technique is to estimate the value of each target variable, using measured values of the source variables that are connected to the target variable. The difference between the estimated and measured values indicates the occurrence of one or more faults. The sensor faults, which correspond to errors the source variables used for the estimation, produce errors in the estimated values. Also, the faults that are connected to the target node result in errors in the measured values. The estimation of our previous study used the PLS model that was built for each decomposed subprocess. By maximizing the covariance between the input data matrix X and the output matrix Y, PLS can predict the values of Y based on X values. The PLS model consists of outer relations (X and Y blocks individually) and an inner relation (linking the two blocks). The scaled and mean-centered X and Y matrixes are decomposed into score and loading vectors, denoted by T and U, and P and Q, respectively. The inner relationship between the two blocks of X and Y is represented as a linear algebraic relation (U ) TB) between their scores:

X ) TPT + E

(1)

Y ) UQT + F

(2)

where E and F are the residual matrixes for X and Y, respectively. In PLS, the loading and score vectors are determined in such a way to maximize the prediction accuracy of Y while describing a large amount of the variation in X. The most common method used to calculate the PLS model parameters is known as noniterative partial least-squares (NIPALS). The prediction of Y from X is done according to the following regression model:

Yˆ ) XBPLS ) XW(PTW)-1BQT

(3)

(4)

The output Y of the PLS model is the estimated value of a target variable, and the input X contains the source variables that are connected to the target variable. To handle the process dynamics accurately, PLS is integrated with ARMAX to become dynamic partial least squares (DPLS). In addition to the past values of the source variables, the resulting input of DPLS for a target variable includes the past values of the target variable as well as the source variables. The input matrix, X, of the DPLS model for the variable i at time t is as follows:

[] xij(t)

l yi(t - ∆ti) xij(t - ∆tj) l Xi ) yi(t - l∆ti) xij(t - l∆tj) l xiNi(t - l∆tNi)

T

(5)

where l is the number of necessary past values, ∆t the sampling rate, and Ni the number of source variables for variable i. The number of past values (time lags, l) and the principal components (PCs) are determined from the learning data. The number l is usually 1 or 2, which indicates the order of the dynamic system. Each DPLS model can be built from the operation data set representing the local relations between the input variables and the output variable of the DPLS models. Therefore, the required data set for each DPLS model can be obtained in the presence of set-point changes or external disturbances, which occur frequently. The hybrid method does not need a faulty case data set, which would otherwise be difficult to obtain. Fault detection is performed by the observation of the residual, which is the difference between the estimated value determined by the DPLS model and the measured value:

ri ) yi - yˆ i

(6)

where ri is the residual of output variable i, and yi and yˆ i are the measured and estimated values of variable i, respectively. A qualitative state, which corresponds to ranges of possible values for the residual, becomes an attribute of the residual. Similar to SDG-based methods, three ranges are used: low, to which the qualitative state is assigned (denoted with the symbol “-”); normal (denoted with the symbol “0”); and high (denoted with a symbol “+”). If a fault occurs, the qualitative state for the residual may be (+) or (-). The abnormal qualitative state for the residual becomes a symptom, which is expressed as a pair of the target variable and the qualitative state of the residual. Those faults causing the abnormalities of each residual are classified along with their symptoms, and the classified faults are stored in a set (called a fault set). The first step of on-line fault diagnosis is the monitoring of the residuals, to detect their qualitative change of state. The monitoring is achieved by CUSUM to ensure a stable detection and diagnosis. The detected residual of a variable becomes an

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 9063

Pulp Mill Process

Figure 2. (a) Measurement delay and (b) snapshot of variables affecting the true values of M1 at time (t - 5).

element in the symptom set. The next step of fault diagnosis is to obtain the minimum set of faults that can explain all of the detected symptoms. Incorporation of Time Delay into the DPLS Model. Time delay can be defined as the time interval between the start of an event at one location and its resulting action at another location.13,14 It is also referred transportation lag, dead time, or distance-velocity lag. In a chemical process, the largest time delay often results from a buffer tank, which is often installed to mitigate disturbances from the upstream process. In the hybrid model, DPLS model represents the instantaneous effects produced from the input variables to an output variable. Therefore, a unit process that has a long time delay between input and output variables of DPLS model can worsen the model accuracy. To increase the accuracy of the estimation, the information on the time delay from the input variables to the output variable should be incorporated into the DPLS model. The input, X of the DPLS models is modified to consider time delay as well as sample delay.15 For example, in Figure 1, M1 has the sampling rate and delay of 2 min, and M2 has the sampling rate of 1 min and no sampling delay. M3 and M4 have the sampling rate and delay of 5 min, and the time delay from M4 to M1 is 10 min (Figure 2a). The minimum diagnosis interval of M1 is the same, with the sampling size of 2 min. The value of M1 at t is estimated by using M2, M3, and M4. At time t, the available true values are the value of M1 up to (t - 2), M2 up to t, and M3 and M4 up to (t - 5). The largest measurement delay is 5 min for M3 and M4, and the input matrix X should be made on the basis of time, (t - 5). As shown in Figure 2a, the true value of M1 at (t - 5) can be measured at (t - 3) because the sampling delay of M1 is 2 min. By the same reason, M2(t 5), M3(t), and M4(t) can give the true values at (t - 5). Because M4 has a time delay of 10 min, the value measured at (t - 10) should be used. As shown in Figure 2b, we should use M1 at (t - 3), M2 at (t - 5), M3 at t, and M4 at (t - 10). The input matrix, X of the DPLS model for the variable i is modified from eq as follows:

[

MAX xij (t + τMD, j - τMD,i - τTD,ij ) l MAX - ∆ti ) yi (t + τMD,i - τMD,i MAX - τTD,ij - ∆tj) xij (t + τMD, j - τMD,i l Xi ) MAX - l∆ti ) yi (t + τMD,i - τMD,i MAX xij (t + τMD, j - τMD,i - τTD,ij - l∆tj ) l MAX - τTD,iNi - l∆tNi ) xiNi(t + τMD,Ni - τMD,i

]

T

(7)

where τMD,i is the measurement delay of i, τTD,ij is the time delay MAX ) max(τMD,j ) 1,Ni,τMD,i). The measured from j to i, and τMD,i value of variable i to calculate the residual i is yi(t + τMD,i MAX ). τMD,i

Process Description. The Kraft pulp mill is a highly integrated process that consists of two major areas: the fiberline and the recovery plant, as depicted in Figure 3.11 The objective of the fiberline area is to produce fibers from wood chips at a desired production rate and quality. Major raw materials of this process are wood chips and chemicals called “white liquor” (WL), which consists primarily of NaOH and NaSH. They are combined in a pressurized impregnation vessel, where wood chips are saturated in the white liquor. They then enter the cook zone of the digester, where lignin in the wood starts to dissolve out. As the materials continue to flow down to the mcc and emcc zones of the digester, white liquor is added and more lignins are dissolved from the wood. The main controlled variable in this unit is the Kappa number, which is a measure of the amount of remaining lignin in the wood. The final zone of the digester is the wash zone, where wash liquor is added and fibers are blown out of the digester. Fibers are further washed in a brown stock drum washing section, to remove chemicals and residual lignin. Fibers then are bleached in several bleaching towers, to further remove lignin and achieve a target brightness coefficient. The bleaching sequence includes postdelignification with oxygen (O) and white liquor, chlorine dioxide (D1), sodium hydroxide (E), and brightening via chlorine dioxide (D2). At the end of each bleaching stage, pulp is washed to removed chemicals and lignin content before going to the next bleaching stage. There are also several recycle loops in this area, including the effluence from the O washer that is used as wash liquor in the brown stock washer and as the emcc dilution water. Furthermore, the effluent from the D2 washer is wash liquor for the E washer. On the other hand, exit chemical streams from the digester and washing stages now has many organic residuals and brown color and, hence, are called “weak black liquors”. To recover chemical components and energy from these streams, the weak black liquors are sent to the recovery plant. The recovery process begins by concentrating the weak black liquor streams in multieffect evaporators until the black liquor solid concentration is high enough to sustain combustion. Subsequently, the strong black liquor stream is fired in the recovery boiler, where the combustion heat is transferred to feedwater to produce highpressure steam and perform the reduction reactions to recover Na2S, Na2CO3, and other sulfur-based salt. This molten salt (which is called smelt) flows out at the bottom of the recovery boiler and dissolves in the weak wash water in the smelt dissolving tank to produce green liquor. Green liquor is further clarified to remove unburned carbon and insoluble components from the liquor. In the slaking/causticizing section, lime is added to the green liquor to produce white liquor. White liquor is then clarified to remove CaCO3 and other insoluble solid before it can be sent to the digester and the oxygen reactor. Lime mud from the white liquor clarifier is burned in the lime kiln to recover lime. Fault Scenario Simulation. Because of the complexity in many unit operations and high interactions among process units, pulp mill represents a very challenging integrated plant for process system engineering community. Recently, Castro and Doyle12 introduced a pulp mill simulator that includes 8200 states, 140 inputs (82 manipulated variables and 58 disturbances), 114 outputs, and 6 operating modes. Major economic variables are wood chip flow rate, steam usage, water usage, chemical (NaOH, oxygen, ClO2, Na2SO4, and CaO) usages, natural gas usage, and pulp production. They also designed the decentralized control scheme to regulate the process

9064


Figure 3. Simplified schematic of the pulp mill process.12

to meet process constraints and operation objectives. This is achieved through the use of 49 regulatory control loops and 13 supervisory control loops. This control structure has been tested under several disturbance and grade-transition scenarios and shown to have satisfactory results, given that all the control elements were functioning perfectly. The pulp mill simulator was used as a testbed for the proposed fault diagnosis method. It can be downloaded from the homepage of Doyle’s group.16 This study modified the simulator to obtain a more-realistic fault diagnosis problem, as follows. At first, white noise was added to the simulator input and output variables. The maximum and minimum values of the noise are (0.1% of the maximum expected errors or changes of the variables, which is specified in the original simulator. For the fault diagnosis, 114 measured output variables and 82 controller outputs are used. For the convenience of programming, the variables of the simulator are renamed by adding the type of

variable and the sequential number. The abbreviation letter to denote variable type is shown in Table 1. The used sequential numbers are 1-40 for measured variables (MVs) of fiber line, 41-78 for MVs of fiber line, 79-152 for measured variables of chemical recovery, and 153-196 for MVs of chemical recovery. For brevity of this paper, we do not present all the variable descriptions of the simulator, because they can be found in the work of Castro and Doyle;12 however, the variables used in fault diagnosis models are described in Table 2. In the table, the two numbers in parentheses represent the sampling rate of more than 5 min and the sampling delay of more than zero, respectively. We investigated cases when there were faults in critical control elements. These faults were imposed on four sensors and four control valves, which were critical to the fiberline operation. These faulty sensors measure important control variables (CVs) of the pulp mill, including the digester Kappa

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 9065 Table 1. Abbreviation Letters Used in the Definition of Faults and Measured Variables Faults

Measured Variables

abbreviation letters

meanings

abbreviation letters

meanings

BR CF D1C KN OH PR WC WL Bias Deg Stuck

final brightness sensor of the pulp in the D2 caustic flow 3 control valve D1 ClO2 flow control valve digester Kappa number sensor [OH-] concentration sensor before the E washer D2 production rate sensor wood chips flow control valve oxygen WL flow control valve sensor or control valve bias sensor degradation control valve stuck

BR C CS DF F KN MV OH P T V

brightness concentration consistency dilution factor flow rate Kappa number manipulated variable [OH-] concentration pressure temperature volume

Table 2. Pulp Mill Variables Used for Fault Diagnosis Model variable name

description

variable name

description

F1 F3 KN4 CD10 T12 T13 T14 T15 T16 F18 KN19 T20 T21 KN22 T23 OH24 T25 BR26 DF31 DF32 CS36 P37 MV41 MV42 MV43 MV44 MV45

digester production rate (10/10) D2 production rate (10/10) digester Kappa no. (10/10) digester upper extract conductivity digester cook zone T digster mcc zone T digster emcc zone T digester WL T digester upper extract T recycle stream number 2 O Kappa number (10/10) OT D1 T E Kappa number (10/10) ET E washer [OH-] concentration (10/10) D2 T D2 brightness washer number 5 DF (5/5) washer number 6 DF (5/5) O washer inlet consistency OP wood chips flow digester cook zone WL flow digester counter zone WL flow digester steam flow 1 digester steam flow 2

MV46 MV52 MV54 MV55 MV56 MV57 MV58 MV59 MV60 MV62 MV63 MV64 MV65 MV66 MV68 MV69 MV76 MV77 MV78 T120 C122 F123 F139 EA148 MV172 MV173 MV186

digester steam flow 3 O caustic flow O steam flow 1 O steam flow 2 O split fraction 1 O steam flow 3 pulp flow to bleach plant D1 water flow D1 ClO2 flow D1 steam flow E caustic flow E water flow E steam flow D2 ClO2 flow D2 wash water flow O split fraction 2 O coolant flow O oxygen flow D2 exit flow WL T WL NaOH concentration (10/10) WL NaOH mass flow (10/10) WL flow estimated wlc EA caustic makeup flow oxygen WL flow O wash water flow

number, the D2 production rate, the final brightness of the pulp in the D2, and the hydroxide concentration before the E tower. These CVs are also under the following control loops: (1) The digester Kappa number is controlled by adjusting the temperature set point of the emcc zone. This cascade controller manipulates the steam flow to the emcc. (2) The D2 production rate is under a proportional integral (PI) controller that adjusts the wood chip flow rate to the digester. (3) The final pulp brightness is controlled by a cascade controller that manipulates the D2 Kappa factor set point. The Kappa factor is controlled by the flow rate of ClO2 to this tower. (4) The hydroxide ion concentration ([OH-]) before the E tower is paired with the set point of the E Kappa factor. This supervisory controller regulates the E Kappa factor to control the E Kappa factor to the set point. Five sensor fault types, including bias (-5% and +20%), precision degradation ((1% and (10%), and drift (0.025% × time) were simulated for these sensors.17 In addition to this, we also simulated three problems of +1% and -5% bias, and sticking for four control valves that manipulate wood chip flow, white liquor flow to the oxygen reactor, ClO2 flow to D1, and the caustic makeup stream that is added to the white liquor in the recovery plant. Deviation of the wood chip flow from the

target will affect the final pulp production rate and pulp quality. In addition, changes in the other three variables will affect the lignin removal efficiency of the process, and thereby affect the brightness of the pulp. Without other disturbances, the variations due to the sticking of control valves may be very small. Therefore, the fault is simulated with the +10% set-point change of the wood chip flow. In this work, we study five disturbance scenarios that are important for pulp mill operations, including variations in wood temperature (WoodTemp) and densities (WoodDens), changes in operating temperature (ECausTemp) and compositions (ECausComp) of the caustic in the E tower, and change in the wash water temperature of the E tower (EbackTemp). The wood chip temperature has an important effect on the reaction rate of the lignin removal process. In this simulator, the wood chips consist of five components (two lignin components and three cellulose components). Thus, variations in wood compositions have a direct impact on the Kappa number of the pulp. In addition, changes in the temperature and compositions of the caustic will affect the bleaching kinetics in the E tower. Furthermore, the wash water temperature will affect the pulp temperature going into the next process stage. Adding the type of fault to the name of the equipment composes the names of the fault. For example, the terminology

9066


Figure 4. Schematic of the digester.12

Figure 5. Schematic of the oxygen reactor.12

KNm5Bias means a -5% bias of the digester Kappa number sensor. The equipments are described in Table 2, and the disturbances are named from their abbreviated meanings, as given previously. The diagnosis interval is 5 min and the total simulation time is 1500 min. Fault occurs at 600 min. The simulation data can be downloaded from the URL http:// chem.chungju.ac.kr/glee.htm. DPLS Model Building. The defined 37 faults (20 sensor faults, 12 control valve faults, and 5 external disturbances) have direct effects on 18 variables among the 114 measured variables.

The variables are F1, F3, KN4, EA8, CD10, T12, T15, and T16 around the digester (Figure 4), K19, T20, T21, KN22, T23, OH24, T25, BR26, and P37 around the bleaching plant (see Figures 5 and 6), and C122 in the recausticizing section. Among these variables, EA8 is excluded, because of too long of a sample rate and a delay of 120 min. In addition, T20 is excepted, because changes in temperature are not sensed from the digester by T20, and there are many variables to replace it. The process is decomposed to center on 16 measured variables. For each decomposed process, the source variables that are connected to


Figure 6. Schematic of the bleach plant.12 Table 3. Input Variables, Number of Principal Components (PCs), and Number of Past Values of the DPLS Model With Dead Time target variable F1 F3 KN4 CD10 T12 T15 T16 KN19 T21 KN22 T23 OH24 T25 BR26 P37 C122

source variables connected to the target variable T12(10), T13(50), T14(5), T15(5), T16(5), F18(5), MV41(0), MV42(10), MV43(10), MV46(10), F123(10), F139(10), MV173(10) MV78(10) F1(0), T12(10), T13(50), T14(10), T15(0), T16(5), F18(5), MV41(10), MV42(10), MV43(10), MV46(10), F123(10), F139(10), MV173(10) T12(45), MV41(5), MV42(50), MV43(50), MV46(10), F123(50), F139(50), MV173(50) T15(0), MV41(0), MV42(5), MV43(5), MV45(5), F139(5), MV173(5) MV42(5), MV43(5), MV44(0), T120(80), F139(5), MV173(5) T12(5), MV41(5), MV42(35), MV43(35), MV46(35), F123(40), F139(40), MV173(40) F1(10), KN4(10), T20(10), P37(5), MV52(10), MV77(10), F123(10), F139(10), MV173(10) F1(245), T20(235), MV56(235), MV57(235), MV58(0), MV59(0), MV60(0), MV69(235), MV76(235), F139(230), MV173(230), MV186(235) F1(330), KN19(320), T21(45), T23(35), DF31(40), CS36(320), MV56(310), MV58(10), MV59(10), MV60(50), MV63(10), MV65(0), F139(310), MV173(310) F1(245), T21(0), DF31(0), CS36(235), MV58(0), MV59(0), MV60(0), MV62(0), MV63(0), F139(230), MV173(230) F1(290), KN19(280), T21(55), T23(35), DF31(40), CS36(290), MV56(300), MV58(40), MV59(40), MV60(40), MV63(40), C122(300), F139(300), MV173(300) F1(245), T23(0), DF32(0), CS36(235), MV58(0), MV59(0), MV60(0), MV63(0), MV64(0), MV65(0), MV66(0), MV68(0), F139(235), MV173(235) F1(390), KN4(380), T21(150), KN22(105), T23(150), OH24(105), T25(105), CS36(390), MV58(145), MV59(140), MV60(140), MV63(140), MV66(140), C122(400), F139(400), MV173(400) F1(0), KN4(0), T20(5), MV52(5), MV54(5), MV555), MV77(5), F123(5), F139(5), MV173(5) F139(0), EA148(0), MV172(0)

the target variables and the sign of the arc from the source variables to the target variables are obtained based on process knowledge and first-principle equations, such as mass and energy balance (Table 3). The time delay from the source variables to the target variables was determined using the open-loop response of the 82 control loops. The time delay shown in the parentheses in Table 3 is the multiple of the diagnosis interval. In the target process, a large time delay can be found in two units of the storage tank and the D2 tower. T21 (D1T), which is located after the storage tank, is controlled by MV57 (O steam flow) located before the

Without Dead Time

number of past values

number of PCs

number of past values

number of PCs

2

13

2

13

2 2

2 13

1 2

2 13

2 1 1 2 2

9 8 6 8 8

1 2 2 1 2

7 10 6 8 9

1

12

2

14

2

18

2

19

2

13

2

14

2

17

2

17

1

15

2

19

2

18

2

17

1

6

2

11

2

2

2

2

storage tank. Figure 7 shows the dynamics of T21 (thick line) and MV66 (thin line) when the T21 set point is increased 10% at 600 min. One can see the time delay of 235 min between two variables. Similarly, the time delay of 140 min is observed between MV66 and BR26, which are located before and after the D2 tower. Because the variable variations due to sensor degradation or control valve sticking are random, the signs of the symptoms can fluctuate between positive (+) and negative (-), which may greatly decrease the diagnosis accuracy. To consider the characteristics of these faults and make a stable diagnosis, the

9068


Figure 7. Dynamics of manipulated and controlled variables: (a) set-point change of T21 and (b) set-point change of BR26. Table 4. Fault Sets of Pulp Mill Process Sign of Arc measured variable F1 F3 KN4

positive (+)

negative (-) WCm5Bias, WL1Bias PRm5Bias KNm5Bias, WCm5Bias, WL1Bias

CD10 T12 T15 T16

WC1Bias, WLm5Bias PR0025Drift, PR20Bias KN0025Drift, KN20Bias, WC1Bias, WLm5Bias WC1Bias, WL1Bias WCm5Bias, WLm5Bias WLm5Bias WCm5Bias

KN19

KLm5Bias, WCm5Bias, WL1Bias

T21

D1Cm5Bias

KN20Bias, KN0025Drift, WC1Bias, WLm5Bias D1C1Bias

KN22

D1Cm5Bias, WL1Bias

D1C1Bias, WLm5Bias

D1C1Bias, OH0025Drift, OH20Bias, WLm5Bias

D1Cm5Bias, Ohm5Bias, WL1Bias

BR26

BR0025Drift, BR20Bias, D1C1Bias, KN0025Drift, KN20Bias, Ohm5Bias

BRm5Bias, D1Cm5Bias, KNm5Bias, OH0025Drift, OH20Bias

P37

KN0025Drift, KN20Bias, WCm5Bias, WL1Bias CF1Bias

KNm5Bias, WC1Bias, WLm5Bias

WCm5Bias, WLm5Bias WC1Bias, WL1Bias WL1Bias WC1Bias

T23 OH24 T25

C122

CFm5Bias

diagnosis strategy is modified for CUSUM to monitor the squared residuals as well as the residuals of each variable, according to the following equation:

ri2 ) (yi - yˆ i)2

(8)

In addition, the sign of the arcs from WoodDens to F1, KN4, etc. cannot be determined because of the complex process behavior. By the same reason, we had difficulties to judge the signs of several arcs from faults to measured variables, and the squared residuals should be used for these arcs. The fault sets for the pulp mill process are obtained as shown in Table 4. The learning data to build the DPLS model is obtained from the closed-loop set-point changes of 82 control loops and the changes are set on the basis of +10%. When the set point of the excess WL split fraction control valve was changed, the process response was opposite to the expectation, and the data set was excluded from the training data sets. The simulation time for each training data set was 1000 min, and the set point was changed after 600 simulation minutes from the run. The total number of samples was 200 × 81 ) 16200. From the learning data, the number of past values l and PCs are deter-

positive (+) or negative (-) WCStuck, WLStuck, WoodDens, WoodTemp PR1Deg, PR10Deg KN1Deg, KN10Deg, WCStuck, WLStuck, WoodDens, WoodTemp WCStuck, WLStuck, WoodDens WCStuck, WLStuck, WoodDens, WoodTemp WLStuck WCStuck, WL1Bias, WLm5Bias, WLStuck, WoodDens, WoodTemp KN1Deg, KN10Deg, WCStuck, WLStuck, WoodDens D1Cstuck, WL1Bias, WLm5Bias, WLStuck, WoodDens D1Cstuck, EbackTemp, ECausComp, WLStuck, WoodDens ECausTemp, WL1Bias, WLm5Bias, WLStuck, WoodDens D1Cstuck, EbackTemp, ECausComp, OH1Deg, OH10Deg, WLStuck, WoodDens EbackTemp, ECausComp, WL1Bias, WLm5Bias, WLStuck, WoodDens BR1Deg, BR10Deg, D1Cstuck, EbackTemp, ECausComp, KN1Deg, KN10Deg, OH1Deg, OH10Deg, WL1Bias, WLm5Bias, WLStuck, WoodDens KN1Deg, KN10Deg, WCStuck, WLStuck, WoodDens CFStuck

mined (see Table 3), and the CUSUM parameters of minimal jump size and threshold size are calculated. The effect size from the source variables to the target variable can be guessed from the regression coefficients of the DPLS model. Among the 16 PLS models, an error was discovered in the model of KN4. Figure 8a shows the regression coefficient of the DPLS model to estimate KN4. The left, center, and right bars for each source variable represent regression coefficients of lag that are equal to 0, 1, and 2, respectively. The summed regression coefficients for each source variable are shown in Figure 8b. In this figure, we can find the sign of the arc from MV173 to KN4 to be positive (+). However, an increase in MV173 means the decrease of the WL flow to the digester, and the arc should have a negative (-) sign. Because the wrong model of KN4 may cause the wrong detection, the DPLS model was removed. The wrong model may be due to high nonlinearity of the digester and theKappa number. When the estimation of KN4 was attempted using a quadratic PLS model, the result was worse. If the unmeasured variables that are connected to the output variable of a PLS model are able to be measured, the estimation result will be better. For more-accurate estimation of KN4, online wood composition analysis will be very helpful.


SPEk ) (xk - tkpTk )2

(9)

In this equation, tk and pk represent the score and loading, respectively.9 The variables with the largest contributions to the SPE become the primary targets for the ensuing fault diagnosis. For fairer comparisons with the proposed method, this study uses dynamic principal component analysis (DPCA) models, which was suggested by Ku et al.18 and integrates PCA with the ARMAX time series model. The data matrix of DPCA includes both the current observation and the previous l observations data. Therefore, the aforementioned equation should be modified that the SPE contribution of each variable is calculated by summation of the SPE contributions of the current and previous observations.

SPEk ) (xk(t) - tk(t)pTk (t))2 + l

Figure 8. Regression coefficients of the PLS model for KN4.

Fault Diagnosis by DPCA. PCA was used to obtain a reference of diagnostic performance, to be compared with the proposed method. PCA finds linear combinations of variables that describe major trends in the data, and it yields a simple and fast fault detection method based on squared prediction error (SPE). The contribution to the SPE from the kth variable is calculated as follows.

Figure 9. Dynamics of measured and estimated variables for WLm5Bias.

(xk(t - i) - tk(t - i)pTk (t - i))2 ∑ i)1

(10)

The number of past values l and PCs are determined by parallel analysis.18 The numbers of PCs for l ) 0, 1, 2, and 3 are 38, 49, 51, and 53, respectively. We chose l ) 2 and 51 PCs for steady-state data of 3000 min. Based on the detection result of steady-state data, the SPE limit is increased from 240.27 to 253.3 (by ∼5%) to remove false detections in the steadystate regime.

9070


Figure 10. Diagnosed variables by PCA for WLm5Bias.

Results and Discussion To measure the diagnostic performance, four parameters were used: accuracy, resolution, wrong detection, and detection delay.7 The accuracy is 1 if the diagnosis is accurate; that is, the true fault is included in the final fault candidates set. Otherwise, the accuracy is 0. The resolution denotes the number of final fault candidates. Wrong detection refers to the number of falsely detected symptoms independent of the true solution. The detection delay refers to the time from fault occurrence to fault diagnosis, and a short detection delay indicates quick detection and diagnosis. This study calculates the average performance parameter from the initial detection time to the last diagnosis time. Also, the worst results of wrong detection and resolution during all diagnostic periods are used to compare the performance. Example: WLm5Bias (-5% Bias of Oxygen at the WL Flow Control Valve). When the fault occurs at 600 min, the oxygen white liquor flow control valve is biased to -5%. As shown in Figure 9, a control loop reacts to increase oxygen WL flow (MV173) and decrease WL flow into the digester. This results in decreases in the digester WL temperature (T15), the digester upper extract EA (EA39), and digester production rate (F1), and it increases the digester Kappa number (KN4). As the WL flow to the oxygen reactor increases, the O Kappa number (KN19) decreases. The thin and thick lines of Figure 9 are the measured and estimated values, respectively.

Figure 11. Residuals of the detected variables for WLm5Bias.

Figure 10 shows the variables that have been diagnosed by DPCA. From 605 min to 695 min, DPCA wrongly diagnosed DF30 or MV77 as having the maximum contributions to the SPE. KN19 (the O Kappa number, from 700 min to 725 min), MV173 (from 730 min to 1175 min), and MV42 (digester cook zone WL flow, from 1180 min to the last diagnosis time) are diagnosed as the variables related to the fault; however, it is another task to find the true fault from these diagnosed variables, as is well-known. Figure 11 shows the residuals of the detected variables, and the bounds of the figure are the minimal jump size of CUSUM. The detection sequence of symptoms are as follows: T152, P37(-), and P372 from 620 min, KN192 from 635 min, T15(+) and KN19(-) from 645 min, CD10(-) from 675 min, T162 from 680 min, T16(-) at 680-825 min, F1(+) from 685 min, and CD102 from 690 min. After the detection of F12 at 695 min, T232, T212, OH242, OH24(+), KN222, KN22(-), BR262, BR26(+), KN22(+), BR26(-), and T21(-) are detected. At 1045-1170 min, the fault candidate is WLStuck, because of the wrong detection of KN22(+), and the true solution of WLm5Bias is missed. Except these periods, the fault candidates are WLm5Bias and WLStuck from 620 min to the last diagnosis time. Therefore, the average accuracy is 0.85 from the detection to the last diagnosis time, and resolution is 2. Consider the diagnosis including the DPLS model of KN4, which has the wrong regression coefficient. KN4(-) is wrongly detected from 745 min, and it causes the true solution to be missed. The accuracy of 0.14 is obtained. In other fault diagnosis cases of KN20Bias, WL1Bias, and WCm5Bias, KN4 is wrongly diagnosed to worsen the diagnosis accuracies (Figure 12). Example: D1Stuck (Sticking of the D1 ClO2 Flow Control Valve). When the +10% set-point change occurred in the wood chips flow at 600 min, the D1 ClO2 flow control valve also stuck at the same time. The increase of wood chips increases the D1 temperature (T21) and E Kappa number (KN22), as shown in Figure 13a. The thin and thick lines of the figure are the measured and estimated values, respectively. PCA-based diagnosis suggests MV41 as one unique solution from 600 min. It represents the set-point change of wood chips flow. PCA has difficulties in regard to diagnosing multiple-


Figure 12. Residual of KN4 for (a) WL1Bias and (b) WCm5Bias.

Figure 13. (a) Dynamics of measured and estimated variables, (b) residual obtained by the DPLS models with time delay, and (c) residuals obtained by the DPLS models without time delay for D1CStuck.

fault, even if the occurred faults are independent of each other. Using the DPLS model, the detection sequence of symptoms is KN22(+) and KN222 from 915 min, T21(+) from 920 min, T212 from 935 min, BR262 from 1185 min, and BR26(-) from 1190 min (see Figure 13b). There is no false detection, and D1Cm5Bias, D1CStuck, WL1Bias, WLStuck, and WoodDens

were obtained as the solution from 920 min. Although the resolution is 5, the operator may easily judge that WL1Bias, WLStuck, and WoodDens are not fault candidates, because there are no detected variables in the digester and oxygen reactor section. Also, the resolution can be reduced using the dynamics of the residuals. Compare the residuals of Figure 13b with the

9072


Figure 14. Residuals of T21 and KN22 for D1Cm5Bias.

Figure 15. Dynamics of the residuals.

residuals for D1Cm5Bias (see Figure 14). The residuals for the control valve sticking increase linearly during 500 min; however, those for the control valve bias show almost constant values after the initial detection time. Also, the residuals dynamics may be obviously different, according to the type of the sensor faults. (Refer the example of BR26 shown in Figure 15a-c.) Except for the three faults of CFStuck, WLStuck, and WoodDens (see Figure 15d), the residuals dynamics can reduce the diagnosis resolution of other faulty cases. This method using the residuals dynamics may be considered as a further study. Figure 13c shows the residual obtained by the previous hybrid model without time delay. KN22(+) is detected from 995 min, KN222 and T21(+) from 1005 min, T212 from 1035 min, BR262 from 1055 min, and BR26(-) from 1105 min. The detection is 80 min later than the proposed method with time delay. Figure 13b and c shows that, evidently, the models with a time delay generate clearer and more accurate residuals that denote fault occurrence than those without a time delay. Diagnosis Result Comparison Table 5 shows the diagnosis result obtained by the proposed method and compares the result with that via DPCA. The performance parameter value given in brackets refers to the worst results during all diagnostic periods. The hybrid1 uses

the DPLS model, considering dead time, and the hybrid2 does not consider dead time. Because the three faults of BR1Deg, CFStuck, and PR1Deg failed via all three methods, they are not shown in the table. The diagnostic performance by hybrid1 for all cases except WoodDens and WoodTemp are much better than that by hybrid2. However, in some cases, hybrid2 has faster detections than hybrid1, because of the wrong detection of T12. In the WoodTemp case, the wood temperature increases and the symptom of T12(+) should be detected. However, hybrid2 wrongly detected T12(-), and the detection was earlier than that of hybrid1. In terms of resolution, this study introduced a new strategy to improve diagnostic performance. The online diagnosis procedure is summarized to obtain the minimum set of faults that can explain all of the detected symptoms. However, it was determined that the resolutions of several sensor faults were over 10, and modification of the strategy was needed to diagnose the sensor faults. For instance, consider the sensor fault of BR10Deg. When BR262 was detected, the final fault candidates were 24 faults, which are elements of the fault set of BR26 (see Table 4). This is due to the fact that, when many sensors were used as the source variables of the DPLS models, the faults occurred on the sensors belonged to the fault sets of the symptoms corresponding to the models. The resolution problems

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 9073 Table 5. Diagnosis Results Comparisons Detection Delay (min) fault BR0025Drift BR10Deg BR20Bias BRm5Bias CF1Bias CFm5Bias D1C1Bias D1CStuck D1Cm5Bias EbackTemp ECausComp ECausTemp KN0025Drift KN10Deg KN1Deg KN20Bias KNm5Bias OH0025Drift OH10Deg OH1Deg OH20Bias OHm5Bias PR0025Drift PR10Deg PR20Bias PRm5Bias WC1Bias WCStuck WCm5Bias WL1Bias WLStuck WLm5Bias WoodDens WoodTemp

Accuracy

PCA

hybrid1

hybrid2

115

215 40 20 20

560 20

65

65

315 20 15 65 15 115 55

395 20 15 65 15 105 45

35 35 65 25 45 25 25 145

25 35 85 25 45 25 25 175

35 45 35 35 15 25 10 10 110 115

35 45 35 35 35 25 10 20 60 70

5 5 10 5 85 5 0 50 0 55 10 45 10 10 50 10 10 10 10 60 20 20 20 20 0 20 85 100 10 85

PCA 0.99 0 1 0.44 0.41 0.49 0.82 0 1 0.33 1 0.17 0.26 0.91 0.47 0.16 0.16 1 0.98 0.44 1 1 1 0.76 0.99 0.71 0.83 1 0.84 0.91 0 0.89 1 0.96

hybrid1 1 1 1 1 0 1 0 1 0.99 0.22 1 1 1 1 0 1 1 1 1 1 0.98 1 1 0 1 1 1 1 0.02 0.91 1 0.85 0.26 1

Wrong Detection hybrid2 0 1 1 0 0 1 0 1 0.36 0.1 0.96 1 1 1 0 0.97 1 1 1 1 1 1 1 0 1 1 1 1 0.0 0.92 1 0.28 0.17 1

Resolution

hybrid1

hybrid2

hybrid1

hybrid2

0 0 0 0

0 0

4(4) 5(5) 4.0(5) 4.0(5)

5(5) 5(5)

0

0

2(2)

2(2)

0 0 0 0 0 0 0

0 0 0 0 0 0 0

5.0(7) 5.1(7) 6(6) 6.2(7) 5(5) 4(4) 4.0(5)

5.0(7) 8.1(9) 6(6) 5.7(7) 5(5) 4(4) 4.0(5)

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

4(4) 3(3) 4(4) 3.5(5) 4.8(5) 4(4) 3(3) 4.1(5)

4(4) 3(3) 4(4) 4.2(5) 4.8(5) 4(4) 3(3) 4.1(5)

1.1(2) 0 0 0 0 0 0 0 0 0

1.3(3) 0 0 0

4(4) 3(3) 6(6) 5.4(6) 4.7(6) 2.3(5) 4.5(5) 2.0(3) 4.6(6) 5.9(6)

4(4) 3(3) 6(6) 4.6(6)

0 0 0 0 0

2.1(5) 2.6(5) 2(2) 4.6(6) 7(8)

occurred when the number of the sensors corresponding to the detected symptoms was 1 or 2. To resolve the problem, the following strategy was used. If the symptom that corresponds to each sensor fault is not detected, the sesnor fault cannot be a fault candidate. For example, the symptom matching to BR10Deg is BR262. In addition, if the number of the sensors corresponding to the detected symptoms is no more than 2, and sensor faults are able to explain all of the detected symptoms, other faults except sensor faults are removed from the final fault candidates set. Although the new strategy greatly improved the diagnostic resolution, as shown in Table 5, the average accuracy of WL1Bias was dropped from 1 to 0.91. In the case, P37(+), P372, KN19(+), and KN192 were detected at 655-730 min. Because the two sensor faults of KN1Deg and KN10Deg could explain all symptoms during these 75 min, the true solution of WL1Bias was removed from the final fault candidates set. Therefore, the strategy should be carefully adopted, considering the resolution as well as the accuracy of the diagnosis.

In comparing the results between DPCA and hybrid1, for most fault cases, DPCA showed faster detection but much worse accuracy than hybrid1. As described previously, DPCA failed to diagnose valve sticking cases, except for WCStuck. In the case of WCStuck, the DPCA result of MV41 can represent the set-point change of wood chips flow, as well as the faults of the wood chips control valve. Therefore, it is difficult to conclude that the DPCA result is accurate in the case. Also, DPCA failed to detect BR10Deg, which was accurately diagnosed by the hybrid methods. Nonetheless, DPCA diagnosed the true solution for the three cases of CF1Bias, D1C1Bias, and KN1Deg, whereas hybrid1 failed. This is due to the fact that the 1% faults of CF1Bias, D1C1Bias, and KN1Deg are too small to be diagnosed by the DPLS models. Although T15 is independent of WCm5Bias and WoodDens, the proposed method wrongly detected T15 (Figure 16) for these faults. Because WLStuck can explain all detected symptoms for these two cases, the accuracy of hybrid1 was worse than DPCA. In

Figure 16. Residual of T15 for WoodDens.

Figure 17. Residual of C122 for CFStuck.

9074


addition, the method failed to diagnose PR10Deg and CFStuck (see Figure 17). These failures mean that the DPLS models of F3, T15, and C122 are not sufficiently accurate, because of noise or modeling errors, although the SDGs of three variables are composed, based on quasi-linear relations such as the control valve equation, mass-balance equation, and heat-transfer equation. Conclusion This study involves the fault diagnosis method for large-scale chemical processes that have a long time delay. It is based on the hybrid method, combining SDG and DPLS approaches, which has the advantage of improving the diagnosis resolution and accuracy and enhancing the ability to diagnose multiple faults. Many commercial chemical processes have unit processes with a time delay, and dead time from the input variables to the output variable, which may worsen fault diagnosis accuracy of the DPLS model. The previous hybrid method was modified to incorporate the information of the dead time into the DPLS models. The modified method was applied to the fault diagnosis of a pulp mill process, which is a realistic and very large benchmark process that includes common characteristics of the industrial processes, such as long measurement and time delay, many components, recycle streams, and high nonlinearity. In addition, the process specific physical properties such as Kappa number and brightness make fault diagnosis even more challenging. Sixteen DPLS models were constructed for the diagnosis of 20 faults in sensors, 12 faults in control valves, and 5 disturbances. The diagnosis results showed the hybrid method using the time delay information can greatly enhance the diagnostic accuracy and detection. Also, the proposed method was compared with DPCA and showed much better accuracy but slower detection than DPCA. The average accuracy of the hybrid method with time delay was the best among the three methods that were compared. When the DPCA, hybrid method with time delay, and hybrid method without time delay were applied for the fault diagnosis of 34 faults, the accuracies that were averaged for 34 cases were 0.67, 0.80, and 0.70, respectively. However, it failed to detect four faults, which were detected by DPCA. This can be resolved by improving the model accuracy. It was also found that the residual dynamics were obviously different, according to the fault types. To better the diagnosis resolution, the fault-type classification that is based on the residual dynamics needs to be studied in the future. Indices such as generalized likelihood ratio and cumulative variance index may be alternatives to classify fault types.17,19 Acknowledgment This work was supported by grant No. (R05-2002-000-000570) from the Basic Research Program of the Korea Science and Engineering Foundation. Nomenclature BPLS ) matrix of PLS regression coefficient E ) residual matrix for X F ) residual matrix for Y l ) number of past values N ) number of source variables P ) loading matrix for X

Q ) loading matrix for Y ri ) residual of variable i T ) score matrix for X U ) score matrix for Y W ) X-weight matrix X ) input matrix Y ) output matrix Yˆ ) predicted output matrix yi ) measured value of variable i yˆ i ) estimated value of variable i Greek Letters ∆t ) sampling rate τMD ) measurement delay τTD,ij ) time delay Literature Cited (1) Venkatasubramanian, V.; Rengaswamy, R.; Yin, K.; Kavuri, S. N. A Review of Process Fault Detection and Diagnosis Part I: Quantitative Model-based Methods. Comput. Chem. Eng. 2003, 27, 293-311. (2) Malthouse, E. C.; Tamhane, A. C.; Mah, R. S. H. Nonlinear Partial Least Squares. Comput. Chem. Eng. 1997, 21, 875-890. (3) Kazantzis, N.; Kravaris, C.; Wright, R. A. Nonlinear Observer Design for Process Monitoring. Ind. Eng. Chem. Res. 2000, 39, 408-419. (4) Han, Z.; Li, W.; Shah, S. L. Fault Detection and Isolation in the Presence of Process Uncertainties. Control Eng. Pract. 2005, 13, 587599. (5) Richard, J. P. Time-Delay Systems: an Overview of Some Recent Advances and Open Problems. Automatica 2003, 39, 1667-1694. (6) Jiang, C.; Zhou, D. H. Fault Detection and Identification for Uncertain Linear Time-Delay Systems. Comput. Chem. Eng. 2005, 30, 228242. (7) Lee, G.; Song, S.-O.; Yoon, E. S. Multiple-Fault Diagnosis Based on System Decomposition and Dynamic PLS. Ind. Eng. Chem. Res. 2003, 42, 6145-6154. (8) Downs, J. J.; Vogel, E. F. A Plant-wide Industrial Process Control Problem. Comput. Chem. Eng. 1993, 17, 245-255. (9) Chiang, L. H.; Russell, E. L.; Braatz, R. D. Fault Detection and Diagnosis in Industrial Systems; Springer: London, 2001. (10) Chen, G.; McAvoy, T. J. Predictive On-Line Monitoring of Continuous Processes. J. Process Control 1998, 8, 409-420. (11) Castro, J. J. Modeling and Control of a Pulp Mill Process, Ph.D. Dissertation, University of Delaware, Newark, DE, 2002. (12) Castro, J. J.; Doyle, F. J., III A Pulp Mill Benchmark Problem for Control: Problem Description. J. Process Control 2004, 14, 17-29. (13) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control; Wiley: New York, 1989. (14) Ahmed, S.; Huang, B.; Shah, S. L. Parameter and Delay Estimation of Continuous-Time Models Using a Linear Filter. J. Process Control 2006, 16, 323-331. (15) Lee, G.; Han, C.; Yoon, E. S. Multiple-Fault Diagnosis of the Tennessee Eastman Process Based on System Decomposition and Dynamic PLS. Ind. Eng. Chem. Res. 2004, 43, 8037-8048. (16) http://www.chemengr.ucsb.edu/∼ceweb/faculty/doyle/docs/benchmarks/mill. (17) Qin, S. J.; Li, W. Detection, Identification, and Reconstruction of Faulty Sensors with Maximized Sensitivity. AIChE J. 1999, 45, 19631976. (18) Ku, W.; Storer, R. H.; Georgakis, C. Disturbance Detection and Isolation by Dynamic Principal Component Analysis. Chemometrics Intell. Lab. Syst. 1995, 30, 179-196. (19) Qin, S. J.; Li, W. Detection and Identification of Faulty Sensors in Dynamic Processes. AIChE J. 2001, 47, 1581-1593.

ReceiVed for reView June 22, 2006 ReVised manuscript receiVed September 22, 2006 Accepted September 25, 2006 IE060793J

Fault Diagnosis Using the Hybrid Method of Signed Digraph and

Recommend Documents