【故障诊断】风力涡轮机自动故障检测和隔离的最优解决方案研究（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客 

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、Simulink、文章讲解

---

💥1 概述

文献来源：

摘要将支持向量机(SVM)用于由三叶片和全变流器组成的变速水平轴风力机的故障检测与隔离。支持向量机方法是基于数据的，因此对处理知识具有鲁棒性。此外，它基于结构风险最小化，增强了泛化能力，并允许使用灵活的内核来处理过程非线性。本文采用径向基函数作为核函数。研究了该工艺的不同部分，包括执行器、传感器和工艺故障。使用重复传感器，我们可以快速检测叶片间距位置、发电机和转子转速的传感器故障(固定值故障为2个采样周期)，但在特定的故障量级限制下。在两个采样周期内可以检测到变矩器的转矩故障(驱动器)。在螺距系统的执行器故障无法检测到。

工艺故障主要是由于传动系统摩擦引起的损坏。在高震级误差约束下，可以检测出其故障。

关键词:故障检测，支持向量机，风力发电机。

摘要:风力发电装机容量在全球范围内急剧增加;这意味着风力涡轮机的可靠性更加重要。这项工作的一部分是改进风力发电机组的故障检测和调节方案。本文提出了一个模拟故障检测和调节方案的基准模型。该基准模型在系统层面上处理风力涡轮机，包括传感器、执行器和桨距系统、传动系统、发电机和转换器系统中的系统故障。

今天，风力涡轮机在世界电力生产中所占的比例越来越大，与此同时，标准涡轮机的尺寸也在增加。目前最常安装的兆瓦级风力涡轮机价格昂贵，因此这些涡轮机的可靠性非常重要。

他们的休息时间应该尽可能少。确保这一点的一个重要部分是在风力涡轮机中引入先进的故障检测、隔离和调节系统。在最先进的工业风力涡轮机故障检测和调节方案是简单的，大多数情况下是保守的。因此，使用先进的故障检测、隔离和调节方法可以提高涡轮机的可靠性，尽管对于某些故障可能会导致电力有限的生产。

或者，状态监测用于监测一些机械部件，如齿轮箱等，参见[Hameeda等人(2009)]。一些工作已经在风力涡轮机上进行了基于模型的故障检测、隔离和调节。[Wei等人(2008)]提出了一种基于观测器的俯仰系统传感器故障检测方案。[Dobrila and Stefansen(2007)]提出了一种基于宇称方程的风力发电机故障检测方案，[Odgaard et al(2009)]提出了一种未知输入观测器用于风力发电机传动系周围传感器故障的检测。

电气转换系统的故障检测可以在[Poure等人(2007)]中找到。

为了测试不同的检测、隔离和调节方案在风力机中的应用，本文提出了风力机系统级的基准模型，包括:传感器、执行器和系统故障。这一基准模型是基于一个现实的通用三叶片水平变速风力涡轮机与全变流器耦合。这种通用涡轮机的额定功率为4.8兆瓦。

在第2节中描述了风力涡轮机的功能，在第3节中介绍了故障场景，在第4节中介绍了风力涡轮机模型，第5节介绍了测试信号。

详细文章讲解见第4部分。

📚2 运行结果

 

 

 

部分代码：

% Fault models
% Fault 9
eta_dt2=.91

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=5;
P_r=4.8e6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sensors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m_vw=1.5;
sigma_vm=0.5;
m_omega_r=0;
sigma_omega_r=0.004*2*pi;
m_omega_g=0;
sigma_omega_g=0.05;
m_tau_g=0;
sigma_tau_g=90;
m_P_g=0;
sigma_P_g=1e3;
m_Beta=0;
sigma_Beta=0.2;

%% Fault models (choose the scenario of fault parameters and magnitude)

% %1st fault scenario
% 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_m2=0.9;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=100;eta_dt2=.92;Constant_Omega_g_m1=140;

