LPV subspace identification for robust fault detection using a set-membership approach: Application to the wind turbine benchmark

Abstract

This paper focuses on robust fault detection for Linear Parameter Varying (LPV) systems using a set-membership approach. Since most of models which represent real systems are subject to modeling errors, standard fault detection (FD) LPV methods should be extended to be robust against model uncertainty. To solve this robust FD problem, a set-membership approach based on an interval predictor is used considering a bounded description of the modeling uncertainty. Satisfactory results of the proposed approach have been obtained using several fault scenarios in the pitch subsystem considered in the wind turbine benchmark introduced in IFAC SAFEPROCESS 2009.

-

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 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. Finally, Section 6 gives some concluding remarks. 2

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.

2.1 Problem formulation

In the model used in identification in [ 9 ], the system matrices depend linearly on the time varying scheduling vector as follows: xk+1 = m ∑ µ (ki)(A(i)xk + B(i)uk + K(i)ek)

yk = Cxk + Duk + ek with xk ∈ Rn, uk ∈ Rr, yk ∈ Rl are the state, input and output vectors and ek denotes the zero mean white innovation process and m is the number of local model or scheduling parameters:

µ k = [1, µ (k2) , ..., µ km]T Eqs.(1) and (2) can be written in the predictor form: xk+1 = m ∑ µ (ki)( A˜(i)xk + B˜(i)uk + K(i)yk) with 2.2

Assumptions and notation

Defining zk = [ukT , yT ]T and using a data window of k length p to define the following vector: (1) (2) (3) with and with and where

⌣ ⌣ ⌣ κ¯pk = [φ p−1,k+1Bk, ..., φ 1,k+p−1Bk+p−2, Bk+p−1] Then, Eq.(3) can be transformed into:

If the system (3) is uniformly exponentially stable the approximation error can be made arbitrarily small then: To calculate the observability matrix Γp times the state X, we first calculate the matrix Γpκp: xk+p ≈ κpN p z¯p

k k

Clp−1 . .

CA(1)lp−1 . .

.

.  Clp  0 Γpκp =  .

 .

0 Then, using the following Singular Value Decomposition (SVD): Γ\pκpZ = [ υ υσ⊥ ] [ ∑

n 0 0 ] [ V ] ∑ V⊥

. pp−1/k+1 . (5) (6) (7)

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

Cmκpin,D ∥Y − CκpZ − DU ∥2F 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:

Y = [yp+1 , ..., yN ] Z = [N1p z¯1p , ..., N Np−p+1 z¯N−p+1 p the controllability matrix can be expressed as: with ZT Z = p ∑(Z1,j )Z1,j j=1

p−j (Zi,j )T Z1,j = ( ∏ µ TN˜+v+j−iµ N˜+v+j−1)(z NT˜+j−iz N˜+j−1) v=0 the state is estimated by: ⌢ X = ∑V n

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 mj j=1

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

Kernel method

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

mαin ∥α∥2F (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 ZT Z and (Zi,j )T (Zi,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.

3 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 ⊕ HBn = { θ0 + Hz : z ∈ Bn} where: θ0 is the nominal model (here obtained with the identification approach, H is matrix uncertainty shape, Bn 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θ ] where θi = θi0 − λi and θi = θi0 + λ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 , ..., θn + θn 2 2 2 )

H = diag(λ1, λ2, ..., λn) For every output, a model can be extracted in the following regressor form: y(k) = φ (k)θ(k) + e(k) (11) (12) (13) (14) (15) • φ (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. where y(k) ∈ [yˆ(k) − σ, yˆ(k) + σ] yˆ(k) = yˆ0(k) − ∥φ (k)H∥1 yˆ(k) = yˆ0(k) + ∥φ (k)H∥1 and yˆ0(k) is the model output prediction with nominal parameters with θ0 =[θ1, θ2, ..., θnθ ]T obtained using the LPV identification algorithm:

yˆ0(k) = φ (k)θ0(k)

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. 4

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. 4.1

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) =

ωn2 s2 + 2ζωns + ωn2

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

x˙1 = x2 x˙2 = −2ξwnx2 − wn2x1 + wn2u (16) (17) (18) (19) (20) (21) with x1 = β, x2 = β˙, u = βr which can be discretised using an Euler approximation. Then, the following system is obtained: with A = B = The pitch parameters wn 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)

Te −2Teξ(P )wn(P ) + 1 ] y(k) = x1(k) = β(k) 4.3

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 = (−Te2wn2 + (2wnξTe − 1)) θ2 = −2wnξTe + 2 θ3 = Te2wn2 5

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 Upper Lower 16 14 12 10 t u tup 8 o x a im 6 n m 4 2 0 deciding if there is a fault. To illustrate the performance of this robust fault detection approach:: ξ ∈[ 0.6 0.63 ] and wn ∈[ 10.34 11.11 ] are considered. Then, a parameter set Θ is bounded by an interval box:

Θ = [θ1, θ1] × [θ2, θ2] × [θ3, θ3] and for i = 1, · · · , 3 λi = ( θi − θi )

2 θi0 = ( θi + θi )

2 using equations (17) and (18), the output bounds are calculated to be used in fault detection test which are given in Fig. 2. yˆ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 (25) (26) (27) 0.635 0.63 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. 6

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 2.2 2.1 2 −0.10

