完整版电力系统潮流程序设计.docx

1、完整版电力系统潮流程序设计潮流程序的设计与开发1 数据结构的设计要求：将设备铭牌参数和有名值作为原始输入数据，潮流结果以有名值的形式输出。支路数据与节点数据分别以一个矩阵的形式输入，矩阵的每一行表示每一个节点或每一条支路，矩阵的每一列表示不同的参数数据。2 变量命名设计变量名称程序变量表示变量名称程序变量表示节点导纳矩阵无功功率电导视在功率电钠雅克比矩阵电压幅值平衡节点编号电压相角有功不平衡量有功功率无功不平衡量电压相角修正量电压幅值修正量3 程序流程图4 程序代码4.1 主程序n=input(请输入节点数:n=);l=input(请输入支路数:l=);%支路数不要包括三绕组变压器sw=inp

2、ut(请输入平衡节点号:sw=);ac=input(请输入误差精度:ac=);SB=input(请输入基准功率:SB=);B1=input(请输入支路参数:B1=)%支路参数节点参数和对地支路参数均以矩阵形式储存%第一列储存支路编号%第二列与第三列分别储存支路的两个端点，分别为p,q%第四列储存支路阻抗%第五列储存支路对地导纳,注意对地导纳不要除以2%第六列储存该支路是否含有变压器，有为1，无为0%第七列储存变压器变比k,p指向q的变压器变比为k:1，且k大于等于1%第八列储存变压器短路损耗%第九列储存变压器短路电压百分数%第十列储存变压器空载损耗%第十一列储存变压器空载电流百分数%第十二列储

