牛顿——拉夫逊潮流计算的matlab程序

%牛顿——拉夫逊潮流计算。 %直角坐标系 %输出值：s_node：节点视在功率 % u_node：节点电压 %输入值：Y： 导纳矩阵 % m： PQ节点数,包含平衡节点 % n： 总结点数 % s： 平衡节点编号，一般定为1号编号 % S_PQ: PQ节点的视在功率[1:m],其中平衡节点的视在功率暂定为0； % P_PV：PV节点的有功功率[m+1:n] % U_PV: PV节点的电压给定值,只有幅值，没有相位角，是实数[m+1:n] % us: 平衡节点的电压给定值，复数。 % Qmax；PV节点无功功率最大值 % Qmin：PV节点无功功率最小值 % epxl: 最大电压变量dumax的最大值，即规定的潮流计算的精度，实数。 %function [s_node,u_node]=newton(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,expl) function [s_node,u_node,Sij,dSij]=newton() %************************************************************************** %以下为实验数据，代替输入数据 Y=[6.250-18.695*j -5.000+15.00*j -1.250+3.750*j 0.000+0.000*j 0.000+0.000*j -5.000+15.000*j 10.833-32.415*j -1.667+5.000*j -1.667+5.000*j -2.500+7.5000*j -1.250+3.750*j -1.667+5.000*j 12.917-38.695*j -10.000+30.000*j 0.000+0.000*j 0.000+0.000*j -1.667+5.000*j -10.000+30.000*j 12.917-38.695*j -1.250+3.750*j 0.000+0.000*j -2.500+7.500*j 0.000+0.000*j -1.250+3.750*j 3.750-11.210*j]; % Y=[7.05-28.21j -5.88+23.5j -1.17+4.71j % -5.88+23.5j 5.88-23.5j 0 % -1.17+4.71j 0 1.17-4.38j]; m=5;%有2个PQ节点，包含一个平衡节点 n=5;%有三个节点 s=1;%第1号节点为平衡节点 S_PQ=[0,0.2+0.2*j,-0.45-0.15*j,-0.4-0.05*j,-0.6-0.1*j]; P_PV=[0.4]; U_PV=[1.1]; us=1.06+0*j; Qmax=50; Qmin=0.002; epxl=0.00001; %以上为实验数据，实际应用时为输入值 %************************************************************************** %以下为牛顿——拉夫逊法的程序 U=DYgaosai(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,0.0001);%利用高斯-塞尔德法得到比较准确的节点电压的初值。 fprintf('牛顿——拉夫逊潮流计算开始\n') dPQV2=zeros(2*(n-1),1);%不平衡量向量； Ykb=zeros(2*(n-1),2*(n-1));%雅克比矩阵； I=zeros(1,n); dU=zeros(2*(n-1),1);%电压不平衡量 dumax=0; k=1000;%设置迭代次数 Ss=0;%平衡节点的视在功率 S=zeros(1,n); Q_PV=ones(1,n-m);%PV节点的无功功率； Q_PV=Q_PV.*0.2; i=1; flag=zeros(1,n-m);%PV节点在计算过程中转化为PQ节点的标志，PV转化为PQ节点时，相应的标志位置1 while(k>=1)%迭代次数控制 %计算节点的不平衡量的dPQV2 for i=1:n %i:节点编号； %******************************* A=0; for a=1:n A=A+(Y(i,a)'*U(a)'); end a=1; A=U(i)*A; %*************************************** if(i<=s-1)%如果i号节点是平衡节点之前的PQ节点 dPQV2(2*i-1)=real(S_PQ(i))-real(A);%PQ节点有功功率不平衡量 dPQV2(2*i)=imag(S_PQ(i))-imag(A);%PQ节点无功功率不平衡量 else if(i>=s+1 && i<=m)%如果i号节点是平衡节点之后的PQ节点 dPQV2(2*(i-1)-1)=real(S_PQ(i))-real(A);%PQ节点有功功率不平衡量 dPQV2(2*(i-1))=imag(S_PQ(i))-imag(A);%PQ节点无功功率不平衡量 elseif(i>=m+1 && i<=n)%如果i号节点是PV节点 dPQV2(2*(i-1)-1)=P_PV(i-m)-real(A);%PV节点有功功率不平衡量 dPQV2(2*(i-1))=U_PV(i-m)^2-abs(U(i))^2;%PV节点电压平方不平衡量 end end end i=1; %节点的不平衡量的dPQV2计算结束 for i=1:n%计算节点注入电流I； if(i<=s-1) for a=1:n I(i)=I(i)+Y(i,a)*U(a); end a=1; elseif(i>=s+1) for a=1:n I(i)=I(i)+Y(i,a)*U(a); end a=1; end end i=1; for i=1:n %计算雅克比矩阵 if(i<=s-1)%如果i号节点是平衡节点之前的PQ节点 for a=1:n if(a~=i) if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点 Ykb(2*i-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数 Ykb(2*i-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数 Ykb(2*i,2*a-1)=-Ykb(2*i-1,2*a);%PQ节点的Jia参数=-Nia Ykb(2*i,2*a)=Ykb(2*i-1,2*a-1);%PQ节点的Lia参数=Hia else if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点,Ykb中的a->a-1;i不变 Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数 Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数 Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1));%PQ节点的Jia参数=-Nia Ykb(2*i,2*(a-1))=Ykb(2*i-1,2*(a-1)-1);%PQ节点的Lia参数=Hia elseif(a>=m+1 && a<=n)%Ykb中的a->a-1;i不变 Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数 Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数 Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1));%PV节点Jia参数=-Nia Ykb(2*i-1,2*(a-1))=Ykb(2*i-1,2*(a-1)-1);%PV节点Lia参数=Hia end end elseif(a==i) if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点 Ykb(2*i-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数 Ykb(2*i-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数 Ykb(2*i,2*a-1)=-Ykb(2*i-1,2*a)-2*real(I(i));%PQ节点的Jia参数 Ykb(2*i,2*a)=Ykb(2*i-1,2*a-1)+2*imag(I(i));%PQ节点的Lia参数 else if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点 Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数 Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数 Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1))-2*real(I(i));%PQ节点的Jia参数=-Nia-2real（I(i)）； Ykb(2*i,2*(a-1))=Ykb(2*i-1,2*(a-1)-1)+2*imag(I(i));%PQ节点的Lia参数=Hia+2imag(I(i)); elseif(a>=m+1 && a<=n) Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数 Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数 Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1))-2*real(I(i));%PQ节点的Jia参数=-Nia-2real（I(i)）； Ykb(2*i-1,2*(a-1))=Ykb(2*i-1,2*(a-1)-1)+2*imag(I(i));%PQ节点的Lia参数=Hia+2imag(I(i)); end end end end a=1; else if(i>=s+1 && i<=m)%如果i号节点是平衡节点之后的PQ节点;Ybk中的i->i-1; for a=1:n if(a~=i) if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点;Ybk中的i->i-1;a不变 Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数 Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数 Ykb(2*(i-1),2*a-1)=-Ykb(2*(i-1)-1,2*a);%PQ节点的Jia参数=-Nia Ykb(2*(i-1),2*a)=Ykb(2*(i-1)-1,2*a-1);%PQ节点的Lia参数Hia else if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点;Ybk中的i->i-1;a->a-1; Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数 Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数 Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1));%PQ节点的Jia参数=-Nia Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1);%PQ节点的Lia参数=Hia elseif(a>=m+1 && a<=