Residual Generator Fuzzy Identification for Wind TurbineBenchmark Fault Diagnosis

Abstract

In order to improve the availability of wind turbines, thus improving their efficiency, it is important to detect and isolate faults in their earlier occurrence. The main problem of model-based fault diagnosis applied to wind turbines is represented by the system complexity, as well as the reliability of the available measurements. In this work, a data-driven strategy relying on fuzzy models is presented, in order to build a fault diagnosis system. Fuzzy theory jointly with the Frisch identification scheme for errors-in-variable models is exploited here, since it allows one to approximate unknown models and manage uncertain data. Moreover, the use of fuzzy models, which are directly identified from the wind turbine measurements, allows the design of the fault detection and isolation module. It is worth noting that, sometimes, the nonlinearity of a wind turbine system could lead to quite complex analytic solutions. However, IF-THEN fuzzy rules provide a simpler solution, important when on-line implementations have to be considered. The wind turbine benchmark is used to validate the achieved performances of the suggested fault detection and isolation scheme. Finally, comparisons of the proposed methodology with respect to different fault diagnosis methods serve to highlight the features of the suggested solution.

1. Introduction

Wind turbines of the megawatt size are expensive, and hence, their overall availability must be high to optimize the energy generation, thus reducing the cost of energy. In the same way, wind turbine downtime must be minimized. This key feature can be achieved by introducing fault detection and isolation (FDI) systems. In the related state-of-the-art, the fault detection schemes for industrial wind turbines can be quite conservative. For example, turbines are simply turned off after minimal faults to wait for service. Consequently, there is a need for an effective FDI for improving wind turbine working conditions, even though it might lead to limited power production in the case of faults. In the last few years, some works have been proposed on wind turbine FDI.

As an example [1], presented a Kalman filter-based diagnosis system for detecting faults in the blade root bending moment sensors. An unknown input observer was designed for the detection of sensor faults around the wind turbine drive train in [2]. On the other hand, in [3], active and passive fault tolerant control schemes are considered. More attention has been drawn on wind turbine electrical conversion systems. Some relevant examples can be found in [4], where observer-based solutions for sensor fault detection are presented. In [5], a fault detection and reconfiguration scheme for a doubly-fed wind turbine converter is shown.

Comparisons of the considered fault diagnosis schemes for wind turbine applications are beneficial to find the most effective methodology. In [6], a wind turbine benchmark was presented, which was proposed for both FDI and FTCsolution comparisons. This benchmark model described a realistic three-blade horizontal variable speed wind turbine with a full-scale converter coupling, and it will be considered in the present study. Several FDI strategies were proposed in the recent literature for this specific benchmark [6]. In particular, the most effective FDI approaches were based on support vector machines with a Gaussian kernel [7], banks of estimators [8], up-down counters for a redundant residual decision [9], combined observer/Kalman filters [10] and automatically generated fault models [11].

This paper proposes a data-driven FDI approach, based on identified fuzzy residual generators, which is applied to the benchmark proposed in [6]. To this aim, the key contributions of the presented study are remarked in as follows. First, the system complexity may not require the design of a sophisticated analytical model of the residual generators. In fact, as shown in this work, a method relying on fuzzy models is proposed, thus obviating the derivation of purely nonlinear mathematical models of the residual functions. Note that the advantages of model-based approaches with respect to data-driven solutions depend on the features of the model under diagnosis, as described, e.g., in [12,13]. For the first time, the suggested methodology is applied to the wind turbine benchmark. Secondly, residual generators in the form of fuzzy models are considered instead of purely nonlinear observer or filter functions. Again, for the first time, the authors have proposed here the design of these residual generators with application to the wind turbine benchmark. Third, by exploiting the failure mode and effect analysis described in the following, both the fuzzy identification and the fault isolation tasks are enhanced, without using complicated unknown input or disturbance decoupling approaches, as addressed, e.g., in [12,13].

This paper suggests the use of fuzzy logic, since it seems the natural tool for handling complicated and uncertain conditions [14,15] that, in general, are not available to the designer. The additional benefits of fuzzy logic in connection with errors-in-variable (EIV) identification techniques include its simplicity and its flexibility [16,17]. Fuzzy logic can model nonlinear functions with arbitrary complexity. Fuzzy models, called fuzzy inference systems (FIS), consist of a number of conditional “IF-THEN” rules. For the designer, these rules are easy to build, and as many rules as necessary can be supplied to describe the system adequately with arbitrary accuracy (although typically, only a moderate number of rules are needed). In fuzzy logic, unlike standard conditional logic, the truth of any statement is a matter of degree. FIS relies on membership functions to explain how to calculate the correct value between zero and one. The degree to which any fuzzy statement is true is denoted by a value between zero and one. Not only do the rule-based approach and flexible membership function scheme make fuzzy systems straightforward to create, but they also simplify the design of systems and ensure that you can easily update and maintain the system over time [15].

As described in the following, the paper suggests the use of the Takagi–Sugeno (TS) model [18], whose parameters are obtained via the identification procedure proposed in [16,17]. This procedure is based on the EIV identification strategy relying on the Frisch scheme, which assumes that the input and output data are affected by noise, thus providing an unbiased parameter estimation. The diagnosis approach being evaluated has also an important implication on the use of on-line diagnosis tools once the wind turbine is under customer operation. Therefore, the further contribution of the study regards the robustness and reliability analysis of the proposed FDI scheme. The effectiveness of the proposed strategy is verified using data sequences acquired from the wind turbine benchmark. Realistic conditions and comparisons with different fault diagnosis approaches have been considered to validate the proposed methodology. Note that the wind turbine simulator exploited in this study was modified by the authors, and the working conditions can be different from the ones shown in [6], since it was oriented toward the assessment of the robustness and the reliability of the considered FDI schemes. Moreover, with respect to previous works by the same authors, this paper aims at describing the comprehensive data-driven approach to wind turbine FDI and its extended comparisons with different FDI strategies applied to the same case study.

Finally, the paper has the following structure. Section 2 briefly recalls the wind turbine benchmark. Section 3 addresses the strategy using the EIV identification approach exploited for obtaining the fuzzy models, which are used as residual generators for the design of the FDI strategy. The proposed FDI methodology is presented in Section 4. The results achieved, which are summarized in Section 5, show the performances of the fault diagnosis scheme, validated on the data directly acquired from the benchmark, and compared also with different methods. Section 6 ends the paper by highlighting the main achievements of the work, providing suggestions for further studies.

2. Wind Turbine Benchmark Description

The three-blade horizontal axis turbine considered in this paper works according to the principle that the wind is acting on the blades and moves the rotor shaft. In order to up-scale the rotational speed to the generator, a gear box is introduced. A more accurate description of the benchmark model can be found in [6].

The rotational speed, and consequently, the generated power, is regulated by means of two controlled inputs: the converter torque ${\tau }_g\left(t\right)$ and the pitch angle ${\beta }_r\left(t\right)$ of the turbine blades. From the wind turbine system, a number of measurements can be acquired. ${\omega }_r\left(t\right)$ is the rotor speed, ${\omega }_g\left(t\right)$ the generator speed and ${\tau }_g\left(t\right)$ the torque of the generator controlled by the converter, which is provided with the torque reference, ${\tau }_r\left(t\right)$. The estimated aerodynamic torque is defined as ${\tau }_{aero}\left(t\right)$. This estimate clearly depends on the wind speed, which is not a very accurate measurement.

The wind turbine model in the continuous-time domain is briefly recalled in this section. In particular, the aerodynamic model is defined as in Equation (1):

$${\tau }_{aero}\left(t\right)=\frac{\rho A{C}_p\left({\beta }_r\left(t\right),\lambda \left(t\right)\right)v^3\left(t\right)}{2{\omega }_r\left(t\right)}$$

(1)

where ρ is the density of the air, A is the area covered by the turbine blades in its rotation, $v\left(t\right)$ is the wind speed, whilst $\lambda \left(t\right)$ is the tip-speed ratio of the blade. $C_p$ represents the power coefficient, here described by means of a two-dimensional map (look-up table). Equation (1) is used to compute the aerodynamic torque ${\tau }_{aero}\left(t\right)$ that depends on $v\left(t\right)$, the measured pitch angle ${\beta }_r\left(t\right)$ and rotor speed ${\omega }_r\left(t\right)$. Due to the uncertainty of the wind speed, the estimate of ${\tau }_{aero}\left(t\right)$ is considered affected by an unknown measurement error, which motivates the approach described in Section 3. Moreover, also the nonlinearity represented by Equation (1) is thus taken into account.

