电力系统稳态分析-潮流计算

1、-. z. PAGE 1电力系统稳态分析 摘 要电力系统潮流计算是研究电力系统稳态运行情况的一种重要的分析计算，它根据给定的运行条件及系统接线情况确定整个电力系统各局部的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗。所以，电力系统潮流计算是进展电力系统故障计算，继电保护整定，平安分析的必要工具。本文介绍了基于MATLAB软件的牛顿-拉夫逊法和P-Q分解法潮流计算的程序，该程序用于计算中小型电力网络的潮流。在本文中，采用的是一个5节点的算例进展分析，并对仿真结果进展比拟，算例的结果验证了程序的正确性和迭代法的有效性。关键词：电力系统潮流计算；MATLAB；牛顿-拉夫逊法；P-Q分解

2、法；-. z.目 次 TOC o 1-2 f h z u HYPERLINK l _Toc3431630271 绪论 PAGEREF _Toc343163027 h 1HYPERLINK l _Toc3431630281.1背景及意义 PAGEREF _Toc343163028 h 1HYPERLINK l _Toc3431630291.2 相关理论 PAGEREF _Toc343163029 h 1HYPERLINK l _Toc3431630301.3 本文的主要工作 PAGEREF _Toc343163030 h 2HYPERLINK l _Toc3431630312 潮流计算的根本理论

3、 PAGEREF _Toc343163031 h 3HYPERLINK l _Toc3431630322.1 节点的分类 PAGEREF _Toc343163032 h 3HYPERLINK l _Toc3431630332.2 根本功率方程式极坐标下 PAGEREF _Toc343163033 h 3HYPERLINK l _Toc3431630342.3 本章小结 PAGEREF _Toc343163034 h 4HYPERLINK l _Toc3431630353 潮流计算的两种算法 PAGEREF _Toc343163035 h 5HYPERLINK l _Toc3431630363.

4、1 牛顿拉夫逊算法 PAGEREF _Toc343163036 h 5HYPERLINK l _Toc3431630373.2 PQ分解算法 PAGEREF _Toc343163037 h 10HYPERLINK l _Toc3431630383.3 本章小结 PAGEREF _Toc343163038 h 15HYPERLINK l _Toc3431630394 算例 PAGEREF _Toc343163039 h 16HYPERLINK l _Toc3431630404.1系统模型 PAGEREF _Toc343163040 h 16HYPERLINK l _Toc3431630414.2

5、 结果分析 PAGEREF _Toc343163041 h 16HYPERLINK l _Toc3431630424.3 本章小结 PAGEREF _Toc343163042 h 19HYPERLINK l _Toc343163043结论 PAGEREF _Toc343163043 h 20HYPERLINK l _Toc343163044参考文献 PAGEREF _Toc343163044 h 21HYPERLINK l _Toc343163045附录 PAGEREF _Toc343163045 h 22-. z.1 绪论1.1背景及意义电力系统稳态分析是研究电力系统运行和规划方案最重要和最

6、根本的手段。电力系统稳态分析根据给定的发电运行方式和系统接线方式来确定系统的稳态运行状态，其中潮流计算针对电力系统的各种正常的运行方式进展稳态分析。潮流计算是根据给定的电网构造、参数和发电机、负荷等元件的运行条件，确定电力系统各局部稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题，其解法都离不开迭代。潮流计算方法的改良过程中，经历了高斯-赛德尔迭代法、阻抗法、分块阻抗法、牛顿-拉夫逊法、

7、改良牛顿法、P-Q分解法等。现在比拟常用的方法就是牛顿-拉夫逊法和P-Q分解法。对潮流计算的要求可以归纳为下面几点：1计算方法的可靠性或收敛性；2对计算机存量的要求；3计算速度；4计算的方便性和灵活性。1.2 相关理论所谓潮流计算，就是电网的接线方式与参数及运行条件，计算电力系统稳态运行各母线电压、各支路电流与功率及网损。对于正在运行的电力系统，通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提

