毕业论文高压输电网潮流的计算机算法程序设计

1、目录第一章课程设计概述11.1 设计目的11.2 设计要求11.3 设计题目11.4 设计内容21.5设计时间2第二章（针对该课程设计题目的）设计思路32.1 潮流计算题目32.2 对课题的分析及求解思路3第三章电力系统潮流计算概述53.1电力系统叙述53.2潮流计算简介53.3潮流计算的意义及其发展6第四章导纳矩阵的原理及计算方法84.1自导纳和互导纳的确定方法84.2潮流计算的基本方程124.3电力系统节点分类154.4潮流计算的约束条件16第五章牛顿拉夫逊法及程序代码185.1牛顿拉夫逊法18牛顿拉夫逊法概要18牛顿拉夫逊法的求解过程19程序框图如下：215.2 潮流计算程序代码21第六

2、章运行结果分析306.1 正常运行结果分析306.1.1 MATLAB运行结果306.1.2 PowerWorld运行结果356. 1. 3 综合分析366.2 非正常运行结果分析376.2.1 若某发电站发电量减半时对系统供电影响376.2.2 若某发电站发电量为零时对系统供电影响396.2.3 若某负荷切掉时对系统供电影响406.2.4 若系统环网解开对系统的影响42第七章总结44参考文献45第一章 课程设计概述1.1 设计目的1. 掌握电力系统潮流计算的基本原理和电力系统运行方式的变化2. 让学生结合课程教学内容了解实际电力系统分析中的某些最基本的计算、分析过程、常规典型计算方法和结果分

3、析等方面知识，加深电力系统分析课程的理解和掌握3. 掌握并能熟练运用一门计算机语言（MATLAB语言或C语言或C+语言）；4. 采用计算机语言对潮流计算进行计算机编程计算。5. 培养和提高学生对电力系统分析与设计能力。 1.2 设计要求1. 程序源代码；2. 选定算例的输入，输出文件；3. 程序说明；4. 选定算例的程序计算过程；5. 选定算例的手算过程（至少迭代2次）（可选）。学生开展课程设计可以7人选择同一题目，但每个人必须独立完成所有给定任务，学生必须按计划开展课程设计, 每一阶段完成情况必须接受教师检查。编程实现部分，学生可以根据自己掌握的编程工具情况选择，如MATLAB、C、VC等，

4、 最后必须提供程序代码，课程设计必须完成对系统的分析与设计任务,编写相应的分析与设计报告。 1.3 设计题目1.高压输电网潮流的计算机算法程序设计（PQ分解法、牛顿-拉夫逊法）2.中压配电网潮流的计算机算法程序设计（前推后代法、同伦延拓法等）3.电力系统短路故障的计算机算法程序设计（要求不限）本设计选择为高压输电网潮流的计算机算法程序设计（PQ分解法、牛顿-拉夫逊法）1.4 设计内容1.根据电力系统网络推导电力网络数学模型，写出节点导纳矩阵；2.赋予各节点电压变量（直角坐标系形式）初值后，求解不平衡量；3.形成雅可比矩阵；4.求解修正量后，重新修改初值，从2开始重新循环计算；5.求解的电压变量

5、达到所要求的精度时，再计算各支路功率分布、功率损耗和平衡节点功率; 6.上机编程调试；连调；7.计算分析给定系统潮流分析。8.准备计算机演示答辩，书写该课程设计说明书（必须计算机打印）。1.5设计时间2013年春季17周至18周第二章（针对该课程设计题目的）设计思路2.1 潮流计算题目G 1 2 3 -10Mvar 180MW 100MvarG 5 4 40MW 30Mvar50MW 30Mvar如图所示，变压器参数为：Zt1=3+j110, Zt2=0.8+j23,变比分别为231/110kV,231/121kV;220kV线路参数为：Z12=5.9+j31.5;110kV线路：Z45=65

6、+j100; 34段，Z=65+j100，所有阻抗都已经归算至220kV侧。2.2 对课题的分析及求解思路等值电路的计算电压是衡量电力系统电能质量的标准之一。电压过高或过低，都将对人身及其用电设备产生重大的影响。保证用户的电压接近额定值是电力系统调度的基本任务之一。当系统的电压偏离允许值时，电力系统必须应用电压调节技术调节系统电压的大小，使其维持在允许值范围内。本文经过手算形成了等值电路图，并编写好了程序得出节点电压标幺值，使其满足所要求的调整范围。我们首先对给定的程序输入部分作了简要的分析，程序开始需要我们确定输入节点数、支路数、平衡母线号、支路参数矩阵、节点参数矩阵。（1）为了保证整个系统

