电力系统潮流程序设计

潮流程序的设计与开发1数据结构的设计要求：将设备铭牌参数和有名值作为原始输入数据， 潮流结果以有名值的形式输出各节点编号节点电压的幅值 节点电压的相角平衡节点节点类型PQ 节点PV 节点输入数据支路编号支路两端节点编号支路数据支路阻抗支路对地导纳变比 短路损耗变压器数据短路电压百分数空载损耗 空载电流百分数支路数据与节点数据分别以一个矩阵的形式输入，矩阵的每一行表示每一个 节点或每一条支路，矩阵的每一列表示不同的参数数据。

各节点电压幅值 各节点电压相角平衡点功率 各支路功率节点数据注入发电功率负荷功率电压输出数据功率2变量命名设计3程序流程图启动输入原始数据形成节点导纳矩阵设节点电压初值计算有功功率和无功功率的不平衡量解修正方程式，得 到各电压的修正量V 判断是否收敛Y计算平衡节点功率 和线路功率4程序代码4.1主程序n=input （'请输入节点数：n='）; l=input （'请输入支路数：1='）; %支路数不要包括三绕组变压器sw=input （'请输入平衡节点号：sw='）; ac=input （'请输入误差精度：ac='）;停止SB=input（'请输入基准功率:SB='）;B仁input（'请输入支路参数:B仁'）%支路参数节点参数和对地支路参数均以矩阵形式储存%第一列储存支路编号%第二列与第三列分别储存支路的两个端点，分别为p,q%第四列储存支路阻抗%第五列储存支路对地导纳，注意对地导纳不要除以2%第六列储存该支路是否含有变压器，有为1，无为0%第七列储存变压器变比k,p指向q的变压器变比为k:1，且k大于等于1 %第八列储存变压器短路损耗%第九列储存变压器短路电压百分数%第十列储存变压器空载损耗%第十一列储存变压器空载电流百分数%第十二列储存变压器低压侧电压%第十三列储存变压器额定功率%第十四列储存归算所取基准电压%注意，将三绕组变压器转换为双绕组变压器输入A1=input（'请输入节点参数:A仁'）；%第一列为节点编号%第二列为注入发电功率%第三列为负荷功率%第四列为节点电压幅值，为方便起见，以标幺值形式表示%第五列为节点电压相角%第六列储存节点对地导纳%第七列为节点的类型，1为PQ节点，2为PV节点，3为平衡节点%%%首先求解双绕组变压器参数for i=1:if B1(i,6)==1Zt(i)=B1(i,8)*B1(i,12F2心000*B1(i,13F2);Xt(i)=B1(i,9)*B1(i,12)A2/(100*B1(i,13));Gt(i)=B1(i,10)/(1000*B1(i,12)A2);Bt(i)=B1(i,11)*B1(i,13)/(100*B1(i,12)A2);endendB2(:,1:7)=B1(:,1:7);for i=1:lB2(i,8)=Zt(i);B2(i,9)=Xt(i);endfor i=1:lif B1(i,6)~=1ZB(i)=B1(i,14)A2/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:lif B1(i,6)==1ZB(i)=B1(i,14)A2/SB;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: nif A1(i,6)~=0丫(i,i)=A1(i,6);endendB3=B2;for i=1:lP=B3(i,2);q=B3(i,3);if B3(i,6)==1 %%%含有变压器支路支路丫(p,p)=Y(p,p)+1./((B3(i,8)+B3(i,9)*(1i))*B3(i,7F2);丫(p,q)=Y(p,q)-1./((B3(i,8)+B3(i,9)*(1i))*B3(i,7));丫(q,p)=Y(p,q);丫(q,q)=Y(q,q)+1./((B3(i,8)+B3(i,9)*(1i)));else%%%无变压器支路丫(p,p)=Y(p,p)+1./B3(i,4)+B3(i,5)/2;丫(p,q)=Y(p,q)-1./B3(i,4);Y(q,p)=Y(p,q);丫(q,q)=Y(q,q)+1./B3(i,4)+B3(i,5)/2;endendA2=A1;A2(:,2)=A1(:,2)/SB;A2(:,3)=A1(:,3)/SB;%功率参数标幺化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: nswitch A2(i,7)case 1delta(i)=delta(i)+ddelta(i-b);V(i)=V(i)+dV(i-a-b);case 2delta(i)=delta(i)+ddelta(i-b);V(i)=V(i);a=a+1;case 3delta(i)=delta(i);V(i)=V(i);b=b+1;endendDeltaS= 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;endVdeltae%%%%下面求平衡节点功率v=V.*cos(delta)+V.*si n(delta)*(1i);for j=1: nyu(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: nfor j=1: nSl(i,j)=v(i)*(conj(v(i))*conj(A2(i,6))+conj(v(i)-v(j))*conj(-Y(i,j))); end endinput('线路功率为');SI%%%%线路上损耗的功率for i=1: nfor j=1: nDertaS1(i,j)=(SI(i,j)+SI(j,i))/2;endendinput('线路上损耗的功率为');DertaSz=sum(sum(DertaS1))4.2计算功率不平衡量程序function DeltaS= bphl(n,sw,A2,V,G,delta,B,P,Q)%计算功率不平衡量for i=1: nif A2(i,7)~=swEP(i)=0;EQ(i)=0;for j=1: nEP(i)=EP(i)+V(i)*V(j)*(G(i,j)*cos(delta(i)- delta(j))+B(i,j)*si n( delta(i)-delta(j)));EQ(i)=EQ(i)+V(i)*V(j)*(G(i,j)*si n( delta(i)-delta(j))- B(i,j)*cos(delta(i)-delta(j)));endP仁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; elsep=2*i-1;DeltaS(p)=P(i)-P1; p=p+1;DeltaS(p)=0;endendendDeltaS(2*sw-1)=[];DeltaS(2*sw-1)=[];for i=1: nif A2(i,7)==2DeltaS(2*i)=[];endendDeltaS=DeltaS';End4.3计算雅克比矩阵程序function [ J ] = jcb( G,B,V,delta,n,B3,A2,sw )%计算雅克比矩阵for i=1:nif A2(i,7)==1 %PQ 节点for j=1: nif j〜=i&j 〜=swH=V(i)*V(j)*(G(i,j)*si n( delta(i)-delta(j))-B(i,j)*cos(delta(i)- delta(j)));J仁-V(i)*V(j)*(G(i,j)*cos(delta(i)- delta(j))+B(i,j)*si n( delta(i)-delta(j)));N=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*si n( delta(i)-delta(j)));L=V(i)*V(j)*(G(i,j)*si n( 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~=swH仁0;for h=1: nH仁H1+(-V(i)*V(h)*(G(i,h)*si n( delta(i)-delta(h))- B(i,h)*cos(delta(i)-delta(h))));endH=H1+V(i)*V(i)*(-B(i,i));J2=0;for h=1: nJ2=J2+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h))+B(i,h)*si n( delta(i)-delta(h))));endJ仁J2-V(i)*V(i)*G(i,i);N仁0;for h=1: nN仁N1+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h))+B(i,h)*si n( delta(i)-delta(h))));endN=N1-V(i)*V(i)*G(i,i)+2*V(i)*V(i)*G(i,i);L仁0;for h=1: nL仁L1+(V(i)*V(h)*(G(i,h)*si n( delta(i)-delta(h))- B(i,h)*cos(delta(i)-delta(h))));endL=L1-V(i)*V(i)*(-B(i,i))-2*V(i)*V(i)*B(i,i);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;endendendendif A2(i,7)==2 %PV 节点for j=1: nif j〜=i&j 〜=swH=V(i)*V(j)*(G(i,j)*si n( delta(i)-delta(j))-B(i,j)*cos(delta(i)- delta(j)));J仁0;N=V(i)*V(j)*(G(i,j)*cos(delta(i)- delta(j))+B(i,j)*si n( delta(i)-delta(j)));L=0;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~=swH仁0;for h=1: nH仁H1+(-V(i)*V(h)*(G(i,h)*si n( delta(i)-delta(h))- B(i,h)*cos(delta(i)-delta(h))));endH=H1+V(i)*V(i)*(-B(i,i));J仁0;N仁0;for h=1: nN仁N1+(V(i)*V(h)*(G(i,h)*cos(delta(i)-delta(h))+B(i,h)*si n( delta(i)-delta(h))));endN=N1-V(i)*V(i)*G(i,i)+2*V(i)*V(i)*G(i,i);L=0;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;endendendendenda=0;for i=1: nswitch A2(i,7)case 1 case 2a=a+1;J(2*i,：)=[]；J(：,2*i)=[]；case 3J(2*i-1-a,:)=[];J(2*i-1-a,:)=[];J(:,2*i-1-a)=[];J(:,2*i-1-a)=[];endendend4.4计算电压修正量的程序function [ ddelta, dV ] = xzl( n,JQeltaS ,A2)%求取电压不平衡量dX=J\DeltaS;a=2;for i=1: nswitch A2(i,7)case 1ddelta(i)=dX(a-1);dV(i)=dX(a);a=a-0;a=a+2;case 2ddelta(i)=dX(a-1);dV(i)=0;a=a-1;a=a+2;case 3ddelta(i)=0;dV(i)=0;a=a-2;a=a+2;endenda=0;for i=1: nswitch A2(i,7)case 1case 2a=a+1; dV(i)=[];case 3ddelta(i)=[]; dV(i-a)=[];endendend5程序调试首先，在MATLAB命令窗口输入主函数名“ shuangraozU运行程序，结果如下：E L I 3.1£LL1azzu幷■■去留料：1=- iE \TF T 5-耳；、TilHA?決壬年應-沈aotn 吐圻■/星电5昌1上“化| /* iRiV 去洱驶穷汀：输入支路参数矩阵如下：[1 2 1 0 0 1 1 665.2 14 189.6 0.299 15.75 246 15.75;2 23 2.12+17.2i 1.176e-4i 0 0 0 0 0 0 0 0 220; 3 2 6 2.65+21.5i 1.47e-4i 0 0 0 0 0 0 0 0 220;4 3 6 1.59+12.9i 8.82e-5i 0 0 0 0 0 0 0 0 220;5 3 7 0 0 1 1 395.65 13.85 74.75 0.0444 10.5 180 10.5;6 47 0 0 1 1 243.25 -0.95 0 0 10.5 180 10.5 ; 7 5 7 0 0 1 1 256.75 9.15 0 0 10.5 180 10.5]随后输入节点参数矩阵如下：[1 150 0 1 0 0 2; 2 0 100+100i 1 0 0 1; 3 0 0 1 0 0 1; 4 0 100+80i 1 0 0 1; 5 0 20+15i 1 0 0 1; 6 0 0 1 0 0 3; 7 0 0 1 0 0 1]程序运行计算可得节点导纳矩阵如下:可见与手算潮流例题中所得的节点导纳矩阵基本一致 计算初次功率不平衡量矩阵如下：DsltaS =1.4M1-1. fl(MO-th 2MO-O.OD34 - 0. 17571OOM + 0.17551 doaoo + o, coon o.om * a. ooovi O.QQ« + Q.COOH O-OOQQ + Q +«IMi 也QQW +山羽叽7 MM + UMOr 084S - 0B <"3I I ■佻帥桓* 0^27721 0. 0D« * 0 OOCWi CL 0QW * 0. CQCIdi + Q,Z217i CL 30-X!會 0. ODOOiCK DOCC + 0. OQMi HL 1342 * 0. 27Tti4 99 IS - G, 77»1 6 帕朝 t 0, DOWi O. VOVQ * 0. OQ«i -a曲黛*①jswi f DOZl * g.32K10. ■ 1)1：o. )00J + a. c-oo: ：山 ooocQ.CDOD + Q.OOQO L -Or 021J + &S9L ：i山owe OrCOfrC + o.ooQai7,帖悯+■ a M 的i -0- «oaiO.CW * 0.000010. 00« 中 0. «00«iQ.OQJI - D* lifiTL o. oon + a OOOIIL -o.raiQrQMQ +h QQQgQ,芮即-a 5BQSi-O^OD J I + Ou igfiTxQ . ooan + m ODaaxa. 26930.OD041 6.. 00041 0« Lif9Hi 0.O L OOOO +D.OOOC 吐 L0M9 + O-rlOWi 血 T U B WKL D/M2 + J. tS7?i OU DOM + D. DOW! + Q. KWi -03«Z - L 357^i计算第一次迭代的雅克比矩阵如下:1 ■- ES4B-i7. ms-4. 3)930 0000◎047 二584 P r?・ tfifli-37-TIBM-3,4]®S0000>aG.33S3-乩4*B9«7.3)003F4tK000Q000-a 7B T1B4-3,416^77.«9€ B. 177S C00◎-12. M3l-0-20120X410S»r,7i*4-8-177#7T.&6M C0000. 2062-12, S93i00a冲2S.41830me. 7L»00a0d亠如4⑻~L«B I47»60尊ISB4?]WQ09Q打0Q-19. «7<■ELM 强Q0QI000D-a MW LB<5«740. »w-IB. 657 Xa-I2r W•丸-Q«20€33B&-TLK-26.41«-山3H6-LI3. ass00Q a QI Q-：05：-]2+W3l135.7155 a me-IB.W74-Wi 9JIL-153+C5S 计算第一次迭代的电压修正量矩阵如下:dde It a =0.0374 QB幽-0. »I55 -d. IO«'OULIP ・0UMT曲=-*. 0042 7 -e. (H74 -0」<1監-0. ](JW通过收敛条件判断电压修正量是否满足精度要求，若满足，则输出结果，若不满足，继续迭代。