【故障诊断分析】基于主成分分析实现三容水箱故障诊断附代码.zip

clear; clc; %%%%%%%%四种故障%%%%% t=1:1:3000 ; %设定时间基准 for i=1:size(t) %设定时间基准 if i>=1500; pianyi(i)=1; else pianyi(i)=0; end if i>=1500&i<=1700; pianyi2(i)=1; else pianyi2(i)=0; end if i>1700; pianyi3(i)=1; else pianyi3(i)=0; end if i<=1500; pianyi4(i)=1; else pianyi4(i)=0; end end %设定误差出现基准 wave=0.01*randn(size(t)); % matlab函数randn：产生正态分布的随机数或矩阵的函数 rootwave=7+wave; PCHA=0.5*pianyi; %设定误差类型 偏差故障 PYI=0.0025*pianyi.*(t-1500); %设定误差类型 漂移故障 JDU=2*pianyi.*normrnd(0,0.05,size(t)) ; %设定误差类型 精度下降故障 WQUAN=0.5*pianyi3+0.0025*pianyi2.*(t-1500)+pianyi4.*wave; %设定误差类型 完全故障 xmin=1450 ; xmax=1750; ymin=6.5; ymax=8; %%%%%%%%%%%水箱PID调节 kp1=30;ki1=0.2;kd1=5; kp2=25;ki2=0.2;kd2=5; h1set=25;h2set=20; c1=0.450289;c2=0.611429;c3=0.461526; noise1=rand(1)*0.01;noise2=rand(1)*0.01;noise3=rand(1)*0.01; %s = rand('state') %rand函数产生由在(0, 1)之间均匀分布的随机数组成的数组。 %Y = rand(n) 返回一个n x n的随机矩阵。如果n不是数量，则返回错误信息。 %Y = rand(m,n) 或 Y = rand([m n]) 返回一个m x n的随机矩阵。 %Y = rand(m,n,p,...) 或 Y = rand([m n p...]) 产生随机数组。 %Y = rand(size(A)) 返回一个和A有相同尺寸的随机矩阵。 g=981; Sp=0.5; u1=0;u2=0; A=154; h1=0;h2=0;h3=0; %%%% Q23min=-80.5;Q23max=80.5; Q13min=-60;Q13max=60; Q20min=0;Q20max=60; h1min=0;h1max=100; h2min=0;h2max=100; h3min=0;h3max=100; u2min=0;u2max=100; u1min=0;u1max=100; error1=0;error2=0; for k=1:3000 Q13(k)=0.5*c1*Sp*sign(h1(k)-h3(k))*sqrt(2*g*abs(h1(k)-h3(k))); Q23(k)=0.5*c3*Sp*sign(h3(k)-h2(k))*sqrt(2*g*abs(h3(k)-h2(k))); Q20(k)=0.5*c2*Sp*sqrt(2*g*h2(k)); Q11(k)=0.5*c1*Sp*sqrt(2*g*h1(k)); Q12(k)=0.5*c2*Sp*sqrt(2*g*h2(k)); if k<=2000 Q11(k)=0;%水箱1发生泄漏 end if k<=2500 Q12(k)=0;%水箱2发生泄漏 end if mod(k,10)==0 error1(k)=h1set-h1(k); else if k==1 error1(k)=error1; else error1(k)=error1(k-1); end end if k==1 u1(k)=kp1*error1+ki1*error1*ones(k,1); else u1(k)=kp1*error1(k)+ki1*error1*ones(k,1)+kd1*(error1(k)-error1(k-1)); %%%u1(k)=kp1*error1(k)+ki1*error1*ones(k,1)+kd1*(error1(k)-error1(k-1))+10*PCHA(k); %u1(k)=kp1*error1(k)+ki1*error1*ones(k,1)+kd1*(error1(k)-error1(k-1))+5*PYI(k); %u1(k)=kp1*error1(k)+ki1*error1*ones(k,1)+kd1*(error1(k)-error1(k-1))+10*JDU(k); % u1(k)=kp1*error1(k)+ki1*error1*ones(k,1)+kd1*(error1(k)-error1(k-1))+10*WQUAN(k); size(PCHA) end if mod(k,5)==0 error2(k)=h2set-h2(k); else if k==1 error2(k)=error2; else error2(k)=error2(k-1); end end if k==1 u2(k)=kp2*error2+ki2*error2*ones(k,1); else u2(k)=kp2*error2(k)+ki2*error2*ones(k,1)+kd2*(error2(k)-error2(k-1)); %u2(k)=kp2*error2(k)+ki2*error2*ones(k,1)+kd2*(error2(k)-error2(k-1))+10*PCHA(k); %PCHA(k) %PYI %JDU %WQUAN end h1(k+1)= h1(k)+0.1*((u1(k)-Q13(k)-Q11(k))/A+noise1); %h1(k+1)= h1(k)+0.1*((u1(k)-Q13(k))/A+noise1); h3(k+1)= h3(k)+0.1*((Q13(k)-Q23(k))/A+noise3); h2(k+1)= h2(k)+0.1*((u2(k)+Q23(k)-Q20(k)-Q12(k))/A+noise2); %h2(k+1)= h2(k)+0.1*((u2(k)+Q23(k)-Q20(k)-Q12(k))/A+noise2); %h1(k+1)=a1(k); %h2(k+1)=a2(k); %h3(k+1)=a3(k); if h1(k)>h1max; h1(k)=h1max; %%%%%%%防止h1,h2,h3饱和输出，满液位后水自动溢出排掉 end %%%%%%%防止h1,h2,h3饱和输出，满液位后水自动溢出排掉 %%%%%%%防止h1,h2,h3饱和输出，满液位后水自动溢出排掉 if h3(k)>h3max; h3(k)=h3max; end if h2(k)>h2max; h2(k)=h2max; end if Q20(k)>Q20max; Q20(k)=Q20max; end if Q13(k)>Q13max; Q13(k)=Q13max; end if Q13(k)<Q13min; Q13(k)=Q13min; end if Q23(k)>Q23max; Q23(k)=Q23max; end if Q23(k)<Q23min; Q23(k)=Q23min; end if u1(k)<u1min; u1(k)=u1min; end if u1(k)>u1max; u1(k)=u1max; end if u2(k)<u2min; u2(k)=u2min; end if u2(k)>u2max; u2(k)=u2max; end end %%%%%%PID调节结束%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%液位高度变化图像%%%%%%%%%%%%% figure(1) subplot(1,3,1) hold on plot(h1,'r'); %axis([0 time/10 0 hmax ]) x1=xlabel('执行次数'); x2=ylabel('液位高速'); %set(x1,'Rotation',0); %set(x2,'Rotation',0); title('液位 H1'); subplot(1,3,2) plot(h3,'r') %axis([0 time/10 0 hmax ]) x1=xlabel('执行次数'); x2=ylabel('液位高速'); %set(x1,'Rotation',0); %set(x2,'Rotation',0); title('液位 H3'); subplot(1,3,3) hold on plot(h2,'r'); %axis([0 time/10 0 hmax ]) x1=xlabel('执行次数'); x2=ylabel('液位高速'); %set(x1,'Rotation',0); %set(x2,'Rotation',0); title('液位 H2'); %%%%%%%%%%%%%%%%%%%%%%%%Q13，Q23，Q20图像%%%%%%%%%%%%%%%%%%%%%%% figure(2) subplot(1,3,1) plot(Q13,'r'); title('Q13'); subplot(1,3,2) plot(Q23,'r'); title('Q23'); subplot(1,3,3) plot(Q20,'r'); title('Q20'); %%%%%%%%%%%%画 出 阀 门 U1 和 U2 开 度 的 图 形%%%%%%%%%%%%%%%%%%% figure(3) subplot(1,2,1) plot(u1,'r'); title('阀门U1流量'); subplot(1,2,2) plot(u2,'r'); title('阀门U2流量'); %%%%%%%%%%%%画 出 阀 门 U1 和 U2 开 度 的 图 形%%%%%%%%%%%%%%%%%%% H1=h1'; H2=h2'; H3=h3'; q13=Q13'; q23=Q23'; q20=Q20'; xlswrite('故障数据.xls',H1,'sheet1','A') xlswrite('故障数据.xls',H2,'sheet1','B') xlswrite('故障数据.xls',H3,'sheet1','C') xlswrite('故障数据.xls',q13,'sheet1','D') xlswrite('故障数据.xls',q23,'sheet1','E') xlswrite('故障数据.xls',q20,'sheet1','F')