Big_homework_of_power_system.zip_Big!_Windows编程_coldpjz_电力系统解析

A1=importdata('power.txt');%读取节点信息 B1=importdata('branch.txt');%读取支路信息 A2=A1.data; A=A2(:,[3,6,7,8,9,11:13]);%取节点信息中有用的列：节点号、节点类型等 B=B1(:,[1,2,7,8,9,15]);%取支路信息中有用的列：首末端节点、电阻等 A=A([16:57,1:15],:);%排列A矩阵 for ei=1:size(B,1) for(j=1:2) if(B(ei,j))>=16 B(ei,j)=B(ei,j)-15; else if (B(ei,j)<16) B(ei,j)=B(ei,j)+42; end end end end%排列B矩阵，将原本16号节点设置为1号 nl0 = B(:,1);%提取支路首节点 nr0 = B(:,2);%提取支路末节点 R0= B(:,3);%提取支路电阻 X0 = B(:,4);%提取支路电抗 Bc10 = 1i*B(:,5);%提取支路左侧对地电纳 GB10 =Bc10;%得到支路左侧的导纳 Z0=R0+1i*X0; %得到支路的阻抗 a0=B(:,6);%提取支路的变比，首节点为高压侧，模型k:1，线路按照1：1考虑 nbr0=length(B(:,1));%得到系统支路的数量 nbus0=max(max(nl0),max(nr0));%得到系统的节点数 y0=ones(nbr0,1)./Z0;%生成支路导纳矩阵 for n=1:nbr0 if a0(n)<=0 a0(n)=1;%如果支路上的变比误填入，采取纠错措施将其改为1：！ end; end; Ybus0=zeros(nbus0,nbus0);% 对节点导纳矩阵进行赋初值为 %支路互导纳 for n=1:80 Ybus0(nl0(n),nr0(n))=Ybus0(nl0(n),nr0(n))-y0(n)/a0(n); Ybus0(nr0(n),nl0(n))=Ybus0(nl0(n),nr0(n)); end; %自导纳 for k=1:nbr0 firstnode= nl0(k); secondnode= nr0(k); if a0(k)~=1 Ybus0(firstnode,firstnode) = Ybus0(firstnode,firstnode)+y0(k)/(a0(k).^2); Ybus0(secondnode,secondnode) = Ybus0(secondnode,secondnode)+ y0(k)+ GB10(k); else Ybus0(firstnode,firstnode) = Ybus0(firstnode,firstnode)+y0(k)/(a0(k).^2)+ GB10(k)/2; Ybus0(secondnode,secondnode) = Ybus0(secondnode,secondnode)+ y0(k)+ GB10(k)/2; end end; %支路节点接地导纳 % for m=1:nbr0 % Ybus0(nl0(m),nl0(m))=Ybus0(nl0(m),nl0(m))+y0(m)/a0(m);%变压器支路的阻抗标幺值视为已经折到左侧（BPA的dat文件中） % % Ybus0(nr0(m),:)=0; % % Ybus0(:,nr0(m))=0; % end; nodey0=Ybus0;%节点导纳矩阵存储于nodey0 swing=find(A(:,1)==3); NPV=find(A(:,1)==2); NPQ=find(A(:,1)==0); kNPQ=length(NPQ); kNPV=length(NPV); k=length(A); P0=zeros(k-1,1); Q0=zeros(kNPQ,1); for n=1:kNPQ P0(n,1)=A(NPQ(n),4)-A(NPQ(n),2); Q0(n,1)=A(NPQ(n),5)-A(NPQ(n),3); end for n=(kNPQ+1):(kNPV+kNPQ) P0(n,1)=A(NPV(n-kNPQ),4)-A(NPV(n-kNPQ),2); end %按照PQ节点有功在前，PV节点有功在后的顺序进行排序 P0=P0/100;%标幺值取100 Q0=Q0/100; for n=1:kNPV Uen(n,1)=A(NPV(n),6); end%取PV节点电压期望值 UPQ=ones(kNPQ,1); Ue=[UPQ;Uen]; e0=ones(nbus0-1,1); f0=zeros(nbus0-1,1); U=e0+1i*f0; [nodeychange1,nodeychange2]=Nodeychange(A,nodey0);%PQ节点在前，PV节点在后排列Y矩阵 U1=[Ue;1.04]; dete=0; detf=0; times=0;%计算迭代次数 error1=1; while error1>0.001%设置误差限位0.001 times=times+1;%迭代次数加一 %求PQU的差值 S1=diag(U1)*conj(nodeychange2)*conj(U1); S1(k,:)=[]; P1=real(S1); Q1=imag(S1); detP=(P0-P1); Q1(kNPQ+1:(k-1),:)=[]; detQ=(Q0-Q1); detU=Ue.*conj(Ue)-U.*conj(U); detU(1:kNPQ,:)=[]; %求PQU的差值 %求取雅克比矩阵 I=nodeychange2*U1; I(k,:)=[]; a=real(I); b=imag(I); H=imag(conj(diag(conj(U))*nodeychange1))+imag(diag(I)); N=real(conj(diag(conj(U))*nodeychange1))+real(diag(I)); J=-1*real(conj(diag(conj(U))*nodeychange1))+real(diag(I)); L=imag(conj(diag(conj(U))*nodeychange1))-imag(diag(I)); J=J([1:kNPQ],:); L=L([1:kNPQ],:); R=2*imag(diag(U)); S=2*real(diag(U)); R=R([kNPQ+1:k-1],:); S=S([kNPQ+1:k-1],:); YKB=[H,N;J,L;R,S]; %求取雅克比矩阵 XXX=[detP;detQ;detU]; det=YKB\[detP;detQ;detU]; error1=max(det); detf=det(1:56,:); det(1:56,:)=[]; dete=det; U=U+dete+1i*detf; U1=[U;1.04]; end S1=diag(U1)*conj(nodeychange2)*conj(U1); S1(k,:)=[]; P1=real(S1); Q1=imag(S1); S1=S1*100; for i=1:57%将直角坐标系转化为极坐标系 [U11(i,2),U11(i,1)]=cart2pol(real(U1(i,1)),imag(U1(i,1))); end U11(:,2)=U11(:,2)*180/pi;%弧度转化为角度