电力系统分析PQ分解法潮流计算终稿

《电力系统分析PQ分解法潮流计算终稿》由会员分享，可在线阅读，更多相关《电力系统分析PQ分解法潮流计算终稿（40页珍藏版）》请在装配图网上搜索。

1、电力系统分析课程设计报告题目：复杂电力系统潮流分析程序设计所在学院 电气工程学院 专业班级 12级电气工程及其自动化一班学生姓名 陈剑秋 学生学号 2 指导教师 房大中 提交日期 2014年 12 月 16日 36 / 40文档可自由编辑打印摘要电力系统潮流计算是电力系统分析课程基本计算的核心部分之一。它既有本身的独立意义，又是电力系统规划设计，运行分析和理论研究的基础。电力系统潮流计算是对复杂电力系统正常和故障条件下稳态运行状态的计算，是电力系统计算分析中的一种最基本的计算。潮流计算是电力系统的各种计算的基础，同时它又是研究电力系统的一项重要分析功能，是进行故障计算，继电保护鉴定，安全分析的

2、工具。电力系统潮流计算是计算系统动态稳定和静态稳定的基础。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用电力系统潮流计算来定量的比较供电方案或运行方式的合理性、可靠性和经济性本课程设计采用P-Q分解法对三机九节点系统进行潮流计算。计算结果的得出是通过MATLAB软件编程实现。关键词：电力系统潮流计算，P-Q分解法，三机九节点系统AbstractPower flow calculation is the core part of the power system analysis course basic computing. It has both the independent

3、 meaning itself, and is the power system planning and design, the basic operation analysis and theoretical study.Power flow calculation is the calculation of the steady state operation of normal and fault conditions of complex power system, is one of the most basic calculation of power system calcul

4、ation and analysis of. Power flow calculation is the basis of all kinds of power system calculation, and it is also an important function of power system analysis study, is the fault calculation, relay protection and identification, safety analysis tool. Power flow calculation is the basis of calcul

5、ation system of dynamic and static stabilization. In the research of operation mode of power system planning and the existing power system, the rationality of power supply, reliability and economy of the program or run mode comparison requires the use of power flow calculation to quantitativeThe cur

6、riculum design using P-Q decomposition method for power flow calculation of the three machine nine bus system. It is concluded that the calculation result is achieved by MATLAB software programming.Keywords: power system power flow calculation, the P-Q decomposition method, three machine nine node s

7、ystem目录摘要IAbstractII一、 绪论11.1 本课题的目的和意义11.2 国内外发展现状11.2.1 高斯-赛德尔迭代法21.2.2 牛顿-拉夫逊法和P-Q分解法21.2.3 基于MATLAB 的电力系统潮流计算发展前景3二、 设计目的4三、 设计要求和设计指标4四、 设计内容44.1 选题内容44.2 基础资料54.3算法原理64.3.1 节点导纳矩阵形成的计算机方法64.3.2 电力网络的潮流方程94.3.3 电力网络极坐标形式的潮流方程104.3.4 极坐标牛顿潮流算法的修正方程雅和可比矩阵104.3.5 PQ分解法潮流计算11五、 程序设计135.1 主函数135.2 子

8、函数16六、 输入与输出结果216.1 输入数据216.2 输出结果226.3 当bus4-6的发生故障被切断后，系统的运行情况256.4 两机五节点网络30七、 结果分析347.1 三机九节点系统正常运行结果分析347.2 三机九节点系统故障切除后结果分析35结论36参考文献37一、 绪论1.1 本课题的目的和意义电力系统潮流计算是对复杂电力系统正常和故障条件下稳态运行状态的计算。其目的是求取电力系统在给定运行方式下的节点电压和功率分布，用以检查系统各元件是否过负荷、各点电压是否满足要求、功率分布和分配是否合理以及功率损耗等，是电力系统计算分析中的一种最基本的计算。潮流计算是电力系统的各种计

9、算的基础，同时它又是研究电力系统的一项重要分析功能，是进行故障计算，继电保护鉴定，安全分析的工具。电力系统潮流计算是计算系统动态稳定和静态稳定的基础。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用电力系统潮流计算来定量的比较供电方案或运行方式的合理性、可靠性和经济性。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置整定计算、电力系统故障计算和稳定计算等提供原始数据。潮流计算的目的在于:确定是电力系统的运行方式；检查系统中的各元件是否过压或过载；为电力系统继电保护的整定提供依据；为电力系统的稳定计算提供初值，为电力系

