两机五节点网络潮流计算方法牛拉法和pq法_电力系统稳态分析课程设计

1、电力系统稳态分析课程设计题目名称两机五节点网络潮流计算方法牛拉法和pq法学生姓名学 号系、专业电气工程系电气工程及其自动化指导教师2012年12月28日第一章原理3摘要简介1.1 对潮流分析地简介1.1.1 潮流计算方法分析比较 .31.2 MATLAB 简介 41.2.1 矩阵地运算51.3 牛顿拉夫逊法计算潮流分布 .6第二章 程序及结果 .102.1 设计资料及参数 .102.1.1 牛顿拉夫逊法地程序框图 132.2 用Matlab设计程序 .142.2.1 程序地编写 142.2.2 程序运行结果 .192.2. 3p_q法程序编写 22总结 3参考文献 .32电力系统稳态分析课程设

2、计1.1 对潮流分析地简介潮流分析是研究电力系统地一种最基本和最重要地计算.最初，电力系统潮流计算是通过人工手算地，后来为了适应电力系统日益发展地需要，采用了交流计算台.随着电子数字计算机地出现，1956年Ward等人编制了实际可行地计算机潮流计算程序.这样，就为日趋复杂地大规模电力系统提供了极其有力地计算手段.经过几十年地时间，电力系统潮流计算已经发展得十分成熟.潮流计算是研究电力系统稳态运行情况地一种计算，是根据给定地运行条件及系统接线情况确定整个电力系统各个部分地运行状态，如各母线地电压、各元件中流过地功率、系统地功率损耗等等.电力系统潮流计算是计算系统动态稳定和静态稳定地基础.在电力系

3、统规划设计和现有电力系统运行方式地研究中，都需要利用电力系统潮流计算来定量地比较供电方案或运行方式地合理性、可靠性和经济性.1.1.1 潮流计算方法分析比较高斯 - 赛德尔潮流计算法原理简单，编程实现容易，特别是对于配网潮流有其独特优势 .但是高斯- 赛德尔潮流计算法在牛顿法以及各种解耦法出现以后似乎成了一种边缘性地方法.牛顿- 拉夫逊法地优点是收敛速度快，若初值选择较好，算法将具有平方收敛特性，一般迭代45 次便可以收敛到一个非常精确地解，而且其迭代次数与所计算地网络规模基本无关.牛顿 - 拉夫逊法也具有良好地收敛可靠性对于呈病态地系统，牛顿-拉夫逊法均能可靠地收敛.牛顿法地缺点是每次迭代地

4、计算量和所需地内存量较大.这是因为雅可比阵元素地数目约为2(n- 1)>2(n- 1)个(直角坐标)，且其数值在迭代过程中不断变化.不过，内存占用量及每次迭代所需地时间与程序设计技巧密切相关.牛顿-拉夫逊法地可靠收敛取决于一个良好地启动初值，如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行地解点上.对于正常运行地系统，各节电电压一般均在额定值附近，偏移不会太大，并且各节电地相角差也不大，所以对各节电可以采用统一地电压初值.P- Q 分解法是为了改进牛顿-拉夫逊法在内存占用量及计算速度方面地不足，P- Q 分解法根据电力系统实际运行状态地物理特点，对极坐标形式地牛顿- 拉夫逊法修

5、正方程式进行了合理地简化它无论在内存占用量还是计算速度方面都比牛顿-拉夫逊法有较大地改进，是目前计算速度最快地潮流算法.1.2 Matlab 地简介MA TLAB 是一种交互式、面向对象地程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大地功能.MATLAB 程序设计语言结构完整，且具有优良地移植性，它地基本数据元素是不需要定义地数组.它可以高效率地解决工业计算问题，特别是关于矩阵和矢量地计算.MATLAB 与 C 语言和 FORTRAN 语言相比更容易被掌握.通过M 语言，可以用类似数学公式地方式来编写算法，大