% %2nd fault scenario
% Constant_Beta_1_m1=6;Gain_Beta_2_m2=1.5;Constant_Beta_3_m1=8;Constant_Omega_r_m1=0.2;Gain_Omega_r_m2=1.2;Gain_Omega_g_m2=0.8;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=100;eta_dt2=.71;Constant_Omega_g_m1=80;

%3rd fault scenario
Constant_Beta_1_m1=4;Gain_Beta_2_m2=1.8;Constant_Beta_3_m1=12;Constant_Omega_r_m1=1.2;Gain_Omega_r_m2=0.7;Gain_Omega_g_m2=1.7;
xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=700;eta_dt2=0.3;Constant_Omega_g_m1=100;

%4th fault scenario
% Constant_Beta_1_m1=2;Gain_Beta_2_m2=3;Constant_Beta_3_m1=15;Constant_Omega_r_m1=1.7;Gain_Omega_r_m2=1.7;Gain_Omega_g_m2=0.7;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=-500;eta_dt2=.6;Constant_Omega_g_m1=130;

% %5th fault scenario
% Constant_Beta_1_m1=-5;Gain_Beta_2_m2=10;Constant_Beta_3_m1=3;Constant_Omega_r_m1=2.5;Gain_Omega_r_m2=0.9;Gain_Omega_g_m2=0.2;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=1500;eta_dt2=.25;Constant_Omega_g_m1=120;
 
%6th fault scenario
% Constant_Beta_1_m1=1;Gain_Beta_2_m2=2;Constant_Beta_3_m1=6;Constant_Omega_r_m1=0.5;Gain_Omega_r_m2=0.2;Gain_Omega_g_m2=2.8;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=2500;eta_dt2=0.1;Constant_Omega_g_m1=150;

% % %7th fault scenario
% Constant_Beta_1_m1=24;Gain_Beta_2_m2=5;Constant_Beta_3_m1=20;Constant_Omega_r_m1=3.5;Gain_Omega_r_m2=3;Gain_Omega_g_m2=5;
% xi2=0.25;omega_n2=8.73;xi3=0.95;omega_n3=1.42;Constant_tau_gc=900;eta_dt2=0.4;Constant_Omega_g_m1=50;

%8th fault scenario
% Constant_Beta_1_m1=-3;Gain_Beta_2_m2=5;Constant_Beta_3_m1=7;Constant_Omega_r_m1=2;Gain_Omega_r_m2=0.5;Gain_Omega_g_m2=1.5;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=-1000;eta_dt2=0.22;Constant_Omega_g_m1=100;

% %9th fault scenario
% Constant_Beta_1_m1=2;Gain_Beta_2_m2=3;Constant_Beta_3_m1=7;Constant_Omega_r_m1=2;Gain_Omega_r_m2=0.5;Gain_Omega_g_m2=1.5;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=800;eta_dt2=0.42;Constant_Omega_g_m1=110;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[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]);
%
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];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fault signals (choose the scenario of fault appearance)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %1st fault scenario
% D_Fault1=[Time' 1-[zeros(1,2000/Ts) ones(1,100/Ts) zeros(1,2300/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_Fault10=[Time' 1-[zeros(1,3600/Ts) ones(1,100/Ts) zeros(1,700/Ts) ]'];

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1]Peter Fogh Odgaard, Jakob Stoustrup, and Michel Kinnaert. Fault tolerant control of wind turbines a benchmark model. In Proc.7th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes, pages 155–160, Barcelona, Spain, 2009.

[2]邓媛媛,刘国奇,蔡磊.基于设备运行数据的风机故障诊断方法[J].沈阳建筑大学学报(自然科学版),2021,37(05):942-948.

[3]刘瑶. 基于时频分析的风机故障智能诊断方法研究[D].浙江大学,2021.DOI:10.27461/d.cnki.gzjdx.2021.001016.

🌈4 Matlab代码、Simulink、文章讲解