8、供原始数据。在运行方式管理中，潮流是确定电网运行方式的根本出发点；在规划领域，需要进展潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了多个在预想操作情况下电网的潮流分布以校验运行可靠性。在电力系统调度运行的多个领域都涉及到电网潮流计算。潮流是确定电力网络运行状态的根本因素，潮流问题是研究电力系统稳态问题的根底和前提。1.3 本文的主要工作本文介绍了电力系统潮流计算方法中的牛顿-拉夫逊法和PQ快速分解法的相关知识及其根本原理，并用MATLAB编写程序，最后通过一个5节点的算例来验证该程序的正确性，并对两种算法的结果进展了分析，比照了两种算法。-. z.2 潮流计算的根本理论2.1

9、节点的分类用一般的电路理论求解网络方程，目的是给出电压源或电流源研究网络的电流或电压分布，作为根底的方程式，一般用线性代数方程式表示。然而在电力系统中，给出发电机或负荷连接母线上电压或电流都是向量的情况是很少的，一般是给出发电机母线上发电机的有功功率和母线电压的幅值，给出负荷母线上负荷消耗的有功功率和无功功率。主要目的是由这些量去求电力系统的各种电气量。所以，根据电力系统中各节点性质的不同，很自然的把节点分成三种类型。1PQ节点对这一类节点，事先给定的是节点有功功率和无功功率P、Q，待求的未知量是节点电压向量(V，)，所以叫PQ节点。通常变电所母线都是PQ节点，当*些发电机的输出功率P、Q给定

10、时，也作为PQ节点。在潮流计算中，系统大局部节点属于PQ节点。(2)PV节点这类节点给出的参数是该节点的有功功率P及电压幅值V，待求量为该节点的无功功率Q及电压向量的相角。这类节点在运行中往往要有一定可调节的无功电源，用于维持给定的电压值。通常选择有一定无功功率储藏的发电机母线或者变电所有无功补偿设备的母线作PV节点处理。(3)平衡节点在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在计算中取该节点电压向量的方向作为参考轴，相当于给定该点电压向量的角度为零。也就是说，对平衡节点给定的运行参数是V和，因此又称为V节点，而待求量是该节点的P，Q，整个系统的功率平衡由这一节点承当。关于

11、平衡节点的选择，一般选择系统中担任调频调压的*一发电厂或发电机，有时也可能按其他原则选择，例如，为提高计算的收敛性，可以选择出线数多或者靠近电网中心的发电厂母线作平衡节点。以上三类节点4个运行参数P、Q、V、中，量都是两个，待求量也是两个，只是类型不同而已。2.2 根本功率方程式极坐标下在潮流计算中任何复杂的电力系统都可以归结为以下元件参数组成：发电机注入电流或功率；负荷注入负的电流或功率；输电线支路电阻、电抗；变压器支路电阻、电抗、变比；母线上对地支路阻抗和导纳；线路上的对地支路一般为线路充电电容导纳。必须指出，如果仅研究稳态情况下的潮流而不涉及暂态过程的计算，则不需要发电机和负荷的阻抗参数

12、，只需要给出发电机和负荷的注入功率和电流，并且规定发电机和负荷的注入功率和电流取正，而负荷取负。在潮流计算中，节点功率可表示为： i=1,2,.n) 2.1假设把电压表示为极坐标的形式，即 2.2 将导纳矩阵中元素表示为 2.3这样，我们可以得到： (i=1,2,.n) 2.4按实部和虚部展开，得到 2.5 上式就是功率的极坐标方程式。该方程组在牛顿法和P-Q分解法中起到了重要作用。2.3 本章小结 在本章里主要介绍了电力系统的潮流计算的根本的理论。首先对电力系统中三种节点进展了详细的阐述；其次，介绍了在极坐标的情况下的电力系统的潮流计算的功率方程，为下文的潮流计算分析打下根底。3 潮流计算的

