AMethod for Designing Fault Diagnosis Filters for LPV Polytopic Systems

The work presented in this paper focuses on the design of robust Fault Detection and Isolation (FDI) filters for dynamic systems characterized by LPV (Linear Parameter Varying) polytopic models. A sufficient condition is established to guarantee sensitivity performance of the residual signal vector to faults. Robustness constraints against model perturbations and disturbances are also taken into account in the design method. A key feature of the proposed method is that the residual structuring matrices are optimized as an integral part of the design, together with the dynamic part (i.e. the filter). The design problem is formulated as a convex optimization problem and solved using LMI (Linear Matrix Inequalities) techniques. The proposed method is illustrated on the secondary circuit of a Nuclear Power Plant.


INTRODUCTION AND MOTIVATIONS
The issue of Fault Detection and Isolation (FDI) in dynamic systems has been an active research area in the last two decades.Model-based FDI techniques use mathematical models of the monitored process and extract features from measured signals, to generate fault indicating signals, that is, the residuals.LTI models have been widely used to solve the problem of FDI.Tools are now available to enhance robustness against small parameter variations and other disturbances (see, [1][2][3] for surveys).The resulting robust FDI problem is generally formulated as a min-max optimization setting to maximize fault sensitivity performance and at the same time, to minimize the influence of unknowns inputs.
More recently, some research works have appeared that consider Linear Parameter Varying (LPV) modeling of the monitored system to take into account wider and more rapid parameters variations.Such models can be used efficiently to represent some nonlinear systems (see, e.g., [4,5]).This motivates some researchers from the FDI community to develop model-based methods using LPV models (see [6][7][8] among others).The two commonly used approaches are fault estimation methods where the fault indicating signal is an es-timate of the fault signal, and residual generation methods where the residuals are synthesized to be robust against modeling errors and unknown inputs, while being sensitive to the faults.In this context, a geometric approach is proposed in [6] to design a LPV observer in a Luenberger form.A procedure is derived to obtain the observer parameters via the construction of a suitable family of invariant subspaces (parameter varying (C, A)-invariant and un-observability subspaces).In [7], a multi-model approach is used to solve the FDI problem for nonlinear systems.The nonlinear system is modeled using polytopic models and a robust polytopic unknown input observer is then synthesized by means of pole assignment.The method uses LMI optimization techniques to synthesize the observer gain.The major limitation of this approach is that sensitivity of the residual signal against faults can only be checked a posteriori.More precisely, if the distribution matrices of the fault model and the effects that faults could have on the decoupled state is not of full column rank, then faults could go undetectable.To overcome this problem, a solution is provided in [8] where the main idea is to build a fault estimate using a LPV filter such that the worst-case gain (i.e., the H ∞ performances measure for LPV systems) from disturbances and faults to the estimation error, is minimized.
In this paper a different approach based on residual generation is considered for LPV systems that can be modeled within a LPV polytopic setting.Robustness against exogenous disturbances and sensitivity against faults are considered in a framework similar to the well-known H ∞ /H − setting for LTI systems.The robustness objectives are expressed in terms of a minimization problem using the H ∞ norm for LPV systems, and the sensitivity requirement is formulated in terms of a maximization constraint using also the H-index for LPV systems.The main difference between this problem and the standard H ∞ problem for LPV systems is that it involves the residual structuring matrices that are unknown.
The paper is organized as follows: In Section 2, the general FDI filter design problem and the corresponding solution are presented.In Section 3, the proposed method is applied to real data set coming from the secondary circuit of a nuclear power plant in France.Finally, some concluding remarks are made in a final section.

