PQ分解法程序.docVIP

PQ分解法程序 [x]=xlsread(潮流参数.xls,A2:A2); [y]=xlsread(潮流参数.xls,B2:B2); e=xlsread(潮流参数.xls,B7:B7); [point]=xlsread(潮流参数.xls,D3:H6); [zhilu]=xlsread(潮流参数.xls,J3:R6); TYPE=zeros(x,1); U=zeros(x,1); a=zeros(x,1); P=zeros(x,1); Q=zeros(x,1); I=zeros(y,1); J=zeros(y,1); Rij=zeros(y,1); Xij=zeros(y,1); Zij=Rij+j*Xij; Y=zeros(x); G=zeros(x); B=zeros(x); B0=zeros(y,1); RT=zeros(y,1); XT=zeros(y,1); ZT=RT+j*XT; KT=zeros(y,1); %------------------------------矩阵赋初值： TYPE=point(:,1); U=point(:,2); a=point(:,3); P=point(:,4); Q=point(:,5); I=zhilu(:,1); J=zhilu(:,2); Rij=zhilu(:,3); Xij=zhilu(:,4); Zij=Rij+j*Xij; B0=zhilu(:,5); RT=zhilu(:,6); XT=zhilu(:,7); ZT=RT+j*XT; KT=zhilu(:,8); W=zhilu(:,9); %------------------------------求节点导纳矩阵Y for m=1:y if KT(m)==0 Y(I(m),J(m))=-1/Zij(m); Y(J(m),I(m))=-1/Zij(m); else Y(I(m),J(m))=-1/(KT(m)*ZT(m)); Y(J(m),I(m))=-1/(KT(m)*ZT(m)); end end for m=1:x for n=1:y if KT(n)==0 if(I(n)==m|J(n)==m) Y(m,m)=Y(m,m)-Y(I(n),J(n))+j*B0(n)/2; end else if W==0 if I(n)==m Y(m,m)=Y(m,m)-Y(I(n),J(n))+(KT(n)-1)/KT(n)*(1/ZT(n)); elseif J(n)==m Y(m,m)=Y(m,m)-Y(I(n),J(n))+(1-KT(n))/(KT(n)^2)*(1/ZT(n)); end else if I(n)==m Y(m,m)=Y(m,m)-Y(I(n),J(n))+(1-KT(n))/(KT(n)^2)*(1/ZT(n)); elseif J(n)==m Y(m,m)=Y(m,m)-Y(I(n),J(n))+(KT(n)-1)/KT(n)*(1/ZT(n)); else Y(m,m)=Y(m,m); end end end end end G=real(Y); %-----------------------求B矩阵及其逆矩阵B1 B=imag(Y); ph=find(TYPE(:,1)==3); BB=B; BB(:,ph)=[]; BB(ph,:)=[]; B1=BB; B1=inv(B1); disp(B矩阵) disp(BB) disp(B逆矩阵) disp(B1) %-----------------------%求B及其逆矩阵B2 phpv=find(TYPE(:,1)1); BB=B; BB(:,phpv)=[]; BB(phpv,:)=[]; B2=BB; B2=inv(B2); disp(BB) disp(B2) %-------------计算各节点有功功率不平衡量deltaPi k=0; kp=1; kq=1; while(((kp~=0)||(kq~=0))k=20) kp=1; kq=1; notph=find(TYP