13、两种算法3.1 牛顿拉夫逊算法牛顿-拉夫逊算法产生于50 年代末期，是一种实用且有竞争力的电力系统潮流计算方法，求解非线性方程式的典型方法。在稀疏矩阵技巧和高斯消去法被应用以后，其真正的价值才表达出来。该方法有较好的收敛性，迭代次数少，在电力系统潮流计算中也得到应用。目前，牛顿法潮流计算是最为广泛、效果最好的一种潮流计算方法。3.1.1 根本原理该方法把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，通常称为逐次线性化过程，经过一步步的迭代，得到最终的结果。设有非线性方程式： 3.1设为该方程式的初值。 3.2为初值的修正量。将式3.2代入式3.1，可以得到 3.3按泰勒级数展开。

14、因为很小，二次以及二次以上的各项均可以略去，可以得到简化的方程式： 3.4上式是对于变量修正量的线性方程式，即修正方程式，解得： 3.5上式求得的不是方程真正的解，我们需要进展反复的迭代，一步步的趋近方程式的解，得到最逼近的值。这样反复下去，就构成了不断求解非线性方程式的逐次线性化过程。第t次迭代时的参数方程为 3.6当时，就满足了原方程式3.1，因而就成为该方程的解。式中是函数 在点的一次导数，也就是曲线在点的斜率，如图3.1所示， 3.7修正量则是由点的切线与横轴的交点来确定，由图3.1可以直观的看出牛顿法的求解过程。图3.1 牛顿法的几何解释现在把牛顿法推广到多变量非线性方程组的情况。设

15、有变量的非线性联立方程组： (3.8)给定各变量初值，假设为其修正量，并使其满足 3.9对以上个方程式分别按泰勒级数展开，当忽略所组成的二次项和高次项时，可以得到 3.10式中：为函数对自变量的偏导数在点处的值。把上式写成矩阵形式： 3.11一般第t次迭代式的修正方程为 (3.12)上式可以简写为 (3.13)同样，也可以写出： (3.14)上式中，为第t次迭代时的雅克比矩阵，为第t次迭代时的修正量向量。这样交替反复求解就可以使趋近方程组的真正解。为了判断收敛情况，可采用一下两个不等式中的一个： (3.15) (3.16)式中，为预先给定的很小正数。3.1.2 算法求解分析将式2.5的功率方程

16、式代入到上面的分析之中，将这两个方程改写成残差的形式，即：3.17其中：为节点i和j之间的电压相角差；和分别为支路电导和电纳；和分别为节点i 和j的电压向量。对式3.16进展泰勒级数展开，取一次项近似，即可得到牛顿法潮流计算的修正方程式，即：3.18式中，电压幅值的修正量采用的形式没有特殊意义，只是为了使雅克比矩阵中各元素具有比拟相似的表达式。其中：和为潮流方程的有功功率和无功功率残差向量；和为母线电压修正量；系数J 为雅可比矩阵对方程式2.9进展变换即可得到变量和的求解公式，即： 3.19雅可比矩阵各元素可表达为： 3.20 3.21 3.22 3.23牛顿法潮流计算的具体步骤大致分为：输入

17、原始数据计算节点导纳矩阵；给出各节点电压初值；将电压初值代入,求出。判断是否满足收敛条件，如果满足，则停顿计算。否则，继续进展下面的步骤；将电压初值代入求出雅可比矩阵J；解潮流残差方程，求出节点电压的修正量；修正节点电压向量；判断是否满足收敛条件，如果满足，则停顿计算，否则，再以 为初值，返回第步进展下一次迭代。= 8 * GB3 计算支路功率分布，PV节点无功功率和平衡节点注入功率。= 9 * GB3 输出结果，并完毕。 流程图如下所示：图3.2 牛顿法流程图使用牛顿-拉夫逊法有以下优点：1收敛速度快，具有平方收敛特性，迭代次数与系统规模根本无关，假设初值选择得较好，一般迭代几次就能收敛；2