Preliminaries
The Euclidean norm is always used for vectors and is written without a subscript; for example x .Similarly in the matrix case, the induced vector norm is used: where σ(A) denotes the maximum singular value of A. Signals, for example w(t) or w, are assumed to be of bounded energy, and their norm is denoted by w 2 , that is, w 2 = ( +∞ −∞ w(t) 2 dt) 1/2 < ∞.LTI models, for example, P(s) or simply P, are assumed to be in RH ∞ , real rational functions with P ∞ = sup ω σ(P( jω)) < ∞.P ∞ , that is, the largest gain of P, is accompanied by the smallest gain of P, inf ω σ(P( jω)), which may be equal to zero for some P (e.g., strictly proper systems), if the frequency range of interest is infinite.This motivated [1,[9][10][11][12] to define the nonzero smallest gain of P, that is, the H-index, as the restriction of inf ω σ(P( jω)) to a finite frequency domain Ω, that is, P − = inf ω∈Ω σ(P( jω)).
In [1,9], an evaluation function • e which is a restriction of the H 2 signal norm to Ω, is defined by the authors as w e = w 2,Ω = ( ω2 ω1 w( jω) 2  2 dω) 1/2 .Then, given P so that y = Pu, it follows that and thus that P − ≤ y e / u e .This motivates the introduction of an evaluation function, denoted • sens , which is defined according to: From (2), it follows that P sens takes the sense of the smallest value of a singular value of P( jω) evaluated on Ω.Then it follows that P − ≤ P sens ≤ P ∞ .
The underlying LPV system is modeled by the following state space representation which is denoted in a compact form as x is the state vector, u is the input vector, y is the output vector and θ(t) is a varying parameter vector.It is assumed that all parameters θ i (t), i = 1, . . ., r are bounded, measurable (or estimated) in real time and take their values in the domain Θ, so that Θ is a convex polytope.
The LPV system (3) admits a (non-conservatism) polytopic model if it is possible to determine a set of matrices M i , i = 1, . . ., N, constituting the vertices of a polytope defined by and such that it corresponds to the image by M of the domain Θ: Then, β i , i = 1, . . ., N define barycentric coordinates of Θ and the following convex decomposition yields: Referring to the LPV system (3), the worst-case RMS gain from u to y which is known as the H ∞ -norm for LPV systems is defined by: Following the definition of the index • sens of a LTI transfer given by (2), we will introduce the following evaluation function, that will be useful in the following to formulate fault sensitivity requirements for LPV fault detection schemes: M(θ) sens is also a generalization of (2) to LPV case.

Problem setting
Consider the general FDI design problem for LPV systems represented on Figure 1.G(θ) is a polytopic LPV model.
3 K is a known controller.d ∈ R qd represents exogenous disturbances and f ∈ R q f represents the faults to be detected.F(θ) is the FDI LPV filter to be designed.z ∈ R qr is an estimation of z = M y y + M u u, a subset of the measurements y ∈ R m and the controlled inputs u ∈ R p .M y ∈ R qr ×m and M u ∈ R qr ×p are the two residual structuring matrices to be designed.
It is assumed that the problem depicted on Figure 1 is well posed and thus the lower fractional transformation F 1 (G(θ), K) always exists.
The FDI design problem we are interested in can then be formulated as follows.
Problem 1. Assume that the faults f are detectable (the interested reader can refer to [13] for a discussion on fault detectability).
The goal is to find the state space matrices m+p) of the (stable) polytopic LPV filter F(θ) and the residual structuring matrices M y ∈ R qr ×m and M u ∈ R qr ×p defining the residual vector (10) such that the residual vector meets the following requirements: T r f (θ) sens > γ 2 (12) where T rd and T r f denote the looped transfers between d and r and f and r, respectively.This problem can be represented by the block diagram illustrated on Figure 2, where P(θ) is derived from F 1 (G(θ), K) so that y u = P(θ) d f .In this formulation, γ 1 and γ 2 are two positive constants referring respectively to the robustness and sensitivity performances levels.
Equation (11) represents the worst-case robustness of the residual to disturbances d for all θ ∈ Θ, in the H ∞ -norm sense.Under plant perturbation, the effect that the exogenous disturbances have on the residuals, can greatly increase and the fault detection performance may then be considerably degraded.A robust fault sensitivity specification is then needed to maintain a detection performance level of the FDI unit.Here the sensitivity measure (9) for LPV fault detection scheme is used to guarantee the worst-case sensitivity of the residual to faults.Of course, the smaller γ 1 and the bigger γ 2 are the better the fault detection performance will be.