7、潮流计算的完整性，我们把凡具有母线及发电机处均选作节点，这样，可以确定有5个节点，节点号见等值电路图。（2）确定完节点及编号后，各条支路也相应确定了，网络中总计有5条支路，我们对各支路参数进行了计算。根据所给实际电路图和题中的已知条件，有以下公式计算各输电线路的阻抗和对地支路电容的标幺值和变压器的阻抗标幺值。选择电压基准值为和功率基准值所以数据的输入该系统中，母线1为平衡节点，保持定值电压Ù=220+j0kV,5节点为PV节点，注入有功功率为40MW，无功功率为 30Mvar，其余节点均为PQ节点，注入功率分别为2=-j10MVA，3=-180-j100MVA,4=-50-j30MV

8、A.首先计算节点数据矩阵如下，%(bus#)(volt) (ang) (p) (q) (bus type)bus= 1 1.00 0.00 0.00 0.00 3; 2 1.00 0.00 0.00 -0.10 1; 3 1.00 0.00 -1.80 -1.00 1; 4 1.00 0.00 -0.50 -0.30 1; 5 1.00 0.00 0.40 0.30 2;然后计算线路数据矩阵如下， % b#1 b#2 (R) (X) (G) (B) (K)line= 1 2 5.9 31.5 0.0 0.0 0; 1 5 3 110 0.0 0.0 1.0; 2 3 0.8 23 0.0 0.

9、0 1.0; 3 4 65 100 0.0 0.0 0; 4 5 65 100 0.0 0.0 0;这就形成了计算潮流的原始数据，把它代入程序中就可以计算潮流分布。第三章 电力系统潮流计算概述3.1电力系统叙述 电力工业发展初期，电能是直接在用户附近的发电站（或称发电厂）中生产的，各发电站孤立运行。随着工农业生产和城市的发展，电能的需要量迅速增加，而热能资源和水能资源丰富的地区又往往远离用电比较集中的城市和工矿区，为了解决这个矛盾，就需要在动力资源丰富的地区建立大型发电站，然后将电能远距离输送给电力用户。同时，为了提高供电的可靠性以及资源利用的综合经济性，又把许多分散的各种形式的发电站，通过送

10、电线路和变电所联系起来。这种由发电机、升压和降压变电所，送电线路以及用电设备有机连接起来的整体，即称为电力系统。现代电力系统提出了“灵活交流输电和新型直流输电”的概念。灵活交流输电技术是指运用固态电子器件与现代自动控制技术对交流电网的电压、相位角、阻抗、功率以及电路的通断进行实时闭环控制，从而提高高压输电线路的诉讼能力和电力系统的稳态水平。新型直流输电技术是指应用现电力电子技术的最新成果，改善和简化变流站的造价等。运营方式管理中，潮流是确定电网运行方式的基本出发点：在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了电网在预想操作预想下的电网的潮流分布以及校验运行

11、的可靠性。在电力系统调度运行的多个领域都涉及到电网潮流计算。潮流是确定电力网咯运行状态的基本因素，潮流问题是研究电力系统稳态问题的基础和前提。3.2潮流计算简介电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压。各元件中流过的功率，系统的功率损耗等等。在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量的分析比较供电方案或运行方式的合理性。可靠性和经济性。此外，电力系统的潮流计算也是计算机系统动态稳定和静态稳定的基础，所以潮流计算是研究电力系统的一种和重要和基础的计算。电力系统潮流计算也

12、分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。利用电子数字计算机进行潮流计算从50年代中期就已经开始了。在这20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的，对潮流计算的要求可以归纳为以下几点：（1） 计算方法的可靠性或收敛性；（2） 对计算机内存量的要求；（3） 计算速度；（4） 计算的方便性和灵活性。3.3潮流计算的意义及其发展电力系统潮流计算是电力系统分析中的一种最基本的计算，是对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标是求取电力系统在给定运行状态

13、的计算，即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有的电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和稳态分析都是以潮流计算为基础。潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算的模型和方法有直接影响。实际电力系统的潮流技术那主要采用牛顿拉夫逊法。运行方式管理中，潮流是确定电网运行方式的基本出发点；在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了多个在预想操作情况下电网的潮流分布以及校验运行可靠性。在电力系统调度运行的多个领域问题是研究电力系

14、统稳态问题的基础和前提。在用数字见算计算机解电力系统潮流问题的开始阶段，普遍采取以节点导纳矩阵为基础的逐次代入法。这个方法的原理比较简单，要求的数字计算机内存量比较差下，适应50年代电子计算机制造水平和当时电力系统理论水平，但它的收敛性较差，当系统规模变大时，迭代次数急剧上升，在计算中往往出现迭代不收敛的情况。这就迫使电力系统的计算人员转向以阻抗矩阵为基础的逐次代入法。阻抗法改善了系统潮流计算问题的收敛性，解决了导纳无法求解的一些系统的潮流计算，在60年代获得了广泛的应用，阻抗法德主要缺点是占用计算机内存大，每次迭代的计算量大。当系统不断扩大时，这些缺点就更加突出，为了克服这些缺点，60年代中

15、期发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需要存储各个地区系统的阻抗矩阵及它们之间联络的阻抗，这样不仅大幅度的节省了内存容量，同时也提高了计算速度。克服阻抗法缺点是另一个途径是采用牛顿-拉夫逊法。这是数学中解决非线性方程式的典型方法，有较好的收敛性。在解决电力系统潮流计算问题时，是以导纳矩阵为基础的，因此，只要我们能在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿法潮流程序的效率。自从60年代中期，牛顿法中利用了最佳顺序消去法以后，牛顿法在收敛性。内存要求。速度方面都超过了阻抗法，成为了60年代末期以后广泛采用的优秀方法。 本

16、次设计是利用牛顿拉夫逊法进行求解，用MATLAB软件编程，可以求解系统潮流分布根据题目的不同要求对参数进行调整，通过调节变压器变比和发电厂的电压，求解出合理的潮流分布，最后用powerworld进行潮流分析，将两者进行比较。第四章 导纳矩阵的原理及计算方法4.1自导纳和互导纳的确定方法 电力网络的节点电压方程： 为节点注入电流列向量，注入电流有正有负，注入网络的电流为正，流出网络的电流为负。根据这一规定，电源节点的注入电流为正，负荷节点为负。既无电源又无负荷的联络节点为零，带有地方负荷的电源节点为二者代数之和。为节点电压列向量，由于节点电压是对称于参考节点而言的，因而需先选定参考节点。在电力系

17、统中一般以地为参考节点。如整个网络无接地支路，则需要选定某一节点为参考。设网络中节点数为（不含参考节点），则，均为n*n列向量。为n*n阶节点导纳矩阵。节电导纳矩阵的节点电压方程： 展开为： ：是一个n*n阶节点导纳矩阵，其阶数就等于网络中除参考节点外的节点数。 节点导纳矩阵的对角元素 (i=1，2，n)成为自导纳。自导纳数值上就等于在i节点施加单位电压，其他节点全部接地时，经节点i注入网络的电流，因此，它可以定义为：节点i的自导纳数值上就等于与节点直接连接的所有支路导纳的总和。节点导纳矩阵的非对角元素 (j=1，2，n;i=1，2，。，n;j=i)称互导纳，由此可得互导纳数值上就等于在节点i

18、施加单位电压，其他节点全部接地时，经节点j注入网络的电流，因此可定义为：节点j，i之间的互导纳数值上就等于连接节点j，i支路到导纳的负值。显然，恒等于。互导纳的这些性质决定了节点导纳矩阵是一个对称稀疏矩阵。而且，由于每个节点所连接的支路数总有一个限度，随着网络中节点数的增加非零元素相对愈来愈少，节点导纳矩阵的稀疏度，即零元素数与总元素的比值就愈来愈高。节点导纳矩阵的性质及意义节点导纳矩阵的性质：（1）为对称矩阵，=。如网络中含有源元件，如移相变压器，则对称性不再成立。（2）对无接地支路的节点，其所在行列的元素之和均为零，即 。对于有接地支路的节点，其所在行列的元素之和等于该点接地支路的导纳。利

19、用这一性质，可以检验所形成节点导纳矩阵的正确性。（3）具有强对角性：对角元素的值不小于同一行或同一列中任一元素（4）为稀疏矩阵，因节点i ，j 之间无支路直接相连时=0，这种情况在实际电力系统中非常普遍。矩阵的稀疏性用稀疏度表示，其定义为矩阵中的零元素与全部元素之比，即 ， 式中Z 为中的零元素。S 随节点数n 的增加而增加：n=50，S可达92%；n=100，S 可达90%；n=500，S可达99%，充分利用节点导纳矩阵的稀疏性可节省计算机内存，加快计算速度，这种技巧称为稀疏技术。节点导纳矩阵的意义：是n*n阶方阵，其对角元素 (i=1，2，-n)称为自导纳，非对角元素(i，j=1，2，n，

20、)称为互导纳。将节点电压方程展开为可见， 表明，自导纳在数值上等于仅在节点i施加单位电压而其余节点电压均为零（即其余节点全部接地）时，经节点i注入网络的电流。其显然等于与节点i直接相连的所有支路的导纳之和。同时可见。表明，互导纳在数值上等于仅在节点j施加单位电压而其余节点电压均为零时，经节点i注入网络的电流，其显然等于()即=。为支路的导纳，负号表示该电流流出网络。如节点ij之间无支路直接相连，则该电流为0，从而=0。注意字母几种不写法的不同意义：粗体黑字表示导纳矩阵，大写字母代矩阵中的第i行第j列元素，即节点i和节点j之间的互导纳。小写字母i，j支路的导纳等于支路阻抗的倒数数，。根据定义直接

21、求取节点导纳矩阵时，注意以下几点：1)。节点导纳矩阵是方阵，其阶数就等于网络中除去参考节点外的节点数。参考节点一般取大地，编号为零。2)。节点导纳矩阵是稀疏矩阵，其各行非零非对角元素就等于与该行相对应节点所连接的不接地支路数。3)。节点导纳矩阵的对角元素就等于各该节点所连接导纳的总和。因此，与没有接地支路的节点对应的行或列中，对角元素为非对角元素之和的负值。4)。节点导纳矩阵的非对角元素等于连接节点i，j支路导纳的负值。因此，一般情况下，节点导纳矩阵的对角元素往往大于非对角元素的负值。5)。节点导纳矩阵一般是对称矩阵，这是网络的互易特性所决定的。从而，一般只要求求取这个矩阵的上三角或下三角部分