10、统规划和经济运行提供分析的基础。因此，电力系统潮流计算是电力系统中一项最基本的计算，既具有一定的独立性，又是研究其他问题的基础。1.2 国内外发展现状利用电子计算机进行潮流计算从20世纪50年代中期就已经开始。此后，潮流计算曾采用了各种不同的方法，这些方法的发展主要是围绕着对潮流计算的一些基本要求进行的。对潮流计算的要求可以归纳为下面几点： （1）算法的可靠性或收敛性 （2）计算速度和内存占用量 （3）计算的方便性和灵活性 电力系统潮流计算属于稳态分析范畴，不涉及系统元件的动态特性和过渡过程。因此其数学模型不包含微分方程，是一组高阶非线性方程。非线性代数方程组的解法离不开迭代，因此，潮流计算方

11、法首先要求它是能可靠的收敛，并给出正确答案。随着电力系统规模的不断扩大，潮流问题的方程式阶数越来越高，目前已达到几千阶甚至上万阶，对这样规模的方程式并不是采用任何数学方法都能保证给出正确答案的。这种情况促使电力系统的研究人员不断寻求新的更可靠的计算方法。 1.2.1 高斯-赛德尔迭代法在用数字计算机求解电力系统潮流问题的开始阶段，人们普遍采用以节点导纳矩阵为基础的高斯-赛德尔迭代法（一下简称导纳法）。这个方法的原理比较简单，要求的数字计算机的内存量也比较小，适应当时的电子数字计算机制作水平和电力系统理论水平，于是电力系统计算人员转向以阻抗矩阵为主的逐次代入法（以下简称阻抗法）。 20世纪60年

12、代初，数字计算机已经发展到第二代，计算机的内存和计算速度发生了很大的飞跃，从而为阻抗法的采用创造了条件。阻抗矩阵是满矩阵，阻抗法要求计算机储存表征系统接线和参数的阻抗矩阵。这就需要较大的内存量。而且阻抗法每迭代一次都要求顺次取阻抗矩阵中的每一个元素进行计算，因此，每次迭代的计算量很大。 阻抗法改善了电力系统潮流计算问题的收敛性，解决了导纳法无法解决的一些系统的潮流计算，在当时获得了广泛的应用，曾为我国电力系统设计、运行和研究做出了很大的贡献。但是，阻抗法的主要缺点就是占用计算机的内存很大，每次迭代的计算量很大。当系统不断扩大时，这些缺点就更加突出。为了克服阻抗法在内存和速度方面的缺点，后来发展

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

14、，成为直到目前仍被广泛采用的方法。 在牛拉法的基础上，根据电力系统的特点，抓住主要矛盾，对纯数学的牛拉法进行了改造，得到了P-Q分解法。P-Q分解法在计算速度方面有显著的提高，迅速得到了推广。 牛拉法的特点是将非线性方程线性化。20世纪70年代后期，有人提出采用更精确的模型，即将泰勒级数的高阶项也包括进来，希望以此提高算法的性能，这便产生了保留非线性的潮流算法。另外，为了解决病态潮流计算，出现了将潮流计算表示为一个无约束非线性规划问题的模型，即非线性规划潮流算法。 近20多年来，潮流算法的研究仍然非常活跃，但是大多数研究都是围绕改进牛拉法和P-Q分解法进行的。此外，随着人工智能理论的发展，遗传

15、算法、人工神经网络、模糊算法也逐渐被引入潮流计算。但是，到目前为止这些新的模型和算法还不能取代牛拉法和P-Q分解法的地位。由于电力系统规模的不断扩大，对计算速度的要求不断提高，计算机的并行计算技术也将在潮流计算中得到广泛的应用，成为重要的研究领域。通过几十年的发展，潮流算法日趋成熟。近几年，对潮流算法的研究仍然是如何改善传统的潮流算法，即高斯-塞德尔法、牛拉法和快速解耦法。牛拉法，由于其在求解非线性潮流方程时采用的是逐次线性化的方法，为了进一步提高算法的收敛性和计算速度，人们考虑采用将泰勒级数的高阶项或非线性项也考虑进来，于是产生了二阶潮流算法。后来又提出了根据直角坐标形式的潮流方程是一个二次

