用matlab电力系统潮流计算

1、题目：潮流计算与mat lab教学单位电气信息学院姓名学号年级专业电气工程及其自动化指导教师职称副教授电力系统稳态分析包括潮流讣算和静态安全分析。本文主要运用的事潮流计算，潮流计算是电力网 络设汁与运行中最基本的运算，对电力网络的各种设汁方案及并种运行方式进行潮流计算，可以得到各种 电网各石点的电压，并求得网络的潮流及网络中的各元件的电力损耗，进而求得电能损耗。本位就是运用 潮流计算具体分析，并有MATLAB仿真。关键词：电力系统潮流计算MATLABAbstractElectr ic power system steady fIow caIcuI at ion and ana lysis of

2、 the stat i c safety ana lysis Th i s paper, by means of the calculation, flow calculation i s the trend of the power network des i gn and operat i on of the most basic operations of electr ic power network, var i ous des i gn scheme and the operat i on ways to tide computation, can get a I I k inds

3、 of each node of the power gr id voltage and seek the trend of the network and the network of the components of the power Ioss, and gett i ng electric power - The standard i s to use the power f I ow ca I cu I at i on and ana lysis, the spec i f i c have MATLAB simulation.Key words: Power system; FI

4、ow calculation; MATLAB simulation1任务提出与方案论证12总体设计21. 1潮流计算等值电路22. 2建立电力系统模型 23. 3模型的调试与运行23详细设计34. 1计算前提35. 2手工计算64设计图及源程序106. 1 MATLAB 仿真 107. 2潮流计算源程序 105总结30参考文献311任务提出与方案论证潮流讣算是在给左电力系统网络结构、参数和决左系统运行状态的边界条件的情况下确左系统稳态运行状态 的一种基本方法，是电力系统规划和运营中不可缺少的一个重要组成部分。可以说，它是电力系统分析中最基本、 最重要的讣算，是系统安全、经济分析和实时控制与调度

5、的基础。常规潮流 计算的任务是根据给左的运行条件和 网路结构确左整个系统的运行状态，如各母线上的电压（幅值及 相角）、网络中的功率分布以及功率损耗等。潮 流计算的结果是电力系统稳左计算和故障分析的基础。在电力系统运行方式和规划方案的研究中，都需要进行潮 流汁算以比较运行方式或规划供电方案的可 行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态， 也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。 在系统规划设讣和 安排系统的运行方式时，采用离线潮流讣算；在电力系统运行状态的实时监控中，则采用在线 潮流计算。是电力系统研究人员长期研究的一

6、个课题。它既是对电力系统规划设计和运行方式的合理性、可靠性 及经济性进行左量分析的依据，又是电力系统静态和暂态稳左计算的基础。潮流il算经历了一个由手工到应用数字电子计算机的发展过程，现在的潮流算法都以计算机的应用为前 提用计算机进行潮流il算主要步骤在于编制计算机程序，这是一项非常复杂的工作。对系统进行潮流分析，本 文利用MATLAB中的SimpowerSystems /具箱设il电力系统，在simul ink环境卜一，不仅可以仿系 统的动态过程，还可以对系统进行稳态潮流分析。2总体设计51叩仙6与丫5©使用$所山而卜环境，可以将该系统中的发电机、变压器，线路等模型联结起来，形 成

7、电力系统仿真模拟图。在加人测量模块，并对各元件的参数进行设置后，用measurement和sink中的仪器 可以观察各元件的电压、电流、功率的大小。2. 1潮流计算等值电路10A/W123nM-25W/!OkVOkV1K平=44刖羽=121 加 /u St=0. 35V,% = 10.5o = 18. 6AW V, =89ArYN. d 1 4x63111 H<A4x63. WW；=0.161Xu = 0. 06:cosp% =0.851 OkV1 OkVYN.dH2x20HWA20. WWA% = 0.5 匕见=10.5r-IOW>0 = 15.7iB-S =73A1V /,%