6、大降低了程序所需地难度并节省了时间，从而可把主要地精力集中在算法地构思而不是编程上.另外， MATLAB 提供了一种特殊地工具：工具箱( TOOLBOXES ) .这些工具箱主要包括：信号处理(SIGNAL PROCESSING )、控制 系 统 ( CONTROL SYSTEMS ) 、 神 经 网络 ( NEURAL NETWORKS ) 、 模 糊 逻 辑(FUZZY LOGIC) 、小波 (WAVELETS) 和模拟(SIMULATION )等等.不同领域、不同层次地用户通过相应工具地学习和应用，可以方便地进行计算、分析及设计工作.MATLAB 设计中，原始数据地填写格式是很关键地一个

7、环节，它与程序使用地方便性和灵活性有着直接地关系. 原始数据输入格式地设计，主要应从使用地角度出发，原则是简单明了，便于修改 .1.2.1 矩阵地运算.但需注意进行数除时，常数通常数与矩阵地运算即是同该矩阵地每一元素进行运算 常只能做除数个：基本函数运算中，矩阵地函数运算是矩阵运算中最实用地部分，常用地主要有以下几det(a)eig(a)求矩阵 a 地行列式 求矩阵 a 地特征值inv(a)或 a A (-1)求矩阵a地逆矩阵rank（a）求矩阵a 地秩trace（a）求矩阵a地迹（对角线元素之和）我们在进行工程计算时常常遇到矩阵对应元素之间地运算.这种运算不同于前面讲地数学运算，为有所区别，

8、我们称之为数组运算.数组地加、减与矩阵地加、减运算完全相同 . 而乘除法运算有相当大地区别，数组地乘除法是指两同维数组对应元素之间地乘除法，它们地运算符为“ .*和 ” “ ./或 ”“. ”. 前面讲过常数与矩阵地除法运算中常数只能做除数.在数组运算中有了“对应关系”地规定，数组与常数之间地除法运算没有任何限制.另外，矩阵地数组运算中还有哥运算（运算符为.a ）、指数运算（exp）、对数运算（log）、和开方运算（sqrt）等.有了对应元素”地规定，数组地运算实质上就是针对数组内部地每个元素进行地 .矩阵地幂运算与数组地幂运算有很大地区别.1.3 牛顿拉夫逊法计算潮流分布一、牛顿拉夫逊法求解

9、过程大致可以分为以下步骤：（ 1 ）形成节点导纳矩阵（ 2）将各节点电压设初值U3）将节点初值代入相关求式，求出修正方程式地常数项向量(4)将节点电压初值代入求式，求出雅可比矩阵元素(5)求解修正方程，求修正向量(6)求取节点电压地新值(7)检查是否收敛，如不收敛，则以各节点电压地新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步(8)计算支路功率分布，PV节点无功功率和平衡节点柱入功率.二、直角坐标系计算(1)牛顿-拉夫逊法潮流方t算地公式.把牛顿法用于潮流计算，采用直角坐标形式.其中电压和支路导纳可表示为：VijfiY =G jBYij GijjBij(1-1)(2)PQ节点地输出

10、有功功率和无功功率是给定地，则第 i节点地给定功率设为R和Qis (称为注入功率).假定系统中地第1、2、m节点为PQ节点，对其中每一个节点地N-R法表达式： nnP =Ps-P =Ps-6(Gjej-Bjfj)-f(GijfjBjej)= 0j 1j 1nnQi-Qis-Qi=Qis(Gjej-Bjfj) e' (GijfjBjej)=0j 1ji = (1、2、m)11-2)(3)PV节点地有功功率和节点电压幅值是给定地.假定系统中地第m+1、m+2、n-1节点为PV节点，则对其中每一 PV节点可以列写方程：nn(1-P = Ps-P = Ps-ei£ (Gjjej-Bj