16、代数方程的特点，提出了采用直角坐标的保留非线性快速潮流算法8。1.2.3 基于MATLAB 的电力系统潮流计算发展前景MATLAB 自1980 年问世以来，以其学习简单、使用方便以及其它高级语言所无可比拟的强大的矩阵处理功能越来越受到世人的关注。目前，它已成为国际控制界最流行、使用最广泛的语言了。它的强大的矩阵处理功能给电力系统的分析、计算带来许多方便。在处理潮流计算时，其计算机软件的速度已无法满足大电网模拟和实时控制的仿真要求，而高效的潮流问题相关软件的研究已成为大规模电力系统仿真计算的关键。随着计算机技术的不断发展和成熟,对MATLAB 潮流计算的研究为快速、详细地解决大电网的计算问题开辟

17、了新思路。MATLAB 语言允许用户以数学形式的语言编写程序, 其比BASIC 语言和FORTRAN 等更为接近书写的数学表达格式, 且程序易于调试。在计算要求相同的情况下, 使用MATLAB 编程, 工作量将会大为减少。基于MATLAB 的电力系统潮流计算使计算机在计算、分析、研究复杂的电力系统潮流分布问题上又前进了一步。矩阵输入、输出格式简单, 与数学书写格式相似; 以双精度类型进行数据的存储和运算, 数据精确度高，能进行潮流计算中的各种矩阵运算, 包括求逆、求积和矩阵L R 分解等, 其程序的编写也因MATLAB 提供了许多功能函数而变得简单易行。另外, MATLAB稀疏矩阵技术的引入,

18、 使电力系统潮流计算由传统方法转变为优化算法成为可能。二、 设计目的通过学习Matlab软件的使用来巩固专业知识，通过设计和仿真复杂系统的潮流分布，进一步理解电力系统，用此软件算出电力系统中的各种矩阵和参数，并能进行校验。三、 设计要求和设计指标电力系统稳态分析课程设计以设计和优化电力系统的潮流分析为重点，提高了我们综合能力为目标，尽可能结合实际工程进行，使我们初步学会运用所学知识解决工程中的实际问题，主要用到了有功功率P，无功功率Q，电压额定值V，和角度等参数。四、 设计内容4.1 选题内容通过课程设计,使学生巩固电力系统三相短路计算的基本原理与方法，掌握短路电流的数值求解方法(节点导纳矩阵

19、，修正方程)，开发系统牛顿拉夫逊法的计算程序。让我们掌握用计算机仿真分析电力系统的方法。同时，通过软件开发，也有助于计算机操作能力和软件开发能力的提高。1、根据电力系统网络推导电力网络数学模型，写出节点导纳矩阵；2、赋予各节点电压变量（直角坐标系形式）初值后，求解不平衡量；3、形成雅可比矩阵；4、求解修正量后，重新修改初值，从2开始重新循环计算；5、求解的电压变量达到所要求的精度时，再计算各支路功率分布、功率损耗和平衡节点功率；6、上机编程调试；7、书写课程报告。4.2 基础资料三机九节点系统又称为WSCC（Western System Coordinating Council），系统接线图如

20、图1。系统的基准功率为100MVA，系统频率为50Hz。图1 三机九节点系统接线图系统参数见下表1-3。表1 母线数据编号基准电压（kV）PL(MW)QL(MVAR)PG(MW)QG(MVAR)V（p.u）q(rad)110.5/1.040.0210.5/163.0/1.03/310.5/85.0/1.03/4230/5230125.050.0/623090.030.0/7230/8230100.035.0/9230/表2 变压器数据首节点编号基准电压（kV）额定电压（kV）末节点编号基准电压（kV）额定电压（kV）R（p.u）X（p.u）110.510.542302200.00.057621

21、0.510.572302200.00.0625310.510.592302200.00.0586表3 输电线线数据首节点编号末节点编号R（p.u）X（p.u）B/2（p.u）890.01190.10080.1045690.0390.170.179460.0170.0920.079450.010.0850.088780.00850.0720.0745570.0320.1610.153潮流计算选择1号母线为参考节点，2、3号母线为PV节点，4号到9号节点为PQ节点。4.3算法原理4.3.1 节点导纳矩阵形成的计算机方法在研究电力网络节点导纳矩阵形成与修改的计算机方法之前，首先要设计一种电力网络参数

