Introduction

The fault diagnosis of industrial processes has become an important topic because of its great influence on the operational control of processes. Reliable diagnosis and early detection of incipient faults avoid harmful consequences. Typically, faults in sensors and actuators and the process itself are considered. In the case of the wind turbine benchmark, a set of pre-defined faults with different locations and types are proposed in [1] where the dynamic change in the pitch system is treated. The procedure of fault detection is based either on the knowledge or on the model of the system [2]. Model-based fault detection is often necessary to obtain a good performance in the diagnosis of faults. The methods used in model-based diagnosis can be classified according if they are using state observers, parity equations and parameter estimation [3]. For linear time invariant systems (LTI), the FD task is largely solved by powerful tools. However, physical systems generally present nonlinear behaviors. Using LTI models in many real applications is not sufficient for high performance design. In order to achieve good performance while using linear like techniques, Linear Parameter Varying systems are recently received considerable attention [4]. Recently, many model-based applications using such systems and the subspace identification method were published [5]. In model-based FD, a residual vector is used to describe the consistency check between the predicted and the real behavior of the monitored system. Ideally, the residuals should only be affected by the faults. However, the presence of disturbances, noises and modeling errors yields the residual to become non zero. To take into account these errors, the fault detection algorithm Figure 1: Fault diagnosis with set estimator schema must be robust. When modeling uncertainty in a deterministic way, there are two robust estimation methods: the first method is the bounded error estimation that assumes the parameters are considered time invariant and there is only an additive error [6]. On the other hand, the second approach is the interval predictor that takes into account the variation of parameters and which considers additive and multiplicative errors [7], [8]. Here, the interval predictor is combined with existing nominal LPV identification presented by [9], allowing to include robustness and minimizing false alarms (see Fig. 1) [10]. Thus, this paper contributes with a new set-membership estimator approach that combines the interval predictor scheme with the LPV identification through subspace methods in one step. To illustrate the methodology proposed in this work, the pitch subsystem of wind turbine system proposed as a benchmark in IFAC SAFEPROCESS 2009 will be used. First, this subsystem is modeled as an LPV model using the hydraulic pressure as the scheduling variable. On the hypothesis that damping ratio and natural frequency have an affine variation with hydraulic pressure, this affine LPV model is estimated by means of the subspace LPV estimation algorithm. Second, the residue is synthesized to take into account the robustness against the uncertainties in the parameters. This work is organized as follows: In Section 2, the LPV subspace estimation method is recalled. In Section 3, the interval predictor approach combined with the LPV subspace method is proposed as tool for robust fault detection. In Section 4, the modeling of the pitch system as a LPV model is introduced. Section 5 deals with simulation experiments that illustrate the implementation and performance of the proposed approach applied to the robust fault detection of wind turbine pitch system. Fi-nally, Section 6 gives some concluding remarks.

LPV Subspace Identification method

In the literature, there are two methods for LPV identification: First, the ones based on global LPV estimation. Second, the ones based on the interpolation of local models [11], However, those approaches could lead to unstable representations of the LPV structure while the original system is stable [12]. That is why in this paper, we propose to use a subspace identification algorithm proposed (see [9] and [13]) to identify LPV systems which does not require interpolation or identification of local models and avoid instability problems.

Problem formulation

In the model used in identification in [9], the system matrices depend linearly on the time varying scheduling vector as follows:

x k+1 = m ∑ i=1 µ (i) k (A (i) x k + B (i) u k + K (i) e k )(1)y k = Cx k + Du k + e k(2)

with

x k ∈ R n , u k ∈ R r , y k ∈ R l areµ k = [ 1, µ(2)

k , ..., µ m k ] T Eqs.( 1) and ( 2) can be written in the predictor form:

x k+1 = m ∑ i=1 µ (i) k ( Ã(i) x k + B(i) u k + K (i) y k )(3)

with

Ã(i) = A (i) − K (i) C B(i) = B (i) − K (i) D

Assumptions and notationDefining z k = [ u T k , y T k

] T and using a data window of length p to define the following vector:

zp k =        z k z k+1 . . . z k+p−1       

and introducing the matrix obtained using the Kronecker product ⊗:

P p/k = µ k+p−1 ⊗ .... ⊗ µ k ⊗ I r+l we can define N p k =      p p/k . . . 0 . p p−1/k+1 . . . . 0 p 1/k+p−1     

Now, by defining the matrices U , Y and Z :

