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

%本程序的功能是用牛拉法进行含风电场的电力系统潮流计算functions1=pf(a,B1,B2)n=a(1);%节点数nl=a(2);%支路数isb=a(3);%平衡节点号pr=a(4);%误差精度fori=1:nforj=1:nG(i,j)=0;B(i,j)=0;endend%求导纳矩阵%B1为支路参数矩阵，其每一行格式为首节点号，末节点号，支路电导，支路电纳，首节点对地电纳，末节点对地电纳fori=1:nlp=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);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);endfori=1:nB(i,i)=B(i,i)+B2(i,5);end%求导纳矩阵%B2为节点参数矩阵，每一行格式为节点注入有功，注入无功，电压实部，电压虚部，对地电纳，节点类型%节点类型：1为平衡节点，2为PQ节点，3为PV节点，4为风电场节点%Bf为风电场参数，格式为：有功功率，定子电抗，转子漏抗，转子电阻，励磁电抗fori=1:nifB2(i,6)==4p(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;elseP(i)=B2(i,1);Q(i)=B2(i,2);ende(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;whileIT2~=0IT2=0;a=a+1;fori=1:nifi~=isbC(i)=0;D(i)=0;forj1=1:nC(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);endP1=C(i)*e(i)+f(i)*D(i);Q1=f(i)*C(i)-D(i)*e(i);V2=e(i)^2+f(i)^2;ifB2(i,6)==2DP=P(i)-P1;DQ=Q(i)-Q1;forj1=1:nifj1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*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;elseifj1==i&j1~=isbX1=-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;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;endendelseifB2(i,6)==3DP=P(i)-P1;DV=V(i)^2-V2;forj1=1:nifj1~=isb&j1~=iX1=-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;elseifj1==i&j1~=isbX1=-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;endendelseDP=P(i)-P1;forj1=1:nifj1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*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;elseifj1==i&j1~=isbX1=-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*a4*V2^2/(sqrt(V2^4*a4-a5)))/(a2*V2^2+a3*(sqrt(V2^4*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;endendendendend%求雅克比矩阵fork=3:N0k1=k+1;N1=N;fork2=k1:N1J(k,k2)=J(k,k2)./J(k,k);endJ(k,k)=1;ifk~=3;k4=k-1;fork3=3:k4fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endifk==N0,break;endfork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endelsefork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endendendfork=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N);k1=k+1;f(L)=f(L)-J(k1,N);endfork=3:N0DET=abs(J(k,N));ifDET=prIT2=IT2+1;endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解w=-J*Vfid1=fopen('out1.txt','wt')fprintf(fid1,‘各节点的实际电压标幺值E和电压大小V为：\n’);fork=1:nV(k)=sqrt(e(k)^2+f(k)^2);fprintf(fid1,'E(%d)=%8.5f+j%8.5f,V(%d)=%8.5f\n',k,e(k),f(k),k,V(k));endforp=1:nifp==isbfprintf(fid1,'平衡节点的功率S为:\n');C(p)=0;D(p)=0;forq=1:nC(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);endSp(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.5f\n',Sp(p),Sq(p));endendfprintf(fid1,'各条支路的首端功率为:\n');fori=1:nlp=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.5f\n',p,q,Sip(i),Siq(i));endfprintf(fid1,'¸各条支路的末端功率为:\n');fori=1:nlp=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,%d)=%8.5f+j%8.5f\n',q,p,Sjp(i),Sjq(i));endfprintf(fid1,'各条支路的功率损耗为:\n');fori=1:nlDSp(i)=Sip(i)+Sjp(i);DSq(i)=Siq(i)+Sjq(i);fprintf(fid1,'DS(%d)=%8.5f+j%8.5f\n',i,DSp(i),DSq(i));endfclose(fid1);