22、的输入方法。假定构成电力网络的元件包括，双绕组变压器，输电线路及母线接地支路。对变压器假定，三绕组变压器已转化为三个双绕组变压器的等值电路，且略去励磁导纳支路，若需考虑励磁导纳支路则将该支路视为接到相应节点上的母线接地支路处理。文中的母线接地支路可以表示母线上的接地电容和接地电感，亦可模拟短路故障等异常情况。本文介绍的电力网络参数的输入方法以电力网络元件为单位如表1.1每一行所示，其中：第一和第三个字段为字符串型数据，表示元件的首节点名和尾节点名。第二和第四个字段为整数型数据，表示元件的首节点和尾节点编号，如果尾节点编号为0表示该行数据对应的元件为接地支路。第五和第六个字段为实数型数据，若该行

23、对应变压器元件，这两个数据分别表示变压器绕组的等值电阻和电抗；若该行对应输电线（或母线接地支路元件），这两个数据分别表示该元件模型的串联等值电阻和电抗。第七和第八个字段为实数型数据，分别表示输电线型等值电路一端的接地电纳和变压器的变比，若第八个字段的数据为空表示该行对应输电线元件，否则为双绕组变压器元件。表1.1所示的输入文件内容为图1-1(b)所示电力网络的输入参数。需要注意的是双绕组变压器的变比统一规定为首节点至尾节点的变比为；另外允许电力网络存在并联元件，例如若图1-1网络增加一条与输电线L1相同的并联输电线，表1.1输入文件中只要再增补一行与第四行相同的数据即可；另外需要指出的是输入文

24、件对每行数据的先后次序无限制。表4 图2所示电力网络形成节点导纳矩阵的输入文件列表首节点名首节点编号i尾节点名尾节点编号jR(p.u.)X(p.u.)B/2(p.u.)k(p.u.)BUS11BUS220.00.1/1.05BUS55BUS330.00.175/1.0BUS22BUS330.030.080.02/BUS22BUS440.020.0530.013/BUS33BUS440.0180.0480.012/BUS44/00.0-20.0/图4(a) 电力系统接线图图4(b) 电力系统采用理想变压器的等值电路图图4(c) 包括节点注入化掉理想变压器的系统等值电路图图4一般情况下，节点电力网

25、络节点相量方程可以简记为： (1)其中节点导纳矩阵为阶复系数矩阵，为维节点电压相量构成的复向量，为维节点注入电流相量构成的复向量。节点导纳矩阵的形成方法对照表4所示的电力网络输入文件格式，形成n节点节点导纳矩阵的方法描述如下：读入每一行数据，且执行以下操作；若该行对应变压器元件，则作： (2)若该行对应输电线元件，则作： (3)若该行对母线接地支路元件，则作： (4)对电力网络输入文件的每一行执行上述操作后，阶复系数节点导纳矩阵便形成了。上述方法的特点是逐个将每一电力网络的元件对节点导纳矩阵的贡献添加到矩阵中。实际电力网络具有高维、稀疏及对称的特点，如何针对上述特点，并结合节点导纳矩阵的实际使

26、用设计节点导纳矩阵的稀疏存储方法，对提高电力系统分析的效率是一件非常有意义的研究工作。4.3.2 电力网络的潮流方程1）节点的功率方程对n节点电力系统，节点i的注入复功率与节点注入电流及电网节点电压的关系如式(5)表示。 (5)式中和表示节点注入电流相量和节点电压相量的共轭复数，表示电力网络节点导纳矩阵元素的共轭复数，表示号后的j节点为网络中和i通过元件相连的节点，并包括i=j的情况。考虑可得 (6)式(6)代入式(5)可得： (7)式中，表示i节点与j节点电压间的相角差。式(7)右端展开并将实部和虚部分开可得如下(8)和(9)电压用极坐标表示的节点功率平衡方程。 (8) 4.3.3 电力网络

