(完整版)电力系统潮流计算课程设计

(完满word版)电力系统潮流计算课程设计(完满word版)电力系统潮流计算课程设计(完满word版)电力系统潮流计算课程设计一、问题重述 课程设计要求1、在读懂程序的基础上画出潮流计算根本流程图2、经过输入数据，进行潮流计算输出结果3、对不同样的负荷变化，解析潮流分布，写出解析说明。4、对不同样的负荷变化，进行潮流的调治控制，并说明调治控制的方法，并列表表示调治控制的参数变化。5、打印利用DDRTS进行潮流解析绘制的系统图，以及潮流分布图。 课程设计题目1、系统图：两个发电厂分别经过变压器和输电线路与四个变电所相连。变电所1变电所2变电所4变电所310kV母线35kV母线10kV母线35kV母线一次侧电压220kV一次侧电压220kV线路长为60km线路长为80km线路长为100km线路长为80km线路长为100km线路长为80km母线1母线2。。。。。。。。。。。。。。。。。。。。。。。。。。母线32*QFQ-50-22*QFS-50-22*TQN-100-22*TQN-100-2电厂一电厂二2、发电厂资料：母线1和2为发电厂高压母线，发电厂一总装机容量为〔400MW〕，母线3为机压母线，机压母线上装机容量为〔100MW〕，最大负荷和最小负荷分别为50MW和30MW；发电厂二总装机容量为〔200MW〕。3、变电所资料：〔1〕变电所1、2、3、4低压母线的电压等级分别为：10KV35KV10KV35KV〔2〕变电所的负荷分别为：50MW 40MW 50MW 60MW1〔3〕每个变电所的功率因数均为 cosφ=0.85；〔4〕变电所2和变电所4分别配有两台容量为75MVA的变压器，短路耗费414KW，短路电压〔%〕=16.7；变电所1和变电所3分别配有两台容量为63MVA的变压器，短路耗费为245KW，短路电压〔%〕=10.5；4、输电线路资料：发电厂和变电所之间的输电线路的电压等级及长度标于图中， 单位长度的电阻为 ，单位长度的电抗为 ，单位长度的电纳为 2.78*10-6S。 课程设计根本内容对给定的网络查找潮流计算所需的各元件等值参数，画出等值电路图。输入各支路数据，各节点数据利用给定的程序进行在变电所在某一负荷情况下的潮流计算，并对计算结果进行解析。随从变电所负荷按必然比率发生变化，进行潮流计算解析。4个变电所的负荷同时以2%的比率增大；4个变电所的负荷同时以2%的比率下降1和4号变电所的负荷同时以2%的比率下降，而2和3号变电所的负荷同时以2%的比率上升；在不同样的负荷情况下，解析潮流计算的结果，若是各母线电压不满足要求，进行电压的调整。〔变电所低压母线电压 10KV要求调整范围在 之间；电压35KV要求调整范围在35-36之间〕轮流断开环网一回线，解析潮流的分布。利用DDRTS软件，进行绘制系统图进行上述各种情况潮流的解析，并进行结果的比较。最后形成课程设计成品说明书。二、问题解析 节点设置及分类依照系统图可知此系统为两端供电网路，将母线1，2设为节点1，10，将变电所1、2、3、4的上下压侧分别设为节点2、3、4、5、6、7、8、9。并且，将节点1设为平衡节点，将节点10设为PV节点，其他节点设为PQ节点。2.2参数求取设定基准值SB100MVA，UBSB484依照题目220KV，因此ZB2UB原始资料，计算发电厂、变压器及线路的参数。〔1〕运用以下公式计算变压器参数：22PkUNRTUK%UN2XT1000S100SN2ZTRTXTZTZT〔2〕计算线路参数ZBZRjX(rjx)L〔3〕变电所负荷分别为：变电所1SL变电所2SL变电所3SL变电所4SL将参数整理，见下表：首端号尾端号阻抗有名值阻抗标幺值电纳有名值电纳标幺值12141623004500466700688900810〔4〕计算变压器分接头变比变压器有5个抽头，电压调治范围为UN2*2.5%，UN对应的分接头开始时设变压器高压侧接主接头 ,降压变压器5个分接头时的非标准变比 k*以备调压时采用k*1220(15%)1011220k*2220(12.5%)1011220k*32201011220k*4220(12.5%)1011220k*5220(15%)1011220对变电所低压母线为 35KV时，非标准变比与 10KV时同样。 计算方法利用牛顿拉夫逊法进行求解，用MATLAB软件编程，能够求解系统潮流分布依照题目的不同样要求对参数进行调整，经过调治变压器变比和发电厂的电压，3求解出合理的潮流分布，最后用 PSAT进行潮流解析，将两者进行比较。 牛顿—拉夫逊法1、牛顿—拉夫逊法大纲第一对一般的牛顿—拉夫逊法作一简单的说明。一个变量 X函数为：f(X) 0(0)到此方程时，由合适的近似值 X 出发，依照：(n1)f(XXX(n)f(X

(n)(n)

)(n 1,2,......))(n)屡次进行计算，当X满足合适的收敛条件就是上面方程的根。这样的方法就是所谓的牛顿—拉夫逊法。(n)这一方法还可以够做下边的讲解，设第n次迭代获取的解语真值之差，即X的误差为时，那么：f(X(n))0把f(X(n))在X(n)周边对用泰勒级数张开2f(X(n))f(X(n))f(X(n))f(X(n))......02!上式省略去 2今后局部f(X(n)) f(X(n)) 0(n)的误差能够近似由上式计算出来。(Xf(X

(n)(n)

))(n)比较两式，能够看出牛顿—拉夫逊法的休整量和 X 的误差的一次项相等。用同样的方法考虑，给出 n个变量的n个方程：f1(X1,X2, ,Xn) 0f2(X1,X2, ,Xn) 0fn(X1,X2, ,Xn) 0对其近似解X1得修正量 X1能够经过解下边的方程来确定：4f1(X1,X2,,Xn)f1f1f1x1x2xnX1f2(X1,X2,,Xn)f2f2f2X2x1x2xnfn(X1,X2,,Xn)fnfnfnXnx1x2xn式中等号右边的矩阵fn都是对于X1,X2,,Xn的值。这一矩阵称为雅可比xn〔JACOBI〕矩阵。按上述获取的修正向量X1,X2,,Xn后，获取以下关系XnXnXn这比X1,X2, ,Xn更凑近真实值。这一步在收敛到希望的值以前重复进行，一般要屡次计算满足maxX1n1X1n1,X2n1X2n1,,Xnn1Xnn1为起初规定的小正数，Xnn1是第n次迭代Xn的近似值。2、牛顿法的框图及求解过程1〕用牛顿法计算潮流时，有以下的步骤：①给这各节点电压初始值e(0),f(0)；②将以上电压初始值代入公式，求修正方程的常数项向量P(0), Q(0),(V2)(0)；③将电压初始值在带入上述公式，求出修正方程中系数矩阵的各元素。④解修正方程式e(0),f(0)；⑤修正各节点电压e(1)e(0)e(0)，f(1)f(0)f(0)；⑥将e(1)，f(1)在带入方程式，求出P(1),Q(1),(V2)(1)；(k) (k)⑦检验可否收敛，即 max Pi , Qi若是收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。若是不收敛，转回②进行下次迭代计算，直到收敛为止。〔2〕程序框图以下：5启动输入原始数据形成节点导纳矩阵分解各节点初始电压的实部和虚部迭代次数 K=0(k)(k)(k)(k)求PQ节点的Pi,Qi，求PV节点的Pi，Ui置节点号 i=0是雅克比矩阵可否形成，i>n否求得雅克比矩阵各元素(k)(k)Jij(k)(k)(k)Sij(k)HijNijLijRij增大节点号i=i+1把雅克比矩阵单位化求解修正方程，得 ei(k)，fi(k)回带各电压新值， K=K+1求解最大修正量| ei(k)|max，| fi(k)|max否可否收敛是计算输出电压大小及相角，节电功率及支路耗费停止三、问题求解3.1等值电路的计算电压是衡量电力系统电能质量的标准之一。电压过高或过低，都将对人身及其用电设备产生重要的影响。保证用户的电压凑近额定值是电力系统调换的根本6任务之一。当系统的电压偏离赞同值时，电力系统必定应用电压调治技术调治系统电压的大小，使其保持在赞同值范围内。本文经过手算形成了等值电路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调整范围。我们第一对给定的程序输入局部作了简要的解析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。〔1〕为了保证整个系统潮流计算的完满性，我们把凡拥有母线及发电机处均选作节点，这样，能够确定有10个节点，节点号见等值电路图。〔2〕确定完节点及编号后，各条支路也相应确定了，网络中总计有13条支路，我们对各支路参数进行了计算。依照所给实质电路图和题中的条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。选择电压基准值为UB220KV和功率基准值SB100MVA因此ZBUB484SB3.2依照题目要求及原始资料画出系统等值电路图 ：35791:k1:k1:k1:k12342468jjj110U=231V发电厂一δ=0P=200MW发电厂二U=231V图1系统等值电路利用psat软件可得计算所需系统图如图 2所示。7图23.3潮流计算3.3.1B1、B2矩阵的形成：依照所求参数，以及B1矩阵的含义，列写B1矩阵以下:B1=[1210;1410;1610;23011;45011;4610;67011;6810;89011;81010];由各个变电所负荷功率能够计算出总功率为200MW，而发电厂一、二的总装机容量分别为400MW和200MW。令发电厂二的功率为200MW，为了减小线路上的耗费，令发电机的电压为额定电压的1.05倍。依照前面表达的节点分类，形成B2矩阵以下：B2=[0023123101;00220002;0220002;800220002;0220002;00220002;0220002;00220002;0220002;200023123103]; 编写程序并运行从结果中能够观察到各个节点电压标幺值分别为 :节点12345678910电压折算到有名值今后，能够观察到节点3、5、7、9的电压都与题目给定的范围对照偏高，因此调治变压器分接头和发电厂电压，经过对系统进行屡次调整，最后获取合理结果，调整结果以下:调治方法：调治方法电厂1电压电厂2电压分接头1分接头2分接头3分接头4未调整1111调整后调治结果：节点号12345678910未调整调整后由上表观察到，进行最后一次调治节气点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费未调整调整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路总耗费由此能够看出，随着变压器分接头的调高有功耗费逐渐增大，吻合实质情况，而最后调治使得电压在规定范围内时，有功耗费增加了，相对与实质耗费小很多，能够认为是合理的。详尽的潮流分布见下表9各条支路的首端功率 Si 各条支路的尾端功率 Sj 各节点的功率S为(节点号从小到大排列)： 获取的各节点的电压标幺值， 电压角度，节点注入无功，支路首端无功如图所示：10值度5幺1角标压压0电电510-551000节点号节点号1功功无无入端注0首点路节支34567891002345678910121图33.3.3 4个变电所负荷同时以 2%的比率增大由于变电所负荷以 2%的比率增大，因此各个变电所负荷变为：10变电所1SL变电所2SL变电所3SL变电所4SL相应地，B2矩阵变为：B2=[0023123101;00220002;0220002;00220002;0220002;00220002;0220002;00220002;0220002;200023123103];将程序中的B2矩阵代替后，运行并进行调治，此系统经过调治后到达要求，对应的调治方法和结果见下表：调治方法：调治方法电厂1电压电厂2电压分接头1分接头2分接头3分接头4调整前1111调整后调治结果：节点电压12345678910调整前调整后由上表观察到，进行最后一次调治节气点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费未调整调整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路总耗费11由电压和功率耗费能够观察到，最后一次的电压和有功耗费均吻合题目要求，能够为是合理的。详尽的潮流分布以下：各条支路的首端功率 Si 各条支路的尾端功率 Sj 各节点的功率S为(节点号从小到大排列)：获取的图像以下：10值度5幺1角标压压0电电510-551000节点号节点号1功功无无入端注0首点路节支34567891002345678910121图43.3.4 4个变电所负荷同时以 2%的比率下降12由于变电所负荷以2%的比率增大，因此各个变电所负荷变为：变电所1SL变电所2SL变电所3SL变电所4SL相应地，B2矩阵变为：B2=[0023123101;00220002;0220002;00220002;0220002;00220002;0220002;00220002;0220002;200023123103];将程序中的B2矩阵代替后，运行并进行调治，对应的调治方法和结果见下表：调治方法：调治方法电厂1电压电厂2电压分接头1分接头2分接头3分接头4调整前1111调整后调治结果见下表：节点号12345678910调整前调整后由上表观察到，进行最后一次调治节气点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费未调整调整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路13(8,10)支路总耗费详尽的潮流分布以下表各条支路的首端功率Si各条支路的尾端功率Sj 各节点的功率S为(节点号从小到大排列)： 获取的各节点的电压标幺值，电压角度，节点注入无功，支路首端无功情况见图7。10值度5幺1角标压压0电电510-551000节点号节点号1功功无无入端注0首点路节支34567891002345678910121图5141、4号变电所负荷同时以2%的比率下降，2、3号变电所负荷同时以2%的比率增加当1、4号变电所负荷同时以2%的比率下降，2、3号变电所负荷同时以2%的比率增加，各个变电所负荷变为：变电所1SL变电所2SL变电所3SL变电所4SL相应地，B2矩阵变为：B2=[0023123101;00220002;0220002;00220002;0220002;00220002;0220002;00220002;0220002;200023123103];将程序中的B2矩阵代替后，运行并进行调治，对应的调治方法和结果见下表：调治方法：调治方法电厂1电压电厂2电压分接头1分接头2分接头3分接头4调整前1111调整后调治结果见下表：节点号12345678910调整前调整后由上表观察到，进行最后一次调治后节点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费未调整调整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路15(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路总耗费由电压和功率耗费能够观察到，最后一次的电压和有功耗费均吻合题目要求，能够为是合理的。详尽的潮流分布以下表各条支路的首端功率 Si 各条支路的尾端功率 Sj 各节点的功率S为(节点号从小到大排列)： 获取的各节点的电压标幺值，电压角度，节点注入无功，支路首端无功情况见图6。1610值度5幺1角标压压0电电510-551000节点号节点号1功功无无入端注0首点路节支567891001234图6 断开1、4支路时的潮流解析当断开1、4支路时，以致支路数变为 9和B1矩阵改变，，其他参数不变。用MATLAB软件进行潮流计算进行调整。B1矩阵变为B1=[1210;1610;2301;4501;4610;6701;6810;8901;81010];节点号12345678910调整前调整后由上表观察到，进行最后一次调治后节点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费 未调整 调整后(1,2)支路 (1,6)支路 (2,3)支路 (4,5)支路 17(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路总耗费由电压和功率耗费能够观察到，最后一次的电压和有功耗费均吻合题目要求，能够为是合理的。详尽的潮流分布以下表各条支路的首端功率 Si 各条支路的尾端功率 Sj各节点的功率S为(节点号从小到大排列)：获取的调治后各节点的电压标幺值，电压角度，节点注入无功，支路首端无功情况见图7。10值5度幺1角标压压0电电510-551000节点号节点号11功功无无入端注0首0点路节支34567891023456789121图718 断开1、6支路时的潮流解析当断开1、6支路时，以致支路数变为 9和B1矩阵改变，，其他参数不变。用MATLAB软件进行潮流计算进行调整。B1矩阵变为B1=[1210;1410;2301;4501;4610;6701;6810;8901;81010];节点号12345678910调整前调整后由上表观察到，进行最后一次调治后节点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费未调整调整后(1,2)支路(1,4)支路(2,3)支路(4,5)支路(4,6)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路总耗费由电压和功率耗费能够观察到，最后一次的电压和有功耗费均吻合题目要求，能够为是合理的。详尽的潮流分布以下表各条支路的首端功率 Si 各条支路的尾端功率 Sj 19 各节点的功率S为(节点号从小到大排列)：获取的调治后各节点的电压标幺值，电压角度，节点注入有功，支路首端无功情况见图8。15值10度幺角5标压压1电电0510-551000节点号节点号1功功无无入端注0首点路0节支34567891023456789121图83.3.8断开4、6支路时的潮流解析当断开4、6支路时，以致支路数变为9和B1矩阵改变，，其他参数不变。用MATLAB软件进行潮流计算进行调整。B1矩阵变为B1=[1210;1410;1610;2301;4501;6701;6810;208901;81010];表17节点号12345678910调整前调整后由上表观察到，进行最后一次调治后节点3，5，7，9的电压均在题目赞同的范围内，对线路耗费进行解析，统计每次调整后各个支路的有功耗费，记录于下表：功率耗费未调整调整后(1,2)支路(1,4)支路(1,6)支路(2,3)支路(4,5)支路(6,7)支路(6,8)支路(8,9)支路(8,10)支路总耗费由电压和功率耗费能够观察到，最后一次的电压和有功耗费均吻合题目要求，能够为是合理的。详尽的潮流分布以下表各条支路的首端功率Si各条支路的尾端功率Sj各节点的功率S为(节点号从小到大排列)：获取的调治后各节点的电压标幺值，电压角度，节点注入无功，支路首端无功情况见图9。2115值10度幺角5标压压1电电0510-551000节点号节点号1功功无无入端注首点0路节支345678910023456789121图9四、误差解析以系统在正常运行的情况下为例。利用PSAT软件，比对MATLAB中的B1、B2矩阵，获取各个节点的电压情况。将PSAT软件运行的结果对MATLAB运行的结果进行校验。同时利用该软件对图形颜色的变化的描绘和数据对结果进行呈现。PSAT软件运行的结果如图 10所示。获取各节点的电压标幺值，并与 MATLAB程序运行获取的结果进行比较。见下表节点12345678910PSATMatlab误差0022VoltageMagnitudeProfile1.u.p[V012345678910Bus#图10经过上述的比较可知:应用DDRTS仿真和MATLAB计算的到的结果误差几乎为0-0.02之间，相差很小，大体同样，能够认为此潮流计算是合理、正确的。23五、心得领悟及总结此次课程设计使我在潮流计算、MATLAB的使用方面均有所提高，但也裸露出了一些问题；理论知识的积蓄还是缺乏，对MATLAB的性能和特点还不能够有一个全面的掌握，对MATLAB中PSAT工具箱使用不够熟悉，相信经过今后的学习能填充这些缺乏，从而到达一个新的层次。潮流计算是电力系统的最根本、最常用的解析计算。用以研究系统规划和运行中提出的各种问题。对规划中的电力系统，经过潮流计算能够检验所提出的电力系统规划方案可否满足各种运行方式的要求；对运行中的电力系统，经过潮流计算能够预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压可否在赞同的范围以内，系统中各种元件(线路、变压器等 )可否会出现过负荷，以及可能出现过负荷时应起初采用哪些预防措施等。实质的情况远比我们计算的情况复杂，这让我深刻认识了潮流计算的重要性。精准的潮流计算不但能够使电网处于牢固且平衡的状态，这对电网的安全运行起到要点性的作用，还可以够实现运行的经济性。因此我认为学好电力系统的要点在于学好潮流计算，电力系统的安全运行不是儿戏，一个很小的失误即可能造成整个电网的崩溃，这给国家和社会带来的危害和损失之大是可想而知的。学好电力系统这门课、学好潮流计算是我们电专业学生的当务之急。没有牢固的基础就不能能建筑出高楼大厦，因此，要想在今后的工作岗位上有所成就，就必定打好基础，努力提高自己的专业涵养。参照文件1、?电力系统解析基础?，韦刚，中国电力初版社 ;2、?电力系统稳态解析?，陈珩，中国电力初版社 ;3、?MATLAB合用指南?，苏金明，王永利，电子工业初版社。24附录一、正常情况下的参数二、MATLAB程序%本程序的功能是用牛顿——拉夫逊法进行潮流计算B1矩阵：1、支路首端号；2、尾端号；3、支路阻抗；4、支路对地电纳5、支路的变比；6、支路首端处于K侧为1，1侧为0B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量6、节点分类标号clear;n=10;%input(' 请输入节点数：n=');nl=10;%input(' 请输入支路数：nl=');isb=1;%input(' 请输入平衡母线节点号： isb=');pr=0.00001;%input(' 请输入误差精度：pr=');B1=[128.5+20.1i0.000556i10;1413.6+32.16i0.0002224i10;1613.6+32.16i0.0002224i10;231.495+40.335i01.051;451.78+53.885i01.0251;4610.2+24.12i0.0001668i10;25671.495+40.335i01.0251;686.8+16.08i0.0004448i10;891.78+53.885i01.0251;8108.5+20.1i0.000556i10];%input('请输入由支路参数形成的矩阵： B1=');B2=[0022922901;00220002;050+30.987i220002;00220002;040+27.79i220002;00220002;050+30.987i220002;00220002;060+37.18i220002;200022922903];%input('请输入各节点参数形成的矩阵：B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);%－－－－－－－更正局部－－－－－－－－－－－－ym=1;SB=100;UB=220;%ym=input('您输入的参数是标么值？〔假设不是那么输入一个不为零的数值〕');ifym~=0%SB=input(' 请输入功率基准值:SB=');%UB=input(' 请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;fori=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1 矩阵B1=');disp(B1)fori=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2 矩阵B2=');disp(B2)end%%%---------------------------------------------------26fori=1:nl %支路数ifB1(i,6)==0 %左节点处于低压侧p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵disp(' 导纳矩阵Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部fori=1:n %给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV节点电压给定模值endfori=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%===================================================================P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;whileIT2~=0IT2=0;a=a+1;fori=1:nifi~=isb %非平衡节点C(i)=0;D(i)=0;forj1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);% Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);% Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);% 节 点 功 率 P 计 算 ei Σ(Gij*ej-Bij*fj)+fi Σ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);% 节 点 功 率 Q 计 算 fi Σ(Gij*ej-Bij*fj)-ei Σ(Gij*fj+Bij*ej)%求P',Q'27V2=e(i)^2+f(i)^2; %电压模平方%=========以下针对非PV节点来求取功率差及Jacobi矩阵元素=========ifB2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; % 节点无功功率差%===============以上为除平衡节点外其他节点的功率计算=================%=================求取Jacobi矩阵===================forj1=1:nifj1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); %dP/df=dQ/deX3=X2; %X2=dp/dfX3=dQ/deX4=-X1; %X1=dP/deX4=dQ/dfp=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~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;% 扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;% 扩展列△PJ(m,q)=X2;endendelse%===============下面是针对 PV节点来求取 Jacobi 矩阵的元素===========DP=P(i)-P1; %PV节点有功误差DV=V(i)^2-V2; %PV节点电压误差forj1=1:nifj1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); %dP/dfX5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;28J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseifj1==i&j1~=isb%非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/dfX5=-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;endendendendend%=========以上为求雅可比矩阵的各个元素 =====================fork=3:N0 %N0=2*n 〔从第三行开始，第一、二行是平衡节点〕k1=k+1;N1=N; % N=N0+1即N=2*n+1扩展列△P、△Qfork2=k1:N1 %扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); %非对角元规格化endJ(k,k)=1; % 对角元规格化ifk~=3 % 不是第三行%============================================================k4=k-1;fork3=3:k4 %用k3行从第三行开始到当前行前的 k4行消去fork2=k1:N1 %k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);% 消去运算endJ(k3,k)=0;endifk==N0break;end%==========================================fork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);% 消去运算endJ(k3,k)=0;29endelsefork3=k1:N0fork2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);% 消去运算endJ(k3,k)=0;endendend%==== 上面是用线性变换方式将 Jacobi矩阵化成单位矩阵=====fork=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); % 更正节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); % 更正节点电压虚部end%------ 更正节点电压-----------fork=3:N0DET=abs(J(k,N));ifDET>=pr % 电压误差量可否满足要求IT2=IT2+1;% 不满足要求的节点数加 1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp(' 迭代次数：');disp(ICT1);disp(' 没有到达精度要求的个数： ');disp(ICT2);fork=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%===============计算各输出量===========================disp(' 各节点的实质电压标幺值 E为(节点号从小到大排列)：');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp(' 各节点的电压大小 V为(节点号从小到大排列)：');30disp(V);VV=V*UB;disp(VV);disp('-----------------------------------------------------');disp(' 各节点的电压相角 sida为(节点号从小到大排列)：');disp(sida);forp=1:nC(p)=0;forq=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp(