电力系统潮流计算详解

1、首先声明一下，这些是从网站上转载的，不是本人上编写的基于MATLAB的电力系统潮流计算%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点%编号，将变压器的串联阻抗置于低压侧处理。%第三列为支路的串列阻抗参数。%第四列为支路的对地导纳参数。%第五烈为含变压器支路的变压器的变比%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点%负荷功率参数；第三列为

2、节点电压参数；第六列为节点类型参数，其中%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地%参数。n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:

3、X=');Y=zeros(n);Times=1; %置迭代次数为初始值%创建节点导纳矩阵for i=1:n1 if B1(i,6)=0 %不含变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4); else %含有变压器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5); Y(q,p)=Y(p,

4、q); Y(p,p)=Y(p,p)+1/B1(i,3); Y(q,q)=Y(q,q)+1/(B1(i,5)2*B1(i,3); endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化%创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理 if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2

5、(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i=isb&B2(i,6)=3 h=h+1; f

6、or j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(

7、i,j)*real(B2(j,3); end endendOrgS%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)=3 t=t+1; PVU(t,1)=B2(i,3); endendPVU%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理 if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);

8、 endendt=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; endendDetaS%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,

9、1)*sqrt(-1)/conj(B2(i,3); endendI%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理 if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j %对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B

10、2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else %非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h

11、-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endendk=0;for i=1:n %对PV节点的处理 if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j %对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(

12、Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else %非对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacb

13、i(2*h,2*k)=0; end if k=(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0; end end end endendJacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理 if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n %对PV节点的处理 if B2(i,6)=3 j=j+1

14、; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendB2%开始循环*while abs(max(DetaU)>prOrgS=zeros(2*n-2,1); %!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,

15、3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(

16、2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end enden

17、dOrgS%创建DetaSh=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); endendt=0;for i=1:n if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag

18、(B2(i,3)2; endendDetaS%创建Ii=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endendI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,

19、3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*

20、real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) k=0; end end end endendk=0;for i=1:n if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Ja

21、cbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-

22、1)=0; Jacbi(2*h,2*k)=0; end if k=(n-1) k=0; end end end endendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);

23、 endendB2Times=Times+1; %迭代次数加1endTimes一个原始数据的例子节点数 5支路数 5平衡节点编号 5精度pr 0.000001B1(支路参数矩阵)1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1B2(节点参数矩阵)0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1X

24、(节点号和对地参数)1 0;2 0;3 0;4 0;5 0 程序2：%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算% n=input('请输入系统节点数目:n=');% nl=input('请输入系统支路数目:nl=');% ph=input('请输入平衡节点号:ph=');% jd=input('请输入误差精度：jd=');% B=input('请输入由支路参数形成的矩阵:B=');% A=input('请输入各节点参数形成的矩阵:A=');clear alln=5;nl=5;ph=5

25、;jd=1E-6;%支路参数矩阵Br=1 2 0.04+0.25*j 0.25j 1; 1 3 0.1+0.35*j 0 1; 2 3 0.08+0.30*j 0.25*j 1; 2 4 0.015*j 0 1.05; 3 5 0.03*j 0 1.05;%节点参数矩阵:2-PQ节点，3-PV节点，4-平衡节点A=-1.6-0.8*j 1 0 2; -2-1*j 1 0 2; -3.7-1.3*j 1 0 2; 5 1.05 1.05 3; 0 1.05 0 4;%雅克比矩阵形成Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);for i=1:

26、nl p=Br(i,1); q=Br(i,2); Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./(Br(i,3)+Br(i,4); Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)2)+Br(i,4);end disp('节点导纳矩阵为：Y=');disp(Y);G=real(Y);B=imag(Y);for i=1:n e(i)=real(A(i,2); f(i)=imag(A(i,2); V(i)=A(i,3);endfor i=1:n S(i)=A(i);endP=r

27、eal(S);Q=imag(S);%雅克比矩阵求取Ci=0;a=1;NO=2*n;N=NO-1;while a=0 a=0; for i=1:n if i=ph C(i)=0; D(i)=0; for p=1:n C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p); D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p); end P1=e(i)*C(i)+f(i)*D(i); Q1=f(i)*C(i)-e(i)*D(i); V2=e(i)2+f(i)2; if A(i,4)=3 DP=P(i)-P1; DQ=Q(i)-Q1; for k=1:n%非平衡节点时，非对

28、角线元素 if k=ph & k=i X1=-G(i,k)*f(i)+B(i,k)*e(i); X2=-G(i,k)*e(i)-B(i,k)*f(i); X3=-X2; X4=X1; p=2*i-1; q=2*k-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X3; J(m,N)=DQ; q=q+1; J(p,q)=X2; J(m,q)=X4; %非平衡节点时对角线元素 elseif k=ph & k=i X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X3=

29、-C(i)+G(i,i)*e(i)+B(i,i)*f(i); X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i); p=2*i-1; q=2*k-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X3; J(m,N)=DQ; q=q+1; J(p,q)=X2; J(m,q)=X4; end end else DP=P(i)-P1; DV=V(i)2-V2; for s=1:n if s=ph & s=i X1=-G(i,s)*f(i)+B(i,s)*e(i); X2=-G(i,s)*e(i)-B(i,s)*f(i); X5=0; X6=0; p=

30、2*i-1; q=2*s-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X5; J(m,N)=DV; q=q+1; J(p,q)=X2; J(m,q)=X6; elseif s=ph & s=i X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i); X5=-2*f(i); X6=-2*e(i); p=2*i-1; q=2*s-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X5; J(m,N)=DV; q=q+1; J(p,q)=X2; J(m

31、,q)=X6; end end end end end disp('系统雅克比矩阵为：J='); disp(J);%对于高压电力系统，最大元素出现在对角线位置，简单计算如下 DJ=J(:,1:NO-2); DP=J(1:NO-2,N); DU=DJDP; disp('*'); disp(' ') disp('第几次的修正值DU？'); disp(Ci); disp('*'); disp(DU); for i=1:NO-2 eps=abs(DP(i); if eps>=jd; a=a+1; end end Ci

32、=Ci+1; for i=1:n-1 e(i)=e(i)-DU(2*i); f(i)=f(i)-DU(2*i-1); end E=e+j*f; disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&'); disp(' '); disp('第几次迭代后的近似计算值'); disp(Ci); disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&am