function[tab2,tab4,tab6,cycle]=newton(filename)
base=100; %容量基值
b=textread('case30.txt'); %读文件
[m1,n1]=size(b);
n=0; %节点数
m=0; %支路数
for k=1:400 %最大节点数400
if b(k,1)~=-999
n=n+1;
else
break;
end
end
c=k+1; %指向第一个-999的后一行数
for k=c:1000 %最大支路数
if b(k,1)~=-999
m=m+1;
else
break
end
end
p=zeros(n,n1+1);%节点矩阵
z=zeros(m,n1+2);%支路矩阵
p(1:n,1:n1)=b(1:n,1:n1);
z(1:m,1:n1)=b(n+2:m1-1,1:n1);
%调整为pq,pv,平衡节点,
x1=0;
for i=1:n
if p(i,5)==1
x1=x1+1;
p(i,n1+1)=x1;
end
end
x2=x1;
for i=1:n
if p(i,5)==2
x2=x2+1;
p(i,n1+1)=x2;
end
end
x3=x2;
for i=1:n
if p(i,5)==3
x3=x3+1;
p(i,n1+1)=x3;
end
end
for i=1:m
for j=1:n
if z(i,1)==p(j,1)
z(i,22)=p(j,22);break;
end
end
for j=1:n
if z(i,2)==p(j,1)
z(i,23)=p(j,22);break;
end
end
end
for k=1:(n-1)%冒泡法根据节点类型把行换顺序
for j=1:(n-k)
if p(j+1,22)<p(j,22)
i=p(j+1,:);
p(j+1,:)=p(j,:);
p(j,:)=i;
end
end
end
%%%%%%%%%%%%%形成节点导纳矩阵Y
B=zeros(n,n);
for i=1:m
z(i,6)=z(i,6);%改变电阻
z(i,7)=z(i,7);%改变电抗
z(i,8)=171*z(i,8);%改变电纳
end
for i=1:m
j=z(i,22);
k=z(i,23);
if z(i,5)==1 %%%%%支路为线路
B(j,j)=B(j,j)+1i*z(i,8)/2+1/(z(i,6)+1i*z(i,7));
B(k,k)=B(k,k)+1i*z(i,8)/2+1/(z(i,6)+1i*z(i,7));
B(j,k)=-1/(z(i,6)+1i*z(i,7));
B(k,j)=-1/(z(i,6)+1i*z(i,7));
else if z(i,5)==2 %%支路为变压器 变压器的形式:i k:1 Zt j 的形式
B(k,k)=B(k,k)+1i*z(i,8);
B(k,k)=B(k,k)+(z(i,14)-1)/(z(i,14)*(z(i,6)+1i*z(i,7)))+1/(z(i,14)*(z(i,6)+1i*z(i,7)));
B(j,j)=B(j,j)+(1-z(i,14))/((z(i,14))^2*(z(i,6)+1i*z(i,7)))+1/(z(i,14)*(z(i,6)+1i*z(i,7)));
B(j,k)=-1/(z(i,14)*(z(i,6)+1i*z(i,7)));
B(k,j)=-1/(z(i,14)*(z(i,6)+1i*z(i,7)));
end
end
end
for k=1:n
B(k,k)=B(k,k)+1i*p(k,17); %加上对地电容
end
Y=B;
G=real(B);
B=imag(B);
q1=0;%pq节点个数
v1=0;%pv节点个数
for k=1:n
if p(k,5)==2
v1=v1+1;
end
end
q1=n-v1-1;
P2=zeros(n-1,1);%真实的P
Q2=zeros(q1,1);
Q3=zeros(n-1,1);
P21=zeros(n-1,1);%P差值
Q21=zeros(q1,1);
U1=zeros(n-1,1);%计算的U
A1=zeros(n-1,1);
U2=zeros(n-1,1);%之前的U
A=zeros(n-1,1);
U21=zeros(q1,1);
A21=zeros(n-1,1);
solve=zeros(n+q1-1,1);
PQ=zeros(n+q1-1,1);
Jac=zeros(n+q1-1,n+q1-1);
U=p(:,13);
A=p(:,7)./(180/pi);
P2=(p(1:n-1,10)-p(1:n-1,8))./base;
Q3=(p(:,11)-p(:,9))./base;
Q2=Q3(1:q1,1);
maxError=0;
Node=0;
for cycle=1:2000
P21=P2;
Q21=Q2;
H=zeros(n-1,n-1);
L=zeros(q1,q1);
M=zeros(q1,n-1);
N=zeros(n-1,q1);
for i=1:n-1
if(i<=q1)
N(i,i)=-2*U(i,1)^2*G(i,i);
L(i,i)=2*U(i,1)^2*B(i,i);
end
for j=1:n
P21(i,1)=P21(i,1)-U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1))));
if(i<=q1)
Q21(i,1)=Q21(i,1)-U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1))));
end
if(i~=j)
H(i,i)=H(i,i)+U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1))));
if(i<=q1)
M(i,i)=M(i,i)-U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1))));
N(i,i)=N(i,i)-U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1))));
L(i,i)=L(i,i)-U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1))));
end
if(j<n)
H(i,j)=-U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1))));
end
if(j<=q1)
N(i,j)=-U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1))));
end
if(i<=q1)&(j<n)
M(i,j)=U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1))));
end
if((i<=q1)&(j<=q1))
L(i,j)=-U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1))));
end
end
end
if P21(i,1)>maxError
maxError=P21(i,1);
Node=i;
end
if(i<=q1)&(Q21>maxError)
maxError=Q21(i,1);
Node=i;
end
end
if(abs(all(P21(:))))<0.0001&(abs(all(Q21(:))))<0.0001
disp('success');
break;
end
Jac=[H N;M L];
PQ=[P21;Q21];
solve=-1*inv(Jac)*PQ;
for k=1:n+q1-1
if (k<n)
A21(k,1)=solve(k,1);
A(k,1)=A(k,1)+A21(k,1);
else
U21(k-n+1,1)=solve(k,1)*U(k-n+1,1);
U(k-n+1,1)=U21(k-n+1,1)+U(k-n+1,1);
end
end
end
P=zeros(n,1);
P(1:n-1,1)=P2.*base;
Q=zeros(n,1);
Q(1:q1,1)=Q2.*base;
for i=q1+1:n
for j=1:n
Q(i,1)=Q(i,1)+base*U(i,1)*U(j,1)*(G(i,j)*sin( (A(i,1)-A(j,1)))-B(i,j)*cos( (A(i,1)-A(j,1))));
if(i==n)
P(i,1)=P(i,1)+base*U(i,1)*U(j,1)*(G(i,j)*cos( (A(i,1)-A(j,1)))+B(i,j)*sin( (A(i,1)-A(j,1))));
end
end
end
for k=1:(n-1)%把各参数数据顺序调回来
for j=1:(n-k)
if p(j+1,1)<p(j,1)
bj=p(j+1,:);
p(j+1,:)=p(j,:);
p(j,:)=bj;
bj=U(j+1,:);
U(j+1,:)=U(j,:);
U(j,:)=bj;
bj=A(j+1,:);
A(j+1,:)=A(j,:);
A(j,:)=bj;
bj=P(j+1,:);
P(j+1,:)=P(j,:);
P(j,:)=bj;
bj=Q(j+1,:);
Q(j+1,:)=Q(j,:);
Q(j,:)=bj;
end
end
end
U0=U.*(cos(A)+sin(A).*1i);
A=A*180/pi;
S1=zeros(m,1);
S2=zeros(m,1);
S21=zeros(m,1);
for k=1:m %求支路损耗
i=z(k,1);
j=z(k,2);
if z(k,5)==1
S1(k,1)=U(i,1)^2*(-1i*z(k,8)/2)+U0(i,1)*(conj(U0(i,1))-conj(U0(j,1)))*conj(1/(z(k,6)+1i*z(k,7)));
S2(k,1)=U(j,1)^2*(-1i*z(k,8)/2)+U0(j,1)*(conj(U0(j,1))-conj(U0(i,1)))*conj(1/(z(k,6)+1i*z(k,7)));
S21(k,1)=S1(k,1)+S2(k,1);
else if z(k,5)==2
S1(k,1)=U(i,1)^2*conj(1i*z(k,8)/2+(1-z(k,14))/((z(k,14))^2*(z(k,6)+1i*z(k,7))))+U0(i,1)*(conj(U0(i,1))-conj(U0(j,1)))*conj(1/(z(k,14)*(z(k,6)+1i*z(k,7))));
S2(k,1)=U(j,1)^2*conj(1i*z(k,8)/2+(z(k,14)-1)/(z(k,14)*(z(k,6)+1i*z(k,7))))+U0(j,1)*(conj(U0(j,1))-conj(U0(i,1)))*conj(1/(z(k,14)*((z(k,6)+1i*z(k,7)))));
S21(k,1)=S1(k,1)+S2(k,1);
end
end
end
S21=S21*100;
Ps=real(S21);
Qs=imag(S21);
tab1=[{'节点'},{'节点电压'},{'节点相角'}];
tab2=[ p(:,1) U A ];
tab3=[{'最大误差'},{'最大误差节点'}];
tab4=[ maxError Node ];
tab5=[{'始节点'},{'终结点'},{'有功损耗'},{'无功损耗'}];
tab6=[ z(:,1) z(:,2) Ps Qs];
xlswrite('newton.xls',tab1,1,'A1');
xlswrite('newton.xls',tab2,1,'A2');
xlswrite('newton.xls',tab3,1,'D1');
xlswrite('newton.xls',tab4,1,'D2');
xlswrite('newton.xls',tab5,1,'F1');
xlswrite('newton.xls',tab6,1,'F2');