8、= 0.5 V/4 = 10.5 2xi6. WHVA 63A/WV- 4 Gly 严44 阳 Y Sp、2kW /.,* = 0. 35V.% = I0. 5 Q10W1x50. W» V塞炸rX2x25A/W:50AW2/>% = 0. 55V/t = 10.5>d = 13. 2nV 甲 =631IV /uSi =0. 55K. « t % = 10.580.WIV 二 6.5 6i2x|o.wi»r4xl5A/W3x|2A/IV=0.136 x：v； =0. 136x： =0.16x. =0. 073cos 仇=0.8.Vu =0. 075co

9、sy、=0. 8A； =0.128x： =0.154 民=0.054 cosp% =0.85=0.128 x：=0.157 X.=0. 0591 cog=0.822建立电力系统模型在Simulink中按照电力系统原型选择元件进行建模。所建立的模型和建立的方法在详细设计中详述。在电力系统模型的建立工程中主要涉及到的是：元器件的选择及其参数的设程；发电机选型；变压器选 择：线路的选择；负荷模型的选择；母线选择。23模型的调试与运行建立系统模型，并设置好参数以后，就可以在Si mu link环境下进行仿真运行。运行的具体结果和分析也在详细设i 1中详述。3.1计算前提3详细设计首先是发电机的参数讣算

10、，先对5个发电厂简化为5台发电机来计算。发电机G1：r = 4x15 = 60MW=60xtan (arcco"). 8) = 45MVt7r发电机G2：P2 =4x63 = 252MW发电机G3：。2 : 252x tan (arcco<0. 85) = 1 56M Var出=3x12 = 36MW发电机G4：。3 = 36x tan (arccosO. 8) = 27 MVar片= 1x50 = 50MW发电机G5：Q = 50 x tan (arccos 0. 85) = 3 MVar其次是变电站的参数讣算，我们还是对7个变电站简化为7台变压器来讣算。变压器T1：R 如m

11、l。'二心小= 3. 450Q (16x10 yXri/"xl0=10-5x,l° 2xl0 = 79. 406Q Sn 16xI03/ %= “()+ /孟 xS” =(0. 0157+变压器T2：(双并联)阳显.346o2S；2 (20xl03)2X” _lxdlxlO 二殳业葺 W7625Q2Sn2 20x103A5o2 =2x (A/?o + j 监 xS ”二(0. 0372+ jO. 2000) MV A变压” 3：（四并联）陷二殳牟鲨xIO'Jx皂啤X1OW. 092Q7- 4 S；4 (63xl05) V %xVXg=A_ ' xIO

12、 =-0, 153G Sn Xg =gx"3IO = 4J 34Q/ %A5o6=2x (A/八 +，益 xIO)二(0.0264+ JO. 1100)M VAp 1 V%xV； “1 10.5x1102 "Xt3 = 一xxIO = -xxIO = 5. 042Q4 Sn463x10'/ %A5o3=4x (A/?o + 7-八X5J = (0. 1760 + jO. SS20)MVA 变压器T4：(双并联)&4= _Lr H=1.7250o2x. /4=l xn =39. 70300二 2ASoi = (0.0314 + J 01600 )MVA 变压器

13、T5：念 5 二4 卷 3 二0. 3680Q Xr5=4Xr3 =20. 168Q AS05=£ A503 =(0.0440+ jO. 2205) MV A变压器T6：(两个三绕组变压器并联) 163x35?,Rz - Rt6-2 二尺 6-3 =yx u () x(丽 *10 = 0. 386OV5i% 二八x 匕卜 2)% + Vv0_3> % 一匕(2_3) % =10. 75 V,2% = |xVi0_2)% + 匕(2 J % - Vi0_3> % = -0. 25 匕 3 % = * x 匕 + 匕 3)% -匕 g %】=6. 751 y %xy2 Xr6