27、极坐标形式的潮流方程当电力网络的节点类型确定后，便可以列写全网的潮流方程。设n节点电网PQ节点的个数是m个，这样对每个PQ和PV节点可列写一个有功功率方程（共有n-1个），对每个PQ节点可列写一个无功功率方程（共有m个）。式(12)和(13)为由式(8)和(9)所得节点功率方程。对PQ、PV节点： (12)对PQ节点： (13)式中和是节点指定的有功和无功注入功率。式(12)和式(13)联立可以表示成如下非线性向量方程的形式。 (14)式中，和分别表示n-1维（次序PQ节点的有功功率方程在前）和m维的非线性向量函数，和分别表示n-1维节点电压相角向量（次序PQ节点的相角在前）和m维节点电压幅值

28、向量。式(14)的方程个数和待求变量的个数皆为n+m-1，称作电力网络极坐标形式的潮流方程。4.3.4 极坐标牛顿潮流算法的修正方程雅和可比矩阵极坐标向量形式的潮流方程式(14)是非线性相量方程，一般采用牛顿法迭代计算其解。牛顿法对式(14)的修正方程如下 (15)式中和分别表示n-1维（次序为PQ节点的不平衡有功功率在前）和m维的节点不平衡功率向量，和分别表示待求的n-1维的节点电压相角修正量向量（次序PQ节点的相角在前）和待求的m维节点电压幅值修正量除以电压幅值的向量，即 (16) (17),和分别为, , 和阶的实系数雅可比子矩阵。 牛顿潮流算法的计算流程如图5所示图5 牛顿潮流算法4.

29、3.5 PQ分解法潮流计算节点的有功功率不平衡量只用于修正电压的相位，节点的无功功率不平衡量只用于修正电压的幅值，将2组的方程组分别进行迭代。在此过程中还要进行最为关键的一部化简就在于把系数矩阵H和L简化成为常数矩阵，在这2个修正方程式中，系数矩阵都由节点导纳矩阵的虚部构成，只是介次不同，矩阵B为n-1阶，不含平衡节点对应的行和列，矩阵B2为m阶，不含平衡节点和PV节点所对应的行和列，由于修正方程的系数矩阵未常数矩阵，只要作一次三角分解，即可反复使用，结合裁采用稀疏技巧，还可进一步的节省机器的内存和计算时间。PQ分解法的特点：1. 以一个n-1阶和一个n-m-1阶线性方程组代替原有的2n-m-

30、1阶线性方程组；2.修正方程的系数矩阵B和B”为对称常数矩阵，且在迭代过程中保持不变；3.P-Q分解法具有线性收敛特性，与牛顿-拉夫逊法相比，当收敛到同样的精度时需要的迭代次数较多；4.P-Q分解法一般只适用于110KV及以上电网的计算。因为35KV及以下电压等级的线路r/x比值很大，不满足上述简化条件，可能出现迭代计算不收敛的情况。4.3.5.1 程序流程五、 程序设计5.1 主函数%PQmatlab-power flowstudy programSbase_MVA=100. %指明基准容量epsilon=0.000001; %指定的潮流计算收敛精度max_no_itearation=20

31、%PQ分解算法fid=fopen(Nodedata.txt); %打开Nodedata.txt文件N=textscan(fid, %s %u %d %f %f %f %f %f %f); %输入母线数据fclose(fid);busnumber=size(N1,1) %数组N的维数即电网母线数%建立Bus结构体for i=1:busnumber Bus(i).name=N1(i); Bus(i).type=N2(i); %规定1平衡节点；3 PV节点；2 PQ节点 Bus(i).no=i; %节点编号自动生成。避免编辑出错 Bus(i).Base_KV=N3(i); Bus(i).PG=N4(

32、i); %母线发电有功（标幺值） Bus(i).QG=N5(i); %母线发电无功（标幺值） Bus(i).PL=N6(i); %母线有功负荷（标幺值） Bus(i).QL=N7(i); %母线无功负荷（标幺值） Bus(i).pb=N8(i); %母线接地电抗（标幺值） Bus(i).V=1.0; %母线电压幅值（标幺值） Bus(i).angle=0; %母线电压辐角（弧度）endfid=fopen(Aclinedata.txt); %打开Nodedata.txt文件A=textscan(fid, %s %s %f %f %f %f); %输入AC线数据fclose(fid);acline

