含风电场的电力系统潮流计算程序

1、精选优质文档-倾情为你奉上%本程序的功能是用牛拉法进行含风电场的电力系统潮流计算 function s1=pf(a,B1,B2) n=a(1);%节点数 nl=a(2);%支路数 isb=a(3);%平衡节点号 pr=a(4);%误差精度for i=1:n for j=1:n G(i,j)=0; B(i,j)=0; endend%求导纳矩阵%B1为支路参数矩阵，其每一行格式为首节点号，末节点号，支路电导，支路电纳，首节点对地电纳，末节点对地电纳for i=1:nl p=B1(i,1);q=B1(i,2) G(p,q)=G(p,q)-B1(i,3); B(p,q)=B(p,q)- B1(i,4)

2、; G(q,p)=G(p,q); B(q,p)=B(p,q); G(p,p)=G(p,p)+ B1(i,3); B(p,p)=B(p,p)+B1(i,4); G(q,q)=G(q,q)+ B1(i,3); B(q,q)=B(q,q)+B1(i,4);endfor i=1:n B(i,i)=B(i,i)+B2(i,5);end%求导纳矩阵%B2为节点参数矩阵，每一行格式为节点注入有功，注入无功，电压实部，电压虚部，对地电纳，节点类型%节点类型：1为平衡节点，2为PQ节点，3为PV节点，4为风电场节点%Bf为风电场参数，格式为：有功功率，定子电抗，转子漏抗，转子电阻，励磁电抗for i=1:n i

3、f B2(i,6)=4 p(i)=Bf(1); a1=2*p(i)2*(Bf(2)+Bf(3)2; a2=-Bf(4)3*Bf(5); a3=Bf(4)2*Bf(5); a4=Bf(4)2; a5=4*p(i)2*(Bf(2)+Bf(3)2*Bf(4)2; else P(i)=B2(i,1); Q(i)=B2(i,2); end e(i)= B2(i,3); f(i)=B2(i,4); V(i)=sqrt(e(i)2+f(i)2);endICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2=0 IT2=0;a=a+1; for i=1:n if i=isb C(i)

4、=0; D(i)=0; for j1=1:n C(i)= C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); D(i)= D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); end P1=C(i)*e(i)+f(i)*D(i); Q1=f(i)*C(i)-D(i)*e(i); V2=e(i)2+f(i)2;if B2(i,6)=2 DP=P(i)-P1; DQ=Q(i)-Q1; for j1=1:n if j1=isb&j1=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X3

5、=X2; X4=-X1; p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ; m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; elseif j1=i&j1=isb X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;

6、m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; end endelseif B2(i,6)=3 DP=P(i)-P1; DV=V(i)2-V2; for j1=1:n if j1=isb&j1=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1)*e(i)-G(i,j1)*f(i); X5=0; X6=0; p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2; el

7、seif j1=i&j1=isb X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X5=-2*e(i); X6=-2*f(i); p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; J(m,q)=X2; end endelse DP=P(i)-P1; for j1=1:n if j1=isb&j1=i X1=-G(i,j1)*e(i)-B(i,j1)*f(i); X2=B(i,j1

8、)*e(i)-G(i,j1)*f(i); X3=X2; X4=-X1; p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ; m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; elseif j1=i&j1=isb X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); x3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); x4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i); x5=2*a1*(a2+a3

9、*a4*V22/(sqrt(V24*a4-a5)/(a2*V22+a3*(sqrt(V24*a4-a5)2; X3=x3-e(i)*x5; X4=x4-f(i)*x5; p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1; J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2; end end end endend%求雅克比矩阵for k=3:N0 k1=k+1;N1=N; for k2=k1:N1 J(k,k2)=J(k,k2)./J(k,k); end J(k,k)=1; if k=3; k4=k-1; for k3=

10、3:k4 for k2=k1:N1 J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end if k=N0,break;end for k3=k1:N0 for k2=k1:N1 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end else for k3=k1:N0 for k2=k1:N1 J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2); end J(k3,k)=0; end endend for k=3:2:N0-1 L=(k+1)./2; e(L)=e(L)

11、-J(k,N); k1=k+1; f(L)=f(L)-J(k1,N);endfor k=3:N0 DET=abs(J(k,N); if DET>=pr IT2=IT2+1; endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解w=-J*Vfid1=fopen('out1.txt','wt')fprintf(fid1,各节点的实际电压标幺值E和电压大小V为：n);for k=1:n V(k)=sqrt(e(k)2+f(k)2); fprintf(fid1,'E(%d)=%8.5f+j%8.5f, V(%d)=%8.5fn

12、',k,e(k),f(k),k,V(k);endfor p=1:n if p=isb fprintf(fid1,'平衡节点的功率S为:n'); C(p)=0;D(p)=0; for q=1:n C(p)=C(p)+G(p,q)*e(q)-B(p,q)*f(q); D(p)=D(p)-G(p,q)*f(q)-B(p,q)*e(q); end Sp(p)=e(p)*C(p)-f(p)*D(p); Sq(p)=(e(p)*D(p)+f(p)*C(p); fprintf(fid1,'S=%8.5f+j%8.5fn',Sp(p),Sq(p); end endfpr

13、intf(fid1,'各条支路的首端功率为:n');for i=1:nl p=B1(i,1);q=B1(i,2); sip(i)=-e(p)*B1(i,5)+(e(q)-e(p)*G(p,q)-(f(q)-f(p)*B(p,q); siq(i)=(f(p)*B1(i,5)+(e(p)-e(q)*B(p,q)-(f(q)-f(p)*G(p,q); Sip(i)=e(p)*sip(i)-f(p)*siq(i); Siq(i)=e(p)*siq(i)+f(p)*sip(i); fprintf(fid1,'S(%d,%d)=%8.5f+j%8.5fn',p,q,Sip(i),Siq(i); endfprintf(fid1,'¸各条支路的末端功率为:n');for i=1:nl p=B1(i,1);q=B1(i,2); sjp(i)=-e(q)*B1(i,6)+(e(p)-e(q)*G(p,q)-(f(p)-f(q)*B(p,q); sjq(i)=(f(q)*B1(i,6)+(e(q)-e(p)*B(p,q)-(f(p)-f(q)*G(p,q); Sjp(i)=e(q)*sjp(i)-f(q)*sjq(i); Sjq(i)=e(q)*sjq(i)+f(q)*sjp(i); fprintf(fid1,'S(%d,