22、。非标准变比变压器等值电路变压器型等值电路更便于计算机反复计算,更适宜于复杂网络的潮流计算.双绕组变压器可用阻抗与一个理想变压器串联的电路表示.理想变压器只是一个参数,那就是变比。现在变压器阻抗按实际变比归算到低压侧为例,推导出变压器型等值电路.a 双绕组变压器原理图b 变压器阻抗归算到低压侧等值模型流入和流出理想变压器的功率相等式中,是理想变压器的变比,和 分别为变压器高,低绕组的实际电压.从图b直接可得：从而可得: 式中,又因节点电流方程应具有如下形式:将式（1-8）与（1-9）比较，得：因此可得各支路导纳为: 由此可得用导纳表示的变压器型等值电路:图 c4.2潮流计算的基本方程在潮流问题

23、中，任何复杂的电力系统都可以归纳为以下元件（参数）组成。（1）发电机（注入电流或功率）（2）负荷（注入负的电流或功率）（3）输电线支路（电阻，电抗）（4）变压器支路（电阻，电抗，变比）（5）母线上的对地支路（阻抗和导纳）（6）线路上的对地支路（一般为线路充电点容导纳）集中了以上各类型的元件的简单网络如图 (a) 潮流计算用的电网结构图(b) 潮流计算等值网络采用导纳矩阵时，节点注入电流和节点电压构成以下线性方程组其中 可展开如下形式 由于实际电网中测量的节点注入量一般不是电流而是功率，因此必须将式中的注入电流用节点注入功率来表示。节点功率与节点电流之间的关系为 式中，因此用导纳矩阵时，PQ节点

