牛顿拉夫逊法潮流计算.zip

clear all;clc format long; ac=0.00001; %设定精度作为收敛条件 nl=6;n=6; %支路数和节点 Vs=[1.05 0 0 0 0]; type=[1 2 2 2 2 3]; % 1--PV节点， 2--PQ节点, 3--平衡节点 %% 满发最大负荷运行方式下 Ps=[8.37 -2.900 -3.100 -2.300 -4.250]; %给定节点有功功率 Qs=[0 -0.350 -0.600 -0.850 -0.800]; %给定节点无功功率 %% 轻载最小负荷运行方式下 %Ps=[2.511 -1.856 -1.984 -1.472 -2.72]; %给定节点有功功率 %Qs=[0 -0.224 -0.384 -0.544 -0.512]; %给定节点无功功率 %% 线路参数 %1(FromBus) 2(ToBus) 3(R+jX) 4(jB/2) B1=[ 2 6 0.001562 + 0.008909i 0.049368i 2 3 0.000837 + 0.004773i 0.02662i 2 4 0.000558 + 0.003182i 0.017424i 1 3 0.00106 + 0.006045i 0.033396i 4 5 0.000446 + 0.002545i 0.014036i 1 5 0.001395 + 0.007955i 0.044095i ]; %% N—1潮流校核用 answer=input('是否进行N-1校核? [y/n]n >> ','s'); if answer == 'Y';answer = 'y';end if answer == 'y' ls=input('请输入N-1校验支路 格式:[i j] i--线路始端, j--线路末端：'); j=0; for i=1:n if ((B1(i,1)==ls(1))&&(B1(i,2)==ls(2)))||((B1(i,1)==ls(2))&&(B1(i,2)==ls(1))) B1(i,3)=2*B1(i,3); B1(i,4)=0.5*B1(i,4); else j=j+1; end end if j==n fprintf('不存在这样的线路！\n'); end end %% 节点导纳矩阵计算 Y=zeros(nl,nl); for i=1:nl p=B1(i,1);q=B1(i,2); Y(p,q)=-1./B1(i,3);Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4); end %% 牛顿拉夫逊法潮流计算 for i=1:n e(i)=1;f(i)=0; end e(1)=1.05;e(nl)=1.04545; %给定电压初值 G=real(Y);B=imag(Y); %提取节点导纳矩阵的实部和虚部 %进入迭代------------------------------------------------ t=0; %迭代次数 S1=0;S2=0; DW=zeros(2*nl-2,1);%节点功率不平衡量矩阵 DV1=zeros(2*nl-2,1);%节点电压不平衡量矩阵 while 1 for i=1:n-1 if type(i)==2 %PQ节点方程 for j=1:n S1=S1+G(i,j)*e(j)-B(i,j)*f(j); S2=S2+G(i,j)*f(j)+B(i,j)*e(j); end DP(i)=Ps(i)-e(i)*S1-f(i)*S2; DQ(i)=Qs(i)-f(i)*S1+e(i)*S2; S1=0;S2=0; else for j=1:n %PV节点方程 S1=S1+G(i,j)*e(j)-B(i,j)*f(j); S2=S2+G(i,j)*f(j)+B(i,j)*e(j); end DP(i)=Ps(i)-e(i)*S1-f(i)*S2; DV1(i)=Vs(i)^2-e(i)^2-f(i)^2; S1=0;S2=0; end end for i=1:n-1 if type(i)==2 DW(2*i-1)=DP(i);DW(2*i)=DQ(i); else DW(2*i-1)=DP(i);DW(2*i)=DV1(i); end end for i=1:2*n-2 M_abs(i)=abs(DW(i)); end DV=zeros(2*nl-2,1); %节点电压修正量矩阵 if max(M_abs)>ac %收敛条件 Ja=Jacobi(G,B,e,f,type,n); %求解雅克比矩阵 DV=-inv(Ja)*DW; %解修正方程式 for i=1:n-1 e(i)=e(i)+DV(2*i-1); f(i)=f(i)+DV(2*i); end %修正各节点的电压 else break; end t=t+1; %记录迭代次数 E(t,:)=e;F(t,:)=f; end %%电压水平、网损计算和载流校验---------------------------------------- V=complex(E(t,:),F(t,:)); %各节点电压精确值 Yl=zeros(nl,nl); %支路阻抗矩阵 yl=zeros(nl,nl); %支路阻抗矩阵 I=zeros(nl,nl); %线路载流矩阵 DS=zeros(nl,nl); %网络损耗矩阵 S=zeros(nl,nl); %线路首端功率 S2=zeros(nl,nl); %线路末端功率 num=zeros(nl,1); for i=1:nl num(i)=i; end for i=1:nl p=B1(i,1);q=B1(i,2); Yl(p,q)=1./B1(i,3);Yl(q,p)=Yl(p,q);%支路阻抗矩阵 yl(p,q)=B1(i,4);yl(q,p)= yl(p,q); %支路电纳矩阵 end for i=1:nl for j=1:nl S(i,j)=abs(V(i))^2*yl(i,j)+V(i)*(conj(V(i))-conj(V(j)))*conj(Yl(i,j));%线路首段功率 Sl(i,j)=V(i)*(conj(V(i))-conj(V(j)))*conj(Yl(i,j)); Pl=real(Sl);Ql=imag(Sl);%网损计算用的功率 % if S(i,j)~=0 % end for k=1:nl if i==B1(k,1)&j==B1(k,2) DS(i,j)=(Pl(i,j)^2+Ql(i,j)^2)/abs(V(i))^2*B1(k,3); %线路损耗 S2(i,j)=-(Sl(i,j)-DS(i,j)-abs(V(j))^2*yl(i,j)); %线路末端功率 I(i,j)=sqrt((Pl(i,j)^2+Ql(i,j)^2)/abs(V(i))^2/3); %线路载流 end end end end %% 下面是Bus Data和Branch Data输出参数设置 v=abs(V'); %幅值 o=angle(V)/pi*180; %幅角 o=o'; DS1=zeros(nl,1);DS2=zeros(nl,1);I1=zeros(nl,1); %提取网损和载流非零量 P1=zeros(nl,1);Q1=zeros(nl,1); %线路首端有功和无功 P2=zeros(nl,1);Q2=zeros(nl,1); %线路末端有功和无功 P3=zeros(nl,1);Q3=zeros(nl,1); %节点注入有功和无功 TotalP=0;TotalQ=0; for i=1:nl j=B1(i,1); k=B1(i,2); P1(i)=real(S(j,k)); Q1(i)=imag(S(j,k)); P2(i)=real(S2(j,k)); Q2(i)=imag(S2(j,k)); DS1(i)=real(DS(j,k)); TotalP=TotalP+DS1(i); DS2(i)=imag(DS(j,k)); TotalQ=TotalQ+DS2(i); I1(i)=I(j,k); if type(i)==2 P3(i)=Ps(i); %给定节点有功功率 Q3(i)=Qs(i); else for m=1:nl P3(i)=P3(i)+real(S(i,m)); Q3(i)=Q3(i)+imag(S(i,m)); end end end fprintf('Power Flow Result: iter = %2d',t); fprintf('\n 每次迭代的电压:'); V=complex(E,F) %记录每次迭代的电压 fprintf('\n ======================================================================================'); fprintf('\n Bus Data:'); fprintf('\n BusNum V(pu) Angle(°) Pe(pu) Qe(pu)'); fprintf('\n %5d %15.8f %15.8f %12.4f %13.4f',[num,v,o,P3,Q3]'); fprintf('\n'); fprintf('\n ======================================================================================'); fprintf('\n Branch Data:'); fprintf('\n From To LineCurrent From Bus Injection To Bus Injection Loss (I^2 * Z) '); fprintf('\n Bus Bus I(pu) P(pu) Q(pu) P(pu) Q(pu) P(pu) Q(pu) '); fprintf('\n --------- ----------- ------------------ ---------------- ---------------'); fprintf('\n %2d %6d %10.5f %11.5f %9.5f %10.5f %9.5f %11.5f %8.5f',[B1(:,1), B1(:,2), I1,P1,Q1,P2,Q2,DS1,DS2]'); fprintf('\n ------- -------'); fprintf('\n Total: %5.5f %8.5f',[TotalP,TotalQ]); fprintf('\n');