33、number=size(A1,1) %数组A的维数即电网AC线条数%建立Acline结构体for i=1:aclinenumber Acline(i).fbname=A1(i); %AC线始节点名 Acline(i).tbname=A2(i); %AC线终节点名 Acline(i).Base_KV=A3(i); %AC线基准电压 Acline(i).R=A4(i); %AC线串联R（标幺值） Acline(i).X=A5(i); %AC线串联X（标幺值） Acline(i).hB=A6(i); %AC线串联0.5*B（标幺值） for k=1:busnumber if strcmp(Aclin

34、e(i).fbname, Bus(k).name) Acline(i).fbno=Bus(k).no; %AC线始节点编号 end if strcmp(Acline(i).tbname, Bus(k).name) Acline(i).tbno=Bus(k).no; %AC线终节点编号 end end endfid=fopen(Transdata.txt); %打开Transdata.txt文件T=textscan(fid, %s %f %f %s %f %f %f %f); %输入变压器数据fclose(fid);transnumber=size(T1,1) %数组T的维数即电网变压器个数（t

35、ansnumber）用变量名tansnumber表示if transnumber=0 Trans(i).fbname=bus0;else for i=1:transnumber Trans(i).fbname=T1(i); %变压器始节点名 Trans(i).fbBase_KV=T2(i); %变压器始节点基准电压 Trans(i).fbrated_KV=T3(i); %变压器始节点额定电压 Trans(i).tbname=T4(i); %变压器终节点名 Trans(i).tbBase_KV=T5(i); %变压器终节点基准电压 Trans(i).tbrated_KV=T6(i); %变压器终

36、节点额定电压 Trans(i).R=T7(i); %变压器线串联R（标幺值） Trans(i).X=T8(i); %变压器线串联X（标幺值） for k=1:busnumber if strcmp(Trans(i).fbname, Bus(k).name) Trans(i).fbno=Bus(k).no; %变压器始节点编号 end if strcmp(Trans(i).tbname, Bus(k).name) Trans(i).tbno=Bus(k).no; %变压器终节点编号 end end Trans(i).k=Trans(i).tbrated_KV*Trans(i).fbBase_KV

37、/Trans(i).fbrated_KV/Trans(i).tbBase_KV; %计算变比 tempx=Trans(i).fbrated_KV2/Trans(i).fbBase_KV2; Trans(i).X=tempx*Trans(i).X; %折算后的变压器线串联X（标幺值） Trans(i).R=tempx*Trans(i).R; %折算后的变压器线串联R（标幺值）endendG,B,B2=FormYmatrix(Bus,busnumber,Acline,aclinenumber,Trans,transnumber); %B:=B;B2:=B% G:节点导纳矩阵Y的实部, B:节点导纳

38、矩阵Y的虚部,B2对应只用支路电抗,变压器电抗形成的B矩阵 JP,JQ=FormJPQmatrix(Bus,B,B2,busnumber);B2 %显示B2矩阵JPiJP=-inv(JP); B %显示B矩阵JQiJQ=-inv(JQ);Bus=solve_flow_eq(Bus,B,G,busnumber,iJP,iJQ,epsilon,max_no_itearation);pause%计算潮流分布Sf_T,St_T,Sf_L,St_L=comput_flow_distribution(Bus,Acline,Trans,aclinenumber,transnumber,Sbase_MVA);

39、if transnumber0Sf_T %输出变压器始节点潮流St_T %输出变压器终节点潮流endSf_L %输出输电线始节点潮流St_L %输出输电线终节点潮流clear5.2 子函数%生成G、B矩阵function G,B,X=FormYmatrix(Bus,busnumber,Acline,aclinenumber,Trans,tansnumber)Y=zeros(busnumber);X=zeros(busnumber);for i=1:busnumberY(i,i)=Y(i,i)+Bus(i).pb*j;endfor i=1:aclinenumber f=Acline(i).fbn

40、o; t=Acline(i).tbno; Y(f,f)=Y(f,f)+Acline(i).hB*j+1/(Acline(i).R+Acline(i).X*j); Y(t,t)=Y(t,t)+Acline(i).hB*j+1/(Acline(i).R+Acline(i).X*j); Y(f,t)=Y(f,t)-1/(Acline(i).R+Acline(i).X*j); Y(t,f)=Y(t,f)-1/(Acline(i).R+Acline(i).X*j); X(f,f)=X(f,f)-1/Acline(i).X; X(t,t)=X(t,t)-1/Acline(i).X; X(f,t)=1/Ac