24、可以表示为把这个关系代入式中 ，得就是电力系统潮流计算的数学模型-潮流方程。它具有如下特点：（1）它是一组代数方程，因而表征的是电力系统的稳定运行特性。（2）它是一组非线性方程，因而只能用迭代方法求其数值解。（3）由于方程中的电压和导纳既可以表为直角坐标，又可表为极坐标，因而潮流方程有多种表达形式-极坐标形式，直角坐标形式和混合坐标形式。a。取 ，得到潮流方程的极坐标形式：b。 取 ， ，得到潮流方程的直角坐标形式：c。取，得到潮流方程的极坐标形式：不同坐标形式的潮流方程适用于不同的迭代解法。例如：利用牛顿-拉夫逊迭代法求解，以直角坐标和极坐标形式的潮流方程为方便；而P-Q解耦法是在极坐标形式

25、的基础上发展而成，本次设计采用牛顿-拉夫逊迭代法求解且采用极坐标形式。（4）它是一组n个复数方程，因而实数方程数为2n个但方程中共含4n个变量：P，Q，U和，i=1，2，n，故必须先指定2n个变量才能求解。4.3电力系统节点分类用一般的电路理论求解网络方程，目的是给出电压源(或电流源)研究网络内的电流(或电压)分布，作为基础的方程式，一般用线性代数方程式表示。然而在电力系统中，给出发电机或负荷连接母线上电压或电流(都是向量)的情况是很少的，一般是给出发电机母线上发电机的有功功率(P)和母线电压的幅值(U)，给出负荷母线上负荷消耗的有功功率(P)和无功功率(Q)。主要目的是由这些已知量去求电力系