18、对于有些病态条件的问题，也能利用该方法求解；3应用了稀疏矩阵技巧，所需计算机存适中。牛顿-拉夫逊法虽是一种广泛使用的方法，但也存在以下缺点。1编程比拟复杂，且收敛速度的快慢和迭代次数与初始值的好坏有很大的关系，如果初始值选择不适宜有可能永远不收敛；2因非对称的雅可比矩阵不是固定的，每次迭代都需要重新计算，大量的求导运算，计算量很大，降低了计算速度。3.2 PQ分解算法针对牛顿-拉夫逊法计算速度方面存在的缺乏和电力系统实现在线控制的要求，在改良牛顿-拉夫逊法的根底上，提出了快速解耦算法。快速解耦算法派生于牛顿-拉夫逊法的极坐标形式，又称为PQ 分解法。其根本思想是：把节点功率表示为电压向量的极坐

19、标方程式，抓住主要矛盾，把有功功率误差作为修正电压向量角度的依据，把无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开进展。它密切地结合了电力系统的固有特点，无论是存占用量还是计算速度方面都比牛顿-拉夫逊法有了较大的改良。当节点功率方程式采取极坐标表达式时，修正方程式为： 3.24这相当于把2(n-1)阶的线性方程组变成了两个n-1阶的线性方程组，将P和Q分开来进展迭代计算，因而大减少了计算工作量。但是H，L在迭代过程中仍然在不断的变化，而且又都是不对称的矩阵。对牛顿法的进一步简化，即把式3.23中的系数矩阵简化为在迭代过程中不变的对称矩阵。在一般情况下，线路两端电压的相角是不大