14、, ,=-x V1A xIO = 6. 584Q 2 Sn灯 0 二殳 403x10 = 6.43gI %,xSn)二(0. 0220+ JO. I 100) A/VA2 IOxIO3必二2x (g + j 希再次是传输线参数汁算，5条传输线的具体讣算如下。根据教材查得心二0 - 210/«加心=0.4。/皿b0 =1 K%xv" 2. 8X| S/kmX.7 =一2 Sn线路L1 ：R/j =AJ)X/i =0. 21x40 = 8. 40X/ 二 X。X/J = 0.4 x40 = 16QBl =2. 8x10-6x40 = 1. 12x10-501 = - _L d

15、M =0. 6776M%" 2线路L2 ：尺八二% x 人=0. 21x130 = 27. 3C也二 X。x & 二 0.4 xl30 = 520.3/2 =/?ox/2 =2. 8x10"x130 = 3. 64x10-$02 = - - Bt2V； = -2. 2022MW/r 2线路L3：(双回路)&3 二 一 xroxl3 - 一 x0. 21 x70 = 7. 35Q 2 2X/3 = x xo x /3 = *x0. 4x70 = 14QBiy =2xbox/3=2x2. 8xI0-6x70 = 3. 92xlO-4SA(2/? =Bz, V；

16、=-2. 3716MV« r 2线路L4：R!A = ? x A = 0. 21 x 60 = 12. 6QX/4 = xo x/4 = 0.4 x60= 24Q 5/4=/?oX/4=2. 8x10-6x60 = 1.68x10 450401 (AMVar2线路L5：（双回路）x/5 C 0.21x20 = 2j 0/ X Xo X /5 =x 0. 4 X 20 = 4Q=2xZ?ox/5 =2x2. 8x10x20 =:xIO*05 二一B,?V； = -0. 0686MV« r23. 2手工计算FLR1：尸IO?/A5尸一(/?0+ jXn) = (3. 450 +

17、 J74. 406)=(0. 0285 + j0. 6562) A/VA 7n 11oS。=+ A5fi + ASoi + JAftj = (10.0442 + J Oj 142)M VA10.04422+0. 11422,S “ 二 (Ru +jXG = (8.4 + JI 6)=(0. 070+ JO. 1334)M 内Vn11UF+G402 + 452,Sc 二一 (/?, 2 + jX r2) =一一 (1.346 + ; 31.7625)=(0. 4032 + y9. 5156) MK411°=5Q-20-ASr2=60 + ;45 - 20 - 0. 4032 - J9.

18、 5156 = (39. 5968 + J35. 4844)MW FLR2Sc = S-S-25-jOLI -(4. 4826 + J35. 9144)M WS/ - 3 二J _ (R门 +评棵(0. 092 + )5. 042)=(0. 6679 + J36. 6024)M<4s； = (4. 4931 + j34. 1048)M<4FLRP4 49312 +34 10482S “二十 JXL2) = -一一步"一(27. 3+ J52)=(2.67 + J5. 0854) MV/1Sd = SG1 + S； -120-A573 - AS03 - jOl2 - ASL

19、2 = (132,9792 + j 149. 229)MV<4 3：AS7.4=-(/?y.4 + jXm)= &(1.725 + J39. 703)=(0. 1091 + )2. 5101)ME4S “ =(133.5955 + ) 149. 9956) MH4P2 +o133 59552 +149 90562S3 =&3 + jX( = " -(,<1 - (7. 35 + ) (4)= (24.51 +J46. 682) M<4Sf =+ Sd - 30 - 25 - A5“ 一 A5W - jU 一 次二(89. 945 + J130. 01

20、51)M<4FLR4：P+d502+312,曲 5 二一(Rts + jXT5) = (0. 368 + J20. 168) = (0. 1052 + J5. 7687) MM匕广110*Se = (92. 7481 + ； 133. 9937) MVA尸92.74872 +133.99372/Su =-(弘 + jX 心二(12.6+)24) = (27.654+ ;52. 674)M<4Sj = SGA + 5； -80- ASr5 - A5o5 - jOL4 - A5L4 =(34. 9449 + yl07. 3469) MVAFLR5:152焰7X(0. 306 + )6.

21、 431) = (0. 0562 + jAS2)MVASh=5 + ASn+AS" + JA (2L5 =(15.0782 + jO. 3422) MK415.07822+0. 3422,x(2. 1 + >4)=(0. 3899 + JO. 743)AfVAS.二 + AS/s + yA<2as +5 =(20. 4945 + )1. 266)MVA » 6-3 152 +37 52= 二 . x (0. 386 + >4. 34)-(0.514 + j5. 7793)MVAS u-2 20 . 65 052+0 . 54512”=-, x (0. 38