26、统内的各种电气量。所以，根据电力系统中各节点性质的不同，很自然地把节点分成三类： PQ节点对这一类点，事先给定的是节点功率(P，Q)，待求的未知量是节点电压向量(U，)，所以叫PQ节点。通常变电所母线都是PQ节点，当某些发电机的输出功率P。Q给定时，也作为PQ节点。PQ节点上的发电机称之为PQ机(或PQ给定型发电机)。在潮流计算中，系统大部分节点属于PQ节点。 PU节点这类节点给出的参数是该节点的有功功率P及电压幅值U，待求量为该节点的无功功率Q及电压向量的相角。这类节点在运行中往往要有一定可调节的无功电源。用以维持给定的电压值。通常选择有一定无功功率储备的发电机母线或者变电所有无功补偿设备的

27、母线做PU节点处理。PU节点上的发电机称为PU机(或PU给定型发电机) 平衡节点在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在计算中取该节点电压向量的方向作为参考轴，相当于给定该点电压向量的角度为零。也就是说，对平衡节点给定的运行参数是U和，因此有城为U节点，而待求量是该节点的P。Q，整个系统的功率平衡由这一节点承担。关于平衡节点的选择，一般选择系统中担任调频调压的某一发电厂(或发电机)，有时也可能按其他原则选择，例如，为提高计算的收敛性。可以选择出线数多或者靠近电网中心的发电厂母线作平衡节点。以上三类节点4个运行参数P。Q。U。中，已知量都是两个，待求量也是两个，只是类型不

28、同而已。4.4潮流计算的约束条件电力系统运行必须满足一定技术和经济上的要求。这些要求够成了潮流问题中某些变量的约束条件，常用的约束条件如下：1. 节点电压应满足2.从保证电能质量和供电安全的要求来看，电力系统的所有电气设备都必须运行在额定电压附近。PU节点电压幅值必须按上述条件给定。因此，这一约束条件对PQ节点而言。3. 节点的有功功率和无功功率应满足4.PQ节点的有功功率和无功功率，以及PU节点的有功功率，在给定是就必须满足上述条件，因此，对平衡节点的P和Q以及PU节点的Q应按上述条件进行检验。5. 节点之间电压的相位差应满足 为了保证系统运行的稳定性，要求某些输电线路两端的电压相位不超过一

29、定的数值。这一约束的主要意义就在于此。 因此，潮流计算可以归结为求解一组非线性方程组，并使其解答满足一定的约束条件。常用的方法是迭代法和牛顿法，在计算过程中，或得出结果之后用约束条件进行检验。如果不能满足要求，则应修改某些变量的给定值，甚至修改系统的运行方式，重新进行计算。第五章牛顿拉夫逊法及程序代码5.1牛顿拉夫逊法牛顿拉夫逊法概要首先对一般的牛顿拉夫逊法作一简单的说明。已知一个变量X函数为：到此方程时，由适当的近似值出发，根据：反复进行计算，当满足适当的收敛条件就是上面方程的根。这样的方法就是所谓的牛顿拉夫逊法。这一方法还可以做下面的解释，设第次迭代得到的解语真值之差，即的误差为时，则：把

