- Email: [email protected]
RESIDUAL GENERATOR IDENTIFICATION AND DESIGN FOR LINEAR MULTIVARIABLE SYSTEMS
RESIDUAL GENERATOR IDENTIFICATION AND DESIGN FOR LINEAR MULTIVARIABLE SYSTEMS Silvio Simani ∗ Roberto Diversi ∗∗ ∗
Dipartimento di Ingegneria, Universit`a di Ferrara Via Saragat 1, 44100 Ferrara, Italy Tel: +39 0532 97 4844, fax: +39 0532 97 4870 e-mail: [email protected] ∗∗ Dipartimento di Elettronica, Informatica e Sistemistica, Universit`a di Bologna Viale del Risorgimento 2, 40136 Bologna, Italy e-mail: [email protected]
Abstract: This paper investigates the design of residual generators in order to realize complete diagnosis schemes in linear multivariable systems with additive faults and disturbances. The use of canonical input–output polynomial forms leads to characterise in a straightforward fashion the basis of the subspace described by all the possible residual generators. These tools show how the mathematical description of these filters can be obtained by following also a black–box identification approach. An example is finally c 2006 IFAC. reported. Copyright Keywords: Fault detection, linear multivariable systems, input–output polynomial models, parameter identification, minimum redundancy.
1. INTRODUCTION
elled in terms of input–output polynomial description, so that the residual generation problem can be reduced to the determination of the null–space of a specific polynomial matrix associated to the process model.
The model–based approach to fault detection and isolation in complex technological systems has received considerable attention during the last 25 years (Patton et al., 1989; Chen and Patton, 1999; Simani et al., 2002).
In particular, the use of canonical input–output forms allows to compute in a straightforward fashion an analytical expression for the basis of such a null– space and upper and lower bounds for the order of the dynamic residual generator. Another advantage of the proposed method can be found in the implementation of black–box identification ideas, that is, the discrete– time residual generator can be obtained without any knowledge of the mathematical model of the process under investigation.
There are a great variety of methods in the literature which can be classified on the basis of the particular mathematical description used for the monitored process. The origins of these methods can be found in parity space methodologies, observed–based approaches and parameter estimation techniques, with many cross connections among the different approaches. This work investigates the residual generation problem for linear multivariable systems with additive faults and disturbances, by following the minimal polynomial approach suggested in (Frisk and Nyberg, 2001). The process under investigation is mod-
890
Let us consider a linear, finite–dimensional, time– invariant, discrete–time system described by the following input–output equation P (z) y(t) = Q(z) u(t) ,
t = 0, 1, . . .
⎡
⎤ c(t) ˜ d (z) Q ˜ f (z) ⎣ d(t) ⎦ , (8) ˜ c (z) Q P˜ (z) y(t) = Q f (t)
2. PROBLEM FORMULATION
where c(t) is the c –dimensional known–input vector, d(t) is the d –dimensional disturbance vector, f (t) is the f –dimensional monitored fault vector and c + d + f = . It is worth noting that, when only additive faults fc (t) on the input sensors of the system are considered, the input vector measurements can be written as
(1)
where z −1 is the unitary delay operator and P (z) and Q(z) are polynomial matrices with dimension (m × m) and (m × ) respectively, with P (z) nonsingular. The terms u(t) and y(t) are, respectively, the – dimensional and m–dimensional input and output vectors of the considered multivariable system. Models of type (1) can be frequently found in practice by applying well–known physical laws to describe the input–output dynamical links of various systems and are a powerful tool in all fields where the system state does not play a direct role, such as identification, decoupling, output controllability, etc. For a generic input–output model P (z), Q(z) , its canonical input–output form is the equivalent rep ˜ resentation P˜ (z), Q(z) satisfying the following properties:
c∗ (t) = c(t) + fc (t)
(9)
˜ c (z)c∗ (t) + and Eq. (8) becomes P˜ (z)y(t) = Q ˜ c (z)fc (t). Analogously, when only ad˜ d (z)d(t) − Q Q ditive faults fo (t) on the output sensors of the system are considered the output vector measurements can be written as y ∗ (t) = y(t) + fo (t) .
(10)
˜ c (z)c(t) + In this case, it results that P˜ (z)y ∗ (t) = Q ˜ ˜ Qd (z)d(t) + P (z)fo (t).
deg p˜ii (z) > deg p˜ji (z) i = j
(2)
3. RESIDUAL GENERATOR DESIGN
deg p˜ii (z) > deg p˜ij (z) j > i
(3)
deg p˜ii (z) ≥ deg p˜ij (z) j < i
(4)
A general linear residual generator for the fault detection process of system (8) is a filter of type
deg p˜ii (z) ≥ deg q˜ij (z)
(5)
∀j,
R(z) r(t) = Sy (z) y(t) + Sc (z) c(t) .
with the polynomials p˜ii (z) monic. A constructive proof of the existenceand uniqueness of a canonical form for a given pair P (z), Q(z) can be found in (Beghelli and Guidorzi, 1976). Note that the integers νi = deg p˜ii (i = 1, . . . , m) equals the corresponding row–degrees of P˜ (z). ˜ The canonical representation P˜ (z), Q(z) directly leads to a correspondent canonical state–space realization of type ˜ u(t) x(t + 1) = A˜ x(t) + B ˜ ˜ y(t) = Cx(t) + Du(t) .
(6)
System (11) processes the known input–output data and generates the residual r(t), i.e. a signal which is “small” (ideally zero) in the fault–free case and is “large” when a fault is acting on the system. Without loss of generality, r(t) can be assumed to be a scalar signal. In such condition R(z) is a polynomial with degree greater than or equal to the row–degree of Sc (z) and Sy (z), in order to guarantee the physical realisability of the filter. Moreover, if R(z) has all roots inside the unit circle filter (11) is asymptotically stable.
(7)
An important aspect of the design of the residual generator consists of the decoupling of the disturbance d(t) in order to produce a correct diagnosis in all operating conditions. In absence of faults, Equation (8) can be rewritten in the form
with order n=
m
νi .
(11)
i=1
˜ d (z) d(t) . ˜ c (z) c(t) = Q P˜ (z) y(t) − Q
The integers νi are the ordered setof Kronecker invari˜ ˜ ants associatedto the pair A, C of every observable realization of P (z), Q(z) (Guidorzi, 1975).
(12)
Premultiplying all the terms in (12) by a row polynomial vector L(z) belonging to the left null–space of ˜ d (z)), we obtain ˜ d (z), N (Q Q
In order to design residual generators for the fault diagnosis system, model (1) can be firstly transformed ˜ into its canonical representation P˜ (z), Q(z) , sat˜ isfying conditions (2)–(5). Then, matrix Q(z) can be decomposed according to the following structure
˜ c (z) c(t) = 0 . L(z) P˜ (z) y(t) − L(z) Q
(13)
Starting from Eq. (13) it is possible to obtain a residual of type (11) by setting:
891
Sy (z) = L(z) P˜ (z) ˜ c (z) Sc (z) = −L(z) Q nf R(z) = z ,
all the residual generators (11) for the fault detection process of system (8) can be obtained by replacing in relation (14) L(z) with B(z), i.e.
(14)
˜ d (z) adj Q ˜ d (z) P˜1 (z) − det Q ˜ d (z) P˜2 (z) Sy (z) = Q 2 1 1
where nf is the maximal row–degree of the pair ˜ c (z) . The choice R(z) = z nf L(z) P˜ (z), L(z) Q leads to an asymptotical stable filter with nf poles equal to zero. It can be observed that Equation (13) involves the input–output samples c(t+i), y(t+i) with i = 0, . . . , nf . Thus, Equation (13) can be rewritten also in the form
˜ d (z) adj Q ˜ d (z) Q ˜ c (z) + det Q ˜ d (z) Q ˜ c (z) Sc (z) = −Q 1 2 2 1 1 n n (21) nf R(z) = diag z f1 z f2 . . . z m−d ,
where nfi (i = 1, . . . , m − d ) is the row–degree of Syi (z), where Syi (z) denotes the i–th row of matrix Sy (z). It can be noted that relation (5) leads to the following inequality: (22) row deg Syi (z) ≥ row deg Sci (z) ,
(15) r(t + nf ) = z nf r(t) ˜ c (z) c(t). = L(z) P˜ (z) y(t) − L(z) Q When a fault is acting on the system the residual generator is governed by the relation ˜ f (z) f (t) r(t + nf ) = −L(z) Q
where Sci (z) denotes the i–th row of matrix Sc (z) so that the residual generator (11), (21) is physically realizable.
(16)
and r(t + nf ) assumes values that are, in general, different from zero. Matrix L(z) in Equation (13) can be designed by taking into account good fault sensitivity properties for r(t+nf ). For example, it can be noted that in absence of disturbances Qd (z) = 0, so ˜ d (z)) coincides with the whole vector space. that N (Q Consequently, a set of residual generators of the fault detection process of system (8) can be expressed as:
Previous considerations can be summarised as follows. A linear generator for the fault detection process of system (8) is a dynamical system (11) whose minimal order nf ∗ is constrained in the following range νmin ≤ n∗f ≤ min (d + 1) νmax , n . (23) The lower bound can be obtained in the no–disturbance case (d = 0) from relation (17) by selecting the row of P˜ (z) associated to the minimal Kronecker invariant.
˜ c (z) c(t),(17) ri (t + νi ) = z νi r(t) = P˜i (z) y(t) − Q i ˜ c (z) ((i = 1, 2, . . . , m)) are the where P˜i (z) and Q i ˜ c (z), respectively, i–th rows of matrices P˜ (z) and Q ˜ ˜ c (z) cannot and νi is the row–degree of Pi (z), since Q i show a greater row–degree. In order to determine all possible residual generators of minimal order it is ˜ d (z)). necessary to compute a minimal basis of N (Q ˜ Under the assumption that matrix Qd (z) is of full ˜ d (z)) has dimension ˜ d (z) = d , N (Q rank, i.e. rank Q m − d and a minimal basis of it can be computed as suggested in (Kailath, 1980).
4. RESIDUAL IDENTIFICATION In this section, the problem of identifying minimal order residual generators for systems of type (8) is investigated. It will be shown that these residuals are insensitive to disturbance signals and can be exploited for both model–based and model–free fault detection purposes. It is assumed that, in absence of faults, a finite sequence of the variables yi (t) (i = 1, . . . , m) and cj (t) (j = 1 . . . , c ) with t = 1, . . . , N is available. If linear relations are present among these variables, they can be described by difference equations of the type
˜ d (z) can be In general, for 1 < d < m matrix Q partitioned in the following way:
˜ ˜ d (z) = Qd1 (z) , (18) Q ˜ d (z) Q 2
mi m
˜ d (z) have dimension ˜ d (z) and Q where matrices Q 1 2 d × d and (m − d ) × d respectively. It can be as˜ d (z) sumed, without loss of generality, that matrix Q 1 is non singular. In this case it can be easily verified ˜ d (z)) is given by the polynomial that a basis of N (Q matrix
i=1 k=0
αik yi (t + k) +
nj c
βjk cj (t + k)
j=1 k=0
= Sy (z) y(t) + Sc (z) c(t) = 0,
(24)
where Sy (z) = sy1 (z) sy2 (z) · · · sym (z) Sc (z) = sc1 (z) sc2 (z) · · · scc (z) ,
˜ d (z) −det Q ˜ d (z) Im− . ˜ d (z) adj Q B(z) = Q 2 1 1 d (19)
(25)
with syi (z) and scj (z) representing polynomials and αi mi and βj nj their coefficients.
˜ c (z) as Q ˜ d (z) in (18) By partitioning P˜ (z) and Q
˜ P˜ (z) ˜ c (z) = Qc1 (z) , (20) Q P˜ (z) = ˜1 ˜ c (z) P2 (z) Q 2
The residual identification problem can be stated as follows. Given the fault–free input–output sequences
892
k¯ = max {mi , (i = 1, . . . , m)
c(·), y(·) generated by a system of type (8) determine the orders mi (i = 1, . . . , m), nj (j = 1, . . . , c ) and the parameters αik , βjk of Equation (24). Define now the following vectors and matrices T yiL (t) = yi (t) . . . yi (t + L − 1) T cL = cj (t) . . . cj (t + L − 1) j (t) (26) XhL (yi ) = yiL (1) . . . yiL (h + 1) L XhL (cj ) = cL j (1) . . . cj (h + 1) ,
nj , (j = 1, . . . , c )} = nf .
¿From Equations (27) and (29) it follows that the integer L must satisfy: L≥
(27)
and compute the sample covariance matrix ΣL (h1 h2 . . . hm+c ) = (28) 1 L = H (h1 h2 . . . hm+c )T H L (h1 h2 . . . hm+c ). L
i = 1, . . . , m
(35)
(29) This residual function is of type (15) and hence insensitive to the disturbance term d(t). Moreover, the test performed on sequence (32) ensures that such residual is of minimal order. Since the number of linearly independent residual generators is m − d , if m > d +1 it is possible to exploit again the algorithm ¯ ¯ . . . , k+1), until above, starting from matrix ΣL (k+1, a second singular matrix is encountered and another residual generator is computed. In the same manner all m − d residuals can be obtained.
(31)
j = 1, . . . , c .
When no faults are present in the process r(t + nf ) is ideally zero, while in presence of faults relation (24) does not longer hold so that r(t + nf ) = 0. The faults affecting the monitored system can thus be detected by comparing the residual function (36) with a fixed threshold ε according to the simple threshold logic: |r(t + nf )| ≤ ε for fault–free case, (37) |r(t + nf )| > ε for faulty cases.
In the following, vector θ, solution of Eq. (29), will be considered with a unitary Euclidean norm. On the basis of these considerations, the previous problem can be solved by means of the algorithm described below. (1) Consider the sequence of symmetrical increasing dimension matrices: ΣL (1, . . . , 1), ΣL (2, . . . , 2), . . .
nj + m + c ,
j=1
r(t + nf ) = z nf r(t) = Sy (z) y(t) + Sc (z) c(t) .(36)
where the entries of θ are the coefficients of the polynomials:
T (30) θ = θyT1 · · · θyTm θcT1 · · · θcTc
T θyi = αi0 αi1 · · · αimi
T θcj = βj0 βi1 · · · βjnj
c
Once that model (24) has been identified, a residual function for system (8) can be immediately obtained as:
If a relation like (24) holds, it is immediate to verify that Σ (m1 . . . mm n1 . . . nc ) θ = 0,
mi +
in order to avoid unwanted linear dependence relations due to limitations in the dimension of the involved vector spaces in relation (29). Note also that, due to (34), the number of available samples N must be ¯ In order to perform greater than or equal to L + k. correctly the independence test in the algorithm above, the Hankel matrix XhLm+1 (c1 ) . . . XhLm+c (cc ) must be of full rank, i.e. the known inputs c1 , c2 , . . . , cc must be persistently exciting of sufficient orders (identifiability conditions). A check of this rank should thus be included in step (1).
XhLm (ym )
L
m i=1
for i = 1, . . . , m, j = 1, . . . , c . Define also the Hankel matrix H L (h1 h2 . . . hm+c ) = [ XhL1 (y1 ) . . . XhLm+1 (c1 ) . . . XhLm+c (cc ) ,
(34)
(32)
There are many ways of defining residual functions and determining thresholds (Chen and Patton, 1999). For example, in (37) such a function has been chosen as a norm of the residual r(t + nf ). The value of the threshold ε depends on the considered system and can be determined from the knowledge of the process operating conditions. It is important to note that Equation (29) can be directly used for fault detection purposes without computing model (24). In fact, it is possible to determine matrix ΣL only on the basis of the input and output sequences y(·), c(·).
and test the linear independence of their columns L ¯ ¯ as long as a singular matrix ΣL ¯ = Σ (k, . . . , k) k L is encountered. The rank of Σk¯ is equal to (m + c )k¯ − 1. (2) Compute a basis for the null space of ΣL ¯: k ker ΣL ¯ | · · · | αm0 · · · αmk ¯| ¯ = α10 · · · α1k k T β10 · · · · · · β1k¯ | · · · | βc 0 · · · βc k¯ . (33) (3) By deleting the redundant null coefficients in (33) (if present) it is possible to determine the orders mi and nj of the mathematical description (24) and the parameters vector θ. Obviously
A residual function may thus be represented by the least eigenvalue λ(ΣL ) of ΣL . The fixed threshold
893
ε is related to the least eigenvalue of ΣL computed in fault–free conditions. Therefore, the fault detection criterion may be performed by monitoring the values of λ(ΣL ) according to the following logic: λ(ΣL ) ≤ ε for fault–free case (38) λ(ΣL ) > ε for faulty cases
u1 (t) u2 (t) u3 (t) u4 (t) u5 (t)
y1 (t) y2 (t) y3 (t)
The direct identification technique of the residual generator described in Section 4 has been applied to the 120MW power plant of Pont–sur–Sambre (Guidorzi and Rossi, 1974).
2
5
H.P.
M.P.
as the residual generator has to be sensitive to the inputs u1 (t), u3 (t) and u5 (t). Therefore, as described in Section 4, the identification of the residual generator is performed by testing the singularity of the symmetrical increasing dimension matrices (32). In this case, the first singular matrix encountered is:
B.P. 6
Ry 3
as the residual generator has to be decoupled from the inputs u2 (t) and u4 (t) and ˜ c (z) = Q ˜ 3 (z), Q ˜ 5 (z) ˜ 1 (z), Q (40) Q
Os Trs
4
9
350 (9, 9, 9, 9, 9, 9), ΣL ¯ =Σ k
7 8
(41)
whose null space computed in the Matlab environment corresponds (under the Matlab floating point relative accuracy) to the normalised coefficients of the polynomial terms of the matrices Sc (z) and Sy (z). Moreover, it is worth noting how these matrices are equal (under the Matlab floating point accuracy) to the ones computed by means of relations (14) with nf = 9 and L(z) given by
Cb Qa
Fig. 1. The structure of the power plant. The block–diagram of the plant is shown in Figure 1 where the numbers refer to: (1) (2) (3) (4) (5) (6) (7) (8) (9)
steam pressure main steam temperature reheat steam temperature
In the first example, the process inputs u2 (t) and u4 (t) have been considered as disturbance signals d1 (t) and d2 (t), respectively. Therefore, using the notation introduced in Section 3, the following matrices can be computed: ˜ d (z) = Q ˜ 4 (z) , ˜ 2 (z), Q (39) Q
It is a double–shaft industrial gas turbine working in parallel with the electrical mains. This industrial process consists mainly of three major components: the reactor, turbine, and condenser. Furthermore, there are several pumps, valves (not highlighted) and one turbine. The boiler boils the water and the steam generated drives the turbine. After the turbine, the condenser cools the steam. In turn, external cooling water cools the condenser. The cooling pumps transport the water from the condenser tank back to the boiler tank.
1
Pv Ts Trs
˜ ˜ A canonical input–output description P(z) , Q(z) of the dynamic process was obtained by means of the identification algorithm for multivariable dynamic systems, where the Kronecker invariants are ν1 = 4, ν2 = 3, ν3 = 2. The input–output data (367 samples) used for the direct identification of the residual generators have been obtained from the simulation described by the identified model of the model ˜ ˜ P(z) , Q(z) under fault–free conditions.
5. APPLICATION EXAMPLE
Ts , Pv
gas flow turbine valves opening super heater spray flow gas dampers air flow
The process output variables (m = 3) are:
for increasing values of L as time t increases.
Qd
Cb Os Qd Ry Qa
super heater (radiation); super heater (convection); super heater; re-heater; dampers; condenser; drum; water pump; burner.
L(z) =
−0.033 + 0.19z − 0.38z 2 + . . . −0.024 + 0.016z − 0.042z 2 + 0.14z 3 + . . . −0.0098 + 0.055z − 0.07z 2 − 0.015z 3 + . . .
. . . + 0.23z 3 + 0.0018z 4 − 0.0056z 5 . . . + 0.5z 4 − 0.69z 5 + 0.046z 6 . . . + 0.078z 4 − 0.023z 5 − 0.043z 6 + 0.029z 7
T (42)
but computed by the Matlab null space (Forney, 1975; ˜ d (z) in the form of EquaChen, 1984) of the matrix Q tion (39).
The process input variables (r = 5) are:
894
In the second example, the process inputs u2 (t) has been considered as disturbance signal d1 (t). Therefore, under this assumption: ˜ d (z) = Q ˜ 2 (z) , (43) Q and ˜ c (z) = Q ˜ 3 (z), Q ˜ 4 (z), Q ˜ 5 (z) , (44) ˜ 1 (z), Q Q
2.5 2 1.5 1 0.5 0 −0.5
in order to obtain a residual generator sensitive to the input signals u1 (t), u3 (t), u4 (t) and u5 (t) but decoupled from the disturbance signal d1 (t) = u2 (t). Therefore, in such a case, the identification of the two residual generators is performed again by testing the singularity of the symmetrical increasing dimension matrices (32). In this situation, the first singular matrix encountered is: 350 (4, 4, 4, 4, 4, 4, 4), ΣL ¯ =Σ k
−1 50
150
200
250
300
Data Samples
350
400
described by all residual generators. Moreover, it is shown that the order and the parameters of these filters can be identified directly from a finite number of input–output samples describing the behaviour of the process in absence of faults. The results obtained in the simulation of real cases are encouraging.
(45)
(46)
whose null spaces computed in the Matlab environment corresponds again to the coefficients of the polynomial terms of the row of the matrices Sc (z) and Sy (z). Again, as described in Section 4, it is worth noting that these matrices are equal to the ones that may be computed by means of relations (14) with nf1 = 4 and nf2 = 5. On the other hand, when one disturbance only has to be decoupled, L(z) has the form: T 2 2 L(z) =
100
Fig. 2. Fault on the first output when one disturbance (ld = 1) has been decoupled.
whilst the second singular matrix: 350 (5, 5, 5, 5, 5, 5, 5) ΣL ¯ =Σ k
r(t)
−12
x 10 3
REFERENCES
−0.029z − 0.19z −0.048z 0.46 + 0.62z + 0.49z 2 −0.21z + 0.44z 2 −0.11z + 0.35z 2 0.42 − 0.69z + 0.32z 2
(47) that can be evaluated by means of the null of the matrix ˜ d (z) given by Equation (43). Q In order to assess the diagnostic capabilities of the suggested techniques, additive faults have been simulated by using the dynamic model of the process in Figure 1. As an example, the residual signals regarding two different faulty situations are reported in Figure 2. Continuous lines represent the fault–free residual functions generated by the disturbance de– coupled dynamic filters, while the dotted ones depict the faulty signals. The faults (ramp functions) have been added to the first output measurement of the considered process, starting at the data sample t = 160. The performance of the residual generator depicted in Figure 2 seems to assess the diagnostic capabilities of the suggested technique. 6. CONCLUSION Residual generation problem for linear multivariable systems with additive faults and disturbances has been analysed in this paper. The use of canonical input– output polynomial representations to model the considered dynamical process leads to a simple characterisation of the polynomial basis of the subspace
895
Beghelli, S. and R. Guidorzi (1976). A new input– otput canonical form for multivariable systems. IEEE Transactions on Automatic Control, AC– 21, 692–696. Chen, C. T. (1984). Linear Systems Theory and Design. Holt, Rinehart and Winston, New York. Chen, J. and R. J. Patton (1999). Robust Model–Based Fault Diagnosis for Dynamic Systems. Kluwer Academic Publishers. Forney Jr., G. D. (1975). Minimal bases of rational vector spaces with applications to multivariable linear systems. SIAM Journal on Control, 13, 493–520. Frisk, E. and M. Nyberg (2001). A minimal polynomial basis solution to residual generation for fault diagnosis in linear systems. Automatica, 37, 1417–1424. Guidorzi, R. (1975). Canonical structures in the identification of multivariable systems. Automatica, 11, 361–374. Guidorzi, R. and R. Rossi (1974). Identification of a power plant from normal operating records. Automatic Control Theory and Applications, 2, 63–67. Kailath, T. (1980). Linear systems. Prentice Hall, Englewood Cliffs, New Jersey. Patton, R. J., Frank, P. M. and R. N. Clark (1989). Fault Diagnosis in Dynamic Systems, Theory and Application. Control Engineering Series, Prentice Hall, London. Simani, S., Fantuzzi, C. and R. J. Patton (2002). Model-based fault diagnosis in dynamic systems using identification techniques. Advances in Industrial Control, Springer–Verlag, London.
Dipartimento di Ingegneria, Universit`a di Ferrara Via Saragat 1, 44100 Ferrara, Italy Tel: +39 0532 97 4844, fax: +39 0532 97 4870 e-mail: [email protected] ∗∗ Dipartimento di Elettronica, Informatica e Sistemistica, Universit`a di Bologna Viale del Risorgimento 2, 40136 Bologna, Italy e-mail: [email protected]
Abstract: This paper investigates the design of residual generators in order to realize complete diagnosis schemes in linear multivariable systems with additive faults and disturbances. The use of canonical input–output polynomial forms leads to characterise in a straightforward fashion the basis of the subspace described by all the possible residual generators. These tools show how the mathematical description of these filters can be obtained by following also a black–box identification approach. An example is finally c 2006 IFAC. reported. Copyright Keywords: Fault detection, linear multivariable systems, input–output polynomial models, parameter identification, minimum redundancy.
1. INTRODUCTION
elled in terms of input–output polynomial description, so that the residual generation problem can be reduced to the determination of the null–space of a specific polynomial matrix associated to the process model.
The model–based approach to fault detection and isolation in complex technological systems has received considerable attention during the last 25 years (Patton et al., 1989; Chen and Patton, 1999; Simani et al., 2002).
In particular, the use of canonical input–output forms allows to compute in a straightforward fashion an analytical expression for the basis of such a null– space and upper and lower bounds for the order of the dynamic residual generator. Another advantage of the proposed method can be found in the implementation of black–box identification ideas, that is, the discrete– time residual generator can be obtained without any knowledge of the mathematical model of the process under investigation.
There are a great variety of methods in the literature which can be classified on the basis of the particular mathematical description used for the monitored process. The origins of these methods can be found in parity space methodologies, observed–based approaches and parameter estimation techniques, with many cross connections among the different approaches. This work investigates the residual generation problem for linear multivariable systems with additive faults and disturbances, by following the minimal polynomial approach suggested in (Frisk and Nyberg, 2001). The process under investigation is mod-
890
Let us consider a linear, finite–dimensional, time– invariant, discrete–time system described by the following input–output equation P (z) y(t) = Q(z) u(t) ,
t = 0, 1, . . .
⎡
⎤ c(t) ˜ d (z) Q ˜ f (z) ⎣ d(t) ⎦ , (8) ˜ c (z) Q P˜ (z) y(t) = Q f (t)
2. PROBLEM FORMULATION
where c(t) is the c –dimensional known–input vector, d(t) is the d –dimensional disturbance vector, f (t) is the f –dimensional monitored fault vector and c + d + f = . It is worth noting that, when only additive faults fc (t) on the input sensors of the system are considered, the input vector measurements can be written as
(1)
where z −1 is the unitary delay operator and P (z) and Q(z) are polynomial matrices with dimension (m × m) and (m × ) respectively, with P (z) nonsingular. The terms u(t) and y(t) are, respectively, the – dimensional and m–dimensional input and output vectors of the considered multivariable system. Models of type (1) can be frequently found in practice by applying well–known physical laws to describe the input–output dynamical links of various systems and are a powerful tool in all fields where the system state does not play a direct role, such as identification, decoupling, output controllability, etc. For a generic input–output model P (z), Q(z) , its canonical input–output form is the equivalent rep ˜ resentation P˜ (z), Q(z) satisfying the following properties:
c∗ (t) = c(t) + fc (t)
(9)
˜ c (z)c∗ (t) + and Eq. (8) becomes P˜ (z)y(t) = Q ˜ c (z)fc (t). Analogously, when only ad˜ d (z)d(t) − Q Q ditive faults fo (t) on the output sensors of the system are considered the output vector measurements can be written as y ∗ (t) = y(t) + fo (t) .
(10)
˜ c (z)c(t) + In this case, it results that P˜ (z)y ∗ (t) = Q ˜ ˜ Qd (z)d(t) + P (z)fo (t).
deg p˜ii (z) > deg p˜ji (z) i = j
(2)
3. RESIDUAL GENERATOR DESIGN
deg p˜ii (z) > deg p˜ij (z) j > i
(3)
deg p˜ii (z) ≥ deg p˜ij (z) j < i
(4)
A general linear residual generator for the fault detection process of system (8) is a filter of type
deg p˜ii (z) ≥ deg q˜ij (z)
(5)
∀j,
R(z) r(t) = Sy (z) y(t) + Sc (z) c(t) .
with the polynomials p˜ii (z) monic. A constructive proof of the existenceand uniqueness of a canonical form for a given pair P (z), Q(z) can be found in (Beghelli and Guidorzi, 1976). Note that the integers νi = deg p˜ii (i = 1, . . . , m) equals the corresponding row–degrees of P˜ (z). ˜ The canonical representation P˜ (z), Q(z) directly leads to a correspondent canonical state–space realization of type ˜ u(t) x(t + 1) = A˜ x(t) + B ˜ ˜ y(t) = Cx(t) + Du(t) .
(6)
System (11) processes the known input–output data and generates the residual r(t), i.e. a signal which is “small” (ideally zero) in the fault–free case and is “large” when a fault is acting on the system. Without loss of generality, r(t) can be assumed to be a scalar signal. In such condition R(z) is a polynomial with degree greater than or equal to the row–degree of Sc (z) and Sy (z), in order to guarantee the physical realisability of the filter. Moreover, if R(z) has all roots inside the unit circle filter (11) is asymptotically stable.
(7)
An important aspect of the design of the residual generator consists of the decoupling of the disturbance d(t) in order to produce a correct diagnosis in all operating conditions. In absence of faults, Equation (8) can be rewritten in the form
with order n=
m
νi .
(11)
i=1
˜ d (z) d(t) . ˜ c (z) c(t) = Q P˜ (z) y(t) − Q
The integers νi are the ordered setof Kronecker invari˜ ˜ ants associatedto the pair A, C of every observable realization of P (z), Q(z) (Guidorzi, 1975).
(12)
Premultiplying all the terms in (12) by a row polynomial vector L(z) belonging to the left null–space of ˜ d (z)), we obtain ˜ d (z), N (Q Q
In order to design residual generators for the fault diagnosis system, model (1) can be firstly transformed ˜ into its canonical representation P˜ (z), Q(z) , sat˜ isfying conditions (2)–(5). Then, matrix Q(z) can be decomposed according to the following structure
˜ c (z) c(t) = 0 . L(z) P˜ (z) y(t) − L(z) Q
(13)
Starting from Eq. (13) it is possible to obtain a residual of type (11) by setting:
891
Sy (z) = L(z) P˜ (z) ˜ c (z) Sc (z) = −L(z) Q nf R(z) = z ,
all the residual generators (11) for the fault detection process of system (8) can be obtained by replacing in relation (14) L(z) with B(z), i.e.
(14)
˜ d (z) adj Q ˜ d (z) P˜1 (z) − det Q ˜ d (z) P˜2 (z) Sy (z) = Q 2 1 1
where nf is the maximal row–degree of the pair ˜ c (z) . The choice R(z) = z nf L(z) P˜ (z), L(z) Q leads to an asymptotical stable filter with nf poles equal to zero. It can be observed that Equation (13) involves the input–output samples c(t+i), y(t+i) with i = 0, . . . , nf . Thus, Equation (13) can be rewritten also in the form
˜ d (z) adj Q ˜ d (z) Q ˜ c (z) + det Q ˜ d (z) Q ˜ c (z) Sc (z) = −Q 1 2 2 1 1 n n (21) nf R(z) = diag z f1 z f2 . . . z m−d ,
where nfi (i = 1, . . . , m − d ) is the row–degree of Syi (z), where Syi (z) denotes the i–th row of matrix Sy (z). It can be noted that relation (5) leads to the following inequality: (22) row deg Syi (z) ≥ row deg Sci (z) ,
(15) r(t + nf ) = z nf r(t) ˜ c (z) c(t). = L(z) P˜ (z) y(t) − L(z) Q When a fault is acting on the system the residual generator is governed by the relation ˜ f (z) f (t) r(t + nf ) = −L(z) Q
where Sci (z) denotes the i–th row of matrix Sc (z) so that the residual generator (11), (21) is physically realizable.
(16)
and r(t + nf ) assumes values that are, in general, different from zero. Matrix L(z) in Equation (13) can be designed by taking into account good fault sensitivity properties for r(t+nf ). For example, it can be noted that in absence of disturbances Qd (z) = 0, so ˜ d (z)) coincides with the whole vector space. that N (Q Consequently, a set of residual generators of the fault detection process of system (8) can be expressed as:
Previous considerations can be summarised as follows. A linear generator for the fault detection process of system (8) is a dynamical system (11) whose minimal order nf ∗ is constrained in the following range νmin ≤ n∗f ≤ min (d + 1) νmax , n . (23) The lower bound can be obtained in the no–disturbance case (d = 0) from relation (17) by selecting the row of P˜ (z) associated to the minimal Kronecker invariant.
˜ c (z) c(t),(17) ri (t + νi ) = z νi r(t) = P˜i (z) y(t) − Q i ˜ c (z) ((i = 1, 2, . . . , m)) are the where P˜i (z) and Q i ˜ c (z), respectively, i–th rows of matrices P˜ (z) and Q ˜ ˜ c (z) cannot and νi is the row–degree of Pi (z), since Q i show a greater row–degree. In order to determine all possible residual generators of minimal order it is ˜ d (z)). necessary to compute a minimal basis of N (Q ˜ Under the assumption that matrix Qd (z) is of full ˜ d (z)) has dimension ˜ d (z) = d , N (Q rank, i.e. rank Q m − d and a minimal basis of it can be computed as suggested in (Kailath, 1980).
4. RESIDUAL IDENTIFICATION In this section, the problem of identifying minimal order residual generators for systems of type (8) is investigated. It will be shown that these residuals are insensitive to disturbance signals and can be exploited for both model–based and model–free fault detection purposes. It is assumed that, in absence of faults, a finite sequence of the variables yi (t) (i = 1, . . . , m) and cj (t) (j = 1 . . . , c ) with t = 1, . . . , N is available. If linear relations are present among these variables, they can be described by difference equations of the type
˜ d (z) can be In general, for 1 < d < m matrix Q partitioned in the following way:
˜ ˜ d (z) = Qd1 (z) , (18) Q ˜ d (z) Q 2
mi m
˜ d (z) have dimension ˜ d (z) and Q where matrices Q 1 2 d × d and (m − d ) × d respectively. It can be as˜ d (z) sumed, without loss of generality, that matrix Q 1 is non singular. In this case it can be easily verified ˜ d (z)) is given by the polynomial that a basis of N (Q matrix
i=1 k=0
αik yi (t + k) +
nj c
βjk cj (t + k)
j=1 k=0
= Sy (z) y(t) + Sc (z) c(t) = 0,
(24)
where Sy (z) = sy1 (z) sy2 (z) · · · sym (z) Sc (z) = sc1 (z) sc2 (z) · · · scc (z) ,
˜ d (z) −det Q ˜ d (z) Im− . ˜ d (z) adj Q B(z) = Q 2 1 1 d (19)
(25)
with syi (z) and scj (z) representing polynomials and αi mi and βj nj their coefficients.
˜ c (z) as Q ˜ d (z) in (18) By partitioning P˜ (z) and Q
˜ P˜ (z) ˜ c (z) = Qc1 (z) , (20) Q P˜ (z) = ˜1 ˜ c (z) P2 (z) Q 2
The residual identification problem can be stated as follows. Given the fault–free input–output sequences
892
k¯ = max {mi , (i = 1, . . . , m)
c(·), y(·) generated by a system of type (8) determine the orders mi (i = 1, . . . , m), nj (j = 1, . . . , c ) and the parameters αik , βjk of Equation (24). Define now the following vectors and matrices T yiL (t) = yi (t) . . . yi (t + L − 1) T cL = cj (t) . . . cj (t + L − 1) j (t) (26) XhL (yi ) = yiL (1) . . . yiL (h + 1) L XhL (cj ) = cL j (1) . . . cj (h + 1) ,
nj , (j = 1, . . . , c )} = nf .
¿From Equations (27) and (29) it follows that the integer L must satisfy: L≥
(27)
and compute the sample covariance matrix ΣL (h1 h2 . . . hm+c ) = (28) 1 L = H (h1 h2 . . . hm+c )T H L (h1 h2 . . . hm+c ). L
i = 1, . . . , m
(35)
(29) This residual function is of type (15) and hence insensitive to the disturbance term d(t). Moreover, the test performed on sequence (32) ensures that such residual is of minimal order. Since the number of linearly independent residual generators is m − d , if m > d +1 it is possible to exploit again the algorithm ¯ ¯ . . . , k+1), until above, starting from matrix ΣL (k+1, a second singular matrix is encountered and another residual generator is computed. In the same manner all m − d residuals can be obtained.
(31)
j = 1, . . . , c .
When no faults are present in the process r(t + nf ) is ideally zero, while in presence of faults relation (24) does not longer hold so that r(t + nf ) = 0. The faults affecting the monitored system can thus be detected by comparing the residual function (36) with a fixed threshold ε according to the simple threshold logic: |r(t + nf )| ≤ ε for fault–free case, (37) |r(t + nf )| > ε for faulty cases.
In the following, vector θ, solution of Eq. (29), will be considered with a unitary Euclidean norm. On the basis of these considerations, the previous problem can be solved by means of the algorithm described below. (1) Consider the sequence of symmetrical increasing dimension matrices: ΣL (1, . . . , 1), ΣL (2, . . . , 2), . . .
nj + m + c ,
j=1
r(t + nf ) = z nf r(t) = Sy (z) y(t) + Sc (z) c(t) .(36)
where the entries of θ are the coefficients of the polynomials:
T (30) θ = θyT1 · · · θyTm θcT1 · · · θcTc
T θyi = αi0 αi1 · · · αimi
T θcj = βj0 βi1 · · · βjnj
c
Once that model (24) has been identified, a residual function for system (8) can be immediately obtained as:
If a relation like (24) holds, it is immediate to verify that Σ (m1 . . . mm n1 . . . nc ) θ = 0,
mi +
in order to avoid unwanted linear dependence relations due to limitations in the dimension of the involved vector spaces in relation (29). Note also that, due to (34), the number of available samples N must be ¯ In order to perform greater than or equal to L + k. correctly the independence test in the algorithm above, the Hankel matrix XhLm+1 (c1 ) . . . XhLm+c (cc ) must be of full rank, i.e. the known inputs c1 , c2 , . . . , cc must be persistently exciting of sufficient orders (identifiability conditions). A check of this rank should thus be included in step (1).
XhLm (ym )
L
m i=1
for i = 1, . . . , m, j = 1, . . . , c . Define also the Hankel matrix H L (h1 h2 . . . hm+c ) = [ XhL1 (y1 ) . . . XhLm+1 (c1 ) . . . XhLm+c (cc ) ,
(34)
(32)
There are many ways of defining residual functions and determining thresholds (Chen and Patton, 1999). For example, in (37) such a function has been chosen as a norm of the residual r(t + nf ). The value of the threshold ε depends on the considered system and can be determined from the knowledge of the process operating conditions. It is important to note that Equation (29) can be directly used for fault detection purposes without computing model (24). In fact, it is possible to determine matrix ΣL only on the basis of the input and output sequences y(·), c(·).
and test the linear independence of their columns L ¯ ¯ as long as a singular matrix ΣL ¯ = Σ (k, . . . , k) k L is encountered. The rank of Σk¯ is equal to (m + c )k¯ − 1. (2) Compute a basis for the null space of ΣL ¯: k ker ΣL ¯ | · · · | αm0 · · · αmk ¯| ¯ = α10 · · · α1k k T β10 · · · · · · β1k¯ | · · · | βc 0 · · · βc k¯ . (33) (3) By deleting the redundant null coefficients in (33) (if present) it is possible to determine the orders mi and nj of the mathematical description (24) and the parameters vector θ. Obviously
A residual function may thus be represented by the least eigenvalue λ(ΣL ) of ΣL . The fixed threshold
893
ε is related to the least eigenvalue of ΣL computed in fault–free conditions. Therefore, the fault detection criterion may be performed by monitoring the values of λ(ΣL ) according to the following logic: λ(ΣL ) ≤ ε for fault–free case (38) λ(ΣL ) > ε for faulty cases
u1 (t) u2 (t) u3 (t) u4 (t) u5 (t)
y1 (t) y2 (t) y3 (t)
The direct identification technique of the residual generator described in Section 4 has been applied to the 120MW power plant of Pont–sur–Sambre (Guidorzi and Rossi, 1974).
2
5
H.P.
M.P.
as the residual generator has to be sensitive to the inputs u1 (t), u3 (t) and u5 (t). Therefore, as described in Section 4, the identification of the residual generator is performed by testing the singularity of the symmetrical increasing dimension matrices (32). In this case, the first singular matrix encountered is:
B.P. 6
Ry 3
as the residual generator has to be decoupled from the inputs u2 (t) and u4 (t) and ˜ c (z) = Q ˜ 3 (z), Q ˜ 5 (z) ˜ 1 (z), Q (40) Q
Os Trs
4
9
350 (9, 9, 9, 9, 9, 9), ΣL ¯ =Σ k
7 8
(41)
whose null space computed in the Matlab environment corresponds (under the Matlab floating point relative accuracy) to the normalised coefficients of the polynomial terms of the matrices Sc (z) and Sy (z). Moreover, it is worth noting how these matrices are equal (under the Matlab floating point accuracy) to the ones computed by means of relations (14) with nf = 9 and L(z) given by
Cb Qa
Fig. 1. The structure of the power plant. The block–diagram of the plant is shown in Figure 1 where the numbers refer to: (1) (2) (3) (4) (5) (6) (7) (8) (9)
steam pressure main steam temperature reheat steam temperature
In the first example, the process inputs u2 (t) and u4 (t) have been considered as disturbance signals d1 (t) and d2 (t), respectively. Therefore, using the notation introduced in Section 3, the following matrices can be computed: ˜ d (z) = Q ˜ 4 (z) , ˜ 2 (z), Q (39) Q
It is a double–shaft industrial gas turbine working in parallel with the electrical mains. This industrial process consists mainly of three major components: the reactor, turbine, and condenser. Furthermore, there are several pumps, valves (not highlighted) and one turbine. The boiler boils the water and the steam generated drives the turbine. After the turbine, the condenser cools the steam. In turn, external cooling water cools the condenser. The cooling pumps transport the water from the condenser tank back to the boiler tank.
1
Pv Ts Trs
˜ ˜ A canonical input–output description P(z) , Q(z) of the dynamic process was obtained by means of the identification algorithm for multivariable dynamic systems, where the Kronecker invariants are ν1 = 4, ν2 = 3, ν3 = 2. The input–output data (367 samples) used for the direct identification of the residual generators have been obtained from the simulation described by the identified model of the model ˜ ˜ P(z) , Q(z) under fault–free conditions.
5. APPLICATION EXAMPLE
Ts , Pv
gas flow turbine valves opening super heater spray flow gas dampers air flow
The process output variables (m = 3) are:
for increasing values of L as time t increases.
Qd
Cb Os Qd Ry Qa
super heater (radiation); super heater (convection); super heater; re-heater; dampers; condenser; drum; water pump; burner.
L(z) =
−0.033 + 0.19z − 0.38z 2 + . . . −0.024 + 0.016z − 0.042z 2 + 0.14z 3 + . . . −0.0098 + 0.055z − 0.07z 2 − 0.015z 3 + . . .
. . . + 0.23z 3 + 0.0018z 4 − 0.0056z 5 . . . + 0.5z 4 − 0.69z 5 + 0.046z 6 . . . + 0.078z 4 − 0.023z 5 − 0.043z 6 + 0.029z 7
T (42)
but computed by the Matlab null space (Forney, 1975; ˜ d (z) in the form of EquaChen, 1984) of the matrix Q tion (39).
The process input variables (r = 5) are:
894
In the second example, the process inputs u2 (t) has been considered as disturbance signal d1 (t). Therefore, under this assumption: ˜ d (z) = Q ˜ 2 (z) , (43) Q and ˜ c (z) = Q ˜ 3 (z), Q ˜ 4 (z), Q ˜ 5 (z) , (44) ˜ 1 (z), Q Q
2.5 2 1.5 1 0.5 0 −0.5
in order to obtain a residual generator sensitive to the input signals u1 (t), u3 (t), u4 (t) and u5 (t) but decoupled from the disturbance signal d1 (t) = u2 (t). Therefore, in such a case, the identification of the two residual generators is performed again by testing the singularity of the symmetrical increasing dimension matrices (32). In this situation, the first singular matrix encountered is: 350 (4, 4, 4, 4, 4, 4, 4), ΣL ¯ =Σ k
−1 50
150
200
250
300
Data Samples
350
400
described by all residual generators. Moreover, it is shown that the order and the parameters of these filters can be identified directly from a finite number of input–output samples describing the behaviour of the process in absence of faults. The results obtained in the simulation of real cases are encouraging.
(45)
(46)
whose null spaces computed in the Matlab environment corresponds again to the coefficients of the polynomial terms of the row of the matrices Sc (z) and Sy (z). Again, as described in Section 4, it is worth noting that these matrices are equal to the ones that may be computed by means of relations (14) with nf1 = 4 and nf2 = 5. On the other hand, when one disturbance only has to be decoupled, L(z) has the form: T 2 2 L(z) =
100
Fig. 2. Fault on the first output when one disturbance (ld = 1) has been decoupled.
whilst the second singular matrix: 350 (5, 5, 5, 5, 5, 5, 5) ΣL ¯ =Σ k
r(t)
−12
x 10 3
REFERENCES
−0.029z − 0.19z −0.048z 0.46 + 0.62z + 0.49z 2 −0.21z + 0.44z 2 −0.11z + 0.35z 2 0.42 − 0.69z + 0.32z 2
(47) that can be evaluated by means of the null of the matrix ˜ d (z) given by Equation (43). Q In order to assess the diagnostic capabilities of the suggested techniques, additive faults have been simulated by using the dynamic model of the process in Figure 1. As an example, the residual signals regarding two different faulty situations are reported in Figure 2. Continuous lines represent the fault–free residual functions generated by the disturbance de– coupled dynamic filters, while the dotted ones depict the faulty signals. The faults (ramp functions) have been added to the first output measurement of the considered process, starting at the data sample t = 160. The performance of the residual generator depicted in Figure 2 seems to assess the diagnostic capabilities of the suggested technique. 6. CONCLUSION Residual generation problem for linear multivariable systems with additive faults and disturbances has been analysed in this paper. The use of canonical input– output polynomial representations to model the considered dynamical process leads to a simple characterisation of the polynomial basis of the subspace
895
Beghelli, S. and R. Guidorzi (1976). A new input– otput canonical form for multivariable systems. IEEE Transactions on Automatic Control, AC– 21, 692–696. Chen, C. T. (1984). Linear Systems Theory and Design. Holt, Rinehart and Winston, New York. Chen, J. and R. J. Patton (1999). Robust Model–Based Fault Diagnosis for Dynamic Systems. Kluwer Academic Publishers. Forney Jr., G. D. (1975). Minimal bases of rational vector spaces with applications to multivariable linear systems. SIAM Journal on Control, 13, 493–520. Frisk, E. and M. Nyberg (2001). A minimal polynomial basis solution to residual generation for fault diagnosis in linear systems. Automatica, 37, 1417–1424. Guidorzi, R. (1975). Canonical structures in the identification of multivariable systems. Automatica, 11, 361–374. Guidorzi, R. and R. Rossi (1974). Identification of a power plant from normal operating records. Automatic Control Theory and Applications, 2, 63–67. Kailath, T. (1980). Linear systems. Prentice Hall, Englewood Cliffs, New Jersey. Patton, R. J., Frank, P. M. and R. N. Clark (1989). Fault Diagnosis in Dynamic Systems, Theory and Application. Control Engineering Series, Prentice Hall, London. Simani, S., Fantuzzi, C. and R. J. Patton (2002). Model-based fault diagnosis in dynamic systems using identification techniques. Advances in Industrial Control, Springer–Verlag, London.