U = [u p+1 , ..., u N ](4)Y = [y p+1 , ..., y N ](5)Z = [N p 1 zp 1 , ..., N p N −p+1 zp N −p+1 ](6)

the controllability matrix can be expressed as:

κ p = [l p , ..., l 1 ] with l 1 = [ B(1) , ..., B(m) ] and l j = [ Ã(1) l j−1 , ..., Ã(m) l j−1 ]

If the matrix [ Z T , U T ] has full row rank, the matrix Cκ p and D can be estimated by solving the following linear regression problem [14]:

min Cκ p ,D ∥Y − Cκ p Z − DU ∥ 2 F (7)

where ∥∥ F represents the Frobenius norm. This problem can be solved by using traditional least square methods as in the case of LTI identification for time varying systems. Moreover, the observability matrix for the first model is calculated as follows:

Γ p =         C C Ã(1) . . . C( Ã(1) ) p−1         with κk p = [φ p−1,k+1 ⌣ B k , ..., φ 1,k+p−1 ⌣ B k+p−2 , ⌣ B k+p−1 ]

and 3) can be transformed into: 3) is uniformly exponentially stable the approximation error can be made arbitrarily small then:

⌣ B k = [ B, K k ] Then, Eq.(x k+p = φ p,k x k + κk p zp k x k+p = φ p,k x k + κ p N p k zp k where φ p,k = ÃK+p−1 ... Ãk+1 Ãk If the system (x k+p ≈ κ p N p k zp k

To calculate the observability matrix Γ p times the state X, we first calculate the matrix Γ p κ p :

Γ p κ p =       Cl p Cl p−1 . . Cl 1 0 CA (1) l p−1 . . C Ã(1) l 1 . . . . 0 C( Ã(1) ) p−1 l 1      

Then, using the following Singular Value Decomposition (SVD):

Γ p κ p Z = [ υ υ σ⊥ ] [ ∑ n 0 0 ∑ ] [ V V ⊥ ]

Proceedings of the 26 th International Workshop on Principles of Diagnosis the state is estimated by:

⌢ X = ∑ n V

Finally, C and D matrix are estimated using output equation (2) and A and B are estimated using the state equation (1). This algorithm can be summarized as follows [9]:

• Create the matrices U , Y and Z using (4),( 5) and ( 6),

• Solve the linear problems given in (7) ,

• Construct Γ p times the state X, • Estimate the state sequence,

• With the estimated state, use the linear relations to obtain the system matrices.

In the case of a very small p, we have in general a biased estimate. However, when the bias is too large, it will be a problem. That is why a large p would be chosen. In the case of a very large p, this method suffers from the curse of dimensionality [13] and the number of rows of Z grows exponentially with the size of the past window. In fact, the number of rows is given by:

ρ Z = (r + ℓ) ∑ p j=1 m j

To overcome this drawback, the kernel method will be introduced in the next subsection [15].

Kernel method

The equation ( 7) has a unique solution if the matrix [ Z T U T ] has full row rank and is given by:

[ Cκ p D ] = Y [ Z T U T ] ( [ Z U ] [ Z T U T ] ) −1

When this is not the case, that will occurs when p is large, the solution is computed by using the SVD of the matrix:

[ Z U ] = [ υ υ ⊥ ] [ ∑ m 0 0 0 ] [ V T V T ⊥ ]

Then, the solution of the minimum norm is given by:

[ Cκ p D ] = Y V ∑ −1 m υ T

To avoid computations in a large dimensional space, the minimum norm results in:

min α ∥α∥ 2 F (8) with Y − α [ Z T Z + U T U ]

= 0 where α are the Lagrange multipliers and

[ Z T Z + U T U ]

is referred as the kernel matrix.

The matrix Γ times the state X can be constructed as follows:

Γκ p Z =               α p ∑ j=1 (Z 1,j ) T Z 1,j α p ∑ j=2 (Z 2,j ) T Z 1,j . . . α p ∑ j=p (Z p,j ) T Z 1,j              (9)

with

(Z i,j ) T Z 1,j = ( p−j ∏ v=0 µ T Ñ +v+j−i µ Ñ +v+j−1 )(z T Ñ +j−i z Ñ +j−1 ) Z T Z = p ∑ j=1 (Z 1,j )Z 1,j(10)

Finally, the estimate sequence is obtained by solving the original SVD problem.

The kernel method can be summarized as follows [9]:

• Create the matrices U T U using (4) and Z T Z and (Z i,j ) T (Z i,j ) using ( 10), • Solve the linear problem given in (8),

• Construct Γ times the state X using ( 9)and ( 10),

• Estimate the state sequence,

• With the estimated state, use the linear relation to obtain the system matrices.

Interval predictor approach

To add robustness to the LPV subspace identification approach presented in the previous section, it will be combined with the interval predictor approach [16]. The interval predictor approach is an extension of classical system identification methods in order to provide the nominal model plus the uncertainty bounds for parameters guaranteeing that all collected data from the system in non-faulty scenarios will be included in the model prediction interval. This approach considers separately the additive and multiplicative uncertainties. Additive uncertainty is taken into account in the additive error term e(k) and modeling uncertainty is considered to be located in the parameters that are represented by a nominal value plus some uncertainty set around. In the literature, there are many approximation of the set uncertain parameter Θ. In our case, this set is described by a zonotope [10] :

Θ = θ 0 ⊕ HB n = {θ 0 + Hz : z ∈ B n }(11)

where: θ 0 is the nominal model (here obtained with the identification approach, H is matrix uncertainty shape, B n is a unitary box composed of n unitary (B = [−1, 1]) interval vectors and ⊕ denotes the Minkowski sum. A particular case of the parameter set is used that corresponds to the case where the parameter set Θ is bounded by an interval box [17]:

Θ = [θ 1 , θ 1 ] × ...[θ i , θ i ] × ...[θ n θ , θ n θ ](12)

where θ i = θ 0 i − λ i and θ i = θ 0 i + λ i with λ i ≥ 0 and i = 1, ..., n θ . In particular, the interval box can be viewed as a zonotope with center θ 0 and H equal to an n θ ×n θ diagonal matrix:

θ 0 = ( θ 1 + θ 1 2 , θ 2 + θ 2 2 , ..., θ n + θ n 2 )(13)H = diag(λ 1 , λ 2 , ..., λ n )(14)

For every output, a model can be extracted in the following regressor form:

y(k) = φ(k)θ(k) + e(k)(15)

where Proceedings of the 26 th International Workshop on Principles of Diagnosis

• φ(k) is the regressor vector of dimension 1× n θ which can contain any function of inputs u(k) and outputs y(k).

• θ(k) ∈ Θ is the parameter vector of dimension n θ ×1.

• Θ is the set that bounds parameter values.

• e(k) is the additive error bounded by a constant where |e(k)| ≤ σ. In the interval predictor approach, the set of uncertain parameters Θ should be obtained such that all measured data in fault-free scenario will be covered by the interval predicted output.

y(k) ∈ [ŷ(k) − σ, ŷ(k) + σ](16)

where

ŷ(k) = ŷ0 (k) − ∥φ(k)H∥ 1 (17) ŷ(k) = ŷ0 (k) + ∥φ(k)H∥ 1 (18)

and ŷ0 (k) is the model output prediction with nominal parameters with θ 0 =[θ 1 , θ 2 , ..., θ n θ ] T obtained using the LPV identification algorithm:

ŷ0 (k) = φ(k)θ 0 (k)(19)

Then, fault detection will be based on checking if ( 16) is satisfied. In case that, it is not satisfied a fault can be indicated. Otherwise, nothing can be said.

Case study: wind turbine benchmark system

In this work, a specific variable speed turbine is considered.

It is a three blade horizontal axis turbine with a full converter. The energy conversion from wind energy to mechanical energy can be controlled by changing the aerodynamics of the turbine by pitching the blades or by controlling the rotational speed of the turbine relative to the wind speed. The mechanical energy is converted to electrical energy by a generator fully coupled to a converter. Between the rotor and the generator, a drive train is used to increase the rotational speed from the rotor to the generator [18]. This model can be decomposed into submodels: Aerodynamic, Pitch, Drive train and Generator [19] [20]. In this paper, we focus on faults in the pitch subsystem as explained in the following subsection.

Pitch system model

In the wind turbine benchmark model, the hydraulic pitch is a piston servo mechanism which can be modeled by a second order transfer function [21] [1]:

β(s) β r (s) = ω 2 n s 2 + 2ζω n s + ω 2 n (20)

Notice that β r refers to reference values of pitch angles. The pitch model can be written in the following state space:

{ ẋ1 = x 2 ẋ2 = −2ξw n x 2 − w n 2 x 1 + w n 2 u(21)

with x 1 = β, x 2 = β, u = β r which can be discretised using an Euler approximation. Then, the following system is obtained:

{ x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k)(22)

with

A = [ 1 T e −T e w 2 n −2T e ξw n + 1 ] B = [ 0 T e w 2 n ] C = [ 1 0 ]

LPV Pitch system model

The pitch parameters w n and ξ are variable with hydraulic pressure P [1] [22]. Then, the pitch model can be written as the following LPV model according to [23] using P as the scheduling variable ϑ :

{ x(k + 1) = A(ϑ)x(k) + B(ϑ)u(k) y(k) = Cx(k)(23)

with

A(ϑ) = [ 1 T e −T e w 2 n (P ) −2T e ξ(P )w n (P ) + 1 ] B = [ 0 T e w 2 n (P ) ] y(k) = x 1 (k) = β(k)

Regressor form pitch system model

The pitch model can be transformed to the following regression form [24]:

y(k) = φ(k)θ(k)(24)

where, φ(k) is the regressor vector which can contain any function of inputs u(k) and outputs y(k). θ(k) ∈ Θ is the parameter vector. Θ is the set that bounds parameter values.

In particular

φ(k) = [ y(k − 2) y(k − 1) u(k − 2)] θ = [ θ 1 θ 2 θ 3 ] T θ 1 = (−T 2 e w 2 n + (2w n ξT e − 1)

)

θ 2 = −2w n ξT e + 2 θ 3 = T 2 e w 2 n

Results

The pitch systems, which in this case are hydraulic, could be affected by faults in any of the three blades. The considered faults in the hydraulic system can result in changed dynamics due to a drop in the main line pressure. This dynamic change induces a change in the system parameters: the damping ratio between 0.6 rad/s and 0.9 rad/s and the frequency between 3.42 rad/s and 11.11 rad/s according to [23]. In this work, a fault detection subspace estimator is designed to determine the presence of a fault. To distinguish between fault and modeling errors, an interval predictor approach is applied and a residual generation is used for

Θ = [θ 1 , θ 1 ] × [θ 2 , θ 2 ] × [θ 3 , θ 3 ](25)

and for i = 1, • • • , 3

λ i = ( θ i − θ i 2 )(26)θ 0 i = ( θ i + θ i 2 )(27)

using equations ( 17) and (18), the output bounds are calculated to be used in fault detection test which are given in Fig. 2. ŷ0 (k) is obtained by the use of the identification approach described in Section 2. To validate this algorithm two cases are used:

-Case 1: In this case, the pressure varies after time 10000s while parameters vary in the interval of parametric uncertainties, that is, damping ratio varies between 0.6 rad/s and 0.63 rad/s and the frequency between 10.34 rad/s and 11.11 rad/s. These parameters are presented respectively in Figures. 3 and 4. The pitch angle in this case is given in Fig. 5 altogether with the prediction intervals. For fault detection, the residual signal, based on the comparison between the measured pitch angle and the estimated one at each sampling instance, is calculated and it is shown in Fig. 6. For fault decision, a fault indicator signal is used and the decision is taken in function of this indicator. If the actual angle is not within the predicted interval given in Eq.( 16), the fault indicator is equal to 1 and the system is faulty. Otherwise, it is equal to 0 and the system is faultfree. The fault indicator signal given in Fig. 7 shows that there is no fault despite the pressure variation. The parameters variation is considered as a modeling error.

-Case 2: In this case, the pressure P varies between time t = 10000s and t = 17000s outside its nominal value. In this time interval, the damping ratio varies between 0.63 rad/s and 0.72 rad/s and the frequency varies between 8.03 rad/s and 10.34 rad/s. On the other hand, the damping ratio varies between 0.6 rad/s and 0.63 rad/s and the natural frequency varies between 10.34 rad/s and 11.11 rad/s outside as shown in Figures 8 and 9. In this case, the pitch angle is given in Fig. 10, while the residual and fault indicator signals are presented in Fig. 11 and Fig. 12, respectively. Fig. 12 shows that the fault indicator signal changes its signature between time 10000s and 17000s which induce that the parameters vary larger than the modeling range due to actuator fault in wind turbine benchmark system between instants t = 10000s and 17000s.

Conclusions

The proposed approach is based on an LPV estimation approach to generate a residual as the difference between the real and the nominal behavior of the monitored system. When a fault occurs, this residual goes out of the interval which represents the uncertainty bounds in non faulty case. These bounds are generated by means of an interval predictor approach that adds robustness to this fault detection method, by means of propagating the parameter uncertainty to the residual or predicted output. The proposed