30、在附近对用泰勒级数展开上式省略去以后部分的误差可以近似由上式计算出来。比较两式，可以看出牛顿拉夫逊法的休整量和的误差的一次项相等。用同样的方法考虑，给出个变量的个方程：对其近似解得修正量可以通过解下边的方程来确定：式中等号右边的矩阵都是对于的值。这一矩阵称为雅可比（JACOBI）矩阵。按上述得到的修正向量后，得到如下关系这比更接近真实值。这一步在收敛到希望的值以前重复进行，一般要反复计算满足为预先规定的小正数，是第n次迭代的近似值。牛顿拉夫逊法的求解过程（1）用牛顿法计算潮流时，有以下的步骤：给这各节点电压初始值；将以上电压初始值代入公式，求修正方程的常数项向量 ；将电压初始值在带入上述公式，

31、求出修正方程中系数矩阵的各元素。解修正方程式；修正各节点电压，；将，在带入方程式，求出；检验是否收敛，即如果收敛，迭代到此结束，进一步计算各线路潮流和平衡节点功率，并打印输出结果。如果不收敛，转回进行下次迭代计算，直到收敛为止。程序框图如下：5.2 潮流计算程序代码主程序：clc;clear allformat long;global nSW nPQ nPV; global nb; global nl; global bus; global line; global Y; global nodenum; global lPQ; global myf; openfile; change; ybu

32、s; NR; PQ; flow; ret; output; 打开bus line的数据文件的子程序：function openfileglobal nb nl;global line;global bus;global myf;dfile,pathname=uigetfile('*.m','Select Data File'); if pathname = 0 error(' you must select a valid data file')else lfile =length(dfile);% strip off .m eval(dfile

33、(1:lfile-2); endnl,ml=size(line); nb,mb=size(bus); outfile=''for I=1:length(dfile)if dfile(I)='.'break;end outfile=strcat(outfile,dfile(I); endoutfile=strcat(outfile,'.txt'); myf=fopen(outfile,'w'); 节点重新编号：function changeglobal nb;global nl;global nPQ;global bus;globa

34、l line;global nodenum;nSW = 0; nPV = 0; nPQ = 0; for I = 1:nb, type= bus(I,6);if type = 3, nSW = nSW + 1; SW(nSW,:)=bus(I,:);elseif type = 2, nPV = nPV +1; PV(nPV,:)=bus(I,:);else nPQ = nPQ + 1; PQ(nPQ,:)=bus(I,:);endendbus=PQ;PV;SW; newbus=1:nb'nodenum=newbus bus(:,1); bus(:,1)=newbus; for I=1:

35、nl for J=1:2for k=1:nbif line(I,J)=nodenum(k,2) line(I,J)=nodenum(k,1);breakendendendend建立节点导纳矩阵：functionY = ybus(bus,line)global nl;global nb;global bus;global line;global Y;global myf;Y=zeros(nb,nb); for k=1:nlI=line(k,1); J=line(k,2); Zt=line(k,3)+j*line(k,4);if J=0 Yt=1/Zt; end%Ym=line(k,5)+j*li

36、ne(k,6); K=line(k,7); if (K=0)&(J=0) Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J);endif (K=0)&(J=0) Y(I,I)=Y(I,I)+Ym;endif K>0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/K/K; Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J);endif K<0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+K*K*Yt;

