1、电力系统分析课程设计报告题目：复杂电力系统潮流分析程序设计所在学院 电气工程学院 专业班级 12 级电气工程及其自动化一班学生姓名 陈剑秋 学生学号 201230088099 指导教师 房大中 提交日期 2014 年 12 月 16 日 华南理工大学广州学院电气工程学院电力系统分析课程设计报告I摘要电力系统潮流计算是电力系统分析课程基本计算的核心部分之一。它既有本身的独立意义，又是电力系统规划设计，运行分析和理论研究的基础。电力系统潮流计算是对复杂电力系统正常和故障条件下稳态运行状态的计算，是电力系统计算分析中的一种最基本的计算。潮流计算是电力系统的各种计算的基础，同时它又是研究电力系统的一项
2、重要分析功能，是进行故障计算，继电保护鉴定，安全分析的工具。电力系统潮流计算是计算系统动态稳定和静态稳定的基础。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用电力系统潮流计算来定量的比较供电方案或运行方式的合理性、可靠性和经济性本课程设计采用 P-Q 分解法对三机九节点系统进行潮流计算。计算结果的得出是通过 MATLAB 软件编程实现。关键词：电力系统潮流计算，P-Q 分解法，三机九节点系统华南理工大学广州学院电气工程学院电力系统分析课程设计报告IIAbstractPower flow calculation is the core part of the power syste
3、m analysis course basic computing. It has both the independent 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 syste
4、m, is one of the most basic calculation of power system calculation 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, safet
5、y analysis tool. Power flow calculation is the basis of calculation 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 req
6、uires the use of power flow calculation to quantitativeThe curriculum 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 calcul
7、ation, the P-Q decomposition method, three machine nine node system华南理工大学广州学院电气工程学院电力系统分析课程设计报告III目录摘要 .IAbstract.II一、 绪论 .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 节
8、点导纳矩阵形成的计算机方法 64.3.2 电力网络的潮流方程 94.3.3 电力网络极坐标形式的潮流方程 104.3.4 极坐标牛顿潮流算法的修正方程雅和可比矩阵 104.3.5 PQ 分解法潮流计算 .11五、 程序设计 .135.1 主函数 .135.2 子函数 .16六、 输入与输出结果 .216.1 输入数据 .216.2 输出结果 .226.3 当 bus4-6 的发生故障被切断后，系统的运行情况 .256.4 两机五节点网络 .30七、 结果分析 .347.1 三机九节点系统正常运行结果分析 .347.2 三机九节点系统故障切除后结果分析 .35结论 .36参考文献 .37华南理工
9、大学广州学院电气工程学院电力系统分析课程设计报告1一、 绪论1.1 本课题的目的和意义电力系统潮流计算是对复杂电力系统正常和故障条件下稳态运行状态的计算。其目的是求取电力系统在给定运行方式下的节点电压和功率分布，用以检查系统各元件是否过负荷、各点电压是否满足要求、功率分布和分配是否合理以及功率损耗等，是电力系统计算分析中的一种最基本的计算。潮流计算是电力系统的各种计算的基础，同时它又是研究电力系统的一项重要分析功能，是进行故障计算，继电保护鉴定，安全分析的工具。电力系统潮流计算是计算系统动态稳定和静态稳定的基础。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用电力系统潮流计算来定量
10、的比较供电方案或运行方式的合理性、可靠性和经济性。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置整定计算、电力系统故障计算和稳定计算等提供原始数据。潮流计算的目的在于:确定是电力系统的运行方式；检查系统中的各元件是否过压或过载；为电力系统继电保护的整定提供依据；为电力系统的稳定计算提供初值，为电力系统规划和经济运行提供分析的基础。因此，电力系统潮流计算是电力系统中一项最基本的计算，既具有一定的独立性，又是研究其他问题的基础。1.2 国内外发展现状利用电子计算机进行潮流计算从 20 世纪 50 年代中期就已经开始。此后，潮流计算
11、曾采用了各种不同的方法，这些方法的发展主要是围绕着对潮流计算的一些基本要求进行的。对潮流计算的要求可以归纳为下面几点：（1）算法的可靠性或收敛性 （2）计算速度和内存占用量 （3）计算的方便性和灵活性 华南理工大学广州学院电气工程学院电力系统分析课程设计报告2电力系统潮流计算属于稳态分析范畴，不涉及系统元件的动态特性和过渡过程。因此其数学模型不包含微分方程，是一组高阶非线性方程。非线性代数方程组的解法离不开迭代，因此，潮流计算方法首先要求它是能可靠的收敛，并给出正确答案。随着电力系统规模的不断扩大，潮流问题的方程式阶数越来越高，目前已达到几千阶甚至上万阶，对这样规模的方程式并不是采用任何数学方
12、法都能保证给出正确答案的。这种情况促使电力系统的研究人员不断寻求新的更可靠的计算方法。 1.2.1 高斯-赛德尔迭代法在用数字计算机求解电力系统潮流问题的开始阶段，人们普遍采用以节点导纳矩阵为基础的高斯-赛德尔迭代法（一下简称导纳法） 。这个方法的原理比较简单，要求的数字计算机的内存量也比较小，适应当时的电子数字计算机制作水平和电力系统理论水平，于是电力系统计算人员转向以阻抗矩阵为主的逐次代入法（以下简称阻抗法） 。 20 世纪 60 年代初，数字计算机已经发展到第二代，计算机的内存和计算速度发生了很大的飞跃，从而为阻抗法的采用创造了条件。阻抗矩阵是满矩阵，阻抗法要求计算机储存表征系统接线和参
13、数的阻抗矩阵。这就需要较大的内存量。而且阻抗法每迭代一次都要求顺次取阻抗矩阵中的每一个元素进行计算，因此，每次迭代的计算量很大。 阻抗法改善了电力系统潮流计算问题的收敛性，解决了导纳法无法解决的一些系统的潮流计算，在当时获得了广泛的应用，曾为我国电力系统设计、运行和研究做出了很大的贡献。但是，阻抗法的主要缺点就是占用计算机的内存很大，每次迭代的计算量很大。当系统不断扩大时，这些缺点就更加突出。为了克服阻抗法在内存和速度方面的缺点，后来发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需存储各个地区系统的阻抗矩阵及它们之间的联络线的阻抗，这样不仅大幅度的
14、节省了内存容量，同时也提高了计算速度。1.2.2 牛顿-拉夫逊法和 P-Q 分解法克服阻抗法缺点的另一途径是采用牛顿-拉夫逊法（以下简称牛拉法） 。牛拉法是数学中求解非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳矩阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流程序的计算效率。自从 20 世纪 60 年代中期采用了最佳顺序消去法以后，牛拉法在收敛性、内存要求、计算速度方面都超过了阻抗法，成为直到目前仍被广泛采用的方法。 华南理工大学广州学院电气工程学院电力系统分析课程设计报告3在牛拉法的基础上，根据电力系统的特点，抓住主要矛盾
15、，对纯数学的牛拉法进行了改造，得到了 P-Q 分解法。P-Q 分解法在计算速度方面有显著的提高，迅速得到了推广。 牛拉法的特点是将非线性方程线性化。20 世纪 70 年代后期，有人提出采用更精确的模型，即将泰勒级数的高阶项也包括进来，希望以此提高算法的性能，这便产生了保留非线性的潮流算法。另外，为了解决病态潮流计算，出现了将潮流计算表示为一个无约束非线性规划问题的模型，即非线性规划潮流算法。 近 20 多年来，潮流算法的研究仍然非常活跃，但是大多数研究都是围绕改进牛拉法和 P-Q 分解法进行的。此外，随着人工智能理论的发展，遗传算法、人工神经网络、模糊算法也逐渐被引入潮流计算。但是，到目前为止
16、这些新的模型和算法还不能取代牛拉法和 P-Q 分解法的地位。由于电力系统规模的不断扩大，对计算速度的要求不断提高，计算机的并行计算技术也将在潮流计算中得到广泛的应用，成为重要的研究领域。通过几十年的发展，潮流算法日趋成熟。近几年，对潮流算法的研究仍然是如何改善传统的潮流算法，即高斯-塞德尔法、牛拉法和快速解耦法。牛拉法，由于其在求解非线性潮流方程时采用的是逐次线性化的方法，为了进一步提高算法的收敛性和计算速度，人们考虑采用将泰勒级数的高阶项或非线性项也考虑进来，于是产生了二阶潮流算法。后来又提出了根据直角坐标形式的潮流方程是一个二次代数方程的特点，提出了采用直角坐标的保留非线性快速潮流算法8。
17、1.2.3 基于 MATLAB 的电力系统潮流计算发展前景MATLAB 自 1980 年问世以来，以其学习简单、使用方便以及其它高级语言所无可比拟的强大的矩阵处理功能越来越受到世人的关注。目前，它已成为国际控制界最流行、使用最广泛的语言了。它的强大的矩阵处理功能给电力系统的分析、计算带来许多方便。在处理潮流计算时，其计算机软件的速度已无法满足大电网模拟和实时控制的仿真要求，而高效的潮流问题相关软件的研究已成为大规模电力系统仿真计算的关键。随着计算机技术的不断发展和成熟,对 MATLAB 潮流计算的研究为快速、详细地解决大电网的计算问题开辟了新思路。MATLAB 语言允许用户以数学形式的语言编写
18、程序, 其比 BASIC 语言和 FORTRAN 等更为接近书写的数学表达格式, 且程序易于调试。在计算要求相同的情况下, 使用 MATLAB 编程, 工作量将会大为减少。基于 MATLAB 的电力系统潮流计算使计算机在计算、分析、研究复杂的电力系统潮流分布问题上又前进了一步。矩阵输入、输出格式简单, 与数华南理工大学广州学院电气工程学院电力系统分析课程设计报告4学书写格式相似; 以双精度类型进行数据的存储和运算 , 数据精确度高，能进行潮流计算中的各种矩阵运算, 包括求逆、求积和矩阵 L R 分解等, 其程序的编写也因 MATLAB 提供了许多功能函数而变得简单易行。另外, MATLAB 稀
19、疏矩阵技术的引入, 使电力系统潮流计算由传统方法转变为优化算法成为可能。二、 设计目的通过学习 Matlab 软件的使用来巩固专业知识，通过设计和仿真复杂系统的潮流分布，进一步理解电力系统，用此软件算出电力系统中的各种矩阵和参数，并能进行校验。三、 设计要求和设计指标电力系统稳态分析课程设计以设计和优化电力系统的潮流分析为重点，提高了我们综合能力为目标，尽可能结合实际工程进行，使我们初步学会运用所学知识解决工程中的实际问题，主要用到了有功功率 P，无功功率Q，电压额定值 V，和角度 等参数。四、 设计内容4.1 选题内容通过课程设计,使学生巩固电力系统三相短路计算的基本原理与方法，掌握短路电流
20、的数值求解方法(节点导纳矩阵，修正方程)，开发系统牛顿拉夫逊法的计算程序。让我们掌握用计算机仿真分析电力系统的方法。同时，通过软件开发，也有助于计算机操作能力和软件开发能力的提高。1、根据电力系统网络推导电力网络数学模型，写出节点导纳矩阵；2、赋予各节点电压变量（直角坐标系形式）初值后，求解不平衡量；3、形成雅可比矩阵；4、求解修正量后，重新修改初值，从 2 开始重新循环计算；华南理工大学广州学院电气工程学院电力系统分析课程设计报告55、求解的电压变量达到所要求的精度时，再计算各支路功率分布、功率损耗和平衡节点功率；6、上机编程调试；7、书写课程报告。4.2 基础资料三机九节点系统又称为 WS
21、CC（Western System Coordinating Council） ，系统接线图如图 1。系统的基准功率为 100MVA，系统频率为 50Hz。279836415 (3)f图 1 三机九节点系统接线图系统参数见下表 1-3。表 1 母线数据编号基准电压（kV）PL(MW)QL(MVAR)PG(MW)QG(MVAR)V（p.u）(rad)1 10.5 / / / / 1.04 0.02 10.5 / / 163.0 / 1.03 /3 10.5 / / 85.0 / 1.03 /4 230 / / / / / /5 230 125.0 50.0 / / / /6 230 90.0 3
22、0.0 /华南理工大学广州学院电气工程学院电力系统分析课程设计报告67 230 / / / / / /8 230 100.0 35.0 / / / /9 230 / / / / / /表 2 变压器数据首节点编号基准电压（kV）额定电压（kV）末节点编号基准电压（kV）额定电压（kV）R（p.u）X（p.u）1 10.5 10.5 4 230 220 0.0 0.05762 10.5 10.5 7 230 220 0.0 0.06253 10.5 10.5 9 230 220 0.0 0.0586表 3 输电线线数据首节点编号 末节点编号 R（p.u） X（p.u） B/2（p.u）8 9 0
23、.0119 0.1008 0.10456 9 0.039 0.17 0.1794 6 0.017 0.092 0.0794 5 0.01 0.085 0.0887 8 0.0085 0.072 0.07455 7 0.032 0.161 0.153潮流计算选择 1号母线为参考节点，2、3 号母线为 PV节点，4 号到 9号节点为 PQ节点。4.3 算法原理4.3.1 节点导纳矩阵形成的计算机方法在研究电力网络节点导纳矩阵形成与修改的计算机方法之前，首先要设计一种电力网络参数的输入方法。假定构成电力网络的元件包括，双绕华南理工大学广州学院电气工程学院电力系统分析课程设计报告7组变压器，输电线路及
24、母线接地支路。对变压器假定，三绕组变压器已转化为三个双绕组变压器的等值电路，且略去励磁导纳支路，若需考虑励磁导纳支路则将该支路视为接到相应节点上的母线接地支路处理。文中的母线接地支路可以表示母线上的接地电容和接地电感，亦可模拟短路故障等异常情况。本文介绍的电力网络参数的输入方法以电力网络元件为单位如表 1.1 每一行所示，其中：第一和第三个字段为字符串型数据，表示元件的首节点名和尾节点名。第二和第四个字段为整数型数据，表示元件的首节点和尾节点编号，如果尾节点编号为 0 表示该行数据对应的元件为接地支路。第五和第六个字段为实数型数据，若该行对应变压器元件，这两个数据分别表示变压器绕组的等值电阻和
25、电抗；若该行对应输电线（或母线接地支路元件） ，这两个数据分别表示该元件模型的串联等值电阻和电抗。第七和第八个字段为实数型数据，分别表示输电线 型等值电路一端的接地电纳和变压器的变比，若第八个字段的数据为空表示该行对应输电线元件，否则为双绕组变压器元件。表 1.1 所示的输入文件内容为图 1-1(b)所示电力网络的输入参数。需要注意的是双绕组变压器的变比统一规定为首节点至尾节点的变比为 ；1/K另外允许电力网络存在并联元件，例如若图 1-1 网络增加一条与输电线 L1相同的并联输电线，表 1.1 输入文件中只要再增补一行与第四行相同的数据即可；另外需要指出的是输入文件对每行数据的先后次序无限制
26、。表 4 图 2 所示电力网络形成节点导纳矩阵的输入文件列表首节点名首节点编号i尾节点名尾节点编号jR(p.u.)X(p.u.)B/2(p.u.)k(p.u.)BUS11 BUS22 0.0 0.1 / 1.05BUS55 BUS33 0.0 0.175/ 1.0华南理工大学广州学院电气工程学院电力系统分析课程设计报告8BUS22 BUS33 0.03 0.08 0.02 /BUS22 BUS44 0.02 0.0530.013/BUS33 BUS44 0.0180.0480.012/BUS44 / 0 0.0 -20.0/ / 1L231T2TCGS GSG4S5图 4(a) 电力系统接线图
27、14230y320yz204 340z2z 51:k3:1kz图 4(b) 电力系统采用理想变压器的等值电路图12y1I 34530y230y200204 430yI III图 4(c) 包括节点注入化掉理想变压器的系统等值电路图图 4一般情况下， 节点电力网络节点相量方程可以简记为：n(1)YU=I其中节点导纳矩阵 为 阶复系数矩阵， 为 维节点电压相量构成的复Ynn向量， 为 维节点注入电流相量构成的复向量。In华南理工大学广州学院电气工程学院电力系统分析课程设计报告9节点导纳矩阵的形成方法对照表 4 所示的电力网络输入文件格式，形成 n 节点节点导纳矩阵的方法描述如下：读入每一行数据，且
28、执行以下操作；若该行对应变压器元件，则作：(2)211;i,i, j,j,i,ji,j j,ij,iYYRjXk(RjX)k() 若该行对应输电线元件，则作：(3)11;22i,i, j,j,i,ji,j j,ij,iBBYjYRXRjXj 若该行对母线接地支路元件，则作：(4)1i,i,YRjX对电力网络输入文件的每一行执行上述操作后， 阶复系数节点导纳n矩阵 便形成了。上述方法的特点是逐个将每一电力网络的元件对节点导Y纳矩阵的贡献添加到矩阵中。实际电力网络具有高维、稀疏及对称的特点，如何针对上述特点，并结合节点导纳矩阵的实际使用设计节点导纳矩阵的稀疏存储方法，对提高电力系统分析的效率是一件
29、非常有意义的研究工作。4.3.2 电力网络的潮流方程1）节点的功率方程对 n 节点电力系统，节点 i 的注入复功率 与节点注入电流及电网节点iS电压的关系如式(5)表示。(5)iiiiiijSPjQUIY式中 和 表示节点注入电流相量 和节点电压相量 的共轭复数， 表iIjUi jUijY示电力网络节点导纳矩阵元素 的共轭复数， 表示 号后的 j 节点为网ijYji络中和 i 通过元件相连的节点，并包括 i=j 的情况。考虑 可得jje(6)jiijijijjIU(GB)e式(6)代入式(5)可得：华南理工大学广州学院电气工程学院电力系统分析课程设计报告10(7)ijiiiijiijijjii
30、jijijijjPQUY(GB)Ue(GB)cosn式中 ，表示 i 节点与 j 节点电压间的相角差。式(7)右端展开并将ijij实部和虚部分开可得如下(8)和(9)电压用极坐标表示的节点功率平衡方程。(8)12iijiijijijjPU(GcosBn),(,n)12iijiijijijjQUsBcos,(,n)4.3.3 电力网络极坐标形式的潮流方程当电力网络的节点类型确定后，便可以列写全网的潮流方程。设 n 节点电网 PQ 节点的个数是 m 个，这样对每个 PQ 和 PV 节点可列写一个有功功率方程（共有 n-1 个） ，对每个 PQ 节点可列写一个无功功率方程（共有 m 个） 。式(12
31、)和(13)为由式(8)和(9)所得节点功率方程。对 PQ、PV 节点：(12)()()0(12)siiijiijijijjP,UGcosBn,i,n对 PQ 节点：(13)()()()siiijiijijijjQ, scs,i,m式中 和 是节点指定的有功和无功注入功率。式(12)和式(13)联立可以siPsi表示成如下非线性向量方程的形式。(14)(),PU0Q式中， 和 分别表示 n-1 维（次序 PQ 节点的有功功率方程在(),PV(),前）和 m 维的非线性向量函数， 和 分别表示 n-1 维节点电压相角向量（次序 PQ 节点的相角在前）和 m 维节点电压幅值向量。式(14)的方程个
32、数和待求变量的个数皆为 n+m-1，称作电力网络极坐标形式的潮流方程。4.3.4 极坐标牛顿潮流算法的修正方程雅和可比矩阵极坐标向量形式的潮流方程式(14)是非线性相量方程，一般采用牛顿法迭代计算其解。牛顿法对式(14)的修正方程如下(15)(k)(k)(k)(k)PHNQMLU式中 和 分别表示 n-1 维（次序为 PQ 节()(k)k)(,PU()(k)k)(,点的不平衡有功功率在前）和 m 维的节点不平衡功率向量， 和(k)华南理工大学广州学院电气工程学院电力系统分析课程设计报告11分别表示待求的 n-1 维的节点电压相角修正量向量（次序 PQ 节点(k)()U的相角在前）和待求的 m
33、维节点电压幅值修正量除以电压幅值的向量，即(16)121T(k)(k)(k)(k)n(17)2T(k)()()()()()(k)()mUU, , 和 分别为 , , 和 阶的(k)H(k)N(k)M()Ln实系数雅可比子矩阵。牛顿潮流算法的计算流程如图 5 所示输入电网及节点注入数据计算节点导纳矩阵参数设置节点电压初值(0)x0k计算 (k)(k)fxU1(k)(k)i iimaxf ?求解修正方程获得 (1)()()kkkxx()kmaxkK?设置 及最大迭代次数 maxK是是否否输出结果结束计算潮流分布12345 - 267891 01 11 21 3计算雅可比矩阵 与B5 - 1 输出潮
34、流不收敛信息华南理工大学广州学院电气工程学院电力系统分析课程设计报告12图 5 牛顿潮流算法4.3.5 PQ 分解法潮流计算节点的有功功率不平衡量只用于修正电压的相位，节点的无功功率不平衡量只用于修正电压的幅值，将 2 组的方程组分别进行迭代。在此过程中还要进行最为关键的一部化简就在于把系数矩阵 H 和 L 简化成为常数矩阵，在这 2 个修正方程式中，系数矩阵都由节点导纳矩阵的虚部构成，只是介次不同，矩阵 B 为 n-1 阶，不含平衡节点对应的行和列，矩阵 B2 为 m阶，不含平衡节点和 PV 节点所对应的行和列，由于修正方程的系数矩阵未常数矩阵，只要作一次三角分解，即可反复使用，结合裁采用稀
35、疏技巧，还可进一步的节省机器的内存和计算时间。PQ 分解法的特点：1. 以一个 n-1 阶和一个 n-m-1 阶线性方程组代替原有的 2n-m-1 阶线性方程组；2.修正方程的系数矩阵 B和 B”为对称常数矩阵，且在迭代过程中保持不变；3.P-Q 分解法具有线性收敛特性，与牛顿-拉夫逊法相比，当收敛到同样的精度时需要的迭代次数较多；4.P-Q 分解法一般只适用于 110KV 及以上电网的计算。因为 35KV 及以下电压等级的线路 r/x 比值很大，不满足上述简化条件，可能出现迭代计算不收敛的情况。华南理工大学广州学院电气工程学院电力系统分析课程设计报告134.3.5.1 程序流程五、 程序设计
36、5.1 主函数%PQmatlab-power flowstudy program华南理工大学广州学院电气工程学院电力系统分析课程设计报告14Sbase_MVA=100. %指明基准容量epsilon=0.000001; %指定的潮流计算收敛精度max_no_itearation=20 %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 的维数即电网母线数%建立 Bu
37、s 结构体for i=1:busnumberBus(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(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; %母线电压幅
38、值（标幺值）Bus(i).angle=0; %母线电压辐角（弧度）endfid=fopen(Aclinedata.txt); %打开 Nodedata.txt 文件A=textscan(fid, %s %s %f %f %f %f); %输入 AC 线数据fclose(fid);aclinenumber=size(A1,1) %数组 A 的维数即电网 AC 线条数%建立 Acline 结构体for i=1:aclinenumberAcline(i).fbname=A1(i); %AC 线始节点名Acline(i).tbname=A2(i); %AC 线终节点名Acline(i).Base_KV
39、=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:busnumberif strcmp(Acline(i).fbname, Bus(k).name)Acline(i).fbno=Bus(k).no; %AC 线始节点编号endif strcmp(Acline(i).tbname, Bus(k).name)Acline(i).tbno=Bus(k).no; %AC 线终节点编号endend
40、endfid=fopen(Transdata.txt); %打开 Transdata.txt 文件T=textscan(fid, %s %f %f %s %f %f %f %f); %输入变压器数据华南理工大学广州学院电气工程学院电力系统分析课程设计报告15fclose(fid);transnumber=size(T1,1) %数组 T 的维数即电网变压器个数（tansnumber）用变量名 tansnumber 表示if transnumber=0Trans(i).fbname=bus0;elsefor i=1:transnumberTrans(i).fbname=T1(i); %变压器始节
41、点名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); %变压器终节点额定电压Trans(i).R=T7(i); %变压器线串联 R（标幺值）Trans(i).X=T8(i); %变压器线串联 X（标幺值）for k=1:busnumberif strcmp(Trans(i).fbname, Bus(k).na
42、me)Trans(i).fbno=Bus(k).no; %变压器始节点编号endif strcmp(Trans(i).tbname, Bus(k).name)Trans(i).tbno=Bus(k).no; %变压器终节点编号endendTrans(i).k=Trans(i).tbrated_KV*Trans(i).fbBase_KV/Trans(i).fbrated_KV/Trans(i).tbBase_KV; %计算变比tempx=Trans(i).fbrated_KV2/Trans(i).fbBase_KV2;Trans(i).X=tempx*Trans(i).X; %折算后的变压器线串
43、联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:节点导纳矩阵 Y 的虚部 ,B2 对应只用支路电抗,变压器电抗形成的 B 矩阵 JP,JQ=FormJPQmatrix(Bus,B,B2,busnumber);B2 %显示 B2 矩阵JPiJP=-inv(JP); B %显示 B 矩阵JQ华南理工大学广州学院电气工程学
44、院电力系统分析课程设计报告16iJQ=-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);if transnumber0Sf_T %输出变压器始节点潮流St_T %输出变压器终节点潮流endSf_L %输出输电线始节点潮流St_L %输出输电线终节点潮流clear5.2
45、子函数%生成 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:aclinenumberf=Acline(i).fbno;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
46、)+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/Acline(i).X;X(t,f)=1/Acline(i).X; endfor i=1:tansnumberf=Trans(i).fbno;t=Trans(i).tbno;Y(f,f)=Y(f,f)
47、+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;华南理工大学广州学院电气工程学院电力系统分析课程设计报告17X(f,f)=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)=
48、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:busnumberif Bus(i).type=1 %该循环将 JP 和 JQ 矩阵对应平衡节点的非对角元素置零,对角元素置 1(相当于去掉了对应平衡节点的 JP 和 JQ 矩阵的行和列)for k=1:busnumberJQ(i,k)=0.;JQ(k,i)=0.;JP(i,k)=0.;JP(k,i)=0.; endJQ(i,i)=1.;JP(i,i)=1.;endif Bus(i).type=3 %该循环将 JQ 矩阵对应 PV 节点的非对角元素置零,对角元素置 1(相当于去掉了对应 PV 节点 JQ 矩阵的行和列)for k=1:busnumberJQ(i,k)=0.;JQ(k,i)=0.; endJQ(i,i)=1.;end endend%解潮流函数function Bus=so