1 0.95 r 0.9 o t a c i d n i ltu0.85 a f 0.8 0.75 Mesaured Max Min 1.6608 1.6608 1.6608 approach is illustrated by implementing a robust fault detection scheme for a pitch subsystem of the wind turbine benchmark. Simulations show satisfactory fault detection performance despite model uncertainties.

References

[1] P. Odgaard , J. Stoustrup , and M. Kinnaert . Fault tolerant control of wind turbines-a benchmark model. In 7th IFAC symposium on fault detection, supervision and safety of technical processes , Barcelona, spain, 2009 .

R. Isermann . Fault diagnosis systems: an introduction from fault detection to fault tolerance . 2006 .

G. Mercere , M. Lovera , and E Laroche . Identification of a flexible robot manipulator using a linear parameter-varying descriptor state-space structure . In Proc. of the IEEE conference on decision and control , Orlando, Florida,USA, 2011 .

[5] J. Dong , B. Kulcsár , and M Verhaegen . Fault detection and estimation based on closed-loop subspace identification for linear parameter varying systems . In DX, Stockholm, 2009 .

[6] J. Bravo , T. Alamo , and E.F. Camacho . Bounded error identification of systems with time-varying parameters . IEEE Transactions on Automatic Control , 51 :1144 U˝ 1150 ., 2006 .

D. Efimov , L. Fridman , T. Raissi , A. Zolghadri , and R. Seydou . Interval estimation for lpv systems applying high order sliding mode techniques . Automatica , 48 : 2365 - 2371 , 2012 .

D. Efimov , T. Raissi , and A. Zolghadri . Control of nonlinear and lpv systems: interval observer-based framework . IEEE Transactions on Automatic Control. , 2013 .

[9] J. Van Willem and M Verhagen . Subspace identification of bilinear and lpv systems for open-and closedloop data . Automatica , 45 : 371 - 381 , 2009 .

[10] J. Blesa , V. Puig , J Romera , and J Saludes . Fault diagnosis of wind turbines using a set-membership approach . In the 18th IFAC world congress, Milano , Italy, 2011 .

[11] H. Tanaka , Y Ohta , and Y Okimura . A local approach to lpv-identification of a twin rotor mimo system . In in proceedings of the 47th IEEE Conference on Decision and Control Cancun, Mexico , 2008 .

[12] R. Toth , F. Felici , P. Heuberger , and P Van den Hof. Discrete time lpv i/o and state-space representations, differences of behavior and pitfalls of interpolation . In in proceedings of the European Control Conference (ECC) , Kos, Greece, 2007 .

[13] J. Van Willem and M Verhagen . Subspace identification of mimo lpv systems: the pbsid approach . In in Proceedings of the 47th IEEE Conference on Decision and Control Cancun, Mexico , 2008 .

[14] P. Gebraad , J. Van Wingerden , G. Van der Veen , and M Verhaegen . Lpv subspace identification using a novel nuclear norm regularization method . In American Control Conference on O'Farrell Street , San Francisco, CA, USA, 2011 .

[15] V. Verdult and M Verhaegen . Kernel methods for subspace identification of multivariable lpv and bilinear systems . Automatica , 41 : 1557 - 1565 , 2005 .

[16] J. Blesa , V. Puig , and J Saludes . Identification for passive robust fault detection using zonotope based set membership appraches . International journal of adaptive control and signal processing , 25 : 788 - 812 , 2011 .

[17] P. Puig , V. Quevedo , T. Escobet , F. Nejjari , and S De las Heras. Passive robust fault detection of dynamic processes using interval models . IEEE Transactions on Control Systems Technology , 16 : 1083 - 1089 , 2008 .

[18] B. Boussaid , C. Aubrun , and M.N Abdelkrim . Setpoint reconfiguration approach for the ftc of wind turbines . In the 18th World Congress of the International Federation of Automatic Control (IFAC) , Milano, Italy, 2011 .

[19] B. Boussaid , C. Aubrun , and M.N Abdelkrim . Twolevel active fault tolerant control approach . In The Eighth International Multi-Conference on Systems, Signals Devices (SSD'11) ,Sousse, Tunisia , 2011 .

[20] B. Boussaid , C. Aubrun , and M.N Abdelkrim . Active fault tolerant approach for wind turbines . In The International Conference on Communications, Computing and Control Applications (CCCA'11) , Hammamet, Tunisia, 2011 .

[21] P. Odgaard , J. Stoustrup , and M Kinnaert . Fault tolerant control of wind turbines-a benchmark model . IEEE Transactions on control systems Technology , 21 : 1168 - 1182 , 2013 .

[22] P. Odgaard and J Stoustrup . Results of a wind turbine fdi competition . In 8th IFAC symposium on fault detection ,supervision and safety of technical processes,Mexico , 2012 .

[23] C. Sloth , T. Esbensen , and J Stoustrup . Robust and fault tolerant linear parameter varying control of wind turbines . Mechatronics , 21 : 645 - 659 , 2011 .

[24] H. Chouiref , B. Boussaid , M.N Abdelkrim , V. Puig , and C Aubrun . Lpv model-based fault detection: Application to wind turbine benchmark . In International conference on electrical sciences and technologies (cistem'14) , Tunis , 2014 .