Code of An Estimation-based Approach of Fault Detection and Isolation of the Wind TurbineBenchmark

by

 

The matlab code of the solution of “FDI Design Competition for Wind Turbines" (winning 2nd place)

FDIBenchMarkData.m
%close all
clear al
 
load winddata.mat
load AeroDynamics.mat


% script containing parameters for the kk-electronic - FDI Bench Mark Model
% By Peter Fogh Odgaard, kk-electronic a/s, Viby J Denmark
% Date 7.11.2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ts=1/100;
Time=Ts:Ts:4400;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Wind data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D_v_wind=[Time', 9+4*sin(0.01*Time)'];

load winddata
D_v_wind=[windtime windspeed];
%% Wind Model
seed1 = 256;
seed2 = 894;
turbulence_seeds=[seed1 seed2];
R=57.5;
H=87;
k =4.7;
a=2.2;
alpha = 0.1;
D_rotor=2*R;
r = D_rotor/2;
m = 1+alpha*(alpha-1)*R^2/(8*H^2);
Length_scale = 600; %[m]
L=Length_scale;
Turbulence_intensity = 12;
turb_int=Turbulence_intensity;
T_sample=0.05;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Pitch and Blade Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

omega_n=11.11;
xi=0.6;
rho=1.225;
R=57.5;
r0 = 1.5;

% cq table
load AeroDynamics

%%Fault models%%
%fault 6
xi2=0.45;
omega_n2=5.73;
%fault 7
xi3=0.9;
omega_n3=3.42;

%transfers to ss models
[Apb,Bpb,Cpb,Dpb]=tf2ss([omega_n^2],[1 2*xi*omega_n omega_n^2]);
[Apb1,Bpb1,Cpb1,Dpb1]=tf2ss([omega_n2^2],[1 2*xi2*omega_n2 omega_n2^2]);
[Apb2,Bpb2,Cpb2,Dpb2]=tf2ss([omega_n3^2],[1 2*xi3*omega_n3 omega_n3^2]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Drive Train Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_dt=775.49;
B_r=7.11;
B_g=45.6;
N_g=95;
K_dt=2.7e9
eta_dt=0.97;
J_g=390
J_r=55e6

Addt=[-(B_dt+B_r)/J_r B_dt/N_g/J_r -K_dt/J_r; eta_dt*B_dt/N_g/J_g -(eta_dt*B_dt/N_g^2+B_g)/J_g eta_dt*K_dt/N_g/J_g; 1 -1/N_g 0];
Bddt=[1/J_r 0; 0 -1/J_g;0 0];
Cddt=[1 0 0;0 1 0];
Dddt=[0 0;0 0];


%fault models
%fault 9
eta_dt2=.92

Addt2=[-(B_dt+B_r)/J_r B_dt/N_g/J_r -K_dt/J_r; eta_dt2*B_dt/N_g/J_g -(eta_dt2*B_dt/N_g^2+B_g)/J_g eta_dt2*K_dt/N_g/J_g; 1 -1/N_g 0];
Bddt2=[1/J_r 0; 0 -1/J_g;0 0];
Cddt2=[1 0 0;0 1 0];
Dddt2=[0 0;0 0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generator & Converter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha_gc=1/20e-3;
eta_gc=0.98;

%Fault models
%fault 8
Constant_tau_gc=100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Controller
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K_opt=.5* rho* pi*R^2*eta_dt*R^3*0.4554/(N_g*8)^3

K_i=1;
K_p=4;
Omega_nom=162;
Omega_delta=15;
P_r=4.8e6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sensors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m_vw=1.5;
sigma_vm=0.5;
Np_vm=0.0071;
m_omega_r=0;
sigma_omega_r=0.004*2*pi;
Np_omega_r=2.3e-4;
m_omega_g=0;
sigma_omega_g=0.05;
Np_omega_g=5e-4;
m_tau_g=0;
sigma_tau_g=90;
Np_tau_g=.9;
m_P_g=0;
sigma_P_g=1e3;
Np_P_g=10;
m_Beta=0;
sigma_Beta=0.2;
Np_beta=1.5e-3;

%fault models
Constant_Beta_1_m1=5
Gain_Beta_2_m2=1.2
Constant_Beta_3_m1=10
Constant_Omega_r_m1=1.4
Gain_Omega_r_m2=1.1
Gain_Omega_g_m1=0.9

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fault signals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

D_Fault1=[Time' 1-[zeros(1,2000/Ts) ones(1,100/Ts) zeros(1,2300/Ts)]'];
%D_Fault1=[Time' 1-[zeros(1,200/Ts) ones(1,100/Ts) zeros(1,4100/Ts)]'];
D_Fault2=[Time' 1-[zeros(1,2300/Ts) ones(1,100/Ts) zeros(1,2000/Ts)]'];
D_Fault3=[Time' 1-[zeros(1,2600/Ts) ones(1,100/Ts) zeros(1,1700/Ts)]'];

D_Fault4=[Time' 1-[zeros(1,1500/Ts) ones(1,100/Ts) zeros(1,2800/Ts)]'];
D_Fault5=[Time' 1-[zeros(1,1000/Ts) ones(1,100/Ts) zeros(1,3300/Ts)]'];
D_Fault6=[Time' 1-[zeros(1,2900/Ts) ones(1,100/Ts) zeros(1,1400/Ts)]'];
D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];

D_Fault8=[Time' 1-[zeros(1,3800/Ts) ones(1,100/Ts) zeros(1,500/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,4000/Ts) ones(1,200/Ts) zeros(1,200/Ts)]'];
D_Fault9=[Time' 1-[zeros(1,200/Ts) ones(1,200/Ts) zeros(1,4000/Ts)]'];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Noise Seeds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

min_seed = 0;
max_seed = 999;
% Generate random seeds for 13 'Random Number' Generators blocks
seed = min_seed + (max_seed-min_seed).*rand(13, 1);
% Round up to integer
seed = ceil(seed);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% FDI design initialization
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% % Np_vm=0;
% % Np_omega_r=0;
% % Np_omega_g=0;
% % Np_tau_g=0;
% % Np_P_g=0;
% % Np_beta=0;
% 
% 
% %pitch system enhanced
% Ppb = [-2; -2.5];
% A_f = .5*Bpb*Cpb;
% Apb_E = Apb+A_f;
% %Apb_E = Apb;
% Lpb = place(Apb_E',Cpb',Ppb);
% Lpb = Lpb';
% Apb_E_bar = Apb_E-Lpb*Cpb;
% 
% %pitch system original
% Lpb_FD = place(Apb',Cpb',Ppb);
% Lpb_FD = Lpb_FD';
% Apb_bar = Apb-Lpb_FD*Cpb;
% 
% %drive train system
% A_DT = Addt(2:3,2:3);
% B_DT = [Bddt(2:3,2) Addt(2:3,1)];
% C_DT = [1 0 ];
% P_DT = [-2; -2.5];
% L_DT = place(A_DT',C_DT',P_DT);
% L_DT = L_DT';
% A_DT_bar = A_DT-L_DT*C_DT;

Contact us