22、6 ； 0153) -(0. 1345 J 0 0533) MVASg 二=26. 336-+98. 7369- %(° 386 + J6. 584) - (3. 2905 - J56. 1256) MVAS/ 一叽 + Sg5 一 “心 一,S N S' -35 二(25. 5114+ j|94. 2)MVA讣算每一个FLR的功率分布和电压分布讣算如下:FLR1 ：PR + QX 40x1.346 + 45x31.7625岭2 V.115N%二115-必二 102 - 1030 如12. 8970Z:VFLR2：功率分布:PR + QX 10. 0442x8.4 + 0.

23、1442x16 八比】二=0. 8489JIV%102. 1030匕匕严101.2541余申HE “Sd+JZ12+Z73Zu+Zd27. 392-J57. 042xVvZn+ZraZ”+Zt3(27.3 )52) x (132.9792 + /49. 229)-1418.67(0.092 J 5042) x (132.9792 + J149.229)-1418. 67= (108. 687-yl22. 62)ME4=(4. 8812-yl3. 8097)MV<4电压分布:Srl =SL2+A5£ 2 = (4. 8812-y 13. 8097) + (2. 67+ J5. 0

24、854) = (7. 5512 - J8. 7243) ME47 5512x27 3 一 8 7243x52二-2. 424ZrVFLR3：102. 1030KfK- yL2 = 102. 103 - (-2. 424)= 104. 527W功率分布:r _Z.4 q _(%3-%)Zu+Zg Zu+Zn= (59. 444-yl6. 846)MV<4(1.725 -)39.703) x (89.945 + /BOOISD-iOS7. 927s u = -1- » + x VvZ £3+Z? 4 Z £3 +Zr 4二(31.811 + J 60 1256)

25、ME4电压分布：(7.35 -)14) x (89.945 + J 130.0151) +1037. 927(0. 368 一 J20. 168) x (34. 9449 + j|07. 3469)-1037. 927Sn =+= (59. 444 - j|9. 846) +83. 954x7. 35 + 26. 836x14刁40f105.5643K=Vrf-AVL3=96. 16”功率分布：Z?3Zla+Zts=(20. 843 + jl 9. 689)(24. 51 + ；46. 682) = (83. 954 + 726. 836) MV<4FLR4：(12.6 -)24) x

26、(34.9449 + J 107.3469) -1037. 927Zl4+Zt*5Zll+Z -= (1.398 + ; 44. 389) AfVA电压分布:» =S/J +ASH =(59.444+ jl 6. 846) + (24. 51 + >46. 682)=(83. 954+ j63. 528)A/VA83 - 954x12 6 + 63 -528x24=24. 4642v-AV/81. 10J105.5643FLR5:这里我们先将千点和发电机G5当做电源，经过Z ”i和Z ”3构成两端供电网络以S点作为运算负荷进行计算。(0.386 - J 4134) x (20

27、2656 + /70. 9293) + (22.0938 -37) x 35(3.900+ ； 25J 75)MVA.(0. 386 - ) 6584) x (2山。羽293) -(22.0938 - 37)S4 X 350. 772-J10. 718二(16. 5061 + ” 7905)MM4电压分布：631 = 763 + A5r63 =(16. 6421 + ;97. 5698) MVA吟 63 =16. 6421x0. 386, 97. 5698x4. 134 =期向=37匕=37 A%63 = 26 - 0814V20. 2656 x 0. 386 + 70. 9293 x (-0

28、 »153)-0. 1 62kV匕二匕一吟 62 = 26. 197626. 081420.4945x2. 1 + 1. 1266x4 = i8i526. 1976% 二匕一 匕 5 二 24. 38264设计图及源程序4. 1 MATLAB 仿真相关的原始数据输入格式如下：1、B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写。2、对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗 卷于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为含变压器支路的变压 器的变比，第六列为变压器是否是否含有变压器的参数