3、存变压器低压侧电压%第十三列储存变压器额定功率%第十四列储存归算所取基准电压%注意，将三绕组变压器转换为双绕组变压器输入A1=input(请输入节点参数:A1=);%第一列为节点编号%第二列为注入发电功率%第三列为负荷功率%第四列为节点电压幅值,为方便起见，以标幺值形式表示%第五列为节点电压相角%第六列储存节点对地导纳%第七列为节点的类型，1为PQ节点，2为PV节点，3为平衡节点%首先求解双绕组变压器参数for i=1:l if B1(i,6)=1 Zt(i)=B1(i,8)*B1(i,12)2/(1000*B1(i,13)2); Xt(i)=B1(i,9)*B1(i,12)2/(100*B1

4、(i,13); Gt(i)=B1(i,10)/(1000*B1(i,12)2); Bt(i)=B1(i,11)*B1(i,13)/(100*B1(i,12)2); endendB2(:,1:7)=B1(:,1:7);for i=1:l B2(i,8)=Zt(i); B2(i,9)=Xt(i);endfor i=1:l if B1(i,6)=1 ZB(i)=B1(i,14)2/SB; YB(i)=1/ZB(i); B2(i,4)=B1(i,4)/ZB(i); B2(i,5)=B1(i,5)/YB(i); endendfor i=1:l if B1(i,6)=1 ZB(i)=B1(i,14)2/S

5、B; YB(i)=1/ZB(i); B2(i,8)=B2(i,8)/ZB(i); B2(i,9)=B2(i,9)/ZB(i); Gt(i)=Gt(i)/YB(i); Bt(i)=Bt(i)/YB(i); A1(B1(i,3),6)=Gt(i)-Bt(i)*(1i); endend%下面求解节点导纳矩阵Y=zeros(n);for i=1:n if A1(i,6)=0 Y(i,i)=A1(i,6); endendB3=B2;for i=1:l p=B3(i,2); q=B3(i,3); if B3(i,6)=1 %含有变压器支路支路 Y(p,p)=Y(p,p)+1./(B3(i,8)+B3(i,

6、9)*(1i)*B3(i,7)2); Y(p,q)=Y(p,q)-1./(B3(i,8)+B3(i,9)*(1i)*B3(i,7); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./(B3(i,8)+B3(i,9)*(1i); else% 无变压器支路 Y(p,p)=Y(p,p)+1./B3(i,4)+B3(i,5)/2; Y(p,q)=Y(p,q)-1./B3(i,4); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./B3(i,4)+B3(i,5)/2; endendYA2=A1;A2(:,2)=A1(:,2)/SB;A2(:,3)=A1(:,3)/SB;%

7、功率参数标幺化G=real(Y);B=imag(Y);S=A2(:,2)-A2(:,3) ;P=real(S);Q=imag(S);V=A2(:,4);delta=A2(:,5);DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q);J = jcb( G,B,V,delta,n,B3,A2,sw); ddelta, dV = xzl( n,J,DeltaS,A2 );e=1; while (max(ddelta)ac | max(dV)ac) a=0; b=0; for i=1:n switch A2(i,7) case 1 delta(i)=delta(i)+ddelt

8、a(i-b); V(i)=V(i)+dV(i-a-b); case 2 delta(i)=delta(i)+ddelta(i-b); V(i)=V(i); a=a+1; case 3 delta(i)=delta(i); V(i)=V(i); b=b+1; end end DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q); J = jcb( G,B,V,delta,n,B3,A2,sw); ddelta, dV = xzl( n,J,DeltaS,A2 ); e=e+1; end V delta e%下面求平衡节点功率 v=V.*cos(delta)+V.*sin(d

9、elta)*(1i);for j=1:n yu(j)=conj(Y(sw,j)*conj(v(j);endS(sw)=sum(yu)*v(sw);input(平衡节点的功率为);S(sw)B2(sw,1)=S(sw);%下面求解线路功率for i=1:n for j=1:n Sl(i,j)=v(i)*(conj(v(i)*conj(A2(i,6)+conj(v(i)-v(j)*conj(-Y(i,j); endendinput(线路功率为);Sl%线路上损耗的功率for i=1:n for j=1:n DertaS1(i,j)=(Sl(i,j)+Sl(j,i)/2; endendinput(线

10、路上损耗的功率为);DertaSz=sum(sum(DertaS1)4.2 计算功率不平衡量程序function DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q)%计算功率不平衡量for i=1:n if A2(i,7)=sw EP(i)=0; EQ(i)=0; for j=1:n EP(i)=EP(i)+V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j)+B(i,j)*sin(delta(i)-delta(j); EQ(i)=EQ(i)+V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j)-B(i,j)*cos(

11、delta(i)-delta(j); end P1=EP(i); Q1=EQ(i); if A2(i,7)=1 %PQ节点 p=2*i-1; DeltaS(p)=P(i)-P1; p=p+1; DeltaS(p)=Q(i)-Q1; else p=2*i-1; DeltaS(p)=P(i)-P1; p=p+1; DeltaS(p)=0; end endendDeltaS(2*sw-1)=;DeltaS(2*sw-1)=;for i=1:n if A2(i,7)=2 DeltaS(2*i)=; endendDeltaS=DeltaS;End4.3 计算雅克比矩阵程序function J = jcb

12、( G,B,V,delta,n,B3,A2,sw )%计算雅克比矩阵for i=1:n if A2(i,7)=1 %PQ节点 for j=1:n if j=i&j=sw H=V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j)-B(i,j)*cos(delta(i)-delta(j); J1=-V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j)+B(i,j)*sin(delta(i)-delta(j); N=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j)+B(i,j)*sin(delta(i)-delta(j

13、); L=V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j)-B(i,j)*cos(delta(i)-delta(j); p=2*i-1; q=2*j-1; J(p,q)=H; m=p+1; J(m,q)=J1; q=q+1; J(p,q)=N; J(m,q)=L; else if j=i&j=sw H1=0; for h=1:n H1=H1+(-V(i)*V(h)*(G(i,h)*sin(delta(i)-delta(h)-B(i,h)*cos(delta(i)-delta(h); end H=H1+V(i)*V(i)*(-B(i,i); J2=0; for h=

14、1:n J2=J2+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h)+B(i,h)*sin(delta(i)-delta(h); end J1=J2-V(i)*V(i)*G(i,i); N1=0; for h=1:n N1=N1+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h)+B(i,h)*sin(delta(i)-delta(h); end N=N1-V(i)*V(i)*G(i,i)+2*V(i)*V(i)*G(i,i); L1=0; for h=1:n L1=L1+(V(i)*V(h)*(G(i,h)*sin(delta(i)-delta(h)-B(i,h)*cos(delta(i)-delta(h); end L=L1-V(i)*V(i)*(-B(i,i)-2*V(i)*V(i)*B(i,i); p=2*i-1; q=2