20、的，因此，我们可认为 3.25此外，与系统各节点无功功率相应的导纳远远小于该节点自导纳的虚部，即 或 3.26 考虑到以上关系，式3.23的系数矩阵中的各元素可表示为 3.27 3.28由于高压电力系统中有功功率潮流主要与各节点电压向量的角度有关，无功功率潮流则主要受各节点电压幅值的影响。因此，将有功功率和无功功率分解开来迭代，可简化为： 3.29将上式中的系数矩阵简化为在迭代过程中不变的对称矩阵。进一步可以把其表示为以下矩阵的乘积： (3.30 3.31将式 (3.29和3.30代入式 3.23中，得到 3.32用和分别左乘以上两式，便得 3.33由此可以将修正方程式变为： 3.34 3.3

21、5 在这两个修正方程式中系数矩阵元素就是系统导纳矩阵的虚部，因而系数矩阵是对称矩阵，且在迭代过程中保持不变，这就大大减少了计算工作量。用极坐标表示的节点功率增量为: 3.36式3.33、3.34、3.35构成了PQ分解法迭代过程的根本方程式。P-Q分解法迭代步骤为：给定各节点电压向量的电压初值、。计算各节点有功功率误差，并求出。解修正方程式，并进而计算各节点电压向量角度的修正量。修正各节点电压向量角度计算各节点无功功率误差，并求出解修正方程式，求出各节点电压幅值的修正量。修正各节点电压幅值。返回进展迭代，知道各节点功率误差及都满足收敛条件。流程图如下所示：输入节点和支路数据数据输入节点和支路数

22、据数据建立节点导纳矩阵形成矩阵和，设PQ节点电压初值，各节点电压相角初值置迭代次数计算相角修正量，求得用公式计算不平衡功率计算计算电压修正量，求得计算平衡节点功率置输出数据是否是是否否置用公式计算不平衡功率计算置是否置图3.3 PQ分解法流程图PQ分解法修正方程式的一个显著的特点就是：用一个n-1阶和一个m阶的系数矩阵和替代原有的n+m-2阶系数矩阵J，提高了计算速度，降低了对存储容量的要求，使得求逆等运算量大幅度的减少。3.3 本章小结本章主要介绍了电力系统潮流计算时两种常用方法。首先，着重介绍了牛顿-拉夫逊法，推到了在极坐标情况下的迭代的方法和公式。该方法具有很好的收敛性，但要求有适宜的初

23、值。其次，介绍了P-Q分解法并推导了相关公式。P-Q分解法是极坐标形式的牛顿-拉夫逊法的简化算法，由于简化只涉及修正方程的系数矩阵，并未改变节点功率平衡方程和收敛判据，不会降低计算结果的精度，但极大提高了计算速度。4 算例4.1系统模型电力系统模型如图4.1所示，系统的参数见图。此系统中共有1个平衡节点、1个PV节点及3个PQ节点。图4.1 电力系统模型4.2 结果分析按照上图中所示的系统数据，得到如下面所示。系统的输入要求以及对于给定的5节点系统的原始数据的输入，如下所示：图4.2 输入说明以及原始数据输入通过程序计算得到的节点导纳矩阵如下：图4.3 节点导纳矩阵通过算法得到牛拉法的计算结果

24、为：图4.4 牛拉法计算结果通过算法得到PQ分解法的计算结果为：图4.5 PQ法计算结果从上面的运行结果可以看出，采用P-Q分解法与牛顿-拉夫逊法进展潮流计算，得到的计算结果是一致的，这充分说明程序的正确性和可行性。P-Q分解法迭代了的10次，就到达了的计算精度，而牛顿-拉夫逊法仅迭代了5次，就到达了精度要求，可见牛顿-拉夫逊法收敛速度快。但是，虽然P-Q分解法迭代次数较多，但是每次迭代的计算量却很小，因此P-Q分解法的计算速度相比牛顿-拉夫逊法有显著的提高，特别是对于大的电力系统而言，可以的显著的表现出来。4.3 本章小结本章介绍了潮流计算的两种经典算法的的编程过程，通过一个5节点的电力系统

25、验证了程序的正确性，并且对两种算法得到的结果进展了分析比拟，两种算法均能满足准确度的要求，只是迭代次数有所不同。当然，本程序还有一些有待改良的地方。本文所用的算例只是一个5节点的系统，对于实际的大型电力系统而言，节点数很多，这就要考虑迭代时间的问题，还有节点输入的编号的问题，这是本程序需要改良的问题。结论潮流计算是电力系统的一项根本计算，它根据给定的运行条件及系统接线情况确定系统的运行状态、母线的电压，各元件中流过的功率，系统的功率损耗等，进而比拟运行方式的合理性、可靠性和经济性。本文首先简单介绍了潮流计算的根本原理和根本功率方程，为后文的潮流计算方程式的推导打下根底。本文接着详细介绍了牛顿-

26、拉夫逊迭代法的思想以及在潮流计算中的应用，再通过简化推导出另一种目前最常用的潮流算法P-Q分解法。介绍了如何从定雅克比潮流迭代公式推导快速分解潮流算法的迭代公式。实际上，有功功率和无功功率的耦合作用已隐含在快速分解法的特殊迭代过程中，因此它具有良好的收敛性。本文在最后通过用MATLAB编写的潮流程序对5节点系统的计算，分别采用了牛顿-拉夫逊迭代法和P-Q分解法。两个迭代所产生的结果一致，证明了该潮流程序的正确性和可行性。-. z.参考文献1 何仰赞等.电力系统分析M.第3版.华中科技大学，2002.2 伯明,寿,严正.高等电力网络分析M.第2版.清华大学,2000.3 祯祥.电力系统分析M .

27、：大学，2009.4王锡凡，方万良，杜正春.现代电力系统分析M .：科学，2003.5 郗忠梅，有安，法起，等. 基于Matlab的电力系统潮流计算J . 农业大学学报，2010，41(2):291-294.6威，志炜，朝纲，马春生.一种基于牛顿拉夫逊的潮流就算方法J.学院学报. 2006.3,25(2):27-30.7徐劲松，宁玉琳，永峰. 基于Matlab的电力系统PQ分解法潮流计算研究J . 电气传动自动化，2011，33(2):10-18.8 祝书萍.电力系统分析课程实际设计与综合实验M. 第1版.中国电力，20009 帆.电力系统潮流计算程序设计J.冶金，2007，1062：42-4

28、4.附 录程序清单如下：clearclcdisplay( 输入说明： );display(B1是支路参数矩阵);display(第一列和第二列是节点编号(节点编号由小到大编写),对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗置于低压侧处理);display(第三列为支路的串列阻抗参数,第四列为支路的对地导纳,第五列为含变压器支路的变压器的变比,第六列为变压器是否含有变压器的参数，其中1为含有变压器，0为不含有变压器);display(B2为节点参数矩阵);display(第一列为节点注入发电功率参数；第二列为节点负荷功率；第三列为节点电压参数；第六列为

29、节点类型参数，其中1为平衡节点，2为PQ节点，3为PV节点参数理);display(*为节点号和对地参数矩阵,第一列为节点编号，第二列为节点对地参数);disp();display( 下面开场输入，进入程序 );n=input(请输入节点数:n=);nl=input(请输入支路数:nl=);isb=input(请输入平衡母线节点号:isb=);pr=input(请输入误差精度:pr=);B1=input(请输入由支路参数形成的矩阵:B1=);%变压器侧为1，否则为0B2=input(请输入各节点参数形成的矩阵:B2=);*=input(请输入由节点号及其对地阻抗形成的矩阵:*=);Y=zero

30、s(n);U=zeros(1,n);cta=zeros(1,n);V=zeros(1,n);O=zeros(1,n);S1=zeros(nl);%求节点导纳矩阵for i=1:n %判断是否有接地导纳if *(i,2)=0;p=*(i,1);Y(p,p)=*(i,2);endendfor i=1:nlif B1(i,6)=0 %不含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,

31、4);else %含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，因此需要调换次序。p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)2*B1(i,3);endenddisp();disp(节点导纳矩阵为:);disp(Y);G=real(Y);B=imag(Y);%别离出导纳矩阵的实部和虚部display( 本程序的功能是用牛顿拉夫逊法进展潮流计算 );for i=1:ncta(i)=angle(

32、B2(i,3);U(i)=abs(B2(i,3);endfor i=1:nS(i)=B2(i,1)-B2(i,2);B(i,i)=B(i,i)+B2(i,5);endP=real(S);Q=imag(S);ICT1=0;IT2=1;while IT2=0IT2=0;t1=1;t2=1;%t1记录PQ、PV节点有功变化数；t2记录PQ节点无功变化数。for i=1:nif i=isbC(i)=0;D(i)=0;for j1=1:n%每个节点的有功和无功功率方程C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-c

33、ta(j1);D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1);endDP(t1)=P(i)-C(i);%PQ、PV节点的有功功率变化量t1=t1+1;if B2(i,6)=2DQ(t2)=Q(i)-D(i);%PQ节点无功功率变化量t2=t2+1;end endendt1=t1-1;t2=t2-1;%判断是否满足精度要求DPQ=DP;DQ; %求DP,DQfor i=1:t1+t2if abs(DPQ(i)pr%不满足精度要求，继续迭代，直到满足IT2=IT2+1;endendH=zeros

34、(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2);%雅克比矩阵初始化for i=1:t1 for j1=1:t1%H矩阵n-1*n-1阶if j1=isb&j1=i%非对角线元素H(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1);elseif j1=isb&j1=i%对角线元素H(i,j1)=U(i)2*B(i,j1)+D(i);endendendfor i=1:t1%N矩阵n-1*m阶for j1=1:t2if j1=isb&j1=i%非

35、对角线元素N(i,j1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1);elseif j1=isb&j1=i%对角线元素N(i,j1)=0-U(i)2*G(i,j1)-C(i);endendendfor i=1:t2%K矩阵m*n-1阶for j1=1:t1if j1=isb&j1=i%非对角线元素K(i,j1)= U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1);elseif j1=isb&j1=i%对角线元素K(i,j1

36、)=U(i)2*G(i,j1)-C(i);endendendfor i=1:t2%L矩阵m*m阶for j1=1:t2if j1=isb&j1=i%非对角线元素L(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1);elseif j1=isb&j1=i%对角线元素L(i,j1)=U(i)2*B(i,j1)-D(i);endendendJ=H,N;K,L;%求雅可比矩阵modify=-JDPQ;Dcta=modify(1:t1,:);%相角变化修正t3=U(:,1:t2);DU=diag(t3,0)*

37、modify(t1+1:t1+t2,:);%电压变化修正t4=1;for i=1:t1if B2(i,6)=1%计算非平衡节点的电压的相角cta(1,i)=cta(1,i)+Dcta(t4,1);t4=t4+1;endendt5=1;for i=1:t2if B2(i,6)=2%计算PQ节点的电压幅值U(1,i)=U(1,i)+DU(t5,1);t5=t5+1;endendICT1=ICT1+1;endfor i=1:n%计算所有节点的电压幅值UU(i)=U(i)*cos(cta(i)+1i*U(i)*sin(cta(i);enddisp();disp(各节点电压U为节点从小到大排列:);di

38、sp(UU);disp();disp(各节点电压相角为节点从小到大排列:);disp(180*angle(UU);disp();disp(按公式计算全部线路功率，结果如下:);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)=-(UU(p)*(conj(UU(p)*conj(B1(i,4)./2)+(conj(UU(p)*B1(i,5)-conj(UU(q)*conj(1./(B1(i,3)*B1(i,5);%各条支路首端功率SiSj(q,p)=-(UU(q)*(conj(UU(q)*co

39、nj(B1(i,4)./2)+(conj(UU(q)./B1(i,5)-conj(UU(p)*conj(1./(B1(i,3)*B1(i,5);%各条支路末端功率Sjf1=p,q,Si(p,q);f=q,p,Sj(q,p);disp(f1);disp(f);enddisp();disp(按公式计算各节点注入功率结果如下:);for p=1:nc(p)=0;for q=1:nc(p)=c(p)+conj(Y(p,q)*conj(UU(q);%计算所有节点的电流ends(p)=-UU(p)*c(p); %计算所有节点节点注入潮流disp( s(p);endSp=0;for i=1:nSp=-(Sp

40、+UU(isb)*conj(Y(isb,i)*conj(UU(i);enddisp();disp(平衡节点的功率:);disp(Sp);disp();disp(迭代次数:);disp(ICT1);disp();display( 本程序的功能是用PQ进展潮流计算 );na=input(请输入PQ节点数na=);for i=1:n%进展电压的实虚局部离e(i)=real(B2(i,3);f(i)=imag(B2(i,3);V(i)=B2(i,4);endfor i=1:nif B2(i,6)=2%求PQ节点的电压V(i)=sqrt(e(i)2+f(i)2);O(i)=atan(f(i)./e(i)

41、;endend%下面进展B和B的计算求取for i=2:n%对B进展三角分解if i=n B(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:n%下三角的处理for j1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendendp=0;q=0;for i=1:n%对PQ节点的处理if B2(i,6)=2p=p+1;k=0;for j1=1:nif B2(j1,6)=2 k=k+1; A(p,k)=B(i

42、,j1); end%将PQ节点导纳矩阵的虚步提出来，按照上三角存放endendendfor i=1:na %对每个虚步进展处理，即是B矩阵if i=na A(i,i)=1./A(i,i);else k=i+1;for j1=k:na%上三角的处理A(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;%进展PQ分解法的迭代while

43、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)+B(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=pr%不满足继续迭代ICT2=ICT2+1;endendend%Np(K)=ICT2;if ICT2=0%进展相角的迭代for i=2:n-1DP(i)=B(i,i)*DP(i);%计算有功功率的变化量if i=nIC1=i+1;for k=IC1:n-1DP(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:n-1O(i)=O(i)-DP