function [bus_res,S_res]=PowerFlow_NR_2
% 序号 横向电压 纵向电压 P Q 节点类型
bus = [ 1 1.06 0 0 0 3;
2 1 0 0.20 0.20 1;
3 1 0 -0.45 -0.15 1;
4 1 0 -0.40 -0.05 1;
5 1 0 -0.60 -0.10 1];
% 结点1 结点2 变比(为i侧电压除以j侧电压,若反向则为副)
line= [ 1 2 0.082 0.427 0 0.028 1;
2 3 0.145 0.581 0 0 1;
3 4 0.104 0.518 0 0.018 1;
5 4 0.031 0.248 0 0 0.95;
4 1 0.163 0.754 0 0.014 1;
2 0 0 0 0 0.04 1];
[nb,mb]=size(bus);
[nl,ml]=size(line);
format long;
EPS=1.0e-5;
YtYm = z_YtYm(line); %线路矩阵导纳化
%bus_Y=y_Y(YtYm,nb); %结点导纳矩阵生成
j=sqrt(-1);
bus_Y = [ 6.250-j*18.75 -5.000+j*15.000 -1.250+j*3.750 0 0;
-5.000+j*15.000 10.834-j*32.500 -1.667+j*5.000 -1.667+j*5.000 -2.500+j*7.500;
-1.250+j*3.750 -1.667+j*5.000 12.917-j*38.750 -10.000+j*30.000 0;
0 -1.667+j*5.000 -10.000+j*30.000 12.917-j*38.750 -1.250+j*3.750;
0 -2.500+j*7.500 0 -1.250+j*3.750 3.750-j*11.25];
for t=1:100;
[dat_PQ,I_bus,count_PQ,count_PV,count_PH]=d_PQ(bus_Y,bus);
Jacobian=Jac(bus_Y,bus,I_bus,count_PQ,count_PV);
Jac_ni=inv(Jacobian);
total=2*(nb-1);
dU=zeros(total,1);
%电压修正量的计算（功率偏差*雅克比矩阵）
for i=1:total;
for j=1:total;
dU(i,1)=dU(i,1)+Jac_ni(i,j)*dat_PQ(j,1);
end
end
%将修正量和初值相加,得到新值
for i=1:(nb-1);
bus(i+1,2)=bus(i+1,2)+dU(2*i);
bus(i+1,3)=bus(i+1,3)+dU(2*i-1);
end
if(max(abs(dU))<EPS)&(max(abs(dat_PQ))<EPS)
break;
end
end
%迭代收敛后计算平衡结点功率
[S_PH,S_load]=Js_S(bus_Y,bus);
bus(1,4)=real(S_PH);
bus(1,5)=imag(S_PH);
myf=fopen('Result.m','w');
fprintf(myf,'----------------------------------by linaolin--------------------------------------------------------\n\n');
fprintf(myf,'---------------------牛顿-拉夫逊潮流计算结果----------------------------------------\n\n');
fprintf(myf,'节点电压及注入功率：\n');
fprintf(myf,'节点号 e f P Q 节点类型\n');
for i=1:nb;
for k=1:6;
if (k==2)||(k==3)||(k==4)||(k==5)
fprintf(myf,'%10.4f ',bus(i,k));
end
if (k==1)||(k==6)
fprintf(myf,'%d ',bus(i,k));
end
end
fprintf(myf,'\n');
end
fprintf(myf,'节点导纳矩阵：\n');
fprintf(myf,'Y=\n');
for i=1:nb;
for k=1:nb;
temp_R=real(bus_Y(i,k));
tenp_I=imag(bus_Y(i,k));
if tenp_I<0;
fprintf(myf,'%8.4f %8.4f i ',temp_R,tenp_I);
end
if tenp_I>=0;
fprintf(myf,'%8.4f+%8.4f i ',temp_R,tenp_I);
end
end
fprintf(myf,'\n');
end
fprintf(myf,'各线路功率：\n');
fprintf(myf,'S=\n');
for i=1:nb;
for k=1:nb;
temp_R=real(S_load(i,k));
tenp_I=imag(S_load(i,k));
if tenp_I<0;
fprintf(myf,'%8.4f %8.4f i ',temp_R,tenp_I);
end
if tenp_I>=0;
fprintf(myf,'%8.4f+%8.4f i ',temp_R,tenp_I);
end
end
fprintf(myf,'\n');
end
fprintf(myf,'迭代次数：\n');
fprintf(myf,'Time= %d：\n',t);
fclose(myf);
%电阻表示的矩阵转化为导纳表示的矩阵
function YtYm=z_YtYm(line)
[nl,ml]=size(line);
YtYm=zeros(nl,5);
YtYm(:,1:2)=line(:,1:2);
j=sqrt(-1);
for k=1:nl;
I=line(k,1);
J=line(k,2);
Zt=line(k,3)+j*line(k,4);
if(Zt~=0)
Yt=1/Zt;
end
Ym=line(k,5)+j*line(k,6);
K=line(k,7);
if(K==1)&(J~=0) %普通线路K=1
YtYm(k,3)=Yt;
YtYm(k,4)=Ym;
YtYm(k,5)=Ym;
end
if(K==1)&(J==0) %对地支路
YtYm(k,4)=Ym;
end
if(K~=1) %变压器线路
if(K>=0)
YtYm(k,3)=Yt/K;
YtYm(k,4)=Yt*(1-K)/K/K;
YtYm(k,5)=Yt*(K-1)/K
end
if(K<0)
fan_K=-K;
YtYm(k,3)=Yt/fan_K;
YtYm(k,4)=Yt*(fan_K-1)/fan_K;
YtYm(k,5)=Yt*(1-fan_K)/fan_K/fan_K;
end
end
end
%结点导纳矩阵的形成，基于线路导纳矩阵
%入口参数：线路导纳矩阵，节点数
function bus_Y=y_Y(Y,nb)
[nl,ml]=size(Y);
bus_Y=zeros(nb,nb);
for k=1:nl;
number1=Y(k,1);
number2=Y(k,2);
%对地支路计算
if (number2==0)
temp=bus_Y(number1,number1);
bus_Y(number1,number1)=temp+Y(k,4);
continue;
end
%自导纳计算*********************************************
temp=bus_Y(number1,number1);
bus_Y(number1,number1)=temp+Y(k,3)+Y(k,4);
temp=bus_Y(number2,number2);
bus_Y(number2,number2)=temp+Y(k,3)+Y(k,5);
%*******************************************************
%互导纳计算：
temp=-Y(k,3);
bus_Y(number1,number2)=temp;
bus_Y(number2,number1)=temp;
end
%功率偏差计算并计算各节点注入电流
function [dat_PQ,I_bus,count_PQ,count_PV,count_PH]=d_PQ(Y,bus);
[nb,mb]=size(bus);
count_PQ=0;
count_PV=0;
count_PH=0;
res_PQ=zeros(nb,mb);
res_PV=zeros(nb,mb);
for k=1:nb;
if(bus(k,6)==1)
count_PQ=count_PQ+1;
res_PQ(count_PQ,:)=bus(k,:);
end
if(bus(k,6)==2)
count_PV=count_PV+1;
res_PV(count_PV,:)=bus(k,:);
end
if(bus(k,6)==3)
count_PH=count_PH+1;
number_PH=bus(k,1);
end
end
total=2*(count_PQ+count_PV);
dat_PQ=zeros(total,1);
new_PQ=zeros(total,1);
temp_line=1;
%针对PQ结点计算偏差
for i=1:count_PQ;
%Pi计算
temp=0;
for j=1:nb;
temp=temp+res_PQ(i,2)*(real(Y(res_PQ(i,1),j))*bus(j,2)-imag(Y(res_PQ(i,1),j))*bus(j,3))+res_PQ(i,3)*(real(Y(res_PQ(i,1),j))*bus(j,3)+imag(Y(res_PQ(i,1),j))*bus(j,2));
end
new_PQ(temp_line,1)=temp;
temp_line=temp_line+1;
%Qi计算
temp=0;
for j=1:nb;
temp=temp+res_PQ(i,3)*(real(Y(res_PQ(i,1),j))*bus(j,2)-imag(Y(res_PQ(i,1),j))*bus(j,3))-res_PQ(i,2)*(real(Y(res_PQ(i,1),j))*bus(j,3)+imag(Y(res_PQ(i,1),j))*bus(j,2));
end
new_PQ(temp_line,1)=temp;
temp_line=temp_line+1;
end
%针对PV结点计算偏差
for i=1:count_PV;
%Pi计算
temp=0;
for j=1:nb;
temp=temp+res_PV(i,2)*(real(Y(res_PV(i,1),j))*bus(j,2)-imag(Y(res_PV(i,1),j))*bus(j,3))+res_PV(i,3)*(real(Y(res_PV(i,1),j))*bus(j,3)+imag(Y(res_PV(i,1),j))*bus(j,2));
end
new_PQ(temp_line,1)=temp;
temp_line=temp_line+1;
%△U
temp=res_PV(i,2)^2+res_PV(i,3)^2
temp_line=temp_line+1;
end
for i=1:count_PQ;
dat_PQ(2*i-1,1)=res_PQ(i,4)-new_PQ(2*i-1,1);
dat_PQ(2*i ,1)=res_PQ(i,5)-new_PQ(2*i ,1);
end
for i=1:count_PV
dat_PQ(count_PQ*2+2*i-1,1)=res_PV(i,4)-new_PQ(count_PQ*2+i,1);
dat_PQ(count_PQ*2+2*i ,1)=res_PV(i,2)^2+res_PV(i,3)^2-new_PQ(count_PQ*2+i,1);
end
I_bus=zeros(nb-1,2);
j=sqrt(-1);
for i=1:count_PQ;
temp=(new_PQ(2*i-1,1)-j*new_PQ(2*i,1))/(res_PQ(i,2)-j*res_PQ(i,3));
I_bus(i,1)=real(temp);
I_bus(i,2)=imag(temp);
end
for i=1:count_PV
temp=(new_PQ(2*count_PQ+2*i-1,1)-j*new_PQ(2*count_PQ+2*i,1))/(res_PV(i,2)-j*res_PV(i,3));
I_bus(count_PQ+i,1)=real(temp);
I_bus(count_PQ+i,2)=imag(temp);
end
%雅克比矩阵的形成
function Jacobian=Jac(Y,bus,I_bus,count_PQ,count_PV)
[nb,mb]=size(bus);
total=2*(nb-1);
Jacobian=zeros(total,total);
res_H=zeros(nb-1,nb-1);
res_N=zeros(nb-1,