11、fj)-f(Gj fj + 8科)=0| j ijW22222_2u2=u；-u；=u；-(e2+ fj).3)i = m m+1、m+2、n-1)(4)形成雅可比矩阵.当j=i时，对角元素为；：P : :e:PQi :e:Qi"U2ef U2fin-(Gjej -Bijfj)-Ghe -fijn，(Gijfj -Bij。) Biie -GiifiHiijn% (Gijfj -Bij0)Biie -Giifi = Liij in-(Gj -Bjfj) GheB"i = Jj i一 2e一 2fi(1-4)当j #i时,矩阵非对角元素为:P丁f Ui2:ejfF Qi:ej=

12、-(Gije + Bij fi) = Nij = -Jij=Bij ei _ Gij fi = H ij = Lij工0:fj(1-5)三、极坐标计算对于潮流计算中待于求出功率地节点功率方程组，在某个近似解附近用泰勒级数展开略去二阶及以上地高阶得到已矩阵表示地修正方程:式中节点n为节点数，m为PV节点，雅阁比矩阵是（2n-2m-2）阶非奇异矩阵，雅可 比矩阵各元素表示如下:H) = WUgj sin % - 得 cos% )C t/ jnv%=俨/ =-"与(与85% +%sm%) UU INi = TUi = Ui Gii - Pi dUi= 5Uj(Gy cosOy + By s

13、in%)de；句=翳巧=-34（的sin的-*户椭）Lii = u； B/i - QidUi牛顿拉夫逊极坐标潮流计算地修正方程地迭代方程为:HgMl«GNnPQAMm-141MnPQZ LnpQaNnPQAnPQ L“tPQ eLrPQ,RPQ2.1设计资料及参数:课程名称电力系统稳态分析课程设计设计题目两机五节点网络潮流计算一牛拉法指导教师刘景霞时间 1周一、教案要求电力系统稳态分析课程设计以设计和优化电力系统地潮流分析为重点，提高学生综合能力为目标，尽可能结合实际工程进行 .设计内容地安排要充分考虑学校现有地设备，设计时间及工程实际需要，并使学生初步学会运用所学知识解决工程中地实

14、际问题二、设计资料及参数（一）设计原始资料1、待设计电气设备系统图2、电力系统网络各元件参数3、电力系统电气元件地使用规范4、电力工程电气设计手册（二）设计参考资料1、电力系统稳态分析，陈琦，中国电力出版社，2007,第三版2、电力系统分析，韩祯祥，浙江大学出版社，2005,第三版3、电力系统分析课程实际设计与综合实验，祝书萍，中国电力出版社，2007,第一版三、设计要求及成果1 .根据给定地参数或H程具体要求，收集和查阅资料；学习相关软件（软件自选）2 .在给定地电力网络上画出等值电路图.3 .运用计算机进行潮流计算.4 .编写设计说明书.基本要求：1.编写潮流计算程序：2 .在计算机上调试

15、通过（？）；3 .运行程序并计算出正确结果（？）；4 .写出课程设计报告（包括以下内容）（1份）（1）程序框图；（2）源程序；（3）符号说明表；（4）算例及计算结果5 .编写计算说明书（1份）.四、进度安排根据给定地参数或工程具体要求，收集和查阅资料（半天）学习软件（MATLAB C语言等）（一天半）编程计算复杂系统潮流计算（三天）编写计算设计书（一天）五、评分标准课程设计成绩评定依据包括以下几点：1）工作态度（占10% ;2）基本技能地掌握程度（占 20% ;3）程序编写是否合理是否有运行结果（40%）;4）课程设计说明书编写水平（占30%）.5）分为优、良、中、合格、不合格五个等级.Q.0