Remark 1.
In Problem 1 formulation, it is assumed that the structuring matrices M y and M u do not depend on θ.If this assumption vanishes, it can be verified that the following theoretical developments still yields.The only difference in such a case is that, if we consider M y (θ) and M u (θ) in Problem 1, a set of structuring matrices M y (Π i ) and M u (Π i ) for each vertex of the polytope Θ would be obtained rather than constant matrices.

Design of the FDI filter
In this section, a solution is provided to compute simultaneously M y , M u and F(θ) so that the requirements ( 11) and ( 12) are satisfied.It is straightforward to verify that the major difficulty in this problem is related to the fault sensitivity requirement ( 12) since ( 11) can be solved using the techniques developed in the robust control community (see, e.g., [14] or [15]).To overcome this problem, a sufficient condition is established in terms of a fictitious H ∞ problem.It is then shown in the following that a solution to this fictitious problem is a solution of the original one.

Standard setup for the filter design problem
To proceed, let Using some algebra manipulations, the filter design problem illustrated on Figure 2, can be re-casted into the setup depicted in Figure 3, where P(θ, M y , M u ) is deduced from P(θ), M y and M u , according to: where I qr and 0 i× j denote respectively the identity matrix of dimension q r and the null matrix of dimension i × j. n is also the order of P(θ, M y , M u ), that is, A(θ) ∈ R n×n .Following the method proposed in [10,11], the requirements ( 11) and ( 12) are now expressed in terms of loop shapes, that is, of desired gain responses for the appropriate closed-loop transfers.These shaping objectives are then turned into uniform bounds by means of the shaping filters.Let W d and W f denote the (dynamical) shaping filters associated with (11) and ( 12) respectively, so that: Assume that W d is invertible (this can be done without loss of generality because it is always possible to add zeros in W d to make it invertible).W d and W f are also defined in order to tune the gain responses for, respectively, T rd (θ) and T r f (θ).Then it is straightforward to verify that the specification (11) yields if the following constraint is satisfied: In this formulation, d ∈ R q d is a fictitious signal generating d through W d (see Figure 4(a) for easy reference) and T r d denotes the looped transfer between r and d.
The following lemma allows the sensitivity constraint (12) to be transformed into a fictitious H ∞ one.

Lemma 1. Consider an invertible transfer matrix
Define the (fictitious) signal r ∈ R qr such that r = r −W F f (see Figure 4(a)).Then a sufficient condition for the specification (12) to hold, is: where T r f denote the looped transfer between r and f .Proof of Lemma 1.Consider the signal r introduced in Figure 4, that is, where W F is define as in Lemma 1. Then it can be verified that the following relation yields: that can be re-written due to the definition of r given by ( 19): Now consider the weighting function W F defined in Lemma 1. Since, W F is supposed to be invertible, we get where W + F ∞,Ω = sup ω∈Ω σ(W + F ( jω))•W + F denotes the inverse of W F which always exists by assumption (see Lemma 1).Then, factorizing the right term of (21) by that can be done since, by definition, W F − / =0.With (22), it then follows that: Now, since by construction W F − > λ, it is straightforward to verify that the following relation yields: Suppose now that inequality (18) yields, that is, and with (24), we get which terminates the proof.
Following (17) and ( 18), the design problem can be re-casted in a fictitious LPV H ∞ -framework as depicted in Figure 4 The residual generation problem can now be formulated in a framework which looks like a standard H ∞ problem for LPV systems, by combining both requirements ( 17) and ( 18) into a single H ∞ constraint.Using Lemma 1, it can be verified that a sufficient condition for ( 17) and ( 18) to hold, is As mentioned above, this equation seems to be similar to a standard H ∞ LPV problem.In fact, this is not the case since the transfer P(θ, M y , M u ) depends on M y and M u that are unknown.In the following section, a procedure is given to overcome this problem.Remark 2. It is clear that a key feature in the proposed formulation is the a priori choice of the shaping filters W d and W f .From a practical point of view, it is required that the residuals r are as "big" and as "fast" as possible, when a fault occurs.Then, if the considered faults manifest themselves in low frequencies, this leads to select W f as a low pass filter with the static gain and the cutting frequency, the highest possible.With regards to the robustness objectives, it is required that the effects of the disturbances on the residuals are as "small" as possible.This implies to choose the gain of W d as "small" as possible in the frequency range where the energy content of the disturbances is located.In other words, it is required a high attenuation level of the disturbances on the residuals in the appropriate frequencies (see Section 3 where a practical case is presented).However, both sensitivity to faults and robustness against disturbances might be not achieved in some cases.Faults having similar frequency characteristics as those of disturbances might go undetected.In such cases, the proposed formulation provides a framework to find a good balance between fault sensitivity and robustness via the construction of the shaping filters W d and W f .Finally, note that the work reported in [16,17] could be an interesting method to select W d and W f .