41、line(i).X; X(t,f)=1/Acline(i).X; endfor i=1:tansnumber f=Trans(i).fbno; t=Trans(i).tbno; Y(f,f)=Y(f,f)+1/(Trans(i).R+Trans(i).X*j); Y(t,t)=Y(t,t)+1/(Trans(i).R+Trans(i).X*j)/Trans(i).k2; Y(f,t)=Y(f,t)-1/(Trans(i).R+Trans(i).X*j)/Trans(i).k; Y(t,f)=Y(t,f)-1/(Trans(i).R+Trans(i).X*j)/Trans(i).k; X(f,f

42、)=X(f,f)-1/Trans(i).X; X(t,t)=X(t,t)-1/Trans(i).X; X(f,t)=1/Trans(i).X; X(t,f)=1/Trans(i).X; endG=real(Y);B=imag(Y);end%生成JP、JQ矩阵function JP,JQ=FormJPQmatrix(Bus,B,B2,busnumber)JP=B2;JQ=B;for i=1:busnumber if Bus(i).type=1 %该循环将JP和JQ矩阵对应平衡节点的非对角元素置零,对角元素置1(相当于去掉了对应平衡节点的JP和JQ矩阵的行和列) for k=1:busnumber

43、 JQ(i,k)=0.; JQ(k,i)=0.; JP(i,k)=0.; JP(k,i)=0.; end JQ(i,i)=1.; JP(i,i)=1.; end if Bus(i).type=3 %该循环将JQ矩阵对应PV节点的非对角元素置零,对角元素置1(相当于去掉了对应PV节点JQ矩阵的行和列) for k=1:busnumber JQ(i,k)=0.; JQ(k,i)=0.; end JQ(i,i)=1.; end endend%解潮流函数function Bus=solve_flow_eq(Bus,B,G,busnumber,iJP,iJQ,epsilon,max_no_itearat

44、ion)for i=1:busnumber NodeV(i)=Bus(i).V; Nodea(i)=Bus(i).angle; dQGQL(i)=Bus(i).QG-Bus(i).QL; dPGPL(i)=Bus(i).PG-Bus(i).PL;endNodeV=NodeV; %节点电压模值列向量初值Nodea=Nodea; %节点电压相角列向量初值dQGQL=dQGQL; %节点注入无功列向量dPGPL=dPGPL; %节点注入无功列向量pausemaxdP=1.; %最大节点不平衡有功初值maxdQ=1.; %最大节点不平衡无功初值epsilon=0.000001; %设置潮流计算精度no

45、iteration=0; %迭代初始次数为零while (maxdPepsilon)&(maxdQepsilon) %条件1：(maxdPepsilon)且条件2：(maxdPepsilon) %都成立时继续循环deltaP,deltaQ,maxdP,maxdQ=FormdPQvector(Bus,NodeV,Nodea,dQGQL,dPGPL,B,G,busnumber);%函数FormdPQvector调用deltaP; %节点不平衡有功向量deltaQ; %节点不平衡有功向量maxdP %最大节点不平衡有功maxdQ %最大节点不平衡无功da=iJP*deltaP; %解修正方程dV=i

46、JQ*deltaQ; %解修正方程Nodea=Nodea+da; %修正节点相角NodeV=NodeV+dV; %修正节点电压noiteration=noiteration+1; %迭代次数累计pause if noiterationmax_no_itearation %if语句，当迭代次数大于最大次数终止迭代 break end endfor i=1:busnumber Bus(i).V=NodeV(i); %计算节点电压结果存入循环体 NodeV(i)=NodeV(i)*Bus(i).Base_KV; %转换成有名值 Bus(i).angle=Nodea(i); %计算节点相角（单位弧度）

47、结果存入循环体 Nodea(i)=Nodea(i)*180/pi; %单位弧度转换成单位度end%计算支路潮流程序在这里加入，利用计算出的节点电压及结构体中的信息noiterationNodea=Nodea %输出母线电压kVNodeV=NodeV %输出母线电压相角/度End%计算偏节点PQ差量function deltaP,deltaQ,maxdP,maxdQ=FormdPQvector(Bus,NodeV,Nodea,dQGQL,dPGPL,B,G,busnumber) deltaQ=dQGQL; %置初值，为无功发电减去无功负荷 deltaP=dPGPL; %置初值，为有功发电减去有功