16、S+j0.240.45-讪 150.4j0.050.OM0.030.04-H0J2Y0.2 jOQ0.6-j0.1系统接线图2.1其中节点1为平衡节点，节点 2、3、.1牛顿拉夫逊法程序框图4、5为PQ节点.电力系统稳态分析课程设计2.2.1 程序地编写%说明：为了使节点按照先PQ,再PV节点，最后平衡节点地次序编号，以便与公式对照，节点1 与节点 5对调 .%节点阻抗矩阵Z=0,0.04+0.12i,0,0.08+0.24i,0 。 0.04+0.12i,0,0.06+0.18i,0.06+0.18i,0.02+0.06i 。0,0.06+0.18i,0,0.01+0.03i,0.08+0.

17、24i 。 0.08+0.24i,0.06+0.18i,0.01+0.03i,0,0 。0,0.02+0.06i,0.08+0.24i,0,0 。%求互导纳for m=1:5for n=1:5if Z(m,n)=0 y(m,n)=0 。elsey(m,n)=1/Z(m,n) 。 endendend%求导纳for m=1:5for n=1:5if m=nY(m,n)=-y(m,n) 。elseY(m,n)=sum(y(m,:) 。end end end %导纳矩阵 G=real(Y) 。 B=imag(Y) 。 %计算各节点功率地不平衡量设U=E+jF ； Y=G+Bj ；YU=1,1,1,1,

18、1.06。E=real(U) 。F=imag(U) 。%设 S=P+Qj；S=-0.60-0.10i。 0.20+0.20i。 -0.45-0.15i。 -0.40-0.05i。 0P=real(S)。 Q=imag(S) 。k=0。 C=1 。while C > 0.00001E(5)=1.06。 F(5)=0。for m=1:4for n=1:5%计算Pi,Qi,设Pi=Pt。Qi=Qt,按照书上地公式：Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n)+F(m)*(G(m,n)*F(n)+B(m,n)*E(n)。Qt(n)=(F(m)*(G(m,n)*E(n)

19、-B(m,n)*F(n)-E(m)*(G(m,n)*F(n)+B(m,n)*E(n)。end %设P,Q地改变量为dP,dQ dP(m)=P(m)-sum(Pt) 。dQ(m)=Q(m)-sum(Qt) 。 end%计算Hii,Nii,Jii,Lii, 由公式 4-41b 左侧公式实现,sum(Ai) ， sum(Bi) 用于实现公式中地sigerma从j到n地求和。for m=1:4for n=1:5Bi(n)=G(m,n)*F(n)+B(m,n)*E(n) 。Ai(n)=G(m,n)*E(n)-B(m,n)*F(n) 。endH(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,

20、m)*F(m)+2*G(m,m)*F(m)。N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)+2*G(m,m)*E(m)。J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)。L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m)。end%设雅可比矩阵为JJ,以下语句用来实现雅可比矩阵中对角线上元素H N J L地排列for m=1:4JJ(2*m-1,2*m-1)=H(m,m) 。JJ(2*m-1,2*m)=N(m,m) 。JJ(2*m,2*m-1)=

21、J(m,m) 。JJ(2*m,2*m)=L(m,m) 。end%以下语句用于实现雅可比矩阵非对角线上元素地排列for m=1:4for n=1:4if m=nelseH(m,n)=-B(m,n)*E(m)+G(m,n)*F(m) 。N(m,n)=G(m,n)*E(m)+B(m,n)*F(m) 。J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m) 。L(m,n)=G(m,n)*F(m)-B(m,n)*E(m) 。JJ(2*m-1,2*n-1)=H(m,n) 。JJ(2*m-1,2*n)=N(m,n) 。JJ(2*m,2*n-1)=J(m,n) 。JJ(2*m,2*n)=L(m,n)

22、。endendendJJ%设由P,Q地改变量组成地8X1矩阵为PQ,由E,F地改变量组成地8X1矩阵为dUfor m=1:4PQ(2*m-1)=dP(m) 。PQ(2*m)=dQ(m) 。enddU=inv(JJ)*PQ' 。 %求逆矩阵C=max(abs(dU)。 %abs绝对值或复数绝对值for n=1:4F(n)=F(n)+dU(2*n-1)E(n)=E(n)+dU(2*n) 。endfor n=1:5 U(n)=E(n)+(F(n)*i 。endk=k+1 。k, dU=dU',PQ,Uend%计算S(5)，也就是题目中地S1,即平衡节点功率for m=1:5I(m)=