A first order model is used to represent the wind turbine rotor and generator dynamics [6]. Thus, the generator torque ${\tau }_g\left(t\right)$ and the reference ${\tau }_r\left(t\right)$ are transformed to the low speed side of the drive-train (rotor side), whilst $p_{gen}$ is the generator power coefficient. The hydraulic pitch model is described as a closed-loop transfer function of the hydraulic pitch system, whilst the drive-train is modeled using the two-mass description [6]. Moreover, the converter dynamics are modeled by a first-order transfer function, and $P_g\left(t\right)$ represents the power produced by the generator, depending on its efficiency. The measurement sensors are modeled by adding the actual variable values with stochastic noise processes. These noise signals are described as Gaussian processes with fixed mean and standard deviation values, depending on the considered measurement sensors [6].

With these assumptions, the complete continuous-time description of the wind turbine benchmark has the form of Equation (2):

$$\left\{\begin{array}{ccc}{\dot{x}}_c\left(t\right)\hfill & =& f_c\left(x_c\left(t\right),u\left(t\right)\right)\hfill \\ & & \\ y\left(t\right)\hfill & =& x_c\left(t\right)\hfill \end{array}\right.$$

(2)

where $u\left(t\right)={[{\beta }_{1m_i}\left(t\right),{\beta }_{2m_i}\left(t\right),{\beta }_{3m_i}\left(t\right),{\tau }_g\left(t\right)]}^T$ and $y\left(t\right)=x_c\left(t\right)=[P_g\left(t\right)$, ${\omega }_{g{m}_i}\left(t\right),{\omega }_{r{m}_i}{\left(t\right)]}^T$ are the control inputs and the monitored output measurements, respectively, measured by the i-th redundant sensor, with $i=1,2$. $f_c\left(·\right)$ represents the continuous-time nonlinear function describing the behavior of the wind turbine benchmark. These measurements will be sampled for obtaining N input-output data, $u\left(k\right)$ and $y\left(k\right)$, with $k=1,2,\dots N$. Regarding the input and output signals, ${\omega }_{g{m}_i}$ is the i-th generator speed measurement, ${\omega }_{g{m}_i}$ the i-th rotor speed measurement, $P_g\left(t\right)$ the generator power measurement and ${\beta }_{j{m}_i}\left(t\right)$ the i-th pitch measurement of the j-th blade. Finally, the model parameters and the map $C_p\left(\beta ,\lambda \right)$ are chosen to represent a realistic wind turbine installation [6].

2.1. Simulated Fault Conditions

The benchmark model implements a number of realistic faults, which are summarized in Table 1.

Table 1. Benchmark fault cases.

Fault	Description
1	Fixed value on Pitch 1, Position Sensor 1
2	Scaling error on Pitch 2, Position Sensor 2
3	Fixed value on Pitch 3, Position Sensor 1
4	Fixed value on Rotor Speed Sensor 1
5	Scaling error on Rotor Speed Sensor 2 and Generator Speed Sensor 2
6	Changed pitch system response, Pitch Actuator 2: high air content in oil
7	Changed pitch system response, Pitch Actuator 3: low pressure
8	Offset in converter torque control
9	Changed dynamics drive train

Table 1. Benchmark fault cases.

Fault	Description
1	Fixed value on Pitch 1, Position Sensor 1
2	Scaling error on Pitch 2, Position Sensor 2
3	Fixed value on Pitch 3, Position Sensor 1
4	Fixed value on Rotor Speed Sensor 1
5	Scaling error on Rotor Speed Sensor 2 and Generator Speed Sensor 2
6	Changed pitch system response, Pitch Actuator 2: high air content in oil
7	Changed pitch system response, Pitch Actuator 3: low pressure
8	Offset in converter torque control
9	Changed dynamics drive train

A more detailed description of these faults, which is beyond the scope of this paper, can be found in [6].

The remainder of this section describes the relations among the fault cases described above and the monitored measurements acquired from the wind turbine process. In this way, it will be shown that both the system identification and the fault isolation tasks can be easily solved. In fact, Table 2 highlights how a single fault affects the measured inputs $u\left(k\right)$ and outputs $y\left(k\right)$. Moreover, the mismatch between each fault-free and faulty measurement is measured by the relative mean squared error (RMSE), computed for the different fault cases of Table 1.

Table 2. Wind turbine failure mode and effect analysis (FMEA) results.

Measurement	${\beta }_{1m_1}\left(t\right)$	${\beta }_{2m_2}\left(t\right)$	${\beta }_{3m_1}\left(t\right)$	${\omega }_{r{m}_1}$	${\omega }_{r{m}_1}$	${\beta }_{2m_1}\left(t\right)$	${\beta }_{3m_2}\left(t\right)$	${\tau }_{gm}$	${\omega }_{g{m}_1}$
Fault	1	2	3	4	5	6	7	8	9
RMSE	11.29	0.98	2.48	1.44	1.45	0.80	0.73	0.84	0.77

Table 2. Wind turbine failure mode and effect analysis (FMEA) results.

Measurement	${\beta }_{1m_1}\left(t\right)$	${\beta }_{2m_2}\left(t\right)$	${\beta }_{3m_1}\left(t\right)$	${\omega }_{r{m}_1}$	${\omega }_{r{m}_1}$	${\beta }_{2m_1}\left(t\right)$	${\beta }_{3m_2}\left(t\right)$	${\tau }_{gm}$	${\omega }_{g{m}_1}$
Fault	1	2	3	4	5	6	7	8	9
RMSE	11.29	0.98	2.48	1.44	1.45	0.80	0.73	0.84	0.77

Note that in Table 2, the variable ${\beta }_{i{m}_j}\left(t\right)$ indicates the i-th blade pitch ($i=1,2,3$) measured by the j-th redundant sensor $j=1,2$. In the same way, the rotor speed is measured by two redundant sensors ${\omega }_{r{m}_j}\left(t\right)$, with $j=1,2$. On the other hand, only one sensor provides the generator torque measurement ${\tau }_{gm}$.

Table 2 was obtained by performing a fault sensitivity analysis, in particular the so-called failure mode and effect analysis (FMEA), as explained in [19]. Table 2 is thus obtained by selecting the most sensitive measurement with respect to the simulated fault conditions. In practice, the monitored fault signals have been injected in the wind turbine simulator described in Section 2. After each fault signal has been added to the corresponding measurement, the relative absolute errors between the fault-free and faulty measured signals have been computed. The measured signal showing the biggest error corresponds also to the most sensitive measurement to the considered fault.

Finally, the analysis summarized in Table 2 enhances the design of the bank of fuzzy estimators that are used for the fault isolation task, as described in Section 4. Moreover, it was assumed that only a single fault may occur in the considered plant. However, on the basis of the fault effect analysis, faults occurring at the same time can be distinguished by analyzing their effects on the monitored measurements.

3. Fuzzy Modeling and Identification

This section describes the design of the estimation and the diagnosis schemes applied to the wind turbine simulator. In particular, the identification method, which is recalled in Section 3.1, leads to the proposed FDI strategy of Section 4.

3.1. Fuzzy Modeling from Data Clustering

In the considered TS model, the rule consequents are crisp functions of the inputs:

$$R_i:\mathrm{IF}x\mathrm{is}{A}_i\mathrm{THEN}{y}_i=f_i\left(x\right)$$

(3)

with $i=1,2,\dots ,K$. x represents the input (antecedent) variable, and $y_i$ is the output (consequent) variable. $R_i$ denotes the i-th rule, and K is the number of rules (or clusters) in the rule base. $A_i$ is the antecedent fuzzy set of the i-th rule, defined by a (multivariate) membership function [20,21].

The consequent $f_i$ are chosen as a suitable parameterized function, whose structure remains equal in all of the rules and only the parameters vary. The parametrization exploited here is the affine form:

$$y_i=a_i^{T}x+b_i,$$

(4)

where $a_i$ is a parameter vector and $b_i$ is a scalar offset. This model is referred to as the affine TS model. The antecedent of each rule defines a (fuzzy) validity region for the corresponding affine consequent model. These TS models are considered here due to their interesting approximation properties [22].

It is worth noting that the data selection for the fuzzy modeling represents an important topic, since representative data should be selected. Under this consideration, a simple and effective method for selecting the most significant input and output variables for building a fuzzy model is proposed and analyzed in Section 2. As already remarked, the FMEA procedure applied to the wind turbine is able to enhance both the system identification and the fault isolation tasks. On the other hand, the optimal number of fuzzy rules K is determined by a proper integration of existing clustering-based methods with the identification approach proposed in [17]. This approach exploits the fuzzy clustering and validates the data partition on the basis of sub-clusters created according to the system identification scheme recalled in Section 3.2. In this way, the most important input variables that independently and significantly influence the model output can be validated. Moreover, the optimal number of fuzzy rules K can be determined separately via the identification approach modified by the authors and proposed for nonlinear system modeling [17]. The simulation results of Section 5 will show that the proposed method can provide also good model structures for the fuzzy modeling oriented toward wind turbine fault diagnosis.

Before the output can be inferred, the degree of fulfillment of the antecedent denoted by ${\lambda }_i\left(x\right)$ must be computed. For rules with multivariate antecedent fuzzy sets given by Equation (3) the degree of fulfillment is simply equal to the membership degree of the given input x, i.e., ${\lambda }_i={\mu }_{A_i}\left(x\right)$. When logical connectives are used, the degree of fulfillment of the antecedent is computed as a combination of the membership degrees of the individual propositions using fuzzy logic operators.

In the TS model, the inference is reduced to a simple algebraic expression, via, e.g., the fuzzy-mean defuzzification formula:

$$y=\frac{{\sum }_{i=1}^K{\lambda }_i\left(x\right)y_i}{{\sum }_{i=1}^K{\lambda }_i\left(x\right)}$$

(5)

where the membership degrees ${\lambda }_i$ are modeled as exponential functions.

In order to introduce dynamics into the model of Equation (3), the consequents are linear ARXmodels, where n is the order of the dynamic system, $x\left(k\right)=\left[y\right(k-1),\cdots ,y(k-n),$ ${u(k-1),\cdots ,u(k-n)]}^T$ and $a_i=[{\alpha }_1^{\left(i\right)},\cdots ,$${\alpha }_n^{\left(i\right)},{\delta }_1^{\left(i\right)},\cdots ,{\delta }_n^{\left(i\right)}]$.

Finally, with reference to the structure of Equation (5), a procedure proposed in [16,17] for the estimation of both parameters $a_i$, $b_i$ and the order n will be summarized in the following. Moreover, the estimation of the membership degrees ${\lambda }_i$ of Equation (5) used for the aggregation of the local affine submodels of Equation (4) has been obtained via a data clustering method available in the literature [15].

Many clustering algorithms have been proposed; see, e.g., [20,21]. In particular, the clustering algorithm exploited in this work is based on optimization of the basic c-means objective function and known as the fuzzy c-means clustering algorithm [23]. Note that the clustering method is not the key issue of the study, since the identification scheme described in [16] requires only that the clustered data are compatible with linear prototypes. However, the clustering algorithm recalled above was proposed here, since it can be easily integrated with the identification approach developed in [16] and already provided as a ready-to-use software program [15].

Finally, it is worth observing that, when considering fuzzy identification via data clustering, an important point concerns the determination of the optimal number of clusters. When clustering real data without any a priori information about the data structure, one usually has to make assumptions about the number of underlying subgroups (clusters) K in the data. The chosen clustering algorithm then searches for K clusters, regardless of whether they are really present in the data or not. The optimization issue was investigated in [16,17].

3.2. Affine Model Identification from Data

This section recalls the method for the construction of the TS models and the procedure for estimating the consequent parameters.

Regarding the antecedent membership functions $A_i$, they are obtained by computing the membership degrees directly in the product space of the antecedent variables. As shown in Section 5, exponential membership functions proved to be suitable for the accurate representation of the cluster shapes. These functions are derived by considering multi-dimensional antecedent membership relations, which are represented analytically by computing an inverse of the distance from the cluster prototypes. The membership degree is computed directly for the entire input vector (without the decomposition). The antecedents of the TS rules in the expression of Equation (3) are expressed as ${\lambda }_i\left(x\right)={\mu }_{A_i}\left(x\right)$.

On the other hand, in general, there are several methods for estimating the consequent parameters $a_i$ and $b_i$ in Equation (4). By assuming that errors are present in both the regressor and the regressand, the set of optimal parameters are estimated by using the procedure recalled in the following. This approach can be seen as the minimization of the prediction errors of the individual local models, solved as a set of K independent problems, as proposed in [16,17]. This scheme, which is usually preferred when the TS model should serve as a predictor, computes the consequent parameters by the so-called Frisch scheme [17].

In order to enhance the readability of the remainder of this section, the meaning of the symbols and variables is briefly summarized in Table 3.

Table 3. Nomenclature.

Variable	Description
$X_n^{\left(i\right)}$	Matrix of the data from the i-th cluster and n regressors
$x_n\left(h\right)$	Vector of the input-output data sequences generated from n regressors
$\tilde{u}\left(k\right)$, $\tilde{y}\left(k\right)$	Additive input-output noise
$u^{*}\left(k\right)$, $y^{*}\left(k\right)$	Noise-free input-output data (unmeasurable)
$u\left(k\right)$, $y\left(k\right)$	Measured input-output data sequences
${\Sigma }_n^{\left(i\right)}$	Sample covariance matrix from the data in the matrix $X_n^{\left(i\right)}$
${\overline{\tilde{\Sigma }}}_n$	Sample covariance matrix of the actual noise affecting the data
${\Sigma }_n^{*\left(i\right)}$	Sample covariance matrix of the (unmeasurable) noise-free data in the i-th cluster
${\tilde{\Sigma }}_n$	Unknown noise sample covariance matrix (estimated by the proposed procedure)
${\overline{\tilde{\sigma }}}_u$, ${\overline{\tilde{\sigma }}}_y$	Actual noise variances
${\tilde{\sigma }}_u$, ${\tilde{\sigma }}_y$	Unknown noise variances (to be identified)
${\Gamma }_n^{\left(i\right)}=0$	2-dimensional function in the unknowns $({\tilde{\sigma }}_u,{\tilde{\sigma }}_y)$
${\overline{\tilde{\sigma }}}_u^{\left(i\right)}$, ${\overline{\tilde{\sigma }}}_y^{\left(i\right)}$	Identified noise values in the i-th cluster

Table 3. Nomenclature.

Variable	Description
$X_n^{\left(i\right)}$	Matrix of the data from the i-th cluster and n regressors
$x_n\left(h\right)$	Vector of the input-output data sequences generated from n regressors
$\tilde{u}\left(k\right)$, $\tilde{y}\left(k\right)$	Additive input-output noise
$u^{*}\left(k\right)$, $y^{*}\left(k\right)$	Noise-free input-output data (unmeasurable)
$u\left(k\right)$, $y\left(k\right)$	Measured input-output data sequences
${\Sigma }_n^{\left(i\right)}$	Sample covariance matrix from the data in the matrix $X_n^{\left(i\right)}$
${\overline{\tilde{\Sigma }}}_n$	Sample covariance matrix of the actual noise affecting the data
${\Sigma }_n^{*\left(i\right)}$	Sample covariance matrix of the (unmeasurable) noise-free data in the i-th cluster
${\tilde{\Sigma }}_n$	Unknown noise sample covariance matrix (estimated by the proposed procedure)
${\overline{\tilde{\sigma }}}_u$, ${\overline{\tilde{\sigma }}}_y$	Actual noise variances
${\tilde{\sigma }}_u$, ${\tilde{\sigma }}_y$	Unknown noise variances (to be identified)
${\Gamma }_n^{\left(i\right)}=0$	2-dimensional function in the unknowns $({\tilde{\sigma }}_u,{\tilde{\sigma }}_y)$
${\overline{\tilde{\sigma }}}_u^{\left(i\right)}$, ${\overline{\tilde{\sigma }}}_y^{\left(i\right)}$	Identified noise values in the i-th cluster

Thus, in order to identify the structure of the TS model of Equation (5) in the i-th cluster with $i=1,\cdots ,K$ and K clusters, the following matrices are defined:

$$X_n^{\left(i\right)}=\left[\begin{array}{ccc}y\left(k\right)& x_n^T\left(k\right)& 1\\ y(k+1)& x_n^T(k+1)& 1\\ ⋮& ⋮\\ y(k+N_i-1)& x_n^T(k+N_i-1)& 1\end{array}\right]$$

(6)

where the subscript n represents the order of the considered local ARX dynamic model (number of regressors), i.e., $x_n\left(h\right)=[y(h-1),\cdots ,y(h-n),$${u(h-1),\cdots ,u(h-n)]}^T$. Therefore:

$${\Sigma }_n^{\left(i\right)}={\left(X_n^{\left(i\right)}\right)}^T{X}_n^{\left(i\right)}.$$

(7)

In order to solve the so-called noise-rejection problem [16,17] in a mathematical framework, it is necessary to follow the assumptions that the noise variables $\tilde{u}\left(k\right)$ and $\tilde{y}\left(k\right)$ are additive on the input-output data $u^{*}\left(k\right)$ and $y^{*}\left(k\right)$ and region independent ($k=1,2,\cdots ,N$). Under this hypothesis, a positive-definite matrix ${\Sigma }_n^{\left(i\right)}$ associated with the sequences belonging to the i-th cluster is expressed as the sum of two terms, i.e., ${\Sigma }_n^{\left(i\right)}={\Sigma }_n^{*\left(i\right)}+{\overline{\tilde{\Sigma }}}_n$, where:

$${\overline{\tilde{\Sigma }}}_n=\mathrm{diag}[{\overline{\tilde{\sigma }}}_y{I}_{n+1},{\overline{\tilde{\sigma }}}_u{I}_n,0]\ge 0.$$

(8)

The solution of the above identification problem requires the computation of the unknown noise variances ${\overline{\tilde{\sigma }}}_u$ and ${\overline{\tilde{\sigma }}}_y$, which can be obtained by solving Equation (9):

$${\Sigma }_n^{*\left(i\right)}={\Sigma }_n^{\left(i\right)}-{\tilde{\Sigma }}_n\ge 0.$$

(9)

in the variables ${\tilde{\sigma }}_u,{\tilde{\sigma }}_y$, where ${\tilde{\Sigma }}_n=\mathrm{diag}[{\tilde{\sigma }}_y{I}_{n+1},{\tilde{\sigma }}_u{I}_n,0]$. It is worth noting that all surfaces ${\Gamma }_n^{\left(i\right)}=0$ determined by the locus of the points satisfying Equation (9) have necessarily at least one common point, i.e., the point $({\overline{\tilde{\sigma }}}_u,{\overline{\tilde{\sigma }}}_y)$ corresponding to the true variances of the noise affecting the input and the output data [16,17].

The search for a solution for the Frisch scheme identification problem can therefore start from the determination of this point in the noise space, if the noise characteristics are common to all of the clusters and all assumptions regarding the Frisch scheme are satisfied [16,17].

However, in real cases, these assumptions have to be relaxed, as no common point can be determined among surfaces ${\Gamma }_n^{\left(i\right)}=0$ in the noise plane. A unique solution to the identification problem cannot be obtained. In this situation, the local model identification is performed by finding the point $({\tilde{\sigma }}_u,{\tilde{\sigma }}_y)\in {\Gamma }_{n+1}^{\left(i\right)}=0$ that makes ${\Sigma }_{n+1}^{*\left(i\right)}$ closer to the double-singular condition. This leads to determining the nearest point of all of the surfaces, even when the assumptions of the Frisch scheme are violated. Moreover, for each i-th cluster, different noise values $({\overline{\tilde{\sigma }}}_u^{(i)}$, ${\overline{\tilde{\sigma }}}_y^{\left(i\right)})$ are considered, and the following relations should be derived:

$${\Sigma }_n^{*\left(i\right)}={\Sigma }_n^{\left(i\right)}-{\tilde{\Sigma }}_n^{\left(i\right)}\ge 0$$

(10)

where ${\tilde{\Sigma }}_n^{\left(i\right)}=\mathrm{diag}[{\overline{\tilde{\sigma }}}_u^{\left(i\right)}{I}_{n+1},{\overline{\tilde{\sigma }}}_y^{\left(i\right)}{I}_n,0]$, whilst $({\overline{\tilde{\sigma }}}_u^{\left(i\right)}$, ${\overline{\tilde{\sigma }}}_y^{\left(i\right)})$ represent the variances of input and output additive noises in the i-th cluster. The considered identification scheme normally assumes that [16,24]:

$$\left\{\begin{array}{ccc}u\left(k\right)& =& u^{*}\left(k\right)+\tilde{u}\left(k\right)\\ y\left(k\right)& =& y^{*}\left(k\right)+\tilde{y}\left(k\right)\end{array}\right.$$

(11)

where $u^{*}\left(k\right)$ and $y^{*}\left(k\right)$ are the noise-free data, noise terms $\tilde{u}\left(k\right)$ and $\tilde{y}\left(k\right)$ are independent of every other term and only $u\left(k\right)$ and $y\left(k\right)$ are known.

Finally, the matrices ${\tilde{\Sigma }}_n^{\left(i\right)}$ can therefore be built and the parameter of the model in each cluster determined by means of the relations:

$$\left({\Sigma }_n^{\left(i\right)}-{\tilde{\Sigma }}_n^{\left(i\right)}\right)a^{\left(i\right)}=0\mathrm{for}i=1,\dots ,K$$

(12)

for a number of K clusters.

4. Fault Diagnosis Scheme Design

It is assumed that the monitored system is described as represented in Figure 1. The prediction (or estimation error) $y\left(k\right)-\widehat{y}\left(k\right)$ in fault-free conditions represents the model-reality mismatch, which accounts for process noise, parameter variations, disturbance and uncertainty.

Figure 1. The monitored system.

Figure 1. The monitored system.

The relations of Equation (11) describe the realistic situation where the variables $u^{*}\left(k\right)$ and $y^{*}\left(k\right)$ are measured by means of sensors affected by both measurement noise and faults.

Neglecting the sensor dynamics, faults acting on the measured input and output signals $u\left(k\right)$ and $y\left(k\right)$ are modeled as:

$$\left\{\begin{array}{ccc}u\left(k\right)\hfill & =\hfill & u^{*}\left(k\right)+f_u\left(k\right)\hfill \\ y\left(k\right)\hfill & =\hfill & y^{*}\left(k\right)+f_y\left(k\right)\hfill \end{array}\right.$$

(13)

where $f_u\left(k\right)$ and $f_y\left(k\right)$ represent additive signals assuming values different from zero only in the presence of faults.

There are different approaches to generate the residuals for fault diagnosis; see, e.g., [12]. In this work, the TS models are used as residual generators for the wind turbine system. Figure 2 shows that the residuals are generated by the comparison of the measured $y\left(k\right)$ and the estimated outputs: $\widehat{y}\left(k\right)$:

$$r\left(k\right)=\widehat{y}\left(k\right)-y\left(k\right)$$

(14)

Figure 2. The residual generation scheme.

Figure 2. The residual generation scheme.

The symptom evaluation refers to a logic device that processes the redundant signals generated by the first block, in order to detect when a fault occurs and to univocally identify the unreliable actuator or sensor.

The fault detection task is performed here by using a simple thresholding logic, even if different strategies are available; see e.g., [12]. It is worth noting that the faults described in Section 2.1 may not be immediately detected, since the delay in the corresponding alarm normally depends on the fault mode. This situation is shown in Figure 3, where suitable fault detection thresholds are fixed according to Equation (15):

$$\left\{\begin{array}{c}\overline{r}-\delta {\sigma }_r\le r\left(k\right)\le \overline{r}+\delta {\sigma }_r\hfill \\ \mathrm{if}\mathrm{fault-free}\hfill \\ \\ r\left(k\right)<\overline{r}-\delta {\sigma }_r\mathrm{or}r\left(k\right)>\overline{r}+\delta {\sigma }_r\hfill \\ \mathrm{if}\mathrm{faulty}\hfill \end{array}\right.$$

(15)

In practice, the residual signal is represented by the random variable $r\left(k\right)$, whose sample mean value and variance values are estimated as follows:

$$\left\{\begin{array}{c}\overline{r}=\frac{1}{N}\sum _{k=1}^{N}r\left(k\right)\hfill \\ \\ {\sigma }_r^2=\frac{1}{N}\sum _{k=1}^N{\left[r\left(k\right)-\overline{r}\right]}^2\hfill \end{array}\right.$$

(16)

$\overline{r}$ and ${\sigma }_r^2$ are the values for the sample mean and variance of the fault-free residual, respectively. N is the number of samples of $r\left(k\right)$. The values of $\overline{r}$ and ${\sigma }_r^2$ depend on the signal $r\left(k\right)$ statistics, usually unknown.

In order to separate normal from faulty behavior, the tolerance parameter δ (normally $\delta \ge 2$) is selected and properly tuned. Hence, the proper choice of this parameter δ leads to a good trade-off between the maximization of the fault detection probability and the minimization of the false alarm rate. This parameter δ could be fixed with empirical rules or, once the values of $\overline{r}$ and ${\sigma }_r^2$ are estimated from the $r\left(k\right)$ signal, using the three-sigma rule. On the other hand, less conservative results are obtained with a procedure that determines via extensive simulations the optimal δ minimizing the false alarm rate and maximizing the detection/isolation probability. This issue will be addressed in Section 5.

Figure 3. Detection thresholds and fault detection and isolation (FDI) delay for incipient faults.

Figure 3. Detection thresholds and fault detection and isolation (FDI) delay for incipient faults.

Moreover, as shown in Figure 3, if a detection delay is tolerable, which depends on the fault severity, the amplitude of the detectable/isolable fault is lower.

Finally, regarding the fault isolation problem, a generalized observer scheme (GOS) is exploited [12]. In particular, as shown in Section 2.1, since different faults $f_u\left(k\right)$ or $f_y\left(k\right)$ can affect the input or output measurements, to uniquely isolate a fault $f_u\left(k\right)$ concerning one of the inputs, under the assumption that the outputs are fault-free, a bank of estimators in the form of Equation (5) is used, as shown in Figure 4.

Figure 4. Fuzzy estimator scheme for actuator fault isolation.

Figure 4. Fuzzy estimator scheme for actuator fault isolation.

The number of these estimators is equal to the number of the faults $f_u\left(t\right)$ that have to be diagnosed. The i-th fuzzy estimator is driven by all but the i-th input (or even more inputs, if required) and all outputs of the system and generates a residual function, which is sensitive to all but the i-th input fault $f_u\left(k\right)$ (or even more inputs, if necessary). The derivation of these fuzzy estimators follows the procedure described in Section 3. In particular, when the fuzzy estimator insensitive to the i-th input has to be designed, the output $y\left(k\right)$ and all but the i-th inputs $u_i\left(k\right)$ are exploited for the identification process.

On the other hand, to uniquely isolate a fault $f_y\left(t\right)$ concerning one of the system outputs, under the hypothesis that inputs are fault-free, a bank of estimators is used again, according to Figure 5.

This observer configuration represents the dedicated observer scheme (DOS) described in [12]. The number of these estimators is equal to the number of faults $f_y\left(t\right)$ that have to be diagnosed, and each device is driven by a single output and all of the inputs of the system. In this case, a fault on the i-th output affects only the residual function of the output observer or filter driven by the i-th output.

In order to summarize the isolation capabilities of the schemes presented, Table 4 shows the ‘fault signatures’ for the case of single fault occurrence. Note that in Table 4, the residual $r_{I_i}$ ($i=1,\cdots ,r$) coincides with the signal $r_i$ of Figure 4. On the other hand, $r_{O_j}$ ($j=1,\cdots ,m$) represents the residual $r_i$ of Figure 5 for output fault isolation.

Figure 5. Fuzzy estimators for sensor fault isolation.

Figure 5. Fuzzy estimators for sensor fault isolation.

Table 4. Fault signatures.

	$u_1$	$u_2$	…	$u_r$	$y_1$	$y_2$	…	$y_m$
$r_{I_1}$	0	1	…	1	1	1	…	1
$r_{I_2}$	1	0	…	1	1	1	…	1
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
$r_{I_r}$	1	1	…	0	1	1	…	1
$r_{O_1}$	1	1	…	1	1	0	…	0
$r_{O_2}$	1	1	…	1	0	1	…	0
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
$r_{O_m}$	1	1	…	1	0	0	…	1

Table 4. Fault signatures.

	$u_1$	$u_2$	…	$u_r$	$y_1$	$y_2$	…	$y_m$
$r_{I_1}$	0	1	…	1	1	1	…	1
$r_{I_2}$	1	0	…	1	1	1	…	1
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
$r_{I_r}$	1	1	…	0	1	1	…	1
$r_{O_1}$	1	1	…	1	1	0	…	0
$r_{O_2}$	1	1	…	1	0	1	…	0
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
$r_{O_m}$	1	1	…	1	0	0	…	1

The residuals affected by input and output faults are described by an entry ‘1’ in the corresponding table entry, while an entry ‘0’ means that the input or output fault does not affect the corresponding residual.

Note how multiple output faults can be isolated, since a fault on the i-th output signal affects only the residual function $r_{O_i}$ of the output estimator driven by the i-th output, but all of the residual functions $r_{I_j}$. On the other hand, multiple faults on the inputs cannot be isolated, since, in general, all residual functions are sensitive to faults regarding different inputs.

Finally, in order to summarize the complete procedure, the different design phases are summarized in Figure 6.

Once the data have been collected and sampled from the wind turbine (Step 1), the FMEA recalled in Section 2.1 is applied. This procedure suggests how to select the measured signals $u_i\left(k\right)$ and $y_j\left(k\right)$ (Step 2) in order to build the fuzzy estimator banks (Step 5) described in Figure 4 and Figure 5. The different fuzzy estimator models have the form of Equation (5) derived using the fuzzy clustering (Step 3) recalled in Section 3.1, followed by the structure identification (Step 4) in Section 3.2. In this way, by means of the threshold test logic of Equation (15), the fault diagnosis is achieved (Step 6).

Figure 6. Sketch of the complete design procedure.

Figure 6. Sketch of the complete design procedure.

5. Simulation Results

The proposed FDI methodology was applied to a sequence of $N=440\times {10}^3$ data samples $u\left(k\right)$ and $y\left(k\right)$ acquired with a sampling rate of 100 Hz. from the wind turbine benchmark.

According to Section 3 and Section 4, the Gustafson–Kessel (GK) clustering method with $K=4$ clusters and a number of shifts $n=3$ was used for the identification of the fuzzy estimator banks of Section 4. These optimal parameters $K=4$ and $n=3$ were obtained as described in [16,17]. After clustering, the parameters $a_i$ and $b_i$, with $i=1,\cdots ,K$, were estimated using the identification method presented in Section 3. Moreover, the membership degrees ${\lambda }_i$ required by the fuzzy estimators of Equation (5) have been modeled as Gaussian functions.

As shown in Figure 4 and Figure 5, the reconstructed output ${\widehat{y}}_i\left(k\right)$ for the FDI task has been generated by a bank of five multiple-input single-output (MISO) predictors of Equation (5). According to Table 2 and Figure 4, this scheme allows the diagnosis of Fault${}_1$, Fault${}_2$, Fault${}_3$, Fault${}_4$ and Fault${}_5$. On the other hand, with reference again to Table 2 and Figure 5, a bank of four output fuzzy estimators for ${\widehat{y}}_i\left(k\right)$ allows the diagnosis of Fault${}_6$, Fault${}_7$, Fault${}_8$ and Fault${}_9$.

For each fault case, by following the FMEA procedure described in Section 2.1 and Table 2, the input and output measurements used for the design of the estimator banks were reported in Table 5.

The approximation capabilities of the fuzzy residual generators can be expressed in fault-free conditions in terms of the so-called variance accounted for (VAF) index [15]. In particular, the VAF values for all identified MISO estimators were always bigger than 99%. Hence, the multiple model scheme seems to approximate the process outputs quite accurately. Note, in fact, that, as described in [16,17], with the choice of parameters $K=4$ and $n=3$, the fuzzy predictors led to the minimization of the reconstruction errors, i.e., the difference between the measured and predicted outputs.

Table 5. Inputs and outputs for the fuzzy residual generator design.

Fault	Inputs	Output
1	$\left[{\beta }_{1m_1}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
2	$\left[{\beta }_{2m_2}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
3	$\left[{\beta }_{3m_1}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
4	$\left[{\beta }_{1m_2}\left(t\right),{\omega }_{g{m}_2}\left(t\right)\right]$	${\omega }_{r{m}_1}\left(t\right)$
5	$\left[{\beta }_{1m_2}\left(t\right),{\omega }_{g{m}_2}\left(t\right)\right]$	${\omega }_{r{m}_2}\left(t\right)$
6	$\left[{\beta }_{2m_1}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
7	$\left[{\beta }_{3m_2}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
8	$\left[{\beta }_{1m_2}\left(t\right),{\tau }_{gm}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
9	$\left[{\beta }_{1m_2}\left(t\right),{\omega }_{g{m}_1}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$

Table 5. Inputs and outputs for the fuzzy residual generator design.

Fault	Inputs	Output
1	$\left[{\beta }_{1m_1}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
2	$\left[{\beta }_{2m_2}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
3	$\left[{\beta }_{3m_1}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
4	$\left[{\beta }_{1m_2}\left(t\right),{\omega }_{g{m}_2}\left(t\right)\right]$	${\omega }_{r{m}_1}\left(t\right)$
5	$\left[{\beta }_{1m_2}\left(t\right),{\omega }_{g{m}_2}\left(t\right)\right]$	${\omega }_{r{m}_2}\left(t\right)$
6	$\left[{\beta }_{2m_1}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
7	$\left[{\beta }_{3m_2}\left(t\right),{\beta }_{1m_2}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
8	$\left[{\beta }_{1m_2}\left(t\right),{\tau }_{gm}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$
9	$\left[{\beta }_{1m_2}\left(t\right),{\omega }_{g{m}_1}\left(t\right)\right]$	${\omega }_{g{m}_2}\left(t\right)$

The rationale of using TS fuzzy models was highlighted in Section 3, whilst their efficacy is analyzed in the following. Figure 7 reports the results of the correlation analysis on the input variables. In fact, in fault-free conditions, the residuals $y\left(k\right)-\widehat{y}\left(k\right)$ associated with the data of each cluster ($K=4$) and the identified affine models should be ideally white and independent of the inputs $u\left(k\right)$. This situation guarantees that the estimators approximate correctly the measurements $y\left(k\right)$ in each cluster.

Figure 7. Residual auto– and cross–correlation examples ($K=4$).

Figure 7. Residual auto– and cross–correlation examples ($K=4$).

Figure 7 shows the estimator residuals $y\left(k\right)-\widehat{y}\left(k\right)$ in fault-free conditions for each cluster ($K=4$), thus highlighting their whiteness and independence. In particular (1) the auto-correlation function of $y\left(k\right)-\widehat{y}\left(k\right)$ and (2) the cross-correlation function between $y\left(k\right)-\widehat{y}\left(k\right)$ and $u\left(k\right)$ are displayed for 20 lags. For these variables, the 99% confidence intervals are also depicted as dotted lines, thus showing that $y\left(k\right)-\widehat{y}\left(k\right)$ is white and independent of $u\left(k\right)$ for each cluster. These results prove experimentally the validity of the TS models proposed in this study.

It is worth noting that the nonlinear benchmark originally developed in [25] was modified by the authors in order to vary the statistical properties of the signals used for modeling process parameter uncertainty and measurement errors. Under this assumption, Table 6 reports the nominal values of the considered wind turbine model parameters with respect to their simulated uncertainty. In this way, a Monte Carlo analysis can be performed for assessing the reliability and the robustness of the considered FDI scheme by modeling the model variables as Gaussian stochastic processes, with zero-mean and standard deviations corresponding to realistic minimal and maximal error values of Table 6.

Table 6. Realistic wind turbine uncertainty.

Variable	Nominal Value	Min Error	Max Error
ρ	$1.225$ kg/m${}^3$	±0.1%	±20%
J	$7.794\times {10}^6$ kg/m${}^2$	±0.1%	±30%
$C_p$	$C_{p0}$	±0.1%	±50%
u	$u_0$	±0.1%	±20%
y	$y_0$	±0.1%	±20%

Table 6. Realistic wind turbine uncertainty.

Variable	Nominal Value	Min Error	Max Error
ρ	$1.225$ kg/m${}^3$	±0.1%	±20%
J	$7.794\times {10}^6$ kg/m${}^2$	±0.1%	±30%
$C_p$	$C_{p0}$	±0.1%	±50%
u	$u_0$	±0.1%	±20%
y	$y_0$	±0.1%	±20%

It is also assumed that the input-output signals u and y and the power coefficient map $C_p$ entries were affected by errors, expressed as percent standard deviations of the corresponding nominal values $u_0$, $y_0$ and $C_{p0}$, also reported in Table 6. Therefore, for the performance evaluation of the FDI methodologies, a sufficient number of Monte Carlo runs was performed.

Note that Table 6 describes the uncertain parameters that have been simulated in order to analyze the reliability and the robustness features of the proposed approach with respect to parameter variations. In fact, the approach was proposed here also for removing the effect of the uncertain wind term $v\left(t\right)$ and not for handling the parameter variations summarized in Table 6.

The simulations of different fault cases have been reported for highlighting the most important features of the proposed approach. In particular, the first example was obtained by considering the fault Case 1, commencing at the instant $t=2,000$ s, and active for 100 s. The considered fault $f_u\left(t\right)$ causes alteration of the signals $u\left(t\right)$ and $y\left(t\right)$ and, therefore, of the residuals $r_{I_i}\left(t\right)$ given by the model of Equation (5). These residuals indicate the fault occurrence according to the logic of Equation (15), whether their values are lower or higher than the thresholds fixed in fault-free conditions. Figure 8 represents the fault-free (grey continuous line) and the faulty (black dashed line) residuals $r_{I_i}\left(t\right)$.

The fault detection thresholds of Equation (15) are represented as dotted constant lines in Figure 8. Their values were properly settled, as described in Section 5.2, in order to minimize the false alarm and missed fault rates, while maximizing the correct detection and isolation rates. In these conditions, the fault is correctly detected and isolated when the corresponding residual signals exceed the thresholds, as indicated in Figure 8.

Figure 8. Residuals $r_{I_i}\left(t\right)$ for the fault Case 1.

Figure 8. Residuals $r_{I_i}\left(t\right)$ for the fault Case 1.

The second example depicted in Figure 9 represents the fault-free (grey continuous line) and the faulty (black dashed line) residuals $r_{O_i}\left(t\right)$ related to the fault case 4, i.e., $f_y\left(t\right)$. It commences at the instant $t=1500$ s, and it is active for 100 s. The fault detection thresholds, represented as dotted constant lines in Figure 9, were optimally fixed, as in the previous case.

The third example was obtained by considering the fault Case 8, i.e., $f_u\left(t\right)$, which is active between 3800 s and 3900 s. Figure 10 represents the fault-free (grey continuous line) and the faulty (black dashed line) residuals $r_{I_i}\left(t\right)$.

Finally, Figure 11 depicts the fault-free (grey continuous line) and the faulty (black dashed line) residuals $r_{O_i}\left(t\right)$ related to the fault Case 9, i.e., $f_y\left(t\right)$, active between 4000 s and 4200 s. Again, the FDI thresholds represented in Figure 11 were optimally fixed, as in the previous case.

Figure 9. Residuals $r_{O_i}\left(t\right)$ for the fault Case 4.

Figure 9. Residuals $r_{O_i}\left(t\right)$ for the fault Case 4.

Figure 10. Residuals $r_{I_i}\left(t\right)$ for the fault Case 8.

Figure 10. Residuals $r_{I_i}\left(t\right)$ for the fault Case 8.

Figure 11. Residuals $r_{O_i}\left(t\right)$ for the fault Case 9.

Figure 11. Residuals $r_{O_i}\left(t\right)$ for the fault Case 9.

5.1. Comparative Studies

This section provides some comparative results with respect to different FDI schemes.

The first alternative approach considered here uses a support vector machine based on a Gaussian kernel (GKSV) developed in [8]. The scheme defines a vector of features for each fault, which contains relevant signals obtained directly from measurements, filtered measurements or their combinations. These vectors are subsequently projected onto the kernel of the support vector machine (SVM), which provides suitable residuals for all of the defined faults. Different kernels have been tested, and it was found that Gaussian kernel with different variance values can be used for all faults. Data with and without faults were used for learning the model for the FDI of the specific faults.

The second scheme consists of an estimation-based (EB) solution shown in [10]. In particular, a fault detection estimator is designed to detect a fault, and an additional bank of estimators is derived to isolate them. The method was designed on the basis of a system linear model and used fixed thresholds, as in Equation (15). Each estimator for fault isolation was computed on the basis of the particular fault scenario under consideration.

The third method relying on up-down counters (UDCs) was addressed in [11]. These tools, borrowed from the aerospace framework, were exploited in the decision logic applied to the FDI residuals. These residuals were obtained using both physical and analytical redundancy schemes, such as parity equations from redundant sensors and Kalman filters. This approach is different from the straightforward thresholding of Equation (15). In fact, the decision to declare the fault occurrence involves discrete-time dynamics and is not simply a function of the residual current value.

The fourth approach combines observer and Kalman filter (COK) methods [7]. It relies on an observer used as a residual generator for diagnosing the faults of the drive-train, in which the wind speed is considered a disturbance. This diagnosis observer was designed to decouple the disturbance and simultaneously achieve optimal residual generation in a statistical sense. For the other two subsystems of the wind turbine, a Kalman filter-based approach was applied. The residual evaluation task used a generalized likelihood ratio test, and cumulative variance indices were applied. For fault isolation purpose, a bank of residual generators was exploited. Sensor and system faults were thus isolated via a decision table.

The fifth method relies on the general fault model (GFM) scheme, which is a method of automatic design [9]. The FDI strategy consists of three main steps. In the first step, a large set of potential residual generators was designed. In the second step, the most suitable residual generators to be included in the final FDI system were selected. In the third step, tests for the selected set of residual generators were performed, which were based on comparisons of the estimated probability distributions of the residuals, evaluated with fault-free and faulty data.

For performance evaluation and comparison of the considered FDI schemes, some indices have been used. They were presented in [26] and here evaluated on $1,000$ Monte Carlo runs. These indices are defined as:

• False Alarm Rate ($r_{fa}$): the number of wrongly detected faults divided by total fault cases;
• Missed Fault Rate ($r_{mf}$): for each fault, the total number of undetected faults, divided by the total number of times that the fault case occurs;
• True Detection/Isolation Rate ($r_{tdi}$): for a particular fault case, the number of times it is correctly detected/isolated, divided by total number of times that the fault case occurs;
• Mean Detection/Isolation Delay (${\tau }_{mdi}$): for a particular fault case, the average detection/isolation delay time.

These criteria are computed for each fault case and for each FDI scheme. Table 7 summarizes also the results obtained by considering the fuzzy predictors as residual generators (FPRG) and with an optimal choice of the threshold parameter δ in Equation (15) that leads to achieving optimal results.

Several comments can be drawn here. GKSV is able to detect and isolate Faults 1, 2, 3, 4, 5 and 8 and some of them with delays bigger than 25 s. For these diagnosable faults, the average detection rate is bigger than 65%, with missed fault and false alarm rates lower than 35%. Moreover, in general, this scheme showed robustness with respect to the working point changes of the wind turbine, but effective only on sensor faults. EB manages the faults quite quickly, apart from Fault 9. For the detectable faults, the detection rate is bigger than 66%, with missed fault and false alarm rates lower than 33%. The diagnosis delays for Faults 2, 6 and 7 can be bigger than 11 s. In particular, for Faults 2 and 9, false alarms can occur. UDC detects and isolates almost all faults, apart from Fault 9, with delay times lower than 69 s. However, false alarm rates bigger than 12% are measured for Faults 2, 3, 4, 5, 6, 7 and 8. In the same way, COK is able to detect and isolate faults, apart from 9, and in general, with a delay time bigger that 10 s. False alarm and missed fault rates bigger that 10% can also occur. The same considerations holds for the GFM solution. The detection delay times can be bigger than 9 s, with false alarm and missed fault rates bigger that 12%. Regarding the proposed FPRG method, it seems to work relatively better than the others, even if optimization stages are required, for example, for the optimal FDI threshold selection. For this method, in general, also for Fault 9, the detection rates are bigger than 83%, with false alarm and missed fault rates lower than 14%. The issue of the optimal threshold selection will be analyzed in Section 5.2. Note finally that the simulator used in this study was modified by the authors, and the working conditions can be slightly different from the ones shown in [6].

5.2. Robustness Evaluation

This section reports further experimental results regarding the performance optimization of the developed FDI scheme with respect to modeling errors and measurement uncertainty. In particular, the simulation of different fault-free and faulty data sequences has been performed by exploiting again the wind turbine simulator and the Monte Carlo method. In fact, the Monte Carlo tool is useful at this stage, since the efficacy of the FDI module depends on both the model approximation and the measurement errors.

Therefore, the indices defined above have been evaluated for each fault case. In particular, Table 8 summarizes the results obtained by considering the fuzzy predictors as residual generators and with a choice of the parameter δ of Equation (15) that leads to optimal performances.

Table 8 shows that the proper selection of the threshold levels of Equation (15) depending on δ allows one to achieve false alarm and missed fault rates of less than 13% and detection/isolation rates larger than 83%, with minimal detection/isolation delay times. The results demonstrate also that Monte Carlo analysis is an effective tool for experimentally tuning and testing the suggested FDI method. In the presence of uncertainty and modeling errors, this latter simulation technique seems to facilitate the assessment of the reliability of the developed FDI methods for application to real test cases.

6. Conclusions

This paper proposed a procedure for the fault detection and isolation of a wind turbine model using fuzzy models identified from uncertain input-output measurements. It was assumed that the process under investigation is nonlinear, and its available measurements were normally not very reliable, due to the wind speed’s uncertain nature. The residual generators considered here for diagnosis purposes have the form of Takagi–Sugeno models. These fuzzy models were derived using fuzzy clustering and dynamic model linear identification. The effectiveness of the proposed approach was tested on the data acquired from a simulated wind turbine benchmark. The detection and isolation of the faults affecting sensors, component and actuators of the process under diagnosis was thus achieved. Future investigations will concern the application of the diagnosis strategy to real wind turbine installations.

Table 7. Comparison of the considered FDI strategies. GKSV, Gaussian kernel support vector; EB, estimation-based; UDC, up-down counter; COK, combined observer and Kalman filter; GFM, general fault model; FPRG, fuzzy predictors as residual generators.

Fault	Indices	GKSV	EB	UDC	COK	GFM	FPRG
1	$r_{fa}$	0.001	0.001	0.001	0.001	0.001	0.001
	$r_{mf}$	0.002	0.003	0.002	0.003	0.002	0.001
	$r_{tdi}$	0.978	0.977	0.987	0.977	0.982	0.999
	${\tau }_{mdi}$	0.03s.	0.03s.	0.04s.	10.32s.	0.05s.	0.02s.
2	$r_{fa}$	0.234	0.224	0.123	0.003	0.235	0.001
	$r_{mf}$	0.343	0.333	0.232	0.029	0.532	0.003
	$r_{tdi}$	0.657	0.667	0.768	0.971	0.468	0.997
	${\tau }_{mdi}$	47.24s.	44.65s.	69.03s.	19.32s.	13.74s.	0.08s.
3	$r_{fa}$	0.004	0.141	0.123	0.056	0.135	0.003
	$r_{mf}$	0.006	0.132	0.241	0.128	0.232	0.008
	$r_{tdi}$	0.974	0.868	0.769	0.872	0.768	0.992
	${\tau }_{mdi}$	0.05s.	0.54s.	0.05s.	19.32s.	0.74s.	0.02s.
4	$r_{fa}$	0.006	0.005	0.123	0.056	0.236	0.004
	$r_{mf}$	0.005	0.006	0.113	0.128	0.333	0.004
	$r_{tdi}$	0.975	0.994	0.887	0.872	0.667	0.996
	${\tau }_{mdi}$	0.15s.	0.33s.	0.04s.	19.32s.	17.64s.	0.02s.
5	$r_{fa}$	0.178	0.004	0.234	0.256	0.236	0.002
	$r_{mf}$	0.223	0.005	0.254	0.329	0.242	0.003
	$r_{tdi}$	0.777	0.995	0.746	0.671	0.758	0.997
	${\tau }_{mdi}$	25.95s.	0.07s.	0.04s.	31.32s.	9.49s.	0.03s.
6	$r_{fa}$	0.897	0.173	0.334	0.156	0.096	0.042
	$r_{mf}$	0.987	0.234	0.257	0.129	0.042	0.033
	$r_{tdi}$	0.013	0.766	0.743	0.871	0.958	0.967
	${\tau }_{mdi}$	95.95s.	11.37s.	12.94s.	34.02s.	9.49s.	3.03s.
7	$r_{fa}$	0.899	0.044	0.134	0.134	0.123	0.047
	$r_{mf}$	0.899	0.035	0.121	0.101	0.098	0.023
	$r_{tdi}$	0.101	0.965	0.879	0.899	0.902	0.977
	${\tau }_{mdi}$	99.95s.	26.17s.	13.93s.	35.01s.	29.79s.	5.07s.
8	$r_{fa}$	0.004	0.045	0.144	0.109	0.099	0.003
	$r_{mf}$	0.007	0.011	0.101	0.032	0.124	0.002
	$r_{tdi}$	0.993	0.989	0.899	0.968	0.876	0.998
	${\tau }_{mdi}$	0.07s.	0.08s.	0.09s.	0.06s.	8.94s.	0.05s.
9	$r_{fa}$	0.778	0.879	0.894	0.956	0.995	0.134
	$r_{mf}$	0.996	0.934	0.947	0.929	0.941	0.165
	$r_{tdi}$	0.004	0.066	0.053	0.071	0.059	0.835
	${\tau }_{mdi}$	95.95s.	91.37s.	92.94s.	94.02s.	99.49s.	0.30s.

Table 7. Comparison of the considered FDI strategies. GKSV, Gaussian kernel support vector; EB, estimation-based; UDC, up-down counter; COK, combined observer and Kalman filter; GFM, general fault model; FPRG, fuzzy predictors as residual generators.

Fault	Indices	GKSV	EB	UDC	COK	GFM	FPRG
1	$r_{fa}$	0.001	0.001	0.001	0.001	0.001	0.001
	$r_{mf}$	0.002	0.003	0.002	0.003	0.002	0.001
	$r_{tdi}$	0.978	0.977	0.987	0.977	0.982	0.999
	${\tau }_{mdi}$	0.03s.	0.03s.	0.04s.	10.32s.	0.05s.	0.02s.
2	$r_{fa}$	0.234	0.224	0.123	0.003	0.235	0.001
	$r_{mf}$	0.343	0.333	0.232	0.029	0.532	0.003
	$r_{tdi}$	0.657	0.667	0.768	0.971	0.468	0.997
	${\tau }_{mdi}$	47.24s.	44.65s.	69.03s.	19.32s.	13.74s.	0.08s.
3	$r_{fa}$	0.004	0.141	0.123	0.056	0.135	0.003
	$r_{mf}$	0.006	0.132	0.241	0.128	0.232	0.008
	$r_{tdi}$	0.974	0.868	0.769	0.872	0.768	0.992
	${\tau }_{mdi}$	0.05s.	0.54s.	0.05s.	19.32s.	0.74s.	0.02s.
4	$r_{fa}$	0.006	0.005	0.123	0.056	0.236	0.004
	$r_{mf}$	0.005	0.006	0.113	0.128	0.333	0.004
	$r_{tdi}$	0.975	0.994	0.887	0.872	0.667	0.996
	${\tau }_{mdi}$	0.15s.	0.33s.	0.04s.	19.32s.	17.64s.	0.02s.
5	$r_{fa}$	0.178	0.004	0.234	0.256	0.236	0.002
	$r_{mf}$	0.223	0.005	0.254	0.329	0.242	0.003
	$r_{tdi}$	0.777	0.995	0.746	0.671	0.758	0.997
	${\tau }_{mdi}$	25.95s.	0.07s.	0.04s.	31.32s.	9.49s.	0.03s.
6	$r_{fa}$	0.897	0.173	0.334	0.156	0.096	0.042
	$r_{mf}$	0.987	0.234	0.257	0.129	0.042	0.033
	$r_{tdi}$	0.013	0.766	0.743	0.871	0.958	0.967
	${\tau }_{mdi}$	95.95s.	11.37s.	12.94s.	34.02s.	9.49s.	3.03s.
7	$r_{fa}$	0.899	0.044	0.134	0.134	0.123	0.047
	$r_{mf}$	0.899	0.035	0.121	0.101	0.098	0.023
	$r_{tdi}$	0.101	0.965	0.879	0.899	0.902	0.977
	${\tau }_{mdi}$	99.95s.	26.17s.	13.93s.	35.01s.	29.79s.	5.07s.
8	$r_{fa}$	0.004	0.045	0.144	0.109	0.099	0.003
	$r_{mf}$	0.007	0.011	0.101	0.032	0.124	0.002
	$r_{tdi}$	0.993	0.989	0.899	0.968	0.876	0.998
	${\tau }_{mdi}$	0.07s.	0.08s.	0.09s.	0.06s.	8.94s.	0.05s.
9	$r_{fa}$	0.778	0.879	0.894	0.956	0.995	0.134
	$r_{mf}$	0.996	0.934	0.947	0.929	0.941	0.165
	$r_{tdi}$	0.004	0.066	0.053	0.071	0.059	0.835
	${\tau }_{mdi}$	95.95s.	91.37s.	92.94s.	94.02s.	99.49s.	0.30s.

Table 8. Monte Carlo analysis with δ of Equation (15).

Fault	$r_{fa}$	$r_{mf}$	$r_{tdi}$	${\tau }_{mdi}$	δ
1	$0.001$	$0.001$	$0.999$	$0.02$s.	$3.9$
2	$0.001$	$0.003$	$0.997$	$0.08$s.	$4.1$
3	$0.003$	$0.008$	$0.992$	$0.02$s.	$3.9$
4	$0.004$	$0.004$	$0.996$	$0.02$s.	$4.3$
5	$0.002$	$0.003$	$0.997$	$0.03$s.	$3.5$
6	$0.042$	$0.033$	$0.967$	$3.03$s.	$4.6$
7	$0.047$	$0.023$	$0.977$	$5.07$s.	$4.5$
8	$0.003$	$0.002$	$0.998$	$0.05$s.	$3.7$
9	$0.134$	$0.165$	$0.835$	$0.30$s.	$2.8$

Table 8. Monte Carlo analysis with δ of Equation (15).

Fault	$r_{fa}$	$r_{mf}$	$r_{tdi}$	${\tau }_{mdi}$	δ
1	$0.001$	$0.001$	$0.999$	$0.02$s.	$3.9$
2	$0.001$	$0.003$	$0.997$	$0.08$s.	$4.1$
3	$0.003$	$0.008$	$0.992$	$0.02$s.	$3.9$
4	$0.004$	$0.004$	$0.996$	$0.02$s.	$4.3$
5	$0.002$	$0.003$	$0.997$	$0.03$s.	$3.5$
6	$0.042$	$0.033$	$0.967$	$3.03$s.	$4.6$
7	$0.047$	$0.023$	$0.977$	$5.07$s.	$4.5$
8	$0.003$	$0.002$	$0.998$	$0.05$s.	$3.7$
9	$0.134$	$0.165$	$0.835$	$0.30$s.	$2.8$