48、负荷 maxdP=0.; %置初值 maxdQ=0.; %置初值 for i=1:busnumber if Bus(i).type=1 %平衡节点 deltaQ(i)=0.; deltaP(i)=0.; end if Bus(i).type=3 deltaQ(i)=0.; %y1=0; %y2=0; y3=0; for k=1:busnumber if (B(i,k)=0|G(i,k)=0) %y1=y1+(G(i,k)*VX(k)-B(i,k)*VY(k); %y2=y2+(G(i,k)*VY(k)+B(i,k)*VX(k); y3=y3+NodeV(k)*(G(i,k)*cos(Nodea

49、(i)-Nodea(k)+B(i,k)*sin(Nodea(i)-Nodea(k); end end deltaP(i)=deltaP(i)-y3*NodeV(i); %deltaP2(i)=(deltaP2(i)-(y1*VX(i)+y2*VY(i)/Bus(i).V; end if Bus(i).type=2 %PQ节点计算有功偏差和无功偏差 %y1=0; %y2=0; y3=0; y4=0; for k=1:busnumber if (B(i,k)=0|G(i,k)=0) %y1=y1+(G(i,k)*VX(k)-B(i,k)*VY(k); %y2=y2+(G(i,k)*VY(k)+B(

50、i,k)*VX(k); y3=y3+NodeV(k)*(G(i,k)*cos(Nodea(i)-Nodea(k)+B(i,k)*sin(Nodea(i)-Nodea(k); y4=y4+NodeV(k)*(G(i,k)*sin(Nodea(i)-Nodea(k)-B(i,k)*cos(Nodea(i)-Nodea(k); end end deltaP(i)=deltaP(i)-y3*NodeV(i); %deltaP2(i)=(deltaP2(i)-(y1*VX(i)+y2*VY(i)/Bus(i).V; deltaQ(i)=deltaQ(i)-y4*NodeV(i); %deltaQ2(i)

51、=(deltaQ2(i)-(y1*VY(i)-y2*VX(i)/Bus(i).V; end %以下程序断找出 deltaP和 deltaQ向量中的绝对值最大元素 if maxdPabs(deltaP(i); maxdP=abs(deltaP(i); end if maxdQ0 for l=1:transnumber Z=Trans(l).R+j*Trans(l).X; %变压器阻抗 f=Trans(l).fbno; % f 始节点编号 Vf=Bus(f).V*cos(Bus(f).angle)+j*Bus(f).V*sin(Bus(f).angle); %始节点电压相量 t=Trans(l).

52、tbno; % t 终节点编号 Vt=(Bus(t).V*cos(Bus(t).angle)+j*Bus(t).V*sin(Bus(t).angle)/Trans(l).k; %终节点理想变压器内侧电压相量 I=(Vf-Vt)/Z; %变压器电流相量 Sf_T(l)=Vf*conj(I)*Sbase_MVA; %变压器始节点复功率/MVar St_T(l)=Vt*conj(I)*Sbase_MVA; %变压器终节点复功率/MVar Trans(l).fS=abs(Sf_T(l); Trans(l).tS=abs(St_T(l); endelse Sf_T(1)=0; St_T(1)=0;end

53、pausefor l=1:aclinenumber Z=Acline(l).R+j*Acline(l).X %线路阻抗 f=Acline(l).fbno; Vf=Bus(f).V*cos(Bus(f).angle)+j*Bus(f).V*sin(Bus(f).angle); %始节点电压相量 t=Acline(l).tbno; Vt=Bus(t).V*cos(Bus(t).angle)+j*Bus(t).V*sin(Bus(t).angle); %终节点电压相量 I=(Vf-Vt)/Z; %线路电流相量 Sf_L(l)=(Vf*conj(I)-j*(Bus(f).V)2*Acline(l).hB)*Sbase_MVA; %线路始节点复功率/MVar St_L(l)=(Vt*conj(I)+j*(Bus(t).V)2*Acline(l).hB)*Sbase_MVA; %线路终节点复功率/MVar Acline(l).fS=abs(Sf_L(l); Acline(l).tS=abs(St_L(l); endend六、 输入与输出结果6.1 输入数据节点数据（Nodedata.txt）bus1 118. 0. 0.0. 0.0.bus