23、Y(5,m)*U(m) 。endS(5)=U(5)*sum(conj(I)3.7500 -11.2500i -2.5000 + 7.5000i0-1.2500 + 3.7500i0-2.5000 + 7.5000i 10.8333 -32.5000i -1.6667 + 5.0000i -1.6667 + 5.0000i -5.0000 +15.0000i0-1.6667 + 5.0000i 12.9167 -38.7500i -10.0000 +30.0000i -1.2500 + 3.7500i-1.2500 + 3.7500i -1.6667 + 5.0000i -10.0000 +30

24、.0000i 12.9167 -38.7500i00-5.0000 +15.0000i -1.2500 + 3.7500i06.2500 -18.7500iS =-0.6000 - 0.1000i0.2000 + 0.2000i-0.4500 - 0.1500i-0.4000 - 0.0500i0JJ =11.2500 3.7500 -7.5000 -2.500000 -3.7500 -1.2500-3.7500 11.2500 2.5000 -7.500000 1.2500 -3.7500-7.5000 -2.5000 33.4000 10.5333 -5.0000 -1.6667 -5.0

25、000 -1.66672.5000 -7.5000 -11.1333 31.6000 1.6667 -5.0000 1.6667 -5.00000 -5.0000 -1.6667 38.9750 12.8417 -30.0000 -10.00000 1.6667 -5.0000 -12.9917 38.5250 10.0000 -30.0000-3.7500 -1.2500 -5.0000 -1.6667 -30.0000 -10.0000 38.7500 12.91671.2500 -3.7500 1.6667 -5.0000 10.0000 -30.0000 -12.9167 38.750

26、0 dU =-0.1076 0.0093 -0.0473 0.0430 -0.0863 0.0154 -0.0922 0.0141PQ =-0.6000 -0.1000 0.5000 1.1000 -0.3750 0.0750 -0.4000 -0.0500U =1.0093 - 0.1076i 1.0430 - 0.0473i 1.0154 - 0.0863i 1.0141 - 0.0922i 1.0600JJ =11.0516 4.3956 -7.3011 -3.330400 -3.6505 -1.6652-5.5956 10.8516 3.3304 -7.301100 1.6652 -3

27、.6505-7.7040-2.962133.159413.0913-5.1360-1.9747-5.1360-1.97472.9621-7.7040-12.580433.60831.9747-5.13601.9747-5.136000 -4.9331 -2.1238 38.3848 16.0299 -29.5988 -12.742700 2.1238 -4.9331 -16.8888 38.0788 12.7427 -29.5988-3.6876-1.6135-4.9168-2.1513-29.5009-12.907938.155316.27271.6135-3.68762.1513-4.91

28、6812.9079-29.5009 -17.072738.0553dU =0.0032 -0.0131 -0.0004 -0.0075 0.0017 -0.0101 0.0021 -0.0108PQ =0.0164 -0.0636 -0.0770 -0.0220 -0.0008 -0.0317 0.0103 -0.0362U =0.9963 - 0.1044i 1.0355 - 0.0477i 1.0053 - 0.0846i 1.0033 - 0.0902i 1.0600JJ =10.9774 4.3253 -7.2110 -3.273800 -3.6055 -1.6369-5.4960 1

29、0.6556 3.2738 -7.211000 1.6369 -3.6055-7.6466-2.946632.933412.9530-5.0977-1.9644-5.0977-1.96442.9466-7.6466-12.584333.33711.9644-5.09771.9644-5.097700 -4.8857 -2.0983 38.0494 15.8299 -29.3142 -12.589800 2.0983 -4.8857 -16.6938 37.6790 12.5898 -29.3142-3.6499-1.5922-4.8665-2.1230-29.1988-12.738037.79

30、9716.06211.5922-3.64992.1230-4.866512.7380-29.1988 -16.844337.6306 dU =1.0e-003 *0.0241 -0.1731 -0.0025 -0.0892 0.0098 -0.1231 0.0128 -0.1327PQ =1.0e-003 *-0.0034 -0.8423 -0.5236 -0.2008 -0.0815 -0.3198 0.0241 -0.3906U =0.9961 - 0.1044i 1.0354 - 0.0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600JJ =10

31、.9764 4.3243 -7.2098 -3.273100 -3.6049 -1.6366-5.4951 10.6529 3.2731 -7.209800 1.6366 -3.6049-7.6459-2.946432.930712.9517-5.0973-1.9643-5.0973-1.96432.9464-7.6459-12.583933.33401.9643-5.09731.9643-5.097300 -4.8851 -2.0980 38.0451 15.8278 -29.3106 -12.588300 2.0980 -4.8851 -16.6919 37.6740 12.5883 -2

32、9.3106-3.6494-1.5920-4.8658-2.1227-29.1950-12.736337.795216.05991.5920-3.64942.1227-4.865812.7363-29.1950-16.842137.6252 dU =1.0e-007 *0.0232 -0.2889 -0.0013 -0.1352 0.0054 -0.1876 0.0087 -0.2038PQ =1.0e-006 *-0.0240 -0.1486 -0.0424 -0.0189 -0.0165 -0.0415 -0.0042 -0.05250.9961 - 0.1044i 1.0354 - 0.

33、0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600-0.6000 - 0.1000i0.2000 + 0.2000i-0.4500 - 0.1500i-0.4000 - 0.0500i1.2982 + 0.2445i223P q法程序编写基于 MATLAB 地 PQ 分解法潮流计算程序MA TLAB 已广泛应用于自动控制、数学运算、信号分析、计算机技术、图像信号处理、财务分析、航天工业和生物医学工程等领域.由于 MATLAB 语言功能强大、人际界面友好、编程效率高、强大而智能化地作图功能，且具有编程语句简洁、灵活、表达和运算能力强等显著特点.程序清单%

34、本程序地功能是用PQ 分解法进行潮流计算n=inPut( 请输入节点数： n= ) 。nl=inPut( 请输入支路数： nl= ) 。isb=inPut( 请输入平衡母线节点号： isb= ) 。Pr=inPut( 请输入误差精度 :Pr= )。B1=inPut( 请输入由支路参数形成地矩阵 :B1= )。B2=inPut( 请输入由节点参数形成地矩阵 :B2= )。X=inPut( 请输入由节点号及其对地阻抗形成地矩阵:X= ) 。na=inPut( 请输入 PQ 节点号 :na= ) 。Y=zeros(n) 。 YI= zeros(n) 。 e= zeros(1,n)。 f= zeros

35、(1,n) 。 V= zeros(1,n) 。 O= zeros(1,n)。 for i=1:nif X(i,2)=0 。P=x(i,1) 。Y(P,P)=1./X(i,2) 。end end for i=1:nlif B1(i,6)= =0p=B1(i,1) 。 q=B1(i,2)。else p= B1(i,2)。 q=B1(i,1) 。endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)。YI(p,q)= YI(p,q) -1./B1(i,3) 。Y(q,p) =Y(p,q) 。YI(q,p) = YI(p,q) 。Y(q,q)= Y(q,q)+1./(B1(i,3)*

36、 B1(i,5)2)+B1(i,4)./2YI(q,q)= YI(q,q) +1./B1(i,3) 。Y(p,p)= Y(p,p) +1./B1(i,3)+ B1(i,4)./2 。YI(p,p)= YI(p,p) +1./B1(i,3) 。end%求导纳矩阵G=real(Y) 。 B=imag(YI) 。 BI=imag(Y) 。for i=1:nS(i)=B2(i,1)-B2(i,2) 。BI(i,i)= BI(i,i)+ B2(i,5) 。endP=real(S)。 Q=imag(S) 。for i=1:ne(i)=real(B2(i,3) 。f(i)=imag(B2(i,3) 。V(i

37、)=B2(i,4) 。endfor i=1:nif B2(i,6)= =2V(i)=sqrt(e(i)2+f(i)2) 。O(i)=atan(f(i)./e(i) 。endendfor i=2:nif i= =nB(i,i)=1./B(i,i) 。else IC1=i+1 。for j1=IC1:nB(i,j1)= B(i,j1)./B(i,i) 。endB(i,i)=1./B(i,i) 。for k=i+1:nfor j1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1) endendendendp=0。 q=0。for i=1:nif B2(i,6)= =2p=p+1

38、 。 k=0 。for j1=1:nif B2(i1,6)= =2k=k+1 。A(p,k)=BI(i,j1) 。endendendendfor i=1:naif i= =naA(i,i)=1./A(i,i) 。else k=i+1 。for j1=k:naA(i,j1)= A(i,j1)./A(i,i) 。endA(i,i)=1. /A(i,i) 。for k=i+1:nafor j1=i+1:naA(k,j1)= A(k,j1)-A(k,i)*A(i,j1)。endendendendICT2=1 。 ICT1=0 。 kp=1,kq=1 。 K=1 。 DET=0。 ICT3=1 。whi

39、le ICT2=0|ICT3=0ICT2=0。 ICT3=0 。for i=1:nif i=isbC(i)=0 。for k=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k)+BI(i,k)*sin(O(i)-O(k) endDP1(i)=P(i)-V(i)*C(i) 。DP(i)=DP1(i)./V(i) 。DET=abs(DP1(i) 。if DET>=prICT2=ICT2+1 。endendendNp(K)=ICT2 。if ICT2=0for i=2:nDP(i)=B(i,i)*DP(i) 。if i =nIC1=i+1 。for k=ic1:nD

40、P(k)=DP(k)-B(k,i)*DP(i) 。endelsefor LZ=3:iL=i+3-LZ 。IC4=L-1 。for MZ=2:IC4I=IC4+2-MZ 。DP(I)=DP(I)-B(I,L)*DP(L) 。endendendendfor i=2:nO(i)=O(i)-DP(i) 。endkq=1 。 L=0 。for i=1:nif B2(i,6)= =2C(i)=0 。 L=L+1 。for k=1:nC(i)= C(i) +V(k)*(G(i,k)*sin(O(i)-O(k)+BI(i,k)*cos(O(i)-O(k) endDQ1(i)=Q(i)-V(i)*C(i) 。D

41、Q(L)=DQ1(i)./V(i) 。DET=abs(DQ1(i) 。if DET>=prICT3=ICT3+1 。endendendendelse kp=0。if kq=0 。L=0。for i=1:nif B2(i,6)= =2C(i)=0 。 L=L+1 。for k=1:nC(i)= C(i) +V(k)*(G(i,k)*sin(O(i)-O(k)-BI(i,k)*cos(O(i)-O(k) endDQ1(i)=Q(i)-V(i)*C(i) 。DQ(L)= DQ1(i)./V(i) 。DET=abs(DQ1(i) 。endendendendNq(K)=ICT3 。if ICT3=

42、0L=0。for i=1:naDQ(i)=A(i,i)*DQ(i) 。if i= =nafor LZ=2:iL=i+2-LZ 。IC4=L-1 。for MZ=1:IC4I=IC4+1-MZ 。DQ(I)=DQ(I)-A(I,L)*DQ(L) 。end endelseIC1=i+1 。for k=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i) 。endendendL=0。for i=1:nif B2(i,6)= =2L=L+1 。V(i)=V(i)-DQ(L) 。endendkp=1 。K=K+1 。elsekq=0。if kp=0K=K+1 。endendfor i=1:nDy

43、(K-1,i)=V(i) 。endenddisp( 迭代次数 ) ；disp( K)；disp( 每次没有达到精度要求地有功功率个数为)；disp( NP)；disp( 每次没有达到精度要求地无功功率个数为)；disp( Nq)。for k=1:nE(k)=V(k)*cos(O(K)+V(k)*sin(O(k)*j。O(k)=O(k)*180./pi 。enddisp( 各节点地电压标么值 E 为(节点号从小到大排)：disp(E)。disp( 各节点地电压 U 大小(节点号从小到大排)： ) 为；disp(U)。disp(各节点地电压相角O (节点号从小到大排)：'力;disp(O)

44、。for p=1:nC(p)=0。for q=1:nC(p)= C(p)+conj(Y(p,q)*conj(E(q) 。 endS(p)=E(p)*C(p) 。enddisp(各节点地功率S (节点号从小到大排)：为；disp( S)；disp( 各条支路地首端功率 Si 为(顺序同您输入B1 时地一样)：)为；for i=1:nlif B1(i,6)= =0p=B1(i,1)。 q=B1(i,2)。else p=B1(i,2) 。 q=B1(i,1) 。 endSi(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5)-conj(E

45、(q)*conj(1./(B1(i,3)*B1(i,5)。disp(Si(p,q) 。enddisp(各条支路地末端功率Sj (顺序同您输入 B1时地一样)：)为；for i=1:nlif B1(i,6)= =0p=B1(i,1)。 q=B1(i,2)。else p=B1(i,2) 。 q=B1(i,1) 。 endSj(q,p)=E(q)*(conj(E(q)*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5)-conj(E(p)*conj(1./(B1(i,3)*B1(i,5)。disp(Sj(q,p) 。enddisp(各条支路地功率损耗DS (顺序同您输入 B1时地

46、一样)：)为；for i=1:nlif B1(i,6)= =0p=B1(i,1)。 q=B1(i,2)。else p=B1(i,2) 。 q=B1(i,1) 。 endDS(i)=Si(p,q)+Sj(q,p) 。disp(DS(i) 。endfor i=1:KCs（i）=i 。for j=1:nDy（K,j）=Dy（K-1,j） 。endenddisp（ 以下是每次迭代后各节点地电压值（如图所示） ） ；plot（Cs,Dy）,xlable（迭代次数 ；）ylable（电压'）,title（电压迭彳弋次数曲线；）设计验证现以一个例子来说明该软件得使用方法:电力系统接线，如4.1 所示

47、,试求潮流分布（1） n为节点数、nl为支路数、isb为平衡母线节点号（固定为1）、pr为误差精度.（ 2）输入由支路参数形成地矩阵B1矩阵 B1 地每行是由下列参数构成地： 某支路地首端号P；末端号Q,且P <Q; 支路地阻抗（R +jX ）； 支路地对地容抗； 支路地变比K ； 折算到哪一侧地标志（如果支路地首端P 处于高压侧则请输入“ 1，”否则请输入 “ 0）”Figure 4.1 schematic（ 3）请输入各节点参数形成地矩阵B2矩阵 B2 地每行是由下列参数构成地： 节点所接发电机地功率SG； 节点负荷地功率SL； 节点电压地初始值； PU 节点电压U 地给定值； 节点所接地无功补偿设备地容量； 节点分类标号igl.输入数据：请输入节点数：n=5请输入支路数：nl=5请输入平衡母线节点号：isb=1请输入误差精度：pr=0.00001请输入由支路参数形成地矩阵：B1=1 2 0.03i 0 1.050 。 2 3 0.08+0.3i 0.5i 1 0 。 2 4 0.1+0.35i0 1 0。 3 4 0.04+0.25i 0.5i 1 0 。 3 5 0.015i 0 1.05 1请输入由节点参数形成地矩阵：B2=0 0 1.05 1.05 0 1 。 0 3.7+1.3i 1.0