Synthesis of the FDI filter
In the following, a solution is derived in terms of a SDP (Semi Definite Programming) problem.To proceed, let {A wd , B wd , C wd , D wd } and {A wF , B wF , C wF , D wF } be the statespace representations of W −1 d and W F respectively, and denote n wd and n wF the associated order, that is, A wd ∈ R nwd×nwd and A wF ∈ R nwF ×nwF .Using some linear algebra manipulations, it can be verified from ( 14) that the matrices A(θ), B(θ), C(θ) and D(θ) of the state-space representation of P(θ, M y , M u ) are defined according to: Having in mind the definition of A(θ), B(θ), C(θ), and D(θ), it can be noted that the H ∞ optimization problem formulated by (30) is non convex since it involves simultaneously the residual structuring matrices M y and M u and the filter state-space matrices A F (θ), B F (θ), C F (θ), and D F (θ).A solution to this problem may consist in chosen heuristically M y and M u .However, as it has been outlined in [10], there is no guarantee to the optimal solution.The following lemma which is an adaptation of Proposition 1 in [8] for our purpose, gives the solution to this problem.The proof is omitted here as it can be found in [8].

. , N be the evaluation of A(θ), B(θ), C(θ), D(θ)
at each vertex Π i of the polytope Θ. Assume that C 2 (θ) and D 21 (θ) do not depend of θ (see Remark 3 for a discussion about this assumption) and let N s = ( C 2 D 21 ) ⊥ .Then, there exists a solution of (31) if and only if there exists γ < 1 and M y ∈ R qr ×m and M u ∈ R qr ×p and two symmetric matrices R ∈ R (n+nwd+nwF )×(n+nwd+nwF ) > 0 and S ∈ R (n+nwd+nwF )×(n+nwd+nwF ) > 0 solving the following SDP problem involving 2N + 1 LMI constraints: Moreover, F(θ) is a full order filter if n F = n + n wd + n wF .The state space realization of the LPV filter F(θ) is then computed using the barycentric coordinates of Θ given by (7), so that: , and D F (Π i ) ∀i = 1, . . ., N are the state space matrices of the N LTI filters F(Π i ) ∀i = 1, . . ., N that are deduced from the unique solution (R, S, M y , M u , γ) following the procedure described in [14].
Remark 3. As it is outlined in Lemma 2, it is required that C 2 and D 21 do not depend on θ.This assumption is also done for N S to be computed.If, by construction, such an assumption is not verified, the solution consists in post filtering (y T u T ) T by a LTI filter at a high cutting frequency.This solution has already been proposed in [14].