29、，其中"1”为含有变压器，“0”为不含有变压器。3、B2为节点参数矩阵，其中第一列为节点注入发电功率参数：第二列为节点负荷功率参数；第三列为节 点电压参数;第六列为节点类型参数，其中“1”为平衡节点，"2”为PQ节点，"3”为PV节点参数。4、X为VF点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地参数。4. 2潮流计算源程序%本程序的功能是用牛顿一一拉夫逊法进行11节点潮流计算clear n=l I ;%input （*W 输入节点数：n='）;nl = l I ;%input （ W输入支路数:nl=*）;isb=l ;%input （

30、9;请输入平衡母线节点号：isb）;pr=O. 0000I ;%input （'请输入误差精度：pr=*）;Bl = l20.03512+0.083061230.0068+0.183751140.05620+0.132891450.00811+0.245491160.05620+0.132891460.04215+0.099671670.0068+0.183751680.02810+0.066451140.05620+0.1328910.1345511001.0238110. 0538211001.0238110. 053821100. 04037110:01.0238110.1076

31、4110:0. 05382110end%i节点无功补偿量P=real;Q=lmag(S) ;%分解出各节点注入的有功和无功功率ICTI=O; IT2=I ;N0=2*n;N=N0+l ;a=0; %迭代次数 ICTI.a；不满足收敛要求的节点数 IT2 while IT2-=0 %加二2"雅可比矢巨阵的阶数；N=N0+1扩展歹IJIT2=0;a=a+l; for i=l:n i f i v i sb%非平衡节点C =0;D =0; for jl=l:nC(i)=C(l)+G(IJI)*eai)-B(IJI)*f(JI) ;% S(GU*eJ-Bi j*f j) D(i)=D(i)+G

32、(iJI)*f(jl)+B(iJI)*e(JI) ； %2(Gi J*f j+Bi J*eJ) endPI=C(i)*e(i)+f(i)*D(i);% 节点功率 P 计算 ei S (GiJ*eJ-BiJ*fj)+fI S (G i J*f j+B i J*e J)QI=C(i)*f(i)-e(i)*D(i);% 节点功率 Q 计算 fiA (Gij*eJ-Bij*f|)-ei S%求i节点有功和无功功率Pl Q，的计算值V2r"+f2%电压模平方%=以下针对非PV节点来求取功率差及Jacobi矩阵元素=ifB2(i,6)=3%#PV 节点DP=P(i)-PI;%节点有功功率差DQ=

