第一篇：2第二章 电力系统潮流计算-2
第二章 电力系统潮流计算
2.1 概 述
2.2 潮流计算问题的数学问题 2.3 潮流计算的牛顿法 2.4 潮流计算的P-Q分解法
2.5 静态安全分析及补偿法
2.5.1 静态安全分析概述
静态安全分析是电力系统规划和调度的常用手段，用以校验输变电设备强迫退出运行后系统的运行状态，回答诸如“假如电网中某一条500kv输电线路开断后，系统运行状态发生什么变化”之类的问题[21,22]。对这个问题的回答可能是系统的潮流和电压都在容许的范围之内，或者出现某些输变电设备过负荷或某些母线电压越界的情况。前者的系统是安全的，后者则是不安全的。因此，静态安全分析是电力系统安全分析的一个重要组成部分，它不涉及电力系统的动态过程的分析，故称为静态安全分析，是以下各节介绍的主要内容。动态安全分析问题的讨论详见第5章及第6章。
利用静态安全分析可以进行事故预想，对一个输电系统规划方案而言，可以校验其承受事故的能力；对运行中的电力系统而言，可以检验其运行方式及接线方式的安全性，进而给出事故前后应采用的防范措施或校正措施。静态安全分析中需要校验的典型事故包括发电机组或输变电设备的强迫停运，也包括短路引起的保护动作致使多个设备同时退出运行的情况。
系统规划设计人员在进行发电系统和输电系统规划时，应利用静态安全分析考虑各种可能的设备开断情况，并评估其后果是否满足安全性的要求。为此，规划设计人员一般需要增加一些冗余的设备或调整计划以减少中断供电的可能性。
在电力系统的运行中，为了避免过负荷和电压越界引起的设备损坏，或由于过负荷设备在系统保护作用下退出运行而导致大面积连锁反应性的停电，在线或实时地进行系统静态安全分析非常重要[23,24]。特别是随着电力市场的进展，电力系统的发输配电各环节由统一管理、统一调度逐步转向双边合同交易和发电厂商的竞价上网，使系统运行出现了诸多不确定因素，对电力系统运行的安全监视和控制提出了更高的要求。
由于不涉及元件动态特性和电力系统的动态过程，静态安全分析实质上是电力系统运行的稳态分析问题，即潮流问题。也就是说，可以根据预想的事故，设想各种可能的设备开断情况，完成相应的潮流计算，即可得出系统是否安全的结论。但是，静态安全分析要求检验的预想事故数量非常大，而在线分析或实时分析又要在短时间内完成这些计算、因此，开发研究了许多专门用于静态安全分析的方法，如补偿法、直流潮流法及灵敏度分析法等，以下将分别介绍这些基本的方法。2.5.2 补偿法
电力系统基本运行方式计算完毕以后，往往还要求系统运行人员或规划设计人员进行一些特殊运行方式的计算，以分析系统中某些支路开断以后系统的运行状态，以下简称断线运行方式。这对于确保电力系统可靠运行，合理安排检修计划都是非常必要的。
发电厂运行状态的变化，如发电厂之间出力的调整和某些发电厂退出运行等情况，在程序中都是比较容易模拟的。因为这时网络结构和网络参数均未发生变化，所以网络的阻抗矩阵、导纳短阵以及P-Q分解法中的因子表都应和基本运行方式一样。因此，我们只要按照新的运行方式给定各发电厂的出力，就可以直接转入迭代程序。应该指出，在这种情况下不必重新送电压初值，利用基本运行方式求得的节点电压作为电压初值可能更有利于收敛。
当系统因故障或检修而开断线路或变压器时，要引起电网参数或局部系统结构发生变化，因此在这种情况下进行潮流计算时，要修改网络的阻抗矩阵或导纳矩阵。
对于牛顿法潮流程序来说，修正导纳矩阵以后，即可转入迭代程序(见图2-5)。
对于P-Q分解法来说，修改导纳矩阵以后，应该先转入形成因子表程序，然后再进行迭代计算(见图2-9)。在程序编制上这样处理比较简单，只需要增加修改导纳矩阵的程序，但是，由于需要重新形成因子表，因此计算速度较慢。为了进一步发挥P-Q分解法的优点，提高计算速度，可以采用补偿法的原理[7]，在原有基本运行方式的因子表的基础上进行开断运行方式的计算。当潮流程序用作在线静态安全监视时，利用补偿法以加速顺序开断方式的检验就显得特别重要。
应该指出，补偿法的概念不仅应用于P-Q分解法潮流程序中，也广泛应用在短路电流、复杂故障以及动态稳定计算程序的网络处理上。以下首先介绍补偿法的基本原理，然后讨论如何利用补偿法进行开断运行方式的计算。
如图2-12所示，设网络N的导纳矩阵已经形成，并对它进行三角分解而得到因子表。现在的问题是，当向网络节点i、j之间追加阻抗zij时．如何根据已知的节点注入电流
利用原电力网络N的因子表，求得新条件下的电压
如果我们能够求得流入原网络N的注入电流向量
&的等效电路 图2-14 求电流Iij那么利用原网络因子表对此I进行消去回代运算就可以得到节点电压向量V。但
&并不知道，因而也就不是在各节点电压求出以前，追加支路Zij上通过的电流Iij能直接利用I求节点电压。
根据迭加原理，可以把图2-12所示网络拆为两个等值网络，如图2-13(a)及(b)所示。节点电压向量V可以表示为
式中：V(0)相当子没有追加支路，或追加支路开路的情况下各节点的电压向量，见图2-13(a)。由于这种情况下各节点的注入电流I已知，因此利用原网络N的因子表不难求得
现在讨论如何求得图2-13(b)中各节点电压V(1)。在这个图中，向原网络注入的电流向量为
图2-13 补偿法原理示意图
&现在暂时还是未知量。但如果假定I&1，则利用原网络因子表就可以求其中Iijij&为单位电流时，网络各节点的电压 得当Iij
&那么由于网络是线性的，这样，如果能求出I就可以按下式求得最终的电压向量： ij
&。为此，需要利用等值发电机原理 因此，现在的关键问题就在于如何求得Iij如上所述，V(0)相当于追加支路Zij开路情况下网络各节点的电压。如果现在把整个系统看成是支路Zij的等值电源，那么这个电源的空载电压就是
电源的等值内阻抗ZT的值应为
(ij)(ij)&&这是因为(V)是在i,j两点分别通入正、负单位电流而在i,j点造成的压Vij降，在数值上应等于从i,j点看进去的输入阻抗。这样，就可以得到图2-14所示的等值电路由该回可以直接求出
&的等效电路 图2-14 求电流Iij
式中：
&代入式(2-90)中．即可得到所需要的节点电压向量V。将式(2-93)求得的I ij以上讨论了补偿法的基本原理。实用上，利用补偿法求解节点电压的过程可按以下步骤进行：
(1)利用原网络的因子表对于单位电流向量
进行消去回代运算．求出Vij。
人(2)利用式(2-92)求等值发电机的内阻抗ZT，并进而根据式(2-94)求Zij(3)利用原网络因子表对节点注入电流向量I进行消去回代运算，求出V0，见图2-13(a)。
&。(4)根据式(2-93)求出流经追加支路Zij的电流Iij(5)利用式(2-90)求出节点电压向量V。
当网络发生一次变化或操作后，需要对不同的节点注入电流I求解节点电压
可以暂时贮存时，步骤(1)及(2)的运其只需要进行一次，把计算结果V、Zijij起来。这样，对不同的I求V时，只需要作步骤(3)~(5)的运算。因此利用补偿法求解网络节点电压和一般用因子表求解网络节点电压相比，在运算量上并没有显著的增加，但是形成一次因子表的运算量约为求解一次网络节点方程运算量的10倍左右，因此，当网络进行一次操作，要求反复求解网络方程的次数小于5次时，用补偿法比重新形成因子表要节约很大的运算量[8]。
补偿法在原理上也可用于网络同时进行两处或多处操作的情况，这时需要递归地套用以上的计算步骤，本书不再详述，有兴趣的读者可以参考文献[3]。
以上介绍了补偿法的原理，下面我们讨论在P-Q分解法潮流程序中如何利用补偿法进行开断运行方式的计算。
对于式(2-81)、式(2-82)所示的修正方程式，可以分别看成是由“导纳矩阵”B及B所描述网络的节点方程式，其注入电流分别为PV及QV，待求的节点电压为V0及V，这样就可以完全套用以上的计算过程。在这种情况下对B及B来说，图2-12中追加支路阻抗应分别为
当开断元件不是线路而是非标准变比的变压器时，式(2-95)的电流表示式应改写为
式中；nT为非标准变比，在j侧。这时式(2-91)~(2-93)相应地变为
式中：
必须注意，以上断线操作在式(2-96)中实际上只考虑了断开电线路和变压器的不接地支路。严格地讲，输电线路对地电容或非标准变比变压器接地支路也应同时断开，但是，这样就成为同时出现3处操作的情况，使计算复杂化。计算实践表明，在利用补偿法进行系统开断运行方式计算时，不计接地支路的影响，给计算带来的误差是很小的，完全可以忽略不计。
2.6 静态安全分析的直流潮流法
直流潮流模型把非线性电力系统潮流问题简化为线性电路问题，从而使分析计算非常方便。直流潮流模型的缺点是精确度差，只能校验过负荷，不能校验电压越界的情况。但直流潮流模型是线性模型，不仅计算快，适合处理断线分析，而且便于形成便于用线性规划求解的优化问题，因此，得到了广泛的应用。2.6.1 直流潮流模型
支路有功潮流可表示为
式中：tij为支路ij的变压器非标准变比；ij为支路ij两端节点电压的相角差；Gij,Bij为节
点导纳矩阵元素的实部与虚部。
式中：rij,xij为支路ij的电阻和电抗，当ij时，将交流潮流很据P-Q分解法的简化条件进行简化，就可以得到如下直流潮流方程
由式(2-104)可知Bij1，但为了以下应用方便起见，我们定义 xij 因此，最后，得到
写成矩阵形式，为
式中：P为节点注入功率向量，其中元素PiPGiPDi，这里PGi和PDi分别为节点i的发电机出力和负荷；为节点电压相角向量；B为节点导纳矩阵的虚部，其元素由式(2-106)和式(2-107)构成。
式(2-109)也可写成另一种形式：
式中：X为B的逆矩阵：
同样，将PQ分解法的简化条件代入支路潮流方程式(2-102)，可以得到
将上式写成矩阵形式，式中；Pl为各支路有功功率潮流构成的向量；为各支路两端相角差向量；Bl为由各支路导纳组成的对角矩阵，设系统的支路数为l，则Bl为l阶方阵。
设网络关联矩阵为A，则有
式(2-109)、式(2-ll0)、式(2-113)均为线性方程，是直流潮流方程的基本形式。当系统运行方式及接线方式给定时，即得到关于的方程(2-109)，通过三角分解或矩阵直接求逆可以由式(2-110)求出状态向量，并进而出式(2-113)求出各支路的有功潮流。
2.6.2 直流潮流的断线模型
由以上讨论可以看出，应用直流潮流模型求解输电系统的状态和支路有功潮流非常简单。现在我们还要指出，由于模型是线性的，故可以快速进行追加和开断线路后的潮流计算。
设原输电系统网络的节点阻抗矩阵为X，支路k两端的节点为i,j。这里的支路是指两节点间各线路的并联，线路是支路中的一个元件。当支路A增加一条电抗为xk的线路(以下称追加线路k)时，形成新的网络。根据1.4节的支路追加原理，新网络的节点阻抗矩阵X应为[见式(1-107)]
式中：
式（2-117）可以简写为
式中：
式中：
式中：Xii,Xjj,Xij均为X中的元素。由式（2-118）可知节点阻抗矩阵的修正量为
根据式(2-121)和线性关系式(2-l05)，在节点注入功率不变的情况下，我们可以直接得到加线路k后的状态向量的增量为
式中：kekT为追加线路前支路k两端电压的相角差。新网络的状态向量为
这样我们就得到了追加线路k后，阻抗矩阵和状态向量的修正公式(2-118)和式(2-123)。当网络去掉或断开支路k时只要将xk换为xk，以上公式同样适用。
应该指出，当网络开断支路k使系统解列时，新的阻抗矩阵X不存在，这时式(2-119)中的k为无穷大，或xkk0。因此，应用直流潮流模型可以方便地找出网络中那些开断后引起系统解列的线路，对于这些线路不能直接进行断线分析。
2.6.3 N1检验与故障排序方法
目前比较常见的网络安全运行要求是满足N1检验，即在全部N条线路中任意开断一条线路后，系统的各项运行指标均应满足给定的要求。在网络规划形成网络结构的初期，最重要的原则是使网络不出现过负荷，即网络能够满足安全的输送电力的要求，为此我们应进行逐条线路开断后的过负荷校验。当任意一条线路开断后能够引起系统其他线路出现过负荷或系统解列时，说明网络没有满足N1检验。
严格的N1检验需要对全部线路进行N次断线分析，计算工作量很大。实际上，网络中有一些线路在开断后并不引起系统过负荷，因此我们可根据各线路开断后引起系统过负荷的可能性进行故障排序，然后按照顺序依次对过负荷可能性较大的线路进行校验。当校验到某条线路开断后不引起过负荷时，则排在其后的线路就可以不再进行校验，从而可以显著地减少计算量，这个过程也称为故障选择。目前国内外已出现了不少故障排序方法[25,26]，这些方法评判系统事故的标准各不相同。本节将介绍一种以是否引起系统过负荷作为标准的故障排序方法。
为了综合反映系统的过负荷情况，定义标量函数PI(Performance lndex)作为系统行为指标：
式中：Pl为线路l的有功潮流；Pl为线路l的传输容量；l为支路l中的并联线路数；l为
线路l的权系数，反映该线路故障对系统的影响；L为网络支路数。
由式(2-124)可以看出，当系统中没有过负荷时，Pl均不大于1，PI指标较Pl小。当系统中有过负荷时，过负荷线路的Pl大于1，正指数将使PI指标变得很Pl大。因此这个指标可以概括地反映系统安全性。为了突出地反映过负荷的情况，甚至可以用高次指数项代替式中的二次项。
通过分析PI指标对各线路导纳变化的灵敏度就可以反映出相应线路故障对系统安全性的影响。当线路k故障时，PI指标的变化量为
式中：Bk即Bk，为线路k的导纳。PIk的值越大，PI值增加越多，说明线路k障引起系统过负荷的可能性越大。
PIk可以用特勒根定理和伴随网络的方法进行计算，有兴趣的读者可参阅参考文献[3]。以下我们将推导一个利用正常情况潮流计算结果的直接计算PIk的公式。
设线路k开断后其他各线路潮流变为Pl(l1,2,L,L;lk)，这时系统行为指标相应地变为
显然
为了便于推导，我们将系统行为指标转化为相角的函数并用矩阵的形式表示。由式(2—113)可知
代入式（2-124）并定义
将式（2-114）代入式（2-129）,后者可进一步表示为
式中：
为一对称矩阵。由上式可知w具有与节点导纳矩阵B相同的结构，相当于以元素2llBl2Pl可取代Bl按形成导纳矩阵的算法直接形成w。这样，PI可表示为
式中：为线路k开断后的节点电压相角向量。
式(2-132)包含了线路k的有关项，但新的系统行为指标PI中不应当包含这一项，因此，将式（2-130）和式（2-133）代入式（2-127）可得
由式（2-123）可知
将以上两式代入式（2-134），有
考虑到矩阵X和w的对称性，令
式中：
将式（2-136）代入式（2-135），后式可简化为
对于开断线路k而言，以上各式中的k应为
将上式代入（2-138）后，可以得到
因为PkBkk，所以
式(2-138)-(2-140)只是表现形式不同，并无本质区别。这些公式中的各量均可由正常情况下的潮流计算数据求得。在已形成矩阵X,w,R,T的情况下，用这些公式计算各条线路开断后的PI值比较方便。
故障排序过程实际上是对所有线路按式(2-38)[或式(2-139)和式(2-140)]计算PI值，并根据PI从大到小排序。在断线分析时，首先对PI值最大的线路进行开断后的潮流计算和检验，直到开断某条线路后不再引起系统过负荷为止，其余PI值较小的线路引起系统过负荷的可能性很小，因而无需做断线分析。但是，采用这种系统行为指标可能存在一定的“退蔽”现象，例如当有个别线路过负荷而其他线路潮流较小时，其PI值可能小于没有过负荷但线路潮流都比较大时的PI值、因而根据这个指标进行故障选择排序可能会出现一定的误差。因此我们建议在实际应用时，应在连续校验几条线路故障都未引起系统过负荷的情况下才终止断线分析。2.7 静态安全分析的灵敏度法 2.7.1 节点功率方程的线性化
第2.6节中介绍的直流潮流模型是一种简单而快速的静态安全分析方法，但这种方法只能进行有功潮流的计算，没有考虑电压和无功问题。采用潮流计算的P—Q分解法和补偿法进行断线分析可以同时给出有功潮流、无功潮流以及节点电压的估计。但为了使计算结果达到一定的精度，要求必须进行反复迭代，否则其计算结果，特别是电压及无功潮流的误羌较大。我们将在本节介绍一种断线分析的灵敏度法[28]。这种方法将线路开断视为正常运行情况的一种扰动，从电力系统潮流方程的泰勒级数展开式出发，导出了灵敏度短阵，以节点注入功率的增量模拟断线的影响，较好地解决了电力系统断线分析计算问题。这种方法简单明了，省去了大量的中间计算过程，显著提高了断线分析的效率。应用本方法既可以提供全面的系统运行指标(包括有功、无功潮流，节点电压、相角)，又具有很高的计算精度和速度，因此是比较实用的静态安全分析方法。
网络断线分析还可以结合故障选择技术(见2.6.3节)，以减少断线分析的次数，进一步提高静态安全的效率。
如前所述，电力系统节点功率方程为[见式(2-9)]
式中：PiS,QiS分别为节点i的有功和无功功率注入量；其余各量的意义与式(2-9)相同
对于正常情况下的系统状态，式(2-141)可概括为
式中；W0为正常情况下节点有功、无功注入功率向量；X0为正常情况下由节点电压、相角组成的状态向量；Y0为正常情况的网络参数。
若系统注入功率发生扰动为W，或网络发生变化Y，状态变量也必然会出现变化，设其变化量为X，并满足方程
将式(2-143)按泰勒级数展开，则有，当扰动及状态改变量不大时，可以忽略X项及高次项，由于fX,Y是Y的线性函数，故fX,Y0。因此式（2-144）可简化为
将式（2-142）代入后，上式成为
由此可求出状态变量与节点功率扰动和网络结构变化的线性关系式为
当不考虑网络结构变化时，Y0式(2-146)成为
式中：
J0为潮流计算这代结束时的雅可比矩阵；S0则称为灵敏度矩阵。因为在潮流计算时J0已经进行了三角分解，所以S0很容易通过回代运算求出。
当不考虑节点注入功率的扰动时，W0式(2-146)变为
或经过变换可改写成如下形式：
式中：I为单位矩阵。
最后，我们得到
与式(2-147)相比，Wy可看作是由于断线而引起的节点注入功率的扰动：
上式中右端各项均可由正常情况的潮流计算结果求出，因此断线分析模拟完全是在正常接线及正常运行方式的基础上进行的。为了校验各种断线时的系统运行情况，只要按式(2-151)求出相应的节点注入功率增量Wy，然后就可利用正常情况下的灵敏度矩阵由式(2-150)直接求出状态变量的修正量。修正后系统的状态变量为
节点状态向量X已知后，即可按下式求出任意支路ij的潮流功率：
式中：tij为支路变比标之值，bij0为支路ij容纳的1/2 2.7.2 断线处节点注人功率增量的计算
断线分析的关键是按式(2-151)求出断线处节点注入功率增量Wy。静态安全校验主要是进行单线开断分析，但也可能涉及到多回线开断的情况，下面我们仅讨论单线开断的情况。对多回线开断情况感兴趣的读者，可参看文献[28]。为叙述方便，暂时假定系统中所有节点均为PQ节点，将式(2-151)简写为
式中：
Wt与断线支路在正常运行情况下的潮流有关。
设系统中总的支路数为b，断线支路两端节点为ij，则在b阶向量Y中只有与支路ij对应的元素为非零元素，即
对于一个节点数为N的网络来说，式(2-156)中的fyX0,Y0为2Nb阶矩阵，由式(2-141)可知，只有节点i和j的注入功率和支路ij的导纳有直接关系，即只有求节点i,j 的注入功率时才用到Gij和Bij。所以该矩阵每列只有4个非零元素。
设支路ij的阻抗角为ij，即
则有
利用以上关系和式（2-141），可以求得
将式（2-153）代入以上两式可得
同理可得到
式(2-158)和式(2-159)中的4个元素即为fYX0,Y0中对应于支路ij的4个非零元素，其他元素为
式中：ki,j表示k不属于节点集i,j。
综合式(2-157)-(2-160)，可得出式(2-156)的简化形式为
X0,Y0是一个2N2Nb阶矩阵，相式(2-155)中的L0为2N2N阶方阵，fxy当于用雅可比矩阵对各支路导纳元素求偏导，每条支路对应一个2N2N阶方阵，其结构如图2-15所示。
X0,Y0的矩阵结构 图2-15 fxy由于当ki,j且mi,j时有
所以对每条支路来说；2N2N阶矩阵中最多只有16个非零元素，它们由雅可比矩阵或由式(2-158)、式(2-159)求出：
同理可对Pj及Qj求出与式(2-163)类似的8个偏导数公式。
以上诸式中，Hij,Nij,Jij,Lij均为雅可比矩阵的元素：
由于Y中只有一个非零元素Yijyij，所以式(2-155)变为
式(2-1 65)中，只有对应于节点i、j两行两列交叉处2i1,2i,2j1,2j有非零元素，其余 元素均为零。
由以上讨论可知，在Wl及L0中只有与断线端点有关的元素才是非零元素，故式(2-154)可以写成更紧凑的形式：
式中：
式中：Sij1,Sij2,Sij3,Sij4等为灵敏度矩阵中行和列都与断线端点有关的元素，且有
式(2-166)中等式左边的向量表示断开线路ij时在节点i、j形成的节点注入功率增量，其他节点的增量为零。据此我们即可由式(2-150)求出各状态变量的修正量。式(2-166)是断线分析的主要公式，式中右端各项均可由牛顿法正常潮流计算结果获得。在形成H矩阵时只需进行两个4阶方阵的运算[见式(2-167)]，因而可以非常简便地求出由于断线引起的注入功率增量，快速进行静态安全分析。2.7.3 快速断线分析计算流程
快速断线分析方法的计算流程如图2-16所示。由图可知，在进行断线分析之前，首先要用牛顿法计算正常运行情况时的潮流，提供断线分析所需的数据。这些数据包括雅可比矩阵J0、灵敏度矩阵S0、正常情况下各节点电压相角和支路潮流等等。
断线分析计算包括3部分(以单线开断为例)：
(1)按式(2-166)求出相应的节点注入功率增量，其中主要的计算是按式(2-167)求出H矩阵。
(2)按式(2-150)求各节点状态变量的改变量，并按式(2-152)求出断线后新的状态变量。
(3)按式(2-153)求出断线后各支路潮流功率。
图2-16 快速断线分析计算流程图
应当指出，当断线使系统分解成两个不相连的子系统时，式(2-167)中H矩阵的逆矩阵不存在，因而不能直接进行断线分析.以上讨论我们曾假定所有节点均为PQ节点。实际上，当与断线相连的节点为PV节点时，在节点功率方程式(2-141)中只有一个与有功功率有关的方程，故断线分析只需计算该节点的有功功率增量，并认为无功功率增量为零，因此式(2-166)和式(2-167)中要除去与无功功率有关的行和列。当断线与系统平衡节点相连时，由于式(2-141)中不包含与平衡节点有关的方程，因此不求平衡节点注入功率的增量。这实际说明，PV节点的无功注入功率和平衡节点的有功及无功注入功率本身就是不定的，所以求它们的增量没有意义。
在静态安全校验中，如果只分析断线对某些关键节点的状态变量和关键支路潮流的影响，那么在图2-16的后两框中可只对这些节点和支路求断线后的数值，从而可进——步减少计算量。
【例2-32】 试对IEEE—14节点系统进行断线分析，并与牛顿—拉弗森法计算结果进行比较。表2-9给出了该系统的原始数据，其中有关数据己化为以100MVA为基准的标么值
【解】 根据断线分析计算流程图2-16，可确定计算步骤如下： 1)用牛顿法计算正常情况下的交流潮流。
当精度为0.0001时，对所给系统迭代3次可以收敛，其节点电压、相角及支路潮流均在表2-10中给出。
2)以断开线路5—6为例说明断线分析计算过程。①计算由于线路5—6开断而引起的节点注入功率增量 首先根据式(2-167)形成H矩阵。
由正常情况潮流计算结果和雅可比矩阵及灵敏度矩阵元素可知[雅可比矩阵和灵敏度矩阵己由潮流计算获知，这里未列出，此外，雅可比矩阵元素也可由式(2-164)算出]
然后由式（2-166）计算断线处的节点注入功率增量为
②根据式(2-150)求各状态变量的改变量 对节点2的相角而言，其改变量2为
同理可求出其他节点状态变量的改变量。
③根据式(2-152)求出各节点断线后新的状态变量。
将正常情况的状态变量与②中求出的状态改变量对应相加即可获得线路5—6开断后各节点新的状态变量。其值如表2-11中的第2、3列所示。表2-11中的第4、5列给出了该线开断情况下直接采用牛顿法计算的结果，表中最后两列为两种方法计算结果之差的绝对值。
由表2-11可以得出电压的平均误差为0.002040，最大误差为0.00517。相角的平均误差为0.00654，最大误差为0．01265。因此应用这种方法可以取得很好的精度，但计算时间却只有牛顿—拉弗森法迭代一次所需时间的1／7。3)断开其他线路时的计算结果。
为全面考察断线分析方法的计算精度，在表2-12中列出了断开其他线路时的计算结果。通过计算可知，总的电压平均误差为0.00478，相角平均误差为0.001199。在计算中可以获知线路5—6开断时的误差最大，这也正是前面选择这条线路为例的缘由。有关支路的情况及误差分析可参看参考文献[28]。
第二篇：电力系统潮流计算
南 京 理 工 大 学
《电力系统稳态分析》
课程报告
姓名
XX
学 号： 5*** 自动化学院 电气工程
基于牛顿-拉夫逊法的潮流计算例题编程报学院(系): 专
业: 题
目: 任课教师 硕士导师 告
杨伟 XX
2015年6月10号
基于牛顿-拉夫逊法的潮流计算例题编程报告
摘要：电力系统潮流计算的目的在于：确定电力系统的运行方式、检查系统中各元件是否过压或者过载、为电力系统继电保护的整定提供依据、为电力系统的稳定计算提供初值、为电力系统规划和经济运行提供分析的基础。潮流计算的计算机算法包含高斯—赛德尔迭代法、牛顿-拉夫逊法和P—Q分解法等，其中牛拉法计算原理较简单、计算过程也不复杂，而且由于人们引入泰勒级数和非线性代数方程等在算法里从而进一步提高了算法的收敛性和计算速度。同时基于MATLAB的计算机算法以双精度类型进行数据的存储和运算, 数据精确度高，能进行潮流计算中的各种矩阵运算，使得传统潮流计算方法更加优化。
一 研究内容
通过一道例题来认真分析牛顿-拉夫逊法的原理和方法（采用极坐标形式的牛拉法），同时掌握潮流计算计算机算法的相关知识，能看懂并初步使用MATLAB软件进行编程，培养自己电力系统潮流计算机算法编程能力。
例题如下：用牛顿-拉夫逊法计算下图所示系统的潮流分布，其中系统中5为平衡节点，节点5电压保持U=1.05为定值，其他四个节点分别为PQ节点，给定的注入功率如图所示。计算精度要求各节点电压修正量不大于10-6。
二 牛顿-拉夫逊法潮流计算 1 基本原理
牛顿法是取近似解x(k)之后，在这个基础上，找到比x(k)更接近的方程的根，一步步地迭代，找到尽可能接近方程根的近似根。牛顿迭代法其最大优点是在方程f(x)=0的单根附近时误差将呈平方减少，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点的电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成节点电压新的初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。2 基本步骤和设计流程图
形成了雅克比矩阵并建立了修正方程式，运用牛顿-拉夫逊法计算潮流的核心问题已经解决，已有可能列出基本计算步骤并编制流程图。由课本总结基本步骤如下：
1)形成节点导纳矩阵Y；
2)设各节点电压的初值，如果是直角坐标的话设电压的实部e和虚部f；如果是极坐标的话则设电压的幅值U和相角a；
3)将各个节点电压的初值代入公式求修正方程中的不平衡量以及修正方程的系数矩阵的雅克比矩阵；
4)解修正方程式，求各节点电压的变化量，即修正量； 5)计算各个节点电压的新值，即修正后的值；
6)利用新值从第（3）步开始进入下一次迭代，直至达到精度退出循环； 7)计算平衡节点的功率和线路功率，输出最后计算结果； ① 公式推导
② 流程图
三
matlab编程代码
clear;
% 如图所示1，2，3，4为PQ节点,5为平衡节点
y=0;
% 输入原始数据，求节点导纳矩阵
y(1,2)=1/(0.07+0.21j);
y(4,5)=0;y(1,3)=1/(0.06+0.18j);
y(1,4)=1/(0.05+0.10j);
y(1,5)=1/(0.04+0.12j);
y(2,3)=1/(0.05+0.10j);
y(2,5)=1/(0.08+0.24j);
y(3,4)=1/(0.06+0.18j);
for i=1:5
for j=i:5
y(j,i)=y(i,j);
end
end
Y=0;
% 求节点导纳矩阵中互导纳
for i=1:5
for j=1:5
if i~=j
Y(i,j)=-y(i,j);
end
end
end
% 求节点导纳矩阵中自导纳
for i=1:5
Y(i,i)=sum(y(i,:));
end
Y
% Y为导纳矩阵
G=real(Y);
B=imag(Y);% 输入原始节点的给定注入功率
S(1)=0.3+0.3j;
S(2)=-0.5-0.15j;
S(3)=-0.6-0.25j;
S(4)=-0.7-0.2j;
S(5)=0;
P=real(S);
Q=imag(S);
% 赋初值，U为节点电压的幅值，a为节点电压的相位角
U=ones(1,5);
U(5)=1.05;
a=zeros(1,5);
x1=ones(8,1);
x2=ones(8,1);
k=0;
while max(x2)>1e-6
for i=1:4
for j=1:4
H(i,j)=0;
N(i,j)=0;
M(i,j)=0;
L(i,j)=0;
oP(i)=0;
oQ(i)=0;
end
end
% 求有功、无功功率不平衡量
for i=1:4
for j=1:5
oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));
end
oP(i)=oP(i)+P(i);
oQ(i)=oQ(i)+Q(i);
end
x2=[oP,oQ]';
% x2为不平衡量列向量
% 求雅克比矩阵
% 当i~=j时，求H,N,M,L
for i=1:4
for j=1:4
if i~=j
H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));
N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
L(i,j)=H(i,j);
M(i,j)=-N(i,j);
end
end
end
% 当i=j时，求H,N,M,L
for i=1:4
for j=1:5
if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))
end
end
N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);
L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);
end
J=[H,N;M,L]
% J为雅克比矩阵
x1=-((inv(J))*x2);
% x1为所求△x的列向量
% 求节点电压新值，准备下一次迭代
for i=1:4
oa(i)=x1(i);
oU(i)=x1(i+4)*U(i);
end
for i=1:4
a(i)=a(i)+oa(i);
U(i)=U(i)+oU(i);
end
k=k+1;
end
k,U,a
% 求节点注入功率
i=5;
for j=1:5
P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i);
Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i);
end
S(5)=P(5)+Q(5)*sqrt(-1);
S
% 求节点注入电流
I=Y*U'
四
运行结果
节点导纳矩阵
经过五次迭代后的雅克比矩阵
迭代次数以及节点电压的幅值和相角（弧度数）
节点注入功率和电流
五 结果分析
在这次学习和实际操作过程里：首先，对电力系统分析中潮流计算的部分特别是潮流计算的计算机算法中的牛顿-拉夫逊法进行深入的研读，弄明白了其原理、计算过程、公式推导以及设计流程。牛顿-拉夫逊法是求解非线性方程的迭代过程，其计算公式为FJX，式中J为所求函数的雅可比矩阵；X为需要求的修正值；F为不平衡的列向量。利用x(*)=x(k+1)+X(k+1)进行多次迭代，通过迭代判据得到所需要的精度值即准确值x(*)。六 结论
通过这个任务，自己在matlab编程，潮流计算，word文档的编辑功能等方面均有提高，但也暴漏出一些问题：理论知识储备不足，对matlab的性能和特点还不能有一个全面的把握，对word软件也不是很熟练，相信通过以后的学习能弥补这些不足，达到一个新的层次。
第三篇：电力系统潮流计算程序设计
电力系统潮流计算程序设计
姓名：韦应顺
学号：2011021052 电力工程学院
牛顿—拉夫逊潮流计算方法具有能够将非线性方程线性化的特点，而使用MATLAB语言是由于MATLAB语言的数学逻辑强，易编译。
【】【】1.MATLAB程序12
Function tisco %这是一个电力系统潮流计算的程序 n=input（‘n请输入节点数：n=’）； m=input(‘请输入支路数：m=’);ph=input(‘n请输入平衡母线的节点号：ph=’)； B1=input（‘n请输入支路信号：B1=’）;%它以矩阵形式存贮支路的情况，每行存贮一条支路 %第一列存贮支路的一个端点 %第二列存贮支路的另一个端点 %第三列存贮支路阻抗
%第四列存贮支路的对地导纳
%第五列存贮变压器的变比，注意支路为1 %第六列存贮支路的序号
B2=input(‘n请输入节点信息：B2=’)； %第一列为电源侧的功率 %第二列为负荷侧的功率 %第三列为该点的电压值
%第四列为该点的类型：1为PQ，2为PV节点，3为平衡节点 A=input（‘n请输入节点号及对地阻抗：A=’）； ip=input（‘n请输入修正值：ip=’）； %ip为修正值);Y=zeros（n）；
Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）；
Y(p,q)=Y(p,q);f=zeros（1，n）；
no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2；
End for i=1：n
G=real(Y);if A（i2）=0
B=imag(Y);p=A(i1)；
Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while
max(abs(De))>ip&max(abs(Dfi)>ip
t=t+1;
e=e+De；
f=f+Df;
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)；
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)；
[De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v);
jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m
a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3));
b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2;
c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%该子程序是用来求取Df for i=1:n
If
i=ph
C(i)=0;
D(i)=0;
For j=i:n
C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end
P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If
B2(i4)=2 p=2*i-1;
Df(p)=P(i)-P1;p=p+1;else p=2*i-1;
Df(p)=P(i)-P1;p=p+1;
Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no=[];end
function [De,Df]=hxf(J,Df,ph,n,no)%该子函数是为求取De Df DX=JDf;DX1=DX;
x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0;
For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else
DX=[DX1,0,0];End k=0;
[x,y]=size(DX);For i=1:2:x K=k+1;
Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i);
X2=G(i,j)*e(i)+B(i,j)*f(i);%该子程序是用来求取jacci矩阵
for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if
J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End
2实例验证 【例题】设有一系统网络结线见图1，各支路阻抗和各节点功率均已以标幺值标示于图1中，其中节点2连接的是发电厂，设节点1电压保持U1＝1.06定值，试计算其中的潮流分布，请输入节点数：n=5 请输入支路数：m=7 请输入平衡母线的节点号：ph=l 请输入支路信息：
BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 请输入节点信息：
B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 请输入节点号及对地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 请输入修正值：ip=0．000 0l
参考文献
[1]陈珩．电力系统稳定分析[M]．北京：中国电力出版社，2002：139—187．
[2]郑阿奇．MATLAB实用教程[M]．北京：电子工业出版社，2005：1-243．
[3] 束洪春，孙士云，等．云电送粤交商流混联系统全过 程动态电压研究[J】．中国电力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn锄ic Voltage stabil时of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4．
[4] 朱新立，汤涌，等．大电网安全分析的全过程动态仿 真技术[J】．电网技术，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn锄ic simulation for the stabilhy a眦lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28．
[5]Roytelm锄I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．发电机励磁系统数学模型及参 数对电网动态稳定性分析结果影响的研究[J】．继电 器，2007，35(21)：22-27．
SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g锄e翰to璐baScd ∞de诅iled excitation system models柚d parameterS t0 power铲id dyn锄ic stabil时【J】．Relay，2007，35(2 1)：22-27．
[7] 方思立，朱方．快速励磁系统对系统稳定的影响[J】．中 国电机工程学报，1986，6(1)：20．28．
FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞
excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28．
[8] 刘取．电力系统稳定性及发电机励磁控制[M】．北京： 中国电力出版社，2007．
LIU Qu．Power system S诅bility锄d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro诅ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996．
[10] 陈利芳，陈天禄．浅谈自并励励磁系统在大容量机组 中的应用【J】．继电器，2007，35(1)：8l培4． CHEN Li-f抽岛CHEN Tian—lIL Application of 辩l仁exci组tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84．
[11] 方思立，刘增煌，孟庆和．大型汽轮发电机自并励励 磁系统的应用条件【J]．中国电力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu锄g，MENG Qin争hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63．
[12]梁小冰，黄方能．利用EMTDC进行长持续时间过程 的仿真研究【J】．电网技术，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan争眦ng．How to cany out simulalion of long dul‘ation processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57．
[13]王卉，陈楷，彭哲，等．数字仿真技术在电力系统中 的应用及常用的几种数字仿真工具【J】．继电器，2004，32(21)：7l一75．
wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眦hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75．
[14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI时studies【s】．
第四篇：电力系统潮流计算程序
电力系统潮流计算c语言程序，两行，大家可以看看，仔细研究，然后在这个基础上修改。谢谢
#include “stdafx.h” #include #include“Complex.h” #include“wanjing.h” #include“gauss.h” using namespace std; int _tmain(int argc, _TCHAR* argv[]){ int i; //i作为整个程序的循环变量 int N=Bus::ScanfBusNo();//输入节点个数 int L=Line::ScanflineNo();//输入支路个数 if((L&&N)==0){return 0;} //如果找不到两个文件中的任意一个，退出 Line *line=new Line[L];//动态分配支路结构体 Line::ScanfLineData(line);//输入支路参数 Line::PrintfLineData(line,L);//输出支路参数 Bus *bus=new Bus[N];//动态分配结点结构体 for(int i=0;i bus[i].Sdelta.real=0; bus[i].Sdelta.image=0;} Bus::ScanfBusData(bus);//输入节点参数 Bus::PrintfBusData(bus,N);//输出结点参数 Complex **X;X=new Complex *[N];for(i=0;i Bus::JisuanNodeDnz(X,line,bus,L,N);//计算节点导纳矩阵 Bus::PrintfNodeDnz(X,N);//输出节点导纳矩阵 int NN=(N-1)*2;double **JacAug;JacAug=new double *[NN];for(i=0;i double *x;x=new double[NN];int count=1; LOOP: Bus::JisuanNodeI(X,bus,N);//计算节点注入电流 Bus::JisuanNodeScal(X,bus,N);//计算节点功率 Bus::JisuanNodeScal(X,bus,N);//计算节点功率 Bus::JisuanNodeSdelta(bus,N);//计算节点功率差值 Bus::PrintfNodeScal(X,bus,N);//输出节点功率差值 int icon=wehcon1(bus,N);//whether converbence看迭代是否结束 if(icon==1){ cout<<“icon=”< Bus::JisuanJacAug(JacAug,X,bus,N);//计算雅可比增广矩阵 // Bus::PrintfJacAug(JacAug,N); gauss::gauss_slove(JacAug,x,NN);//解方程组求出电压差值 Bus::ReviseNodeV(bus,x,N);//修正节点电压 // Bus::PrintfNodeV(bus,N); count++; goto LOOP;} else { for(i=0;i { int statemp,endtemp; Complex aa,bb,cc,dd,B; B.real=0; B.image=-line[i].B; statemp=line[i].start; endtemp=line[i].end; aa=Complex::productComplex(Complex::getconj(bus[statemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].stoe=Complex::productComplex(bus[statemp-1].V,dd); aa=Complex::productComplex(Complex::getconj(bus[endtemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].etos=Complex::productComplex(bus[endtemp-1].V,dd); } cout<<“icon=”< Bus::JisuanNodeScal(X,bus,N);//计算节点功率 for(i=0;i { bus[i].Scal.real = bus[i].Scal.real + bus[i].Load.real;//发电机功率=注入功率+负荷功率 bus[i].Scal.image= bus[i].Scal.image+ bus[i].Load.image; bus[i].V=Complex::Rec2Polar(bus[i].V); } cout<<“====节点电压===============发电机发出功率======”< for(i=0;i { cout<<“节点”<<(i+1)<<'t'; Complex::PrintfComplex(bus[i].V); coutt(bus[i].Scal.real); coutt(bus[i].Scal.image); cout< } cout<<“======线路传输功率==========”< for(i=0;i { int statemp,endtemp; statemp=line[i].start; endtemp=line[i].end; cout< Complex::PrintfComplex(Complex::ComDivRea(line[i].stoe,0.01)); Complex::PrintfComplex(Complex::ComDivRea(line[i].etos,0.01)); cout< } } return 0;} #include “stdafx.h” #include #include“Complex.h” #include“wanjing.h” #include“gauss.h” using namespace std; int _tmain(int argc, _TCHAR* argv[]){ int i; //i作为整个程序的循环变量 int N=Bus::ScanfBusNo();//输入节点个数 int L=Line::ScanflineNo();//输入支路个数 if((L&&N)==0){return 0;} //如果找不到两个文件中的任意一个，退出 Line *line=new Line[L];//动态分配支路结构体 Line::ScanfLineData(line);//输入支路参数 Line::PrintfLineData(line,L);//输出支路参数 Bus *bus=new Bus[N];//动态分配结点结构体 for(int i=0;i bus[i].Sdelta.real=0; bus[i].Sdelta.image=0;} Bus::ScanfBusData(bus);//输入节点参数 Bus::PrintfBusData(bus,N);//输出结点参数 Complex **X;X=new Complex *[N];for(i=0;i Bus::JisuanNodeDnz(X,line,bus,L,N);//计算节点导纳矩阵 Bus::PrintfNodeDnz(X,N);//输出节点导纳矩阵 int NN=(N-1)*2;double **JacAug;JacAug=new double *[NN];for(i=0;i double *x;x=new double[NN];int count=1; LOOP: Bus::JisuanNodeI(X,bus,N);//计算节点注入电流 Bus::JisuanNodeScal(X,bus,N);//计算节点功率 Bus::JisuanNodeScal(X,bus,N);//计算节点功率 Bus::JisuanNodeSdelta(bus,N);//计算节点功率差值 Bus::PrintfNodeScal(X,bus,N);//输出节点功率差值 int icon=wehcon1(bus,N);//whether converbence看迭代是否结束 if(icon==1){ cout<<“icon=”< Bus::JisuanJacAug(JacAug,X,bus,N);//计算雅可比增广矩阵 // Bus::PrintfJacAug(JacAug,N); gauss::gauss_slove(JacAug,x,NN);//解方程组求出电压差值 Bus::ReviseNodeV(bus,x,N);//修正节点电压 // Bus::PrintfNodeV(bus,N); count++; goto LOOP;} else { for(i=0;i { int statemp,endtemp; Complex aa,bb,cc,dd,B; B.real=0; B.image=-line[i].B; statemp=line[i].start; endtemp=line[i].end; aa=Complex::productComplex(Complex::getconj(bus[statemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].stoe=Complex::productComplex(bus[statemp-1].V,dd); aa=Complex::productComplex(Complex::getconj(bus[endtemp-1].V), B); bb=Complex::subComplex (Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V)); cc=Complex::productComplex(bb , Complex::getconj(line[i].Y)); dd=Complex::CaddC(aa,cc); line[i].etos=Complex::productComplex(bus[endtemp-1].V,dd); } cout<<“icon=”< Bus::JisuanNodeScal(X,bus,N);//计算节点功率 for(i=0;i { bus[i].Scal.real = bus[i].Scal.real + bus[i].Load.real;//发电机功率=注入功率+负荷功率 bus[i].Scal.image= bus[i].Scal.image+ bus[i].Load.image; bus[i].V=Complex::Rec2Polar(bus[i].V); } cout<<“====节点电压===============发电机发出功率======”< for(i=0;i { cout<<“节点”<<(i+1)<<'t'; Complex::PrintfComplex(bus[i].V); coutt(bus[i].Scal.real); coutt(bus[i].Scal.image); cout< } cout<<“======线路传输功率==========”< for(i=0;i { int statemp,endtemp; statemp=line[i].start; endtemp=line[i].end; cout< Complex::PrintfComplex(Complex::ComDivRea(line[i].stoe,0.01)); Complex::PrintfComplex(Complex::ComDivRea(line[i].etos,0.01)); cout< } } return 0;} #include class Complex//定义复数类 { public: double real;double image;int RecPolar;//0表示直角坐标，1表示极坐标 static Complex CaddC(Complex c1,Complex c2);//求两个复数和 static Complex subComplex(Complex c1,Complex c2);//求两个复数差 static Complex productComplex(Complex c1,Complex c2);//求两个复数积 static Complex divideComplex(Complex c1,Complex c2);//求两个复数商 static Complex ComDivRea(Complex c1,double r2);//除数 static Complex getconj(Complex c1);//求一个复数共轭 static Complex getinverse(Complex c1);//取倒数 static double getComplexReal(Complex c1);//求一个复数实部 static double getCompleximage(Complex c1);//求一个复数虚部 static void PrintfComplex(Complex c1);//显示一个复数 static void PrintfmultiComplex(Complex C,int N);//显示多个复数 static void zeroComplex(Complex c1);//将复数复零 static Complex Rec2Polar(Complex c1);//取极坐标 Complex(){ RecPolar=0;} }; Complex Complex::Rec2Polar(Complex c1)//极坐标表示 { Complex Node;Node.real=sqrt(c1.real*c1.real+c1.image*c1.image);Node.image=atan2(c1.image,c1.real)*180/3.1415926;Node.RecPolar=1;return Node;} Complex Complex::CaddC(Complex c1,Complex c2)//复数加法 { Complex Node; Node.real=c1.real+c2.real; Node.image=c1.image+c2.image; return Node;} Complex Complex::subComplex(Complex c1,Complex c2)//复数减法 { Complex Node; Node.real=c1.real-c2.real; Node.image=c1.image-c2.image; return Node;} Complex Complex::productComplex(Complex c1,Complex c2)//复数乘法 { Complex Node; Node.real=c1.real*c2.real-c1.image*c2.image; Node.image=c1.image*c2.real+c2.image*c1.real; return Node;} Complex Complex::divideComplex(Complex c1,Complex c2)//复数除法 { Complex Node; Node.real=(c1.real*c2.real+c1.image*c2.image)/(pow(c2.real,2)+pow(c2.image,2));Node.image=(c1.image*c2.real-c1.real*c2.image)/(pow(c2.real,2)+pow(c2.image,2));return Node;} Complex Complex::ComDivRea(Complex c1,double r1)//复数除数 { Complex Node;Node.real=c1.real/(r1);Node.image=c1.image/(r1);return Node;} Complex Complex::getconj(Complex c1)//取共轭 { Complex Node; Node.real=c1.real;Node.image=-c1.image; return Node;} Complex Complex::getinverse(Complex c1)//取倒数 { Complex Node;Node.real=1;Node.image=0;Node=(Complex::divideComplex(Node,c1));return Node;} double Complex::getComplexReal(Complex c1)//取实部 { return c1.real;} double Complex::getCompleximage(Complex c1)//取虚部 { return c1.image;} void Complex::PrintfComplex(Complex c1)//按直角坐标输出 { if(c1.RecPolar==0){ cout.precision(6); cout.width(8); cout.setf(ios::right); cout< ”; cout.precision(6); cout.width(8); cout.setf(ios::left); cout< ”;} else { cout< Complex::zeroComplex(Complex c1)//清零 { c1.real=0;c1.image=0;} class gauss { public: static void gauss_slove(double **a,double *x,int NN);static void gauss_output();}; void gauss::gauss_slove(double **a,double *x,int NN){ int n,i,j,k,*pivrow,**pivcol;double eps,pivot,sum,aik,al; n=NN;pivrow=new int[n];pivcol=new int *[n]; for(i=0;i pivot= fabs(a[k][k]); pivrow[k]=k;//行 pivcol[k][0]=k;pivcol[k][1]=k;//列n*2矩阵 for(i=k;i { for(j=k;j { if(pivot { pivot=fabs(a[i][j]); pivrow[k]=i;//行 pivcol[k][1]=j;//列 } } } if(pivot { cout<<“error”< getchar(); exit(0); } if(pivrow[k]!=k)//行变换 { for(j=k;j<(n+1);j++) { al=a[pivrow[k]][j]; a[pivrow[k]][j]=a[k][j]; a[k][j]=al; } } if(pivcol[k][1]!=k)//列变换 { for(i=0;i { al=a[i][pivcol[k][1]]; a[i][pivcol[k][1]]=a[i][k]; a[i][k]=al; } } if(k!=(n-1))//将矩阵化为上三角形 式 { for(i=(k+1);i { aik=a[i][k]; for(j=k;j<(n+1);j++) { a[i][j]-=aik*a[k][j]/a[k][k]; } } } } x[n-1]=a[n-1][n]/a[n-1][n-1];//解方程 for(i=(n-2);i>=0;i--){ sum=0; for(j=(i+1);j { sum +=a[i][j]*x[j];0.182709 0.016894-0.0310701 -0.0402051 0.156702 -0.0355909-0.0668055 -0.00703229-0.0886481 -0.0129814-0.0390805 -0.0135062-0.1023 -0.0460568 -0.0342827 -0.00382402-0.102896 -0.0184062 } x[i]=(a[i][n]-sum)/a[i][i];} for(k=(n-2);k>=0;k--){ al=x[pivcol[k][1]]; x[pivcol[k][1]]=x[pivcol[k][0]]; x[pivcol[k][0]]=al;} cout<<“节点电压修正量”< cout< } ====节点功率计算值==== 0.935261 -0.159048 0.573909 0.0789973-0.00289889 -0.00796623-0.0791247 -0.0168362-0.436255 -0.0580392 0.0359139 -0.0106592-0.229118 -0.0885419-0.136179 -0.148207 0.0446243 0.0111298-0.0223764 -0.00695775-0.0237482 -0.198318 -5.24266e-015 -0.0354071 -0.0925078 -1.05629e-015 -0.0391348 0.014529 0.00158644 -0.0258771 -0.109514 icon=1进行第2次迭代 节点电压修正量 =================-0.00164889-0.000540034-0.00261067-0.00532027-0.00235315-0.00600971-0.00189677-0.00643874-0.0023631-0.00650659-0.00170949-0.0074907-0.00164545-0.00485415-0.00493977-0.0119042-0.00331285-0.0175611-0.00207908 -0.00347744-0.0869347-9.48909e-015-0.0110778-0.0538236-7.53784e-016-0.0168097 7.049e-005-0.00146487-0.00458276 0.00251645 -0.00336375-0.00530645-0.0147816-0.000326161-0.00640487-0.00251701-0.0169829-0.00175286-0.0174333-0.0239063 -0.0119192-0.076014 -0.0160104-0.441997 -0.0750285 0.000250012 3.72542e-005-0.228052 -0.108844-0.100078 -0.105634 0.000410707 0.000378067-0.057497 -0.0195879 0.200039 0.0582563-0.00307326-0.0163809-0.00232773-0.0175806 8.74293e-005-0.0192018 0.000558996-0.0197776-0.000247851-0.0193784-0.00115346-0.0185848-0.00127275-0.0186244-0.00010108-0.0188966 0.000553585-0.0200901-3.76315e-005-0.0208303 0.00308341-0.0219386-0.00195916-0.0205356-0.00184757-0.0076401 0.00197593-0.0245534 0.00434657-0.027534 ====节点功率计算值==== 0.98623 -0.134163 0.583136 0.166278-0.111173 0.199792 -0.0621041 -0.0821379 -0.0350785 -0.0902383 -0.0320461 -0.0951562 -0.0220362 -0.175458 4.72557e-015 -0.0320661 -0.0871134 -7.03489e-017 -0.0350769 0.000273455 1.51804e-005 -0.0240417 -0.10604 icon=1进行第3次迭代 节点电压修正量 =================-2.67079e-005-2.30128e-006-2.20543e-005-6.00686e-005-2.33043e-005-6.85601e-005-3.22294e-005-2.61107e-005-2.80198e-005-6.6167e-005-2.34528e-005 -0.0739846 0.0227868-0.0158709-0.0248173-0.0179447-0.0578368-0.00890719-0.0337091-0.00693706-0.111601 1.21429e-014-0.0159145-0.0667319 9.24355e-016-0.0228592 7.10354e-005-6.6188e-006-0.00889343-0.0184098 -5.66132e-005-4.4646e-005-1.74668e-005-4.50947e-005-0.000181763-3.81763e-006-0.000286581-6.68993e-005-1.28441e-005-5.17172e-005-0.000223284-4.54717e-005-2.47586e-005 4.32335e-007-0.000258494 1.82635e-005-0.000272051-6.95195e-006-0.000251969 1.11318e-005-0.000279418 5.74737e-005-0.000307368 6.86998e-005-0.000320274 5.38112e-005-0.00031447 3.59531e-005-0.00030494 3.37607e-005-0.000307449 5.26532e-005-0.000310721 6.92761e-005-0.000350373 5.60942e-005-0.00040977 0.000123641-0.000440259 1.36149e-005-0.000426973-1.70227e-005-9.37794e-005 0.000113675-0.000544011 0.000176034-0.000636202 ====节点功率计算值==== 0.986878 -0.133979 0.583 0.167193-0.024 -0.012-0.076 -0.016-0.442 -0.0748606 1.43501e-008 1.07366e-008-0.228 -0.109 -0.0999999 -0.104049 4.51318e-008 8.98835e-008-0.0579999 -0.0199999 0.2 0.0591018-0.112 -0.0749997 0.2 0.0242519-0.062 -0.016-0.082 -0.025-0.035 -0.018 -0.0900001 -0.058-0.032 -0.00899997-0.095 -0.0339999-0.022 -0.00699998-0.175 -0.112 -6.07156e-015 -1.19217e-014-0.032 -0.016-0.087 -0.0669999 7.03078e-017 -9.23979e-016-0.035 -0.0229999 1.09492e-007 4.45699e-008 1.54958e-009 -2.01531e-010-0.024 -0.00899994-0.106 -0.0189996 icon=0,迭代结束。 ====节点电压===============发电机发出功率====== 节点1 1.05 0。 98.6878-13.3979 节点2 1.045 -1.846。 29.4193 节点3 1.02384-3.83352。 0 节 点25 1.01216-9.68486。 0 0 0 节点4 1.01637-4.55698。 0 节 点26 0.994393 -10.1089。 0 0 0 节点5 1.01 -6.48617。 节 点27 1.02012-9.42025。 0 11.5139 0 节点6 1.01332-5.38073。 0 节 点28 1.00992-5.86244。 0 0 0 节点7 1.00489-6.38368。 0 节 点29 1.00022-10.6579。 0 0 节点8 19.5951 节点9 0 节点10 0 节点11 5.91018 节点12 0 节点13 2.42519 节点14 0 节点15 0 节点16 0 节点17 0 节点18 0 节点19 0 节点20 0 节点21 0 节点22 0 节点23 0 节点24 0 1.01 -5.62974。 1.03905-6.78143。 1.03595-8.69362。 -4.5962。 1.04711-7.80323。 1.05 -6.34392。 1.03242-8.7401。 1.02788-8.86784。 1.03458-8.45044。 1.03051-8.83678。 1.01845-9.5141。 1.01604-9.70326。 1.02022-9.50938。 1.0237-9.17478。 1.02432-9.17024。 1.01802-9.36719。 1.01339-9.68362。 0 20 节 点30 0.988705 -11.5464。 0 0 0 ====== 线路传输功率========== 2to1 -57.7373 5.41674i 58.3454 0 -15.1827i 3to1 -39.659 -7.75964i 40.3424 1.78481i 4to2 -30.87 -9.74186i 31.4153 0 3.58352i 4to3 -37.0772 -7.78596i 37.259 6.55964i 5to2 -44.3717 -9.78456i 45.2968 0 4.84242i 6to2 -38.4766 -8.22625i 39.3252 0 2.87667i 6to4 -34.946 1.92384i 35.0885 0 -3.28202i 7to5 -0.16304 -6.41767i 0.171702 0 2.2985i 7to6 -22.637 -4.48233i 22.7745 0 1.44238i 8to6 -11.8939 -5.48098i 11.913 0 3.70557i 6to9 12.3737 -12.3826i -12.3737 0 13.0033i 6to10 10.9107 -3.80907i -10.9107 0 4.53223i 11to9 5.91018i 0 -5.08963i 10to9 -32.652 -2.3712i 32.652 0 3.46974i 4to12 23.5411 -11.5375i -23.5411 0 13.2407i 13to12 2.42519i 1.05 -1.90978i 1.66484i 14to12 -7.9019 -2.06732i 7.97894 30to29 -3.6702 -0.542564i 3.70398 2.22749i 0.606393i 15to12 -18.254 -5.74885i 18.4835 28to8 -1.89152 -3.79982i 1.89395 6.20089i-4.9239i 16to12-7.53872 -2.90237i 7.59633 28to6 -14.7868 -2.82565i 14.8234 3.02352i 0.294601i 15to14-1.69544 -0.461488i 1.70189 请按任意键继续...0.467323i 17to16-4.03014 1.10238i 18to15-6.08074 1.46028i 19to18-2.87549 0.478389i 20to19 6.6418-2.93222i 20to10 -8.8418 3.85077i 17to10-4.96987 4.76656i 21to10-16.1562 9.42843i 22to10-7.87782 4.21401i 22to21 1.34443-2.01837i 23to15-5.59369 2.25006i 24to22-6.48186 2.08163i 24to23-2.38596 0.579814i 25to24-0.167617 0.281364i 26to25 -3.5 2.3674i 27to25 3.39433-2.08638i 28to27 16.1446 3.13006i 29to27-6.10398 1.67047i 30to27-6.92979-1.07089i-1.37839i-0.467767i 2.96679i-3.66679i-4.72911i-9.18162i-4.10132i 2.01969i-2.17981i-2.00141i-0.56401i -0.28102i-2.29999i 2.11848i-2.10093i-1.50639i -1.3574i 4.03872 6.12096 2.88074 -6.62452 8.9242 4.98423 16.2709 7.93248 -1.34378 5.62846 6.53339 2.39369 0.167814 3.54513 -3.37751 -16.1446 6.19083 7.09313 高等电力系统分析 IEEE30节点潮流程序 班级：电研114班 姓名：王大伟 学号：2201100151 电力系统潮流计算发展史 对潮流计算的要求可以归纳为下面几点： （1）算法的可靠性或收敛性（2）计算速度和内存占用量（3）计算的方便性和灵活性 电力系统潮流计算属于稳态分析范畴，不涉及系统元件的动态特性和过渡过程。因此其数学模型不包含微分方程，是一组高阶非线性方程。非线性代数方程组的解法离不开迭代，因此，潮流计算方法首先要求它是能可靠的收敛，并给出正确答案。随着电力系统规模的不断扩大，潮流问题的方程式阶数越来越高，目前已达到几千阶甚至上万阶，对这样规模的方程式并不是采用任何数学方法都能保证给出正确答案的。这种情况促使电力系统的研究人员不断寻求新的更可靠的计算方法。 在用数字计算机求解电力系统潮流问题的开始阶段，人们普遍采用以节点导纳矩阵为基础的高斯-赛德尔迭代法（一下简称导纳法）。这个方法的原理比较简单，要求的数字计算机的内存量也比较小，适应当时的电子数字计算机制作水平和电力系统理论水平，于是电力系统计算人员转向以阻抗矩阵为主的逐次代入法（以下简称阻抗法）。 20世纪60年代初，数字计算机已经发展到第二代，计算机的内存和计算速度发生了很大的飞跃，从而为阻抗法的采用创造了条件。阻抗矩阵是满矩阵，阻抗法要求计算机储存表征系统接线和参数的阻抗矩阵。这就需要较大的内存量。而且阻抗法每迭代一次都要求顺次取阻抗矩阵中的每一个元素进行计算，因此，每次迭代的计算量很大。 阻抗法改善了电力系统潮流计算问题的收敛性，解决了导纳法无法解决的一些系统的潮流计算，在当时获得了广泛的应用，曾为我国电力系统设计、运行和研究作出了很大的贡献。但是，阻抗法的主要缺点就是占用计算机的内存很大，每次迭代的计算量很大。当系统不断扩大时，这些缺点就更加突出。为了克服阻抗法在内存和速度方面的缺点，后来发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需存储各个地区系统的阻抗矩阵及它们之间的联络线的阻抗，这样不仅大幅度的节省了内存容量，同时也提高了节省速度。 克服阻抗法缺点的另一途径是采用牛顿-拉夫逊法（以下简称牛顿法）。牛顿法是数学中求解非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳矩阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流程序的计算效率。自从20世纪60年代中期采用了最佳顺序消去法以后，牛顿法在收敛性、内存要求、计算速度方面都超过了阻抗法，成为直到目前仍被广泛采用的方法。 在牛顿法的基础上，根据电力系统的特点，抓住主要矛盾，对纯数学的牛顿法进行了改造，得到了P-Q分解法。P-Q分解法在计算速度方面有显著的提高，迅速得到了推广。 牛顿法的特点是将非线性方程线性化。20世纪70年代后期，有人提出采用更精确的模型，即将泰勒级数的高阶项也包括进来，希望以此提高算法的性能，这便产生了保留非线性的潮流算法。另外，为了解决病态潮流计算，出现了将潮流计算表示为一个无约束非线性规划问题的模型，即非线性规划潮流算法。 近20多年来，潮流算法的研究仍然非常活跃，但是大多数研究都是围绕改进牛顿法和P-Q分解法进行的。此外，随着人工智能理论的发展，遗传算法、人工神经网络、模糊算法也逐渐被引入潮流计算。但是，到目前为止这些新的模型和算法还不能取代牛顿法和P-Q分解法的地位。由于电力系统规模的不断扩大，对计算速度的要求不断提高，计算机的并行计算技术也将在潮流计算中得到广泛的应用，成为重要的研究领域。第五篇：电力系统潮流计算发展史
点击咨询 