APPLICATION: THE SECONDARY CIRCUIT OF A NPP
To illustrate the benefits of the proposed approach, experimental results obtained from the secondary circuit of a Nuclear Power Plant (NPP) are provided in this section.In a NPP, the secondary circuit erosion can occur in the steam generators, releasing radio nuclides into the secondary coolant.This problem is now well understood and has been the subject of some studies initiated by EDF (Eléctricité De France) for its pressurized water reactors (PWR) to overcome and master the steam generator corrosion problems.The main degradation process is to be controlled by careful optimization of the secondary water chemistry.There is a need to ensure that the optimum secondary chemistry regime is selected and maintained.So, the process of erosion-corrosion of steel piping and other components is of critical importance during operation of a NPP.Feed water pH is adjusted by hydrazine, so that the pH is maintained within the limits specified by the nuclear authority (norm ISO-14253-1).
The NPP under consideration is a 900 MW pressurized water reactors (PWR), located in France.During the win- ter 2002, and thanks to a mandatory period of maintenance operations, it has been possible to measure and record a set of experimental data on the secondary circuit.The aim was to optimally control hydrazine and pH through an adaptive LQG control scheme.The designed controller has been successfully tested under real operational conditions [18].
The dynamics of hydrazine and ammoniac concentrations behavior in the secondary circuit of the NPP can be expressed as follows: V is the circuit's water volume, Q ext the water extraction flow rate of the condenser and β is a parameter depending of the NPP operating conditions.[N 2 H 4 ], [NH 3 ], and [O 2 ] represent hydrazine, ammoniac and oxygen concentrations respectively.u is the flow rates of the pumps used to inject hydrazine in the circuit.Moreover, the system present a time delay (τ = 560 s) corresponding to the chemical reaction after the introduction of hydrazine in the circuit.It is assumed that the pH is measured, that is, where n pH denotes the measurement noise, [Mo] also denotes the morpholine concentration, K b and K Mo are the basicity constants of the ammoniac and morpholine respectively.
Taken into account the dynamical equations (39), it follows that where ) T represents the state vector.n N2H4 represents the measurement noise of the hydrazine sensor.n NH3 is the image of the measurement noise of the pH sensor via the relation (40).Then, since (40) is static we will consider that an ammoniac concentration measure is available through the pH sensor.Characteristics of n NH3 are then deduced from the following equation, which is the inverse of (40): Figure 5 presents the behavior of the time varying parameter Q ext (t) during a period of 3 days.We assume here that Q ext (t) is not affected by the considered faults.As it can be seen on Figure 5, Q ext (t) varies between Q min and Q max with: The considered faults are hydrazine sensor faults and pH sensor faults.Note that monitoring of pH sensors is a key feature for well functioning the overall system.
The relation between the ammoniac measurement and the pH measurement is only algebraic (see (41)), then the pH sensor fault is directly transmitted to the ammoniac measurement.Therefore, we can consider an ammoniac sensor fault in place of a pH sensor fault.Consequently, the following state space representation derived from (42) models the failing behavior of the secondary circuit (the notations are chosen to be consistent with those used in Section 2): f NH3 represent ammoniac sensor faults and f N2H4 hydrazine sensor faults.In this model, the time delay due to actuators is approximated using a first order Pade approximation.The model ( 44) corresponds to the model described in Figure 1, where 878 ≤ θ(t) = Q ext (t) ≤ 1097.It follows that the considered polytope Θ = {θ : 878 ≤ θ ≤ 1097} becomes a simple segment since dim (θ) = 1.
For the FDI purpose, two filters F 1 (θ) and F 2 (θ) are computed such that the two residuals r 1 (t) and r 2 (t) satisfy the following requirements: (i) r 1 (t) is sensitive to pH sensor faults through the (fictitious) ammoniac sensor and robust to hydrazine sensor faults.(ii) r 2 (t) is sensitive to hydrazine sensor faults and robust to ammoniac sensor faults.
This method allows to isolate sensor faults uniquely.Therefore, we consider a different model for each filter to be synthesized.For the design of F 1 (θ), the following model is used where the disturbances vector d 1 includes the oxygen concentration, the measurement noises and the hydrazine sensor fault, that is, For the design of F 2 (θ), the following model is retained where the augmented disturbances vector takes into account the ammoniac sensor faults; According to the methodology developed in the Section 2, two polytopic models P 1 (θ) and P 2 (θ) are built as illustrated on Figure 2. To save place and for a better understanding, the different steps are only detailed for F 1 (θ).Here, because the system is placed in an open-loop control law, it follows from G 1 (θ) that (for clarity the index "1" is ignored from now): where Following the developments in Section 2, the fault detector design problem turns out to be the design of (F(θ), M y , M u ) satisfying the following objectives: Figure 6 gives an illustration of this problem.Finally, following the method describes in Section 2, the problem is recasted into the setup depicted in Figure 3 where the model P(θ, M y , M u ) is defined according to: M y1 and M y2 are the two components of the structuring matrix M y .
As it is outlined in Remark 2 an important step in the proposed method is the choice of the shaping filters W d and W f .Similarly to the developments presented in Section 2, W d is refereed to the robustness objectives against d and W f to the sensitivity requirements against f NH3 .Here, due to the definition of d, it is natural to choose W d such as: The weighting functions W u , W and W f allow to manage separately the robustness objectives against u, [O 2 ], n NH3 , n N2H4 , f N2H4 and f NH3 respectively.The interested reader can refer to [10] or [19] if necessary.These weighting functions are deduced from an off line spectral analysis procedure of available measurements according to: The parameters γ u , γ [O2] , γ nNH 3 , γ nN 2 H 4 , γ [N2H4] and γ 2 allow to manage the gain of the weighting functions separately.They are optimized by performing an iterative refinement.Remember that the goal is to minimize the effects of disturbances on the residual r(t) and maximize the effects of faults on r(t).The numerical values of them have been fixed to The method described in Section 2.2.2 is then used to synthesize the filter F(θ), and the structuring matrices M y and M u .
For the SDP optimization problem computation, the SDPT3 solver is used.
To analyze the computed solution, the principal gains σ(T k dr ( jω)) and σ(T f NH 3 r ( jω)) of the closed loop transfers T k dr ( jω) and T f NH 3 r ( jω) are plotted versus the objectives W k d and W f for some θ ∈ Θ (see Figure 7).The notation "k" is introduced to outline that the analysis is performed with respect to each component of d.As it can be seen on the figures, for each synthesis, σ(T k dr ( jω)) < σ(W k d ( jω)) ∀ω and σ(T f NH 3 r ( jω)) > σ(W f ( jω)) ∀ω ∈ Ω ≈ [0; 10 −4 [rd/s for all considered values of θ(t).This indicates that the requirements (48) are satisfied for the considered values of θ and by virtue of Lemma 2, we know that it still yields for all values of θ.Finally, a sequential Wald decision test is also implemented within the simulator to make a final decision about the faults.The probabilities of non-detection and false alarms have been fixed to 0.1%.The results are presented in Figures 8, 9 and 10.As it can be seen, all faults are successfully detected and isolated.

CONCLUSION
The problem which is addressed in this paper is that of designing FDI filters for dynamic systems that can be described by LPV polytopic models.The method can be seen as a generalization of the well known H ∞ /H − setting for LTI systems.The H ∞ norm for LPV systems is used to formulate the ro-  bustness specifications and a new index, that is, the H sens index, which is deduced from the H-norm for LTI systems, is introduced for fault sensitivity specifications.As a result, various design goals and trades-off can be formulated and managed in a systematic way by means of some high level design parameters formulated in terms of dynamic weighting functions.A key feature of the proposed technique is that the remaining control and measurement canals are optimally merged to build the fault indicating signals.The resulting static matrices are also optimized via LMI together with the dynamic FDI filter.The proposed technique is also appropriate for fault diagnosis in nonlinear systems which can be approximated efficiently by LPV models to cover a wider range of operating, and to cope with rapid parameter variations.The method has been successfully applied to experimental data set coming from the secondary circuit of a nuclear power plant in France.

Figure 1 :
Figure 1: The general FDI filter design problem.

Figure 2 :
Figure 2: General setup for FDI/LPV filter design problem.
(a), where d and r are two fictitious signals, so that d = W d d and r = r − W F f .Then, including γ 1 , λ, W F and W −1 d into the model P(θ, M y , M u ) leads to the equivalent block diagram of Figure 4(b), where P(θ, M y , M u ) is defined according to:

Figure 4 :
Figure 4: Fictitious quadratic H ∞ formulation for the filter design problem.

Figure 8 :
Figure 8: Behavior of r i (t), i = 1, 2 and the decision test-fault-free situation.