含分布式电源的配电网潮流计算和分析的范文

...wd......wd......wd...毕业设论题目含分布式电源的配电网潮流计算与分析含分布式电源的配电网潮流计算与分析摘要电力系统迅速开展，分布式电源的应用也越来越多。分布式电源在并入电网以后，会对配电网的电压、网络损耗等产生一定的影响，所以对分布式电源并网后的潮流计算研究显得尤为重要。本文中，首先对分布式电源的定义和国内外的形势做了介绍。对于用到的数学算法进展了介绍，然后总结了常用的潮流算法。针对PQ算法用于配电网的缺乏提出了改良的方法〔BX法〕，并且在此根基上参加了二范数。对于5、9、14以及33节点的网络进展了算例分析，在5、9节点的网络中改善的算法能够使PQ算法很好的收敛，但是在14、33节点的网络中效果不是很明显，该算法还有待提高。总的来说BX是能够改良PQ算法的，但是还存在缺乏，需要对这个算法进展改善，一定会改善PQ算法的，使该算法能够完全的适用于所有的配电网。关键词：分布式电源、PQ算法、配电网、BX法AbstractWithrapiddevelopmentofthepowersystem,theapplicationofdistributedgeneratorsupplybecomeordinary.Powergridwithdistributedgeneratorwillinfluencethevoltage,networklossofthepowerdistributionnetwork.Sotheresearchofpowergridwithdistributedgenerator‘spowerflowcalculationisparticularlyimportant.Inthisarticle,thedefinitionofdistributedgeneratorsupplyandthesituationofbothathomeandabroadisproducedfirstly.Usedformathematicalalgorithmareintroduced,andthensummarizesthecommontrendofthealgorithmAimingattheshortcomingsofthePQalgorithmusedindistributionnetworkisputforwardtoimprovethemethod,themethod(BX).For5,9,14and33nodesofnetworkhascarriedontheexampleanalysis,in5and9nodesinthenetworktoimprovethealgorithmcanmakethePQgoodconvergence,butat14,33nodesinthenetworkeffectisnotobvious,sothealgorithmneedstobeimproved.TheBXmethodisabletoperfectthePQalgorithm,butthereareinsufficient,itisnecessarytoimprovethealgorithm,IthinkBXwillimprovePQalgorithm,thealgorithmcanbecompletelyappliedtoallofdistributionnetwork.Keywords:distributedgenerator,thePQalgorithm,theBXmethod,powergrid目录TOC\o"1-3"\h\u17482第一章绪论 ⑷对变化后的矩阵用矩阵中的元素进展回代运算，(3-9)(3-10)根据以上〔3-9〕、〔3-10〕那么可求出未知数,即。然后进展有功功率不平衡量的判断，看其是否满足精度要求，满足那么进展其他配电网中的其他的参数的运算，不满足那么形成新的电压相角，并利用新的电压相角参数来进展下一次的迭代计算。在无功功率不平衡量的计算中与其原理是一样的。以上的所有过程就是因子表法的在本次毕业设计重的具体应用，算法中的主程序也是这样进展编写的。3.2二范数的介绍及应3.2.1二范数的介绍范数为在向量空间内对所有的向量赋一个非零的大小。其在数值代数中研究数值的误差和收敛性方面显得很重要。由于本次的设计正是有关算法的收敛性问题，因此才参加了二范数的概念，在本文的算法验证中也表达了二范数对于结果的收敛性有所改善。在向量范数中，描述了向量与实数之间的关系。在不同的情况下可以根据定义设置不同的范数。在常用的向量范数中常用的有1-范数、2-范数和-范数。以下是三者的具体运算公式：1-范数：设有向量X=〔〕,；2-范数：设有向量X=〔〕,；-范数：设有向量X=〔〕，。在这些范数中通过类似的变化，就可以得其通用的式子，即p-范数，公式为在这些范数中存在以下的大小关系：，那也就是说明，所以在本次设计中运用了二范数来代替了无穷范数，而且在运行结果中也可以看出加快了收敛特性，减少了运算时间。3.2.2二范数的应用本次算法中是讨论利用PQ算法进展配电网的潮流计算的时候，所得到的结果是否收敛，判断的依据是,原有的计算中是在一组解中找出最大的然后进展判断，这样做固然可以使所有的满足精度要求，但是误差太大。因此利用二范数来处理所得到的，利用二范数小于无穷范数来减小迭代次数，加快潮流计算的时间并且取到了较好的效果。3.3BX法的介绍及应用在本次的毕业设计中是利用PQ算法来对配电网进展计算，但是一般说来PQ算法并不适应于配电网的潮流计算，在前面的1.3中已经说明不适用的原因，故在此不再进展说明。下来的具体问题就是若何解决PQ算法比照值的敏感性的问题。在此设计中采用的BX法其实是指在PQ算法运算的最初阶段形成系数矩阵和的时候对它的矩阵参数进展改良【8】。原本的两个系数矩阵的参数是电网中的导纳矩阵的虚部。而在本次的设计中分别对这个系数矩阵进展了不同的处理。例如，因此形成的的对角元素就按如下公式计算：(3-11)(3-12)在形成系数矩阵没有考虑对地并联的导纳，忽略了支路电阻，也不考虑变压器变比的影响，因此形成的的对角元素就按如下的式子进展计算：(3-13)(3-14)其中，为线路的电抗值，是为节点的接地导纳的数值。通过以上的对系数矩阵在形成的过程中的简化方法，来改变PQ算法比值的敏感性问题。经过验证，在此种方法下，对有些小的配电网的系统的收敛性有所改善，使其能够收敛。第四章潮流计算结果的分析4.1基于P-Q算法的潮流算法的程序流程图开场开场读取数据读取数据形成系数矩阵形成系数矩阵数据标幺化数据标幺化形成因子表形成因子表计算有有功不平衡量计算有有功不平衡量Y有功满足精度要求Y有功满足精度要求NN修正电压相角修正电压相角计算无功不平衡量计算无功不平衡量N无功满足要求N无功满足要求YY计算电压、网损等计算电压、网损等完毕完毕图4-1PQ算法的总流程图4.2算例结果分析4.2.1五节点的算例结果分析在5节点配电网的计算中，节点1位平衡节点，其余的节点都是PQ节点。其中分布式电源加在节点2上面。五节点网络拓扑构造见附录。对于五节点的网络分别进展了配电网的一般的PQ潮流计算、用BX法进展改良的潮流计算以及加二范数的潮流计算。使用一般的PQ法进展潮流计算，结果不收敛。下面是用BX法进展潮流计算的结果，其中基准电压，，精度要求为0.00001.表4-1五节点电压节点编号电压值电压相角11020.966858-1.5369630.921383-1.5407240.89052-1.53714650.874982-1.536913表4-2五节点网损线路首端线路末端线路功率120.002133+0.002251*i230.001620+0.001615*i340.001044+0.000509*i450.000509+0.000509*i21-0.002120-0.002120*i32-0.001544-0.001539*i43-0.001009-0.001013*i54-0.0005-0.0005*i线路总损耗0.000133+0.000251*i在算法中参加二范数之后，其结果如下表所示：表4-35节点电压节点编号电压值电压相角11020.966858-1.53695930.921455-1.54451540.890528-1.53757950.874982-1.536938表4-4五节点网损线路首端线路末端线路功率120.002133+0.002251*i230.001619+0.001611*i340.001045+0.001052*i450.000509+0.000510*i21-0.002120-0.002120*i32-0.001543-0.001535*i43-0.001010-0.001017*i54-0.000500-0.000501*i线路总损耗0.000133+0.000251*i下表是利用不同的算法时的结果收敛情况：表4-5不同方法的运行结果算法类型是否收敛迭代次数时间(s)一般PQ法否含BX法是382.734BX法+二范数是341.906从以上的表格中可以看出，当在进展配电网的潮流计算在使用PQ时候，计算的结果不收敛，这就是之前所说的PQ法对电阻与电抗比值的敏感性问题。而从表4-3中可以看出当使用BX法对PQ算法进展改良的时候，算法的结果收敛，然后在此根基上参加二范数的时候，迭代的次数有所减少，并且计算的速度也加快了。说明在五节点的配电网中，BX法对于算法的改良是很适用的。而且在二范数的参加PQ算法的时候，对比两组数据明显可以看出所得出的结果几近相似，而且所得的总网损的结果是完全一样，因此二范数的参加并没有结果的精度。因此可以得出在形成系数矩阵的时候的时候，忽略支路的电阻对于PQ算法有较好的改善。4.2.2九节点算例结果分析在9节点的配电网络系统中，节点1为平衡节点，其余的节点都是PQ节点，其中分布式电源加在节点2上面。九节点网络拓扑构造见附录。下面是对其电网的线路损耗和电压以及运算结果的情况。下表为用BX法进展的P-Q算法的结果：表4-6九节点电压节点编号电压值电压相角11020.982881-0.79798430.947173-0.81181940.916606-0.79902350.8907-0.79799760.869776-0.79794270.853971-0.79794880.84338-0.79795490.838067-0.797956表4-7九节点线路功率线路首端线路末端线路功率120.034270+0.035238*i230.039667+0.039162*i340.032297+0.032773*i450.026551+0.026684*i560.020928+0.015438*i670.015438+0.015438*i780.010158+0.010158*i890.005032+0.005032*i21-0.034162-0.034162*i32-0.038235-0.037730*i43-0.031248-0.031683*i54-0.025897-0.025930*i65-0.020437-0.020438*i76-0.015158-0.015158*i87-0.010032-0.010032*i98-0.005000-0.005000*i表4-8九节点网络的损耗线路编号支路损耗10.000108+0.001076*i20.001432+0.014132*i30.026551+0.026684*i40.001049+0.001049*i50.000754+0.000754*i60.000281+0.000281*i70.000126+0.000126*i80.000032+0.000032*i总损耗0.004272+0.005241*i以下的表格是在BX法的根基上参加了二范数的运行结果：表4-9九节点电压节点编号电压值电压相角11.0000000.00000020.982881-0.79798530.948027-0.76014040.916713-0.79227050.890708-0.79741660.869774-0.79806070.853968-0.79811880.843377-0.79810990.838064-0.798101表4-10九节点线路功率线路首端线路末端线路功率120.034270+0.035238*i230.037781+0.039163*i340.033886+0.032792*i450.026855+0.026690*i560.020950+0.020930*i670.015440+0.015438*i780.010157+0.010158*i890.005031+0.005032*i21-0.034162-0.034163*i32-0.036416-0.037799*i43-0.032785-0.031690*i54-0.026095-0.025930*i65-0.020458-0.020438*i76-0.015159-0.015158*i87-0.010031-0.010032*i98-0.005000-0.005000*i表4-11九节点网络的损耗线路编号支路损耗10.000108+0.001076*i20.001365+0.001365*i30.001102+0.001102*i40.000759+0.000759*i50.000492+0.000492*i60.000281+0.000281*i70.000126+0.000126*i80.000032+0.000032*i总损耗0.004264+0.005232*i下表为在9节点的配电网中分别用一般的PQ法、用BX进展改良的PQ法以及在该方法的根基上参加了二范数所得到的运行结果：表4-12九节点网络的运行结果算法类型是否收敛迭代次数时间(s)一般PQ法否含BX法是442.469BX法+二范数是351.905从表4-12可以看出一般的PQ法在九节点的配电网的潮流计算中结果不收敛，即不能直接用PQ潮流计算。发现在用BX法进展改良的算法中结果收敛，在此根基上参加二范数时计算的速度和迭代次数都有所改善。而且比照在这两种方法的使用情况下所得到的电压值、网络的节点功率、线路损耗以及网络的总损耗等数据情况，可以看出二范数虽然加快了计算速度，也减少了迭代次数，但是所得的数据都基本相等，所造成的误差在允许的范围之内，说明BX法加上二范数可以对九节点的网络潮流计算有所改善。4.2.3十四节点算例结果分析在十四节点的网络中，其系统的基准电压，系统的基准容量为,精度要求为0.00001。表4-13十四节点网络的运行结果算法类型是否收敛迭代次数时间(s)一般PQ法否含BX法否BX法+二范数否从表4-13可以看出，对于十四节点的配电网络，在参加了BX法和二范数的情况下结果均不收敛，说明改良之后的PQ算法不适合14节点的配电网络。由于时间的原因，在程序中也就没能够参加新的算法，所以十四节点结果的网络未能够得到改善。本次设计也用IEEE三十三节点标准配电网络进展测试，结果也不是很理想。将本章中所有表格之间的比照，可看出如改良PQ算法，配网中所得潮流计算结果不收敛。但由于其计算的速度快以及需要容量少的优点，在小网络中对PQ法进展改良取得的效果还是挺好的，对于大网络，效果不明显。第五章总结与展望5.1总结随着智能电网的开展，分布式电源必然会在配电网的应用中越来越广泛，它的应用有利于改善环境、节约能源，符合可持续开展的国策。最重要的就是其在电网发生故障的时候可以为那些重要的负荷供电，减少电力事故带来的损失。本次设计的主要目的就上改善P-Q算法，使该算法能够应用于配电网中，由于该算法对于电阻与电抗比值的敏感性，所以它并不适用于配电网的潮流计算中，但由于该算法的快速性以及减少容量要求，所以就意图改善它的缺乏，使其能够应用于配电网的潮流计算中。本次设计用了BX法对P-Q算法进展了改良，又利用二范数来处理数据收敛条件的，加快了计算速度，减少了迭代次数。但是。在十四节点的网络中这种方法并没有使结果收敛，因此本次设计还有缺乏，有待提高。5.2展望在本次的毕业设计中，改良的P-Q算法虽然在五节点和九节点的配电网络中收敛的效果对比好，但是在十四节点以及更大的配电网络中仍然会出现不收敛的问题。因此，需要在算法方面做更多的研究，使P-Q算法能够完全的适用于配电网的计算。可以从以下几个方面入手：=1\*GB3①对于BX法可以在形成系数矩阵的时候，对考虑电阻和忽略电阻进展不同的组合，找出收敛情况最好的组合：=2\*GB3②可以使用支路分解法来解决PQ算法的收敛性问题【15】；=3\*GB3③可以更多的参考国内外的参考文献，研究出新的算法来改善-Q算法的缺乏。致谢参考文献【1】李琼慧，黄碧斌，蒋莉萍.国内外分布式电源定义及开展现况比照分析[J].中国能源,34(8):2012:31-34【2】杨建华.PQ分解法潮流计算的收敛性的影响因素[C].1998年全国高等学校电力系统及其自动化专业第14届学术年会论文，1999.17-20【3】陈海焱,陈金富,段献忠.含分布式电源的配电网潮流计算[J].电力系统自动化,30〔1〕：2006：35-40【4】马超,李春兰，岳勇，石砦.基于PQ分解法的含分布式电源的配电网潮流计算研究[N].沈阳农业大学学报，44〔3〕：2013：327-332【5】王锡凡.现代电力系统分析［M］.科学出版社，2004：17-21,55-88【6】周全仁，张清益.电网计算与程序设计[M].湖南科学技术出版社，2007：1-15【7】杨建华.潮流计算PQ分解法的最正确方式探讨[J].电网技术，20〔1〕：1996：30-32【8】李志强.含分布式电源的配电网潮流计算[D].郑州：郑州大学，2012：24-30【9】刘庆.含分布式电源配电网潮流计算及应用［D］.西安：西安科技大学，2011：3-18【10】YangWenyu，YangXuying，DuanJiandongWanxiaozhworksContainingDistributedGeneration[J].ElectricityDong,FanYue.PowerFlowCalculationinDistributionNetistribution,2008.CICED2008.ChinaInternationalConferenceon.IEEEConferences,2008:1【11】JasonKennedy,BrendanFox,DJohnMorrow.DistributedDieselGenerationtoFacilitateWindPowerIntegration[J].PowerTech,2007461-466【12】刘建，毕鹏翔，杨文宇，程红丽.配电网及应用[M].中国水利水电出版社，2007：112-122【13】张雪松，刘焯.配电网潮流算法对比研究[J].电网技术.22〔4〕：1998：1-5【14】李新，彭怡，赵晶晶，任亚英，张从力.分布式电源的并网及潮流计算[J].电力系统保护与控制.37(17):2009:1-5【15】李俊，肖宏，唐荣平，龙江.一种潮流计算的PQ改良算法研究[J].群众科技，16〔182〕：2014：103-105附录一.五节点配电网络:附图25节点配电网络附表25节点配电网支路参数支路编号起始节点到达节点R(p,μ)X(p,μ)1120.122321.22322231.22321.22323341.22321.22324451.22321.22325————0.122321.2232附表35节点功率节点编号P(MW)Q(MVAR)C(MVAR)10002-50-500350500450500550500二.九节点配电网络:附图49节点配电网络附表59节点配电网络参数支路首节点末节点电阻〔Ω〕电抗〔Ω〕有功功率〔KW〕无功功率〔Kvar〕10、10.16031.60280021、21.60281.6028-50-5032、31.60281.6028505043、41.60281.6028505054、51.60281.6028505065、61.60281.6028505076、71.60281.6028505087、81.60281.60285050三.程序clc[n]=xlsread('pqinput.xls','A2:A2');%n为节点数[nl]=xlsread('pqinput.xls','B2:B2');%nl为支路数[isb]=xlsread('pqinput.xls','A4:A4');%isb为平衡节点编号[pr]=xlsread('pqinput.xls','B4:B4');%pr为准确度[fname,pname]=uigetfile('pqinput.xls');%把矩阵数值输到Excel里，然后从Matlab中调出该Excel里的矩阵数值file=[pnamefname];x(:,1)=xlsread('pqinput.xls','D3:D1000');%起始节点编号(fbus)x(:,2)=xlsread('pqinput.xls','E3:E1000');%终止节点编号(tbus)x(:,3)=xlsread('pqinput.xls','F3:F1000');%支路电阻(r)x(:,4)=xlsread('pqinput.xls','G3:G1000');%支路电抗(x)x(:,5)=xlsread('pqinput.xls','H3:H1000');%支路对地导纳(b)x(:,6)=xlsread('pqinput.xls','I3:I1000');%变压器变比(ratio)x(:,7)=xlsread('pqinput.xls','J3:J1000');B1=[x(:,1),x(:,2),(x(:,3)+x(:,4)*i),x(:,5)*i,x(:,6),x(:,7)];%由支路参数形成的矩阵y(:,1)=xlsread('pqinput.xls','N3:N1000');%节点所接发电机的有功功率(Pg)y(:,2)=xlsread('pqinput.xls','O3:O1000');%节点所接发电机的无功功率(Qg)y(:,3)=xlsread('pqinput.xls','P3:P1000');%节电负荷的有功功率(Pd)y(:,4)=xlsread('pqinput.xls','Q3:Q1000');%节电负荷的无功功率(Qd)y(:,5)=xlsread('pqinput.xls','R3:R1000');%节点电压的初始值(Vm)y(:,6)=xlsread('pqinput.xls','S3:S1000');%PV节点电压给定值y(:,7)=xlsread('pqinput.xls','T3:T1000');%节点类型〔1-平衡节点，2-PQ节点，〕B2=[(y(:,1)+y(:,2)*i),(y(:,3)+y(:,4)*i),y(:,5),y(:,6),y(:,7)];%由节点参数形成的矩阵%z(:,1)=xlsread('pqinput.xls','A9:A1000');%节点编号%z(:,2)=xlsread('pqinput.xls','B9:B1000');%对地阻抗值%X=[z(:,1),z(:,2)];%由节点号及其对地阻抗形成的矩阵[na]=xlsread('pqinput.xls','A6:A6');%na为PQ节点数%************************输入参数**************************UB=6;SB=10;ZB=UB^2/SB;B1(:,3)=B1(:,3)/ZB;B2(:,2)=B2(:,2)/SB/1000;tic;t1=clock;%*********************求节点导纳矩阵Y**********************Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);%fori=1:n%ifX(i,2)~=0;%p=X(i,1);%Y(p,p)=1./X(i,2);%end%endfori=1:nlifB1(i,6)==0p=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矩阵互导纳YI(p,q)=YI(p,q)-1./B1(i,3);Y(q,p)=Y(p,q);%求Y矩阵互导纳YI(q,p)=YI(p,q);%求YI矩阵互导纳Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对变压器支路求自导纳YI(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%求导纳矩阵%*********************求节点导纳矩阵Y**********************%G=real(Y);B=imag(YI);BI=-YII/2;%Y的实部虚部，YI的虚部。(YI是不考虑变压器变比和对地导纳的节点导纳矩阵的虚部，理解为节点导纳矩阵的理想电纳。G=real(Y);B=imag(YI);BI=imag(Y);%因为求第一因子表B′和第二因子表B″时，为了加速收敛，使它们更快修正，可以考虑尽量去掉那些与形成B′和B″无关或影响较小的因素，参考?现代电力系统?P83〕fori=1:nS(i)=B2(i,1)-B2(i,2);%发电机功率-节点负荷功率=线路功率endP=real(S);Q=imag(S);fori=1:ne(i)=real(B2(i,3));%节点电压初始值的实部f(i)=imag(B2(i,3));%节点电压初始值的虚部V(i)=B2(i,4);%设V(i)为PV节点电压的给定值幅值endfori=1:nifB2(i,5)~=3%当节点为PQ节点时V(i)=sqrt(e(i)^2+f(i)^2);%该节点的幅值=节点电压初始值的实部和虚部的算术平方根；因为此节点不是PV节点而是PQ节点，所以第75行求Vi的式子求不出电压来，所以需此行求电压O(i)=atan(f(i)./e(i));%相角endend%*****************以上是求PQ节点电压初值，各节点电压相角初值********fori=2:nifi==nB(i,i)=1./B(i,i);%去倒数并赋值elseIC1=i+1;%IC1=3:nforj1=IC1:nB(i,j1)=B(i,j1)./B(i,i);%B〔i,j1〕是从第二行起的上三角矩阵，j1>iendB(i,i)=1./B(i,i);%去倒数并赋值fork=i+1:nforj1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);%消去运算endendendend%**********以上是用YI矩阵形成因子表，并求逆阵〔计算有功功率〕**********p=0;q=0;fori=1:nifB2(i,5)==2;%PQ节点时p=p+1;k=0;forj1=1:nifB2(j1,5)==2;%J1为PQ节点时k=k+1;%PQ节点个数A(p,k)=BI(i,j1);%这里的A为系数矩阵B″endendendend%*********************以上是将Y矩阵中的元素存在A矩阵中**************fori=1:naifi==naA(i,i)=1./A(i,i);elsek=i+1;forj1=k:naA(i,j1)=A(i,j1)./A(i,i);endA(i,i)=1./A(i,i);fork=i+1:naforj1=i+1:naA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendend%*********************以上是用Y矩阵形成因子表，并求逆阵〔计算无功功率〕**********************ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;DET1=0;ICT3=1;%ICT2为不满足有功功率精度值的个数;ICT3为不满足无功功率精度值的个数;%*********************迭代次数K置零，Kp=1，Kq=1开场*****************while(ICT2~=0|ICT3~=0)ICT2=0;ICT3=0;%进展下一次迭代fori=1:nifi~=isb%非平衡节点行C(i)=0;fork=1:n%所有列〔包括PV节点及平衡节点列〕C(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);%节点功率P计算DP(i)=DP1(i)./V(i);%节点功率差除以该点电压%DE(i)=abs(DP(i));endend%DET(K)=max(DE);DET(K)=sqrt((sum(DP1.^2))/n);ifDET(K)>=pr%节点功率误差是否满足要求ICT2=ICT2+1;%本次迭代有功不满足的节点个数end%===============用因子表求解修正方程的修正量角度sida===============Np(K)=ICT2;ifICT2~=0fori=2:nDP(i)=B(i,i)*DP(i);%规格化ifi~=nIC1=i+1;fork=IC1:nDP(k)=DP(k)-B(k,i)*DP(i);%用因子表作前代消去;这里的DP为ΔθiendelseforLZ=3:iL=i+3-LZ;IC4=L-1;forMZ=2:IC4I=IC4+2-MZ;DP(I)=DP(I)-B(I,L)*DP(L);%按列回代;这里的DP为Δθiendendendend%===========================对角度进展修正=================fori=2:nO(i)=O(i)-DP(i);end%=================对无功Q及电压U的关系进展解算===========kq=1;L=0;fori=1:nifB2(i,5)==2%i节点为PQ节点C(i)=0;L=L+1;fork=1:n%所有列〔包括PV节点及平衡节点列〕C(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);%节点功率Q计算,ΔQDQ(L)=DQ1(i)./V(i);%节点功率差除以该点电压,ΔQ/V%DE1(L)=abs(DQ(L));endend%DET1(K)=max(DE1);DET1(K)=sqrt((sum(DQ1.^2))/n);ifDET1(K)>=pr%节点功率误差是否满足要求ICT3=ICT3+1;%本次迭代无功不满足的节点个数endelsekp=0;end%===============用因子表求解修正方程的修正量电压模===============Nq(K)=ICT3;ifICT3~=0L=0;fori=1:naDQ(i)=A(i,i)*DQ(i);%规格化,求ΔViifi==naforLZ=2:iL=i+2-LZ;IC4=L-1;forMZ=1:IC4I=IC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);%回代endendelseIC1=i+1;fork=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i);%前代消去endendend%===================对PQ节点电压模进展修正=======================L=0;fori=1:nifB2(i,5)==2%i节点为PQ节点L=L+1;V(i)=V(i)-DQ(L);%修正电压模值endendkp=1;%有功功率存在不满足的节点K=K+1;elsekq=0;%无功功率都满足要求ifkp~=0%有功功率存在不满足的节点K=K+1;endendend%以下是输出数据，注意涉及k，k=1:n时，把fprintf放在循环中，否那么只会读出k=n时的数据tic;t2=clock;pause(3*rand)disp(['etime程序总运行时间:',num2str(etime(clock,t1))]);fid=fopen('pqoutput.txt','wt');fprintf(fid,'************************