37、Y(I,J)=Y(I,J)+K*Yt; Y(J,I)=Y(I,J);endend牛拉法解方程组：function NRglobal nb;global nPQ;global bus;global line;global Y;global myf;max1=100; eps1=1.0e-10; eps2=1.0e-10;fprintf(myf, '-节点导纳矩阵Y-n'); for I=1:nbfor J=1:nb fprintf(myf, '%10f+j*(%10f) ', real(Y(I,J),imag(Y(I,J);end fprintf(myf, 

38、9;n');endfor i=1:max1 angl(:,1)=bus(1:nb-1,3); u(:,1)=bus(1:nPQ,2); x=angl;u; Jac=form_jac(bus,Y); del=dPQ(Y,bus); dx=Jacdel; fprintf(myf,'-第%d次迭代结果-n',i); fprintf(myf,'-第%d次迭代的雅比矩阵J-n',i);for I=1:nb+nPQ-1for J=1:nb+nPQ-1 fprintf(myf, '%10f ', Jac(I,J); end fprintf(myf, &

39、#39;n');end fprintf(myf,'-第%d次迭代的功率偏差dP和dQ-n',i);for I=1:nb+nPQ-1 fprintf(myf, '%10e ', del(I,1); fprintf(myf, 'n');end fprintf(myf,'-第%d次迭代的节点相角和电压的偏差dx-n',i);for I=1:nb+nPQ-1 fprintf(myf, '%10e ', dx(I,1); fprintf(myf, 'n');endfor I=nb:nb+nPQ-1 d

40、x(I,1)=dx(I,1)*x(I,1); end x=x-dx; fprintf(myf,'-第%d次迭代的节点相角delta和电压U-n',i); angl=x(1:nb-1,1); u=x(nb:nb+nPQ-1,1);for I=1:nb-1 bus(I,3)=angl(I,1); fprintf(myf, 'ang%d %10f n', I,angl(I,1);endfor I=1:nPQ % bus(I,2)=u(I,1); fprintf(myf, 'U%d %10f n', I,u(I,1);endif(max(abs(dx)&

41、lt;eps1)&(max(abs(dPQ(Y,bus)<eps2) breakendendif i=max1 error('超过最大迭代次数，不收敛停机！');end计算节点注入功率：global nb;global nPQ;global bus;global line;global Y;for I=nPQ+1:nb if bus(I,6)=3 sum=0; for J=1:nb ang=bus(I,3)-bus(J,3); A=real(Y(I,J)*cos(ang)+imag(Y(I,J)*sin(ang); sum=sum+bus(I,2)*bus(J,2

42、)*(A);end bus(I,4)=sum; end sum=0;for J=1:nb ang=bus(I,3)-bus(J,3); B=real(Y(I,J)*sin(ang)-imag(Y(I,J)*cos(ang); sum=sum+bus(I,2)*bus(J,2)*B;endbus(I,5)=sum; end计算线路功率及损耗：global nl;global lPQ;global bus;global line;for k=1:nlI=line(k,1); J=line(k,2); lPQ(k,1)=I; lPQ(k,2)=J;if J=0 Zt=line(k,3)+j*line

43、(k,4); Yt=1/Zt;end Ym=line(k,5)+j*line(k,6); K=line(k,7); Ui=bus(I,2)*(cos(bus(I,3)+j*sin(bus(I,3); if J=0 Uj=bus(J,2)*(cos(bus(J,3)+j*sin(bus(J,3); endif (K=0)&(J=0) Iij=Ui*(Yt+Ym)-Uj*Yt; Iji=Uj*(Yt+Ym)-Ui*Yt;endif (K=0)&(J=0) Iij=Ui*Ym; Iji=0;endif K>0 Iij=(Ui-Uj)*Yt/K+Ui*(Ym+Yt*(K-1)/K

44、); Iji=(Uj-Ui)*Yt/K+Uj*Yt*(1-K)/K2;endif K<0 K=-K; Iij=(Ui-Uj)*Yt*K+Ui*(Ym+Yt*(1-K); Iji=(Uj-Ui)*Yt*K+Uj*Yt*K*(K-1);end Sij=Ui*conj(Iij); Sji=Uj*conj(Iji);delS=Sij+Sji; lPQ(k,3:5)=Sij Sji delS; end恢复原节点编号：global nl;global nb;global bus;global line;global Y;global nodenum;global lPQ;for I=1:nb for

45、 J=1:nb if nodenum(J,2)=Ibreak;endendtem(I,:)=bus(nodenum(J,1),:); tem(I,1)=I;endbus=tem;for I=1:nl for J=1:2for k=1:nbif lPQ(I,J)=nodenum(k,1) lPQ(I,J)=nodenum(k,2);breakendendendend输出：function output global nb;global nl;global bus;global lPQ;global myf;fprintf(myf, 'n-牛顿拉夫逊法潮流计算结果-n');fprintf(myf, ' 节点计算结果：n'); fprintf(myf, '节点 节点电压 节点相角（角度） 节点注入功率n');for I=1:nb, fprintf(myf, ' %2d %10f %10f %10f+j%10f n', bus(I,1),bus(