配电网故障诊断（Matlab代码实现）

close all clear all clc % Loading the data load('sensor'); %% Q.1 ACP MODEL %1. The temporal visualation of some pollutants %==> Creating time Interval Sampling_Time = 15; %[Minutes] Size_Samples = length(X1(1,:)); T_min=0:Sampling_Time:(Size_Samples-1)*Sampling_Time; T_hours = T_min/60 ; % 1hour = 60 minutes %==> The plots figure() plot(T_hours,X1(1,:),'r',T_hours,X1(2,:),'b',T_hours,X1(3,:),'g'); xlabel('Time_{[Hour]}'); set(gca,'FontSize',16); ylabel('Sensors Measurements'); set(gca,'FontSize',16); legend('Ozon Measurements','Nitrogen Oxyde Measurements','Dioxyde Measurements'); set(gca,'FontSize',16); title('The Temporal Evoluation Of Pollutants in Location 1'); set(gca,'FontSize',16); % figure() % plot(T_min,X1(1,:),'r',T_min,X1(2,:),'b',T_min,X1(3,:),'g'); % xlabel('Time_{[Minute]}'); % ylabel('Sensors Measurements'); % legend('Ozon Measurements','Nitrogen Oxyde Measurements','Dioxyde Measurements'); % title('The Temporal Evoluation Of Pollutants'); l1=1:3:18; l2=2:3:18; l3=3:3:18; figure () plot(1:length(X1),X1(l1,:)); %Ozon measurements across 6 locations xlabel('Time_{[Hour]}'); set(gca,'FontSize',16); ylabel('Sensors Measurements'); set(gca,'FontSize',16); legend('Location 1','Location 2','Location 3','Location 4','Location 5','Location 6'); set(gca,'FontSize',16); title('The Temporal Evoluation Of Ozon across 6 locations'); set(gca,'FontSize',16); figure () plot(1:length(X1),X1(l2,:)); %Nitrogen oxyde measurements across 6 locations xlabel('Time_{[Hour]}'); set(gca,'FontSize',16); ylabel('Sensors Measurements'); set(gca,'FontSize',16); legend('Location 1','Location 2','Location 3','Location 4','Location 5','Location 6'); set(gca,'FontSize',16); title('The Temporal Evoluation Of Nitrogen Oxyde across 6 locations'); set(gca,'FontSize',16); figure () plot(1:length(X1),X1(l3,:)); %Dioxyde measurements across 6 locations xlabel('Time_{[Hour]}'); set(gca,'FontSize',16); ylabel('Sensors Measurements'); set(gca,'FontSize',16); legend('Location 1','Location 2','Location 3','Location 4','Location 5','Location 6'); set(gca,'FontSize',16); title('The Temporal Evoluation Of Dioxyde across 6 locations'); set(gca,'FontSize',16); %2. Centering And Normalizing The Data Mean_X1 = mean(X1')'; Std_X1 = std(X1')'; for i=1:Size_Samples Norm_X1(:,i) = (X1(:,i)-Mean_X1)./(Std_X1); end %3. Computing The Covaraince Of X1 centered and Normalized Cov_X1 = (Norm_X1*Norm_X1')/(length(Norm_X1(1,:)-1)); %4. The EigenVectors And The EigenValue of The Cov Matrix [P,lambda,P] = svd(Cov_X1); %5. Plotting The cumulative Sum Of The Eigen Values figure() imagesc(diag(Cov_X1)) colormap hsv colorbar xlabel('Proportion Of Variance Matrix Diag'); set(gca,'FontSize',16); ylabel('Components(Axis)'); set(gca,'FontSize',16); title('Proportion of Diag Matrix varianceFor Each Component'); set(gca,'FontSize',16); figure() bar(cumsum(diag(lambda))/length(Norm_X1(:,1))) xlabel('Axis'); set(gca,'FontSize',16); ylabel('Proportion of variance expressed'); set(gca,'FontSize',16); title('Proportion of variance expressed by each axis'); set(gca,'FontSize',16); %6. The Projection Of Each Vector In X1 on 12 Components Frame New_P = P(:,1:12); %P(q) Pr_Comp_X1 = New_P'*Norm_X1; %12 Principas components C(q) Mean_P = mean(Pr_Comp_X1'); %Mean of C(q) Cov_P = cov(Pr_Comp_X1'); %Covariance of C(q) M_dist = (Pr_Comp_X1)'*inv(lambda(1:12,1:12))*(Pr_Comp_X1); %Mahalanobis distance M_dist_mean = mean(diag(M_dist)); M_dist_std = std(diag(M_dist)); % Ma_dist = sqrt((Pr_Comp_X1-Mean_P)*inv(Cov_P)*(Pr_Comp_X1-Mean_P)'); %figure() %plot(T_hours,diag(M_dist)); %7. The Estimation of X1,X^1 Estm_X1 = New_P * Pr_Comp_X1; Eucl_dist = (Estm_X1-Norm_X1)'*(Estm_X1-Norm_X1); Eucl_dist_mean = mean(diag(Eucl_dist)); Eucl_dist_std = std(diag(Eucl_dist)); figure() subplot(3,1,1); plot(T_hours,Estm_X1(1,:),T_hours,Norm_X1(1,:)); xlabel('Time (hours)'); %set(gca,'FontSize',16); ylabel('O3 concentration'); %set(gca,'FontSize',16); title('Reconstruction and true measurements of ozon concentration in location 1'); %set(gca,'FontSize',16); legend('Estimated O3 concentration','True O3 concentration'); %set(gca,'FontSize',16); subplot(3,1,2); plot(T_hours,Estm_X1(2,:),T_hours,Norm_X1(2,:)); xlabel('Time (hours)'); %set(gca,'FontSize',16); ylabel('NO2 concentration'); %set(gca,'FontSize',16); title('Reconstruction and true measurements of NO2 concentration in location 1'); %set(gca,'FontSize',16); legend('Estimated NO2 concentration','True NO2 concentration'); %set(gca,'FontSize',16); subplot(3,1,3); plot(T_hours,Estm_X1(3,:),T_hours,Norm_X1(3,:)); xlabel('Time (hours)'); %set(gca,'FontSize',16); ylabel('CO2 concentration'); %set(gca,'FontSize',16); title('Reconstruction and true measurements of CO2 concentration in location 1'); %set(gca,'FontSize',16); legend('Estimated CO2 concentration','True CO2 concentration'); %set(gca,'FontSize',16); %8. Square Norm (Squared Prediction Error) Q = (Norm_X1-Estm_X1)'*(Norm_X1-Estm_X1); %%Part.2 Fault Detection %Q.1 Applying the ACP model On X2 %% Creating The Time Interval Size_Samples_X2 = length(X2(1,:)); T_min_X2=0:Sampling_Time:(Size_Samples_X2-1)*Sampling_Time; T_hours_X2 = T_min_X2/60 ; % 1hour = 60 minutes for i=1:length(T_hours_X2) Norm_X2(:,i) = (X2(:,i)-Mean_X1)./Std_X1; Pr_Comp_X2(:,i) = New_P'*Norm_X2(:,i); %% Malalobinos Distance Checking M_dist_X2(i) = (Pr_Comp_X2(:,i))'*inv(lambda(1:12,1:12))*(Pr_Comp_X2(:,i)); if M_dist_X2(i)> 12 X2_Decision(i) = 1; else X2_Decision(i) = 0; end %Decision based onquared estimation error Estm_X2(:,i) = New_P * Pr_Comp_X2(:,i); %Slide27 %%Eucliedien Checking Eucludien_Dist2(i) = (Estm_X2(:,i)-Norm_X2(:,i))'*(Estm_X2(:,i)-Norm_X2(:,i)); if Eucludien_Dist2(i)> 1 Eu_decision(i) = 1; else Eu_decision(i) = 0; end Q_X2(i)=(Estm_X2(:,i)-Norm_X2(:,i))'*(Estm_X2(:,i)-Norm_X2(:,i)); if Q_X2(i)> (0.99*Mean_X1/(2*Std_X1)) X2_Decision