fangzhen.m

能够利用偏最小二乘回归模型分析电力系统谐波发射水平

load upcc3 ;
load ipcc3;%%%从模型中的倒出数据PCC点
%在数据段分段成1000段变成1000*1000矩阵
for t=1:1000
    c(:,t)=upcc((1000*t-999+50000):(1000*t+50000),2);
    d(:,t)=ipcc((1000*t-999+50000):(1000*t+50000),2);
end
%%1000个点做一次傅立叶变化得到各次谐波的幅值和相位信息；
q=fft(c)/500;w=fft(d)/500;%注意：在做变换时候要除以采样点个数的一般从而得到原始的幅值和相位信息！！！
q=q(4,:);w=w(4,:);%提取3次谐波的复数
u1=(abs(q))';
u2=(angle(q))';
i1=(abs(w))';
i2=(angle(w))';
% ux=u1.*cos(u2);
%uy=u1.*sin(u2);
 ux=u1.*cos(u2);
uy=u1.*sin(u2);
%U=complex(ux,uy);
%i3=ones(1440,1)*0; 
%ix=i1.*cos(i2);
%iy=i1.*sin(i2);
ix=i1.*cos(i2);
iy=i1.*sin(i2);
%I=complex(ix,iy);
%y=ux.*iy+ix.*uy;
%x=[iy,ix,-(ix.*ix+iy.*iy)];
%{采用模型Vpccx==vsx+zsy*ipccy-zsxipccx}
   x1_60=[iy,-ix];
   y1_60=ux;
   m=20;%(从1到20每一个一次，可以算的个值，最后取平均就可以)
     x1_60=x1_60((50*m-49):(50*m),:);
     y1_60=y1_60((50*m-49):(50*m),:);
 %对50个点进行标准化处理  
 for i=1:50;
     for j=1:2
         E(i,j)=[x1_60(i,j)-mean(x1_60(:,j))]/sqrt(var(x1_60(:,j),1));
         F(i,1)=[y1_60(i,1)-mean(y1_60(:,1))]/sqrt(var(y1_60(:,1),1));
     end 
  
 end
 %执行偏最小二乘算法  
 
 for n=1:2;%% n是根据我具体的残差矩阵来判断循环个数的
    
               Q=sqrt((F.')*E*(E.')*F);
               W1=((E.')*F)/Q;%求第一个特征向量
        %[W1,M1,L1]=RegrPLS(E,F,2)
           T1=E*W1;%第一个主成分
           Q1=(T1.')*T1;
           P1=((E.')*T1)./Q1;
           E1=E-T1*(P1.');%第一个残差矩阵
           T(:,n)=T1;
           W2(:,n)=W1;
           P2(:,n)=P1;
           E=E1;
 end

 %求取F对E的回归系数
   W3(:,1)=W2(:,1);
   W3(:,2)=(eye(2)-W2(:,1)*(P2(:,1).'))*W2(:,2);
   %W3(:,3)=(eye(2)-W2(:,1)*(P2(:,1).'))*(eye(2)-W2(:,2)*(P2(:,2).'))*W2(:,3)
   %  W3(:,m)=c\ T(:,m);
     

   R=regress(F,T);
   A=R(1,1)*W3(:,1)+R(2,1)*W3(:,2);%+R(3,1)*W3(:,3);
    % A(1,1)=R(1,1)*W3(1,1)+R(2,1)*W3(2,1)+R(3,1)*W3(3,1);
    % A(1,2)=R(1,1)*W3(1,2)+R(2,1)*W3(2,2)+R(3,1)*W3(3,2);
   % A(1,3)=R(1,1)*W3(1,3)+R(2,1)*W3(2,3)+R(3,1)*W3(3,3);
%%%还原y对x的回归系数
   zy(m,1)=sqrt(var(y1_60(:,1),1))* A(1,1)/sqrt(var(x1_60(:,1),1));
   zx(m,1)=sqrt(var(y1_60(:,1),1))* A(2,1)/sqrt(var(x1_60(:,2),1));  
   %B(1,3)=sqrt(var(y(:,1),1))* A(3,1)/sqrt(var(x(:,3),1));  
 
   %for n=1:500;
   %   a(n,:)=i(2,(256*n-255):256*n);
   %end
   Vsx(m,1)=[(A(1,1)*mean(x1_60(:,1))/sqrt(var(x1_60(:,1),1))+A(2,1)*mean(x1_60(:,2))/sqrt(var(x1_60(:,2),1)))]*sqrt(var(y1_60(:,1),1));%+A(3,1)*mean(x1_60(:,3))/sqrt(var(x1_60(:,3),1)))]*sqrt(var(y1_60(:,1),1));
   Vsx(m,1) =mean(y1_60)-Vsx(m,1);

   
end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%五次谐波仿真
load upcc5 ;
load ipcc5;

for t=1:1000
    c(:,t)=u((1000*t-999+50000):(1000*t+50000),2);
    d(:,t)=i((1000*t-999+50000):(1000*t+50000),2);
end
q=fft(c)/500;w=fft(d)/500;
q=q(6,:);w=w(6,:);
u1=(abs(q))';
u2=(angle(q))';
i1=(abs(w))';
i2=(angle(w))';
% ux=u1.*cos(u2);
%uy=u1.*sin(u2);
 ux=u1.*cos(u2);
uy=u1.*sin(u2);
%U=complex(ux,uy);
%i3=ones(1440,1)*0; 
%ix=i1.*cos(i2);
%iy=i1.*sin(i2);
ix=i1.*cos(i2);
iy=i1.*sin(i2);
 x1_60=[iy,-ix];
   y1_60=ux;
x1_60=[iy,-ix];
   y1_60=ux;
   m=20;%(从1到20每一个一次，可以算的个值，最后取平均就可以)
     x1_60=x1_60((50*m-49):(50*m),:);
     y1_60=y1_60((50*m-49):(50*m),:);
   
 for i=1:50;
     for j=1:2
         E(i,j)=[x1_60(i,j)-mean(x1_60(:,j))]/sqrt(var(x1_60(:,j),1));
         F(i,1)=[y1_60(i,1)-mean(y1_60(:,1))]/sqrt(var(y1_60(:,1),1));
     end 
  
 end
        for n=1:2;
    
               Q=sqrt((F.')*E*(E.')*F);
                  W1=((E.')*F)/Q;
        %[W1,M1,L1]=RegrPLS(E,F,2)
           T1=E*W1;Q1=(T1.')*T1;
           P1=((E.')*T1)./Q1;
           E1=E-T1*(P1.');
            T(:,n)=T1;
           W2(:,n)=W1;
            P2(:,n)=P1;
            E=E1;
        end

 
   W3(:,1)=W2(:,1);
   W3(:,2)=(eye(2)-W2(:,1)*(P2(:,1).'))*W2(:,2);
   %W3(:,3)=(eye(2)-W2(:,1)*(P2(:,1).'))*(eye(2)-W2(:,2)*(P2(:,2).'))*W2(:,3)
   %  W3(:,m)=c\ T(:,m);
     

   R=regress(F,T);
   A=R(1,1)*W3(:,1)+R(2,1)*W3(:,2);%+R(3,1)*W3(:,3);
    % A(1,1)=R(1,1)*W3(1,1)+R(2,1)*W3(2,1)+R(3,1)*W3(3,1);
    % A(1,2)=R(1,1)*W3(1,2)+R(2,1)*W3(2,2)+R(3,1)*W3(3,2);
   % A(1,3)=R(1,1)*W3(1,3)+R(2,1)*W3(2,3)+R(3,1)*W3(3,3);

   zy(m,1)=sqrt(var(y1_60(:,1),1))* A(1,1)/sqrt(var(x1_60(:,1),1));
   zx(m,1)=sqrt(var(y1_60(:,1),1))* A(2,1)/sqrt(var(x1_60(:,2),1));  
   %B(1,3)=sqrt(var(y(:,1),1))* A(3,1)/sqrt(var(x(:,3),1));  
 
   %for n=1:500;
   %   a(n,:)=i(2,(256*n-255):256*n);
   %end
   Vsx(m,1)=[(A(1,1)*mean(x1_60(:,1))/sqrt(var(x1_60(:,1),1))+A(2,1)*mean(x1_60(:,2))/sqrt(var(x1_60(:,2),1)))]*sqrt(var(y1_60(:,1),1));%+A(3,1)*mean(x1_60(:,3))/sqrt(var(x1_60(:,3),1)))]*sqrt(var(y1_60(:,1),1));
   Vsx(m,1) =mean(y1_60)-Vsx(m,1);

   
end