33、Q(i).QI;%节点无功功率差%=以上为除平衡节点外其它节点的功率计算% =求取Jacob i 矢巨阵二二二二二二二二二二二二二二二for Jl=l:nifjH=isb&JI-=i%非平衡节点&非对角元Xl=-G(ijl)*e(i)-B(ijl)*f(i); %(!P/de=-dQ/df X2=B(iJl)*e(i)-G(iJI)*f；% dP/df=dQ/deX3=X2;% X2=dp/df X3=dQ/deX4=-X1;% XI=dP/de X4=dQ/dfp=2*i-l;q=2*JI-l;J(P, q) =X3;J(p, N)=DQ;m=p+l;% X3=dQ/deJ(

34、P, N)=DQ节点无功功率差J(m, q)=XI;J(m, N) =DP;q=q+l;% XI =dP/deJ(m, N)=DP节点有功功率差J(P, q)=X4;J (m, q)=X2;% X4=dQ/df X2=dp/(Ifelseifjl = i&jl-=isb%非平衡节点&对角元X 1=-C(i)-G(i, i)*e(i)-B(i, i)*f(i); % dP/de X2=-D(i)+B(i, i)*e(i)-G(i, ； % dP/df X3=D(i)+B(iJ)*e(i)-G(l, i)*f(i); % dQ/de X4=-C(i)+G(i, i)*e(i)+B(

35、i, i)*f (i) ； % dQ/df P=2*i-I； q=2*JI-I ;J(p,q)=X3;J(prN)=DQ; % 扩展列 g m=p+l;J(m, q) =XI;q=q+1; J(p, q) =X4; J(m, N) =DP;% 扩展列 ZkP J(in, q)=X2;endend下面是针对PV节点来求取Jacobi矩阵的元素DP=P(i)-PI;% PV节点有功误差DV=V(i)A2-V2;% PV节点电压误差for j1=1:n If j |-=jsb&j |-=i%非平衡节点&非对角兀XI=-G(iJI)*e(D-B(iJI)*f(i) ；% dP/deX2

36、=B(iJI)*e(i)-G(iJI)*f(i) ;% dP/dfX5=0;X6=0; p=2*i-1 ；q=2*j I. I; J (p, q) =X5; J (p, N) =DV;% PV 节点电压误m=p+l; J (m, q) =XI; J (m, N) =DP;q=q+1; J (p, q) =X6;% PV 节点有功误差J(m, q)=X2;elseifj 1 = i&JI-=isb%非平衡节点&对角元，丑上cb HTQ% PV节点有功Xl=-C(i)-G(i, l)*e(i)-B(l, i)*f(i); % dP/de X2=-D+B(i, i)*e(i)-G(i

37、, i)*f ；% dP/df X5=-2*e(i); X6=-2*f(i); / 一 / m=p+l;for k=3:NO 节点)kl=k+l;NI=N;AUfor k2=kl:NIJ(m, q) =XI; J(m, N)=DP;q=q+l; J(p, q) =X6; 误差 J(m, q)=X2; encl end end end end 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差% N0=2*n (从第三行开始，第一、二行是平衡% N=N0+1 即 N=2*n+1 扩展列 ZkP、AQ 或%从k+1列的Jacobi元素到扩展列的AQJ(k,k2)=J(k,k2”(J(k,k)+e

38、ps);%用K彳亍K列对角元素去除K彳亍K列后 的非对角 元素进行规格化endJ(k, k) = l;%对角元规格化K行K列对角元素赋1%=回代运算ifk二3%不是第三行k>3k4=k-l;for k3=3:k4%用k3行从第三行开始到当前行的前一行R4行消去for k2=kl:NI%k3行后各行上三角元素J(k3, k2)=J(k3, k2)-J(k3, k)*J(k, k2) ;%消去运算(当前行k列元素消 为0) end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素 J(k3,k)=0; %当前行第k列元素已消为0 end if k=N0%若已到最后一行break;e

39、nd%=前 代 运 算for k3=k I: NO % 从 k+1 行到 2*n 最后一行for k2=kl:NI %从k+1列到扩展列消去k+1行后各行下三角元素 J(k3, k2)=J(k3, k2) - J (k3, k)*J(k, k2);% ；肖去运算end %用当前行K2列元素减去当前行k列元素乘以第k行K2列 元素、J(k3, k)=0; %当前行第k列元素己消为0end else%是第二行k=3% =第三行k=3的前代运算for k3=kl:N0 %从第四行到2n行(最后一行)for k2=kl:NI %从第四列到2n+l列 (即扩展列)J(k3, k2)=J(k3, k2)-

40、 J(k3, k)*J(k,k2);%消去运算(当前行 3 列元素 消为0)end %用当前行K2列元素减去当前行3列元素乘以第三行K2列 元素 J(k3, k)=0; %当前行第3列元素已消为0 end%=上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵 for k=3:2:N0-lL=(k+l)J2;%修改节点电压实部e(L)=e(L)-J(k, N);kl=k+l; f(L)=f(L)-J(kl,N);%修改节点电压虚部end%修改节点电压for k=3:N0DET=abs (J (k, N);if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足

41、要求的节点数加1 end endICT2 (a) = IT2;%不满足要求的节点数ICT1 = ICT1+11%迭代次数end%用高斯消去法解* w=-J*V dispC迭代次数：);disp(ICTI);dispr没有达到精度要求的个数：，);disp(ICT2); for k=l:n V(k)=sqrt(e(k)2+f(k)"2);%计算各节点电压的模值s ida(k) =atan (f (k). /e (k) *180yp i;% 计算各节点电压的角度E (k) =e (k) +f (k) * J;%将各节点电压用复数表示end%=计算各输出量=d i spr各节点的实际电压

42、标幺值E为(节点号从小到大排列)； disp(E);%显示各节点的实际电压标幺值E用复数 表示dispC');dispr各节点的电压大小V为(节点号从小到大排列)； disp(V);%显示各节点的电压大小V的模值dispC *);disp(,各节点的电压相角sida为(节点号从小到大排列)D ; disp(sida); %显示各节点的 电压相角for p=l:nC (p) =0 ； for q=l:nC(p) =C(p) +conJ(Y(p, q)*conj (E(q); %计算各节点的注入电流的共辄值end S(p)=E(p)*C(p) ;%计算各节点的功率S二电压X注入电流的共轨值

43、end dispr各节点的功率S为(节点号从小到大排列)：J； disp；显示各节点的注入功率 dispC');dispr各条支路的首端功率S1为(顺序同您输入B1时一致)：，)；for i = l:nl p=BI(i, l)；q=BI(i,2);if Bl (i, 6)=0S i(pr q) =E (p) *(conj (E (p) *conJ (Bl (i,4)J2) + (conj (E(p)*BI (ir5). -conJ(E(q)*conJ(IV(BI (i, 3)*BI (i, 5) ; Siz(i)=Si (p, q); elseSi (p, q) =E (p) * (c

44、on j (E(p)*conj (Bl (i, 4). /2) + (conj (E(p) JBI (ir 5). conj (E(q)*conJ(l./(BI (ir 3)*BI (i, 5); Siz(i)=Si(pr q);enddisp(Si (p, q);SSi (p, q) =S i (p, q);ZF= S(, nuin2str (p),1, num2str(q),*)=*, num2str (SSi (p, q); disp(ZF); dispC *);enddispr各条支路的末端功率Sj为(顺序同您输入B1时一致):J； for i=l:nlp=BI(if l);q=BI(

45、i,2);if Bl (i,6)=0Sj (q, p) =E(q)*(conj (E(q)*conj (Bl (i,4)./2) + (conJ(E(q)./BI (i,5). conj (E(p)*conj(IJ(BI (i, 3)*BI (i, 5);SJy =SJ(q, p);elseSJ (q, p) =E (q) * (conj (E (q) *conJ (BI(L4)J2) + (conj (E(q)*BI (i). "conJ(E(p)*conJ(IJ(BI (13)*BI (t,5);SJy =SJ(q, p);enddisp(SJ (q, p)；SSJ (q, p)

46、 =SJ (q, p);ZF= S(nuni2str (q),f,f, num2str (p)/) =num2str (SSJ(q, p); disp(ZF); dispr *);encldispr各条支路的功率损耗DS为(顺序同您输入B1时一致)； for i=l:nlp=BI(i, l);q=BI(l,2);DS(i)=Si(p, q)+SJ(q, p);disp(DS(i);DDS =DS(i);ZF= ' DS (nuin2str (p)/5num2str (q)/) =ni i in2str (DDS(i); d i sp (ZF);dispC');end%本程序的功

47、能是用牛顿一一拉夫逊法进行1。节点潮流计算%本程序的功能是用牛顿一一拉夫逊法进行潮流计算clear;n=10;%inputC 请输入节点数：n=')；nl=10;%inputCW 输入支路数：nl i);isb=l;%inputr请输入平衡母线节点#:isb=r);pr=O. 00001; % i nput。请输入误差精度:pr=f)；BI = I0.03512+0. 0830610.134551230.0068+0.183751140.05620+0.132891450.00811+0.245491160.05620+0.13289100. 05382100. 053821460.0

48、4215+0.099671670.0068+0.183751680.02810+0.066451810 0. 00811+0. 245491890.035X2+0. 083061B2= 001,100100.343+0.21256100100.204+0.12638100100.306+0J8962I0010.5000. 343+0. 212561a0. 04037100.10764100.134551L 1.10100100I 00II L111.0238111. 02381111.0238111100000000000；1:0;1；0；0；1:0；1*01;2;2;2;2;2;2;2;3;

49、2诊Bl矩阵：1、支路首端号；2、末端号；支路阻抗；4.支路对地电纳%5、支路的变比；6.支路首端处于K侧为1.1侧为。%B2矩阵：4该节点发电机功率；2该节点负荷功率；3.节点电压初始值% 4、PV节点电压V的给定值；5.节点所接的无功补偿设备的容量6、节点分类标号为平衡节点(应为1号节点)；2为PQ节点;%3为PV节点；%inpu如请输入各节点参数形成的矩阵：B2=*);Y=zeros(n) ;e=zeros(lf n) ;f=zeros(I, n) ;V=zeros(lr n) ;sida=zeros(lF n) ;SI=zeros(nl); % % %-for i=l:n,%支路数if

50、 Bl (i,6)=0 p=BI(i, I) ;q=BI (i, 2);%左节点处于 1 侧else p=BI (ir2) ;q=BI (i, I);6nd%左节点处干K侧Y(p, q)=Y(pr q)-l./(BI(ir3)*BI(i,5); %非对角元Y(q,p)=Y (p,q)；% 非对角元Y(qfq)=Y(q, q) + l /(Bl (i, 3)*BI (i, 5) A 2)+BI (i, 4) /2; %对角元 K 侧Y(p, p)=Y(pr p) + l./BI(i3)+BI(i4)./2；% 对角元 1 侧%求导纳矩阵 dispr导纳矩阵Y-)；d i sp (Y) %G=re

51、al(Y );B=imag(Y); for i = l :n e(i)=real (B2(i, 3)； f (i) = i inag(B2(i, 3);V(i)=B2(l,4);end for i=l:nS(i)=B2(l, l)-B2(i,2); B(i, i)=B(i, i)+B2(ir5);%分解出导纳阵的实部和虚部%给定各节点初始电压的实部和虚部% PV节点电压给定模值%给定各节点注入功率%i 节点注入功率SG-SL %1节 点无功补偿量endP=real;Q=imag(S);%分解出各节点注入的有功和无功功率I CT I =0; I T2= I; N0=2*n; N=N0+1; a=

52、0; %迭代次数 ICTI.a；不满足收敛要求的节点数 IT2 while IT2-=0 % N0=2*n雅可比矩阵的阶数；N=N0+1扩展列IT2=0;a=a+l; for i=l:n ifi=isb%非平衡节点5)=0;D(l)=0; for Jl=l:n C(i)=C(i)+G(iJI)*e(JI)-B(iJI)*f(jl);%E(Gij*ej-Bij*fj)D(i)=D(i)+G(iJI)*f(jl)+B(iJI)*e(JI) ； %S(Gi J*f j+Bi J*eJ) end PI=C(i)*e(i)+f(i)*D(i) ;% 节点功率 P 计算 ei S (Gi J*ej-Bi

53、j*fJ)+f I SQI=C(i)*f(i)-e(i)*D(i);% 节点功率 Q 计算 fl W (Gij*fj)-ei E%求i节点有功和无功功率，,Q 口的计算值V2=e(i)A2+f人2;%电压模平方%=二=以下针对非PV节点来求取功率差及Jacobi矩阵元素ifB2(i, 6)=3 %#PV 节点DP=P(i)-PI;%节点有功功率差DQ=Q.QI;%节点无功功率差%=以上为除平衡节点外其它节点的功率计算=%二二二二二二二二二二二二二二二二二求取 Jacob i 失巨阵二二二二二二二二二二二二二二二二二二二for jl=l:nifjl-=isb&JI-=i%非平衡节点&非对角元Xl= G(iJI)*e(i)-B(iJI)*f(i) ; % dP/de=-dQ/df X2=B(iJl)*e(i)-G(iJI)*f(i);% <IP/df=dQ/deX3=X2;% X2=dp/df