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;%弧度转化为角度
