摘要
电力系统是现代社会不可或缺的基础设施之一，它的正常运行关系到整个社会的发展和稳定。电力系统潮流计算是电力系统运行分析的基础，它能够计算出电力系统各节点的电流、电压、功率等参数，为电力系统的稳定运行提供重要的参考依据。本文基于牛顿-拉夫逊法和PQ分解法，利用matlab编程实现电力系统潮流上机计算。本文的设计意义在于探索不同的潮流计算方法，提高电力系统潮流计算的准确性和效率。对比了牛顿-拉夫逊法和PQ分解法两种潮流计算方法，分析了它们的特点。牛顿-拉夫逊法适用于复杂的电力系统，可以处理任意类型的负荷和发电机，但需要较长的计算时间和较大的存储器。PQ分解法则适用于较小的电力系统，计算时间和存储空间要求较低，但对于非线性系统和不平衡负荷处理不够灵活。本文以陈珩教授的《电力系统稳态分析（第四板）》中的例4-4为例，解析了题目的给定条件，并通过矩阵计算得到了系统的节点电压和各支路潮流。程序设计中，本文给出了详细的程序流程图，并对结果进行了分析。通过对比牛顿-拉夫逊法和PQ分解法的结果，发现两种方法得出的潮流值基本一致，但牛顿-拉夫逊法耗时更长。
综上，本文基于牛顿-拉夫逊法和PQ分解法，利用matlab编程实现了电力系统潮流上机计算。本文的研究对电力系统的潮流计算具有一定的参考价值，为电力系统潮流计算提供了新的思路和方法。希望通过该研究，可以为电力系统的运行和管理提供更科学、更有效的决策支持，进而促进电力系统的发展和优化。
关键词：电力系统潮流上级计算；MATLAB；牛顿-拉夫逊法；PQ分解法
Power system power flow computer calculation
Abstract
Author's name: Zhang Yuchao Zip code: 071000
The power system is one of the indispensable infrastructures of modern society, and its normal operation is related to the development and stability of the whole society. Power system power flow calculation is the basis of power system operation analysis. It can calculate the current, voltage, power and other parameters of each node of the power system, and provide an important reference for the stable operation of the power system. Based on the Newton-Raphson method and PQ decomposition method, this paper uses matlab programming to realize the power system power flow calculation on the computer. The design significance of this paper is to explore different power flow calculation methods to improve the accuracy and efficiency of power system power flow calculation. Two power flow calculation methods, Newton-Raphson method and PQ decomposition method are compared, and their characteristics are analyzed. The Newton-Raphson method is suitable for complex power systems and can handle any type of loads and generators, but requires a long calculation time and a large memory. The PQ decomposition method is suitable for smaller power systems, and requires less computing time and storage space, but it is not flexible enough for nonlinear systems and unbalanced loads. This paper takes example 4-4 in Professor Chen Heng's "Steady State Analysis of Electric Power System (Fourth Board)" as an example, analyzes the given conditions of the topic, and obtains the node voltage of the system and the power flow of each branch through matrix calculation . In program design, this paper gives a detailed program flow chart and analyzes the results. By comparing the results of the Newton-Raphson method and the PQ decomposition method, it is found that the power flow values obtained by the two methods are basically the same, but the Newton-Raphson method takes longer.
In summary, based on the Newton-Raphson method and PQ decomposition method, this paper realizes the power system power flow calculation on the computer by using matlab programming. The research in this paper has a certain reference value for the power flow calculation of the power system, and provides a new idea and method for the power system power flow calculation. It is hoped that through this research, more scientific and effective decision support can be provided for the operation and management of the power system, and then the development and optimization of the power system can be promoted.
Key words: superordinate calculation of power system power flow; MATLAB; Newton-Raphson method; PQ decomposition method
目 录
三、题目解析（题目电气接线图，题目的给定条件，题目的矩阵计算等）. 8
一、概述
1.1 设计背景
电力系统是指由发电机、输电线路、变电站、配电网和终端用户组成的电力生产、输送、分配和使用系统。电力系统的安全运行和合理经济的运行是保障人们生活和经济发展的重要保障之一。为了保障电力系统的安全运行和合理经济的运行，电力系统潮流计算技术应运而生。
电力系统潮流计算是指在给定系统负荷和发电机输出的条件下，计算各节点的电压、功率、线路电流等电气量的一种电力系统分析方法。计算的结果可以用于电力系统的规划、设计、运行与管理等方面，为电力系统的安全、稳定和经济运行提供重要支持。
我国电力系统潮流计算技术的发展始于20世纪60年代末期，迄今已有50多年的历史。在改革开放以来的30多年间，我国电力系统取得了快速发展，电网规模和复杂度不断提高，电力负荷和电力需求不断增长，电力市场化程度不断提高，电网安全稳定运行的要求也越来越高。在这样一个背景下，电力系统潮流计算技术得到了广泛的应用和发展。目前，我国电力系统潮流计算技术已经进入到了智能化、高效化、多功能化的阶段。深度融合了计算机技术、通信技术、控制技术、数据挖掘技术、人工智能等新技术，实现了电力系统智能化运行、自动化控制、多功能数据管理和分析等功能，推动了电力系统的科技进步和发展。未来，随着政策、技术、市场等多重因素的影响，我国电力系统潮流计算技术将会继续向着智能化、高效化、协同化的方向发展。
1.2 设计意义
电力系统潮流计算是电网计算的核心内容之一，它主要是研究电力系统运行状态下各节点的电压相角以及电流大小，从而确定电力系统的电力负荷分配和电力功率分配。潮流计算是电力系统设计、运行和故障分析的基础，对于确保电力系统的安全、稳定、经济运行具有重要意义。
1.2.1帮助电力系统规划处实现电力系统规划的最优化
通过潮流计算，可以确定电力系统的潮流分布、电压大小和电流距离等参数，从而为电力系统规划提供准确的数据和参考。通过对潮流分布的分析，可以确定电网容量、电力负荷分配和电压控制方案，进而实现电力系统规划的最优化。
1.2.2帮助电力系统运行部门实现电力系统的稳定运行
随着电力系统规模的扩大和负荷不断增加，电力系统的安全稳定性越来越受到关注。通过潮流计算，可以预测电力系统各节点的电压和电流变化情况，提前发现电力系统可能存在的问题，并采取相应的措施进行调整和优化。
1.2.3帮助电力系统故障分析和排除
在电力系统运行过程中，可能会出现各种故障，如线路短路、设备故障等。通过潮流计算，可以快速准确地定位故障点，并及时采取措施进行排除，确保电力系统的正常运行。
1.3设计要求
（1）采用matlab或者C语言编程计算
（2）输入原始的数据单独编写一个输入文件
（3）输出各个节点的迭代次数的数据和迭代次数
（4）自己设计编程结果分析曲线
（5）文献查阅，根据关键词在中国知网查阅
（6）附录中放完整的程序
（7）按照报告格式完成各个部分内容
（8）用visio软件画网络接线图、等值网络图、流程图等
二、计算方法
2.1牛顿拉夫逊法潮流计算
牛顿-拉夫逊法是求解潮流的最常用的方法。其核心在于修正方程的建立及求解。注意的是，修正方程的雅各比矩阵不是对称矩阵，但是稀疏矩阵；由于雅各比矩阵的元素与电压大小和相位有关，因此在每次迭代过程中都要重新形成雅各比矩阵，这是限制牛顿-拉夫逊法速度的最大因素。牛顿-拉夫逊法的收敛速度比高斯-塞德尔法快很多。
2.2 PQ分解法潮流计算
P-Q分解法由牛顿-拉夫逊法的节点电压以极坐标表示时发展而来。主要是根据电力网络的特性对牛-拉法的雅各比矩阵进行简化，变成常系数矩阵，因此在每次迭代过程中都不用重新形成系数矩阵。而且P-Q分解法的系数矩阵阶数较牛-拉法的低，还是对称矩阵。因此其收敛速度较牛-拉法快（其迭代次数比牛-拉法多，但其每次迭代的耗时少）。虽然P-Q分解法是在一定简化的基础上发展得到的，但由于其功率不平衡量的求解与牛-拉法完全一样（即P-Q分解法只对雅各比矩阵简化，不对功率不平衡量简化），而且收敛要求都一样的，因此最终得到的结果跟牛-拉法完全一样。注意的是在运用P-Q分解法时是有限制的，必须在电力网络符合简化要求情况下才能运用。相比而言，牛-拉法没有限制。
三、题目解析（题目电气接线图，题目的给定条件，题目的矩阵计算等）
3.1题目接线图与给定条件
发电厂F母线Ⅱ上所连发电机发出给定功率40+j30MVA，其余功率由母线Ⅰ所连发电机供给。连接母线Ⅰ、Ⅱ的联络变压器容量为60MVA，R_T=3Ω，X_T=110Ω；线路末端降压变压器总容量为240MVA，R_T=0.8Ω，X_T=23Ω；220kV线路，R_l=5.9Ω，X_l=31.5Ω；110kV线路，xb段，R_l=65Ω，X_l=100Ω；bⅡ段，R_l=65Ω，X_l=100Ω。所有阻抗均按线路额定电压的比值归算至220kV侧。
图3.1 网络接线图
依据题中给出的阻抗数据，我们可以将网络接线图改绘为图2 以阻抗表示的等值电路，并将相关线路的阻抗值标注在图上，需要注意的是，由于电压等级的不同，我们需要进行电阻的归算；节点导纳矩阵的形成还需考虑对地支路对自阻抗、互阻抗的影响。
3.2节点导纳矩阵的形成
以阻抗表达的等值电路如下：
图3.2 以阻抗表示的等值电路
导纳的计算
y_12 = 1/z_12 = 1/(5.9+j31.5) = 0.005745-j0.030670 (S)
y_23 = 1/z_23 = 1/(0.8+j23) = 0.001510-j0.043426 (S)
y_34 = 1/z_34 = 1/(65+j100) = 0.004569-j0.007030 (S)
y_45 = 1/z_45 = 1/(65+j100) = 0.004569-j0.007030 (S)
y_51 = 1/z_51 = 1/(3+j110) = 0.000248-j0.009084 (S)
变压器的修正
k_(1*)=(U_ⅡN U_Ⅰ)/(U_ⅠN U_Ⅱ ) = (110×231)/(220×110) = 1.050000
k_(2*)=(U_ⅡN U_Ⅰ)/(U_ⅠN U_Ⅱ ) = (110×231)/(220×121) = 0.954545
由此可得以导纳和理想变压器表示的等值线路如下：
图3.3 以导纳和理想变压器表示的等值网络
由于变压器变比不匹配，事实上无法按照实际变比归算网络参数[3]。可选用变压器Π型等值电路，它的优点是不必进行参数和变量的归算，变压器Π型电路如下：
将变压器改为Π型后等值电路如下：
图3.4 变压器以Π型等值电路表示时的等值网络
此时节点导纳参数如下：
y_10 = Y_T ((1-k_(1*)))/(k_(1*)^2 ) = -0.000011+j0.000412
y_50 = Y_T ((k_(1*)-1))/k_(1*) =0.000012-j0.000433
y_15 = Y_T/k_(1*) = 0.000236-j0.008651
y_20 = Y_T (1ⓜ┤-ⓜk_(2*) )/(k_(2*)^2 ) = 0.000075-j0.002166
y_30 = Y_T (k_(2*)ⓜ┤-ⓜ1)/k_(2*) = -0.000072+j0.002068
y_23 = Y_T/k_(2*) = 0.001582-j0.045494
由此，节点导纳矩阵各元素可求取如下：
Y_11 = y_10+y_12+y_15 = 0.005970-j0.038909
Y_12 = Y_21 = -y_12 = 0.005745+j0.030670
Y_15 = Y_51 = -y_15 = -0.000236+j0.008651
以此类推，形成节点导纳矩阵如下：
Y_B 0.005967-j0.038909 -0.005745+j0.03067 0 0 -0.000236+j0.008651
-0.005745+j0.03067 0.007402-j0.07833 -0.001582+j0.045494 0 0
0 -0.001582+j0.045494 0.006079-j0.050456 -0.004569+j0.00703 0
0 0 -0.004569+j0.00703 0.009138-j0.01406 -0.004569+j0.00703
-0.000236+j0.008651 0 0 -0.004569+j0.00703 0.004817-j0.016114
节点导纳矩阵有如下特点：
①节点导纳矩阵是方阵，它的阶数等于网络中除参考节点外的节点数n。
②节点导纳矩阵是稀疏矩阵；由于网络的互易特性，一般情况下该矩阵为对称矩阵
③节点导纳矩阵的对角元元素等于该节点所连接导纳的总和。
④节点导纳矩阵的非对角元等于连接节点 支路导纳的负值。
注：可利用编程实现节点导纳矩阵的求解；这一部分为简单网络人工求解提供思路。
3.3牛顿-拉夫逊法原理
由于本题节点导纳矩阵具有明显的对称性、稀疏性、初值确定等特点，采用牛顿-拉夫逊法求解较为简便；
牛顿-拉夫逊法是常用的解非线性方程组的方法，也是当前广泛采用的计算潮流的方法。下面简述牛顿拉夫逊法原理：
设有非线性方程组
(3-1)
其近似解为 ，设近似解与精确解分别相差 ，则有如下关系成立
(3-2)
上式均可按照泰勒级数展开。以第一式为例：
(3-3)
由于近似解 与精确解相差不大，故而可以略去 的高次方，从而麦克劳林余式 也可略去，将方程组改写为矩阵方程，整理如下：
(3-4)
上式简写为
(3-5)
式中 称为雅可比矩阵； 为由 组成的列向量； 称为不平衡量的列向量。
将 代入，可得 中的各元素。然后解线性方程组，求得 ，从而得到第一次迭代后的新值 ；
再将求得的 代入，可求得 中各元素的新值。求解下一个线性方程组得 ，从而得到第二次迭代后的新值 ；
以此循环，最终得到的 已无限接近精确解，迭代退出的条件也就是迭代收敛条件满足 。显然，由于泰勒展开略去了高次项，这意味着迭代收敛取决于 的初值接近精确解[4]。
3.4牛顿-拉夫逊法在潮流计算的应用
由电力系统稳态分析中可得网络的功率方程为：
(3-6)
令 ，令 ，从而将上式改写为：
(3-7)
将实部与虚部分列，可得：
(3-8)
建立类似(3-4)的修正方程如下：
(3-9)
式中的 分别为注入功率和节点电压平方的不平衡量，它们分别为
(3-10)
式中雅可比矩阵各个元素分别为：
(3-11)
由此循环迭代求得电压幅值、相角、线路功率、输入输出功率的最终迭代结果。
3.5PQ分解法在潮流计算的应用
1. 对于P节点（平衡节点），潮流方程为：
∑(V_i * (G_ij * cosθ_ij + B_ij * sinθ_ij)) - P_i + P_si = 0
V_i表示节点i的电压幅值；G_ij和B_ij分别为节点i和节点j之间的导纳的实部和虚部；θ_ij表示节点i和节点j之间的相角差；P_i表示节点i的实际注入功率；P_si表示节点i的静态注入功率。
2. 对于Q节点（非平衡节点），潮流方程为：
∑(V_i * (G_ij * sinθ_ij - B_ij * cosθ_ij)) - Q_i + Q_si = 0
Q_i表示节点的无功注入功率；Q_si表示节点的静态无功注入功率。
3.考虑节点的电压平衡条件，即：
V_i^2 = V_i^A + V_i^Q
V_i^A为节点i的有功注入功率引起的电压，V_i^Q为节点i的无功注入功率引起的电压。
过对以上方程进行迭代，可以逐步求解得到系统中各个节点的电压和相角，从而获得节点的潮流信息。
四、程序设计（包括程序流程图，结果分析等）
4.1牛顿拉夫逊法潮流计算
4.1.1程序流程图
|
结束 |
|
输入原始数据 |
|
给定初始电压幅值、相角 |
|
K=0 |
|
计算功率偏移量 |
|
求雅可比矩阵各元素的值 |
|
求解修正方程得修正量的值 |
|
对各节点电压进行修正 |
|
|
|
计算节点功率及全部线路功率 |
|
开始 |
|
K=K+1 |
图4.1 牛顿拉夫逊法潮流计算程序流程图
4.1.2结果分析
（1）输入
表4.1 节点参数
| 节点编号 | 节点电压 | 节点相角（弧度制） | 有功注入 | 无功注入 | 节点类型 |
| 1 | 242 | 0.00 | 0 | 0 | 3 |
| 2 | 220 | 0.00 | 0 | 10 | 1 |
| 3 | 220 | 0.00 | -180 | -100 | 1 |
| 4 | 220 | 0.00 | -50 | -30 | 1 |
| 5 | 220 | 0.00 | 40 | 30 | 1 |
表4.2 支路参数
| 节点i | 节点j | 线路电阻 | 电抗 | 电导 | 电纳 | 变压器变比（普通线路为零） |
| 1 | 2 | 5.9 | 31.5 | 0 | 0 | 0 |
| 3 | 2 | 0.8 | 23 | 0 | 0 | 0.954545 |
| 3 | 4 | 65 | 100 | 0 | 0 | 0 |
| 4 | 5 | 65 | 100 | 0 | 0 | 0 |
| 5 | 1 | 3 | 110 | 0 | 0 | 1.05 |
（2）输出
---------------------- 第1次迭代算，U最大误差为18.4675047132，角度最大误差为0.1671410078 -----------------------
表4.3 第1次迭代节点计算结果
| 节点计算结果： |
| 节点 节点电压 节点相角（角度） 节点注入功率 |
| 1 242.000000 0.000000 196.084682 + j 113.830126 |
| 2 221.547359 5.175770 0.000000 + j 10.000000 |
| 3 219.451852 9.576474 180.000000 + j100.000000 |
| 4 214.755133 7.655453 50.000000 + j 30.000000 |
| 5 238.467505 2.123851 40.000000 + j 30.000000 |
表4.4 第1次迭代线路计算结果
| 线路计算结果： |
| 节点i 节点j 线路功率Sij 线路功率Sji 线路损耗△Sij |
| 1 2 178.028659 + j 130.723039 173.114069 + j104.484128 4.914590 + j 26.238911 |
| 3 2 173.683224 + j108.085492 174.378392 + j 128.071568 0.695168 + j 19.986076 |
| 3 4 6.275278 + j 14.650836 6.618135 + j 14.123364 0.342857 + j 0.527473 |
| 4 5 56.882914 + j 11.564787 61.631673 + j 18.870570 4.748759 + j 7.305783 |
| 5 1 18.021494 + j 18.158987 18.056023 + j 16.892913 0.034529 + j 1.266074 |
---------------------- 第2次迭代算，U最大误差为6.1929384291，角度最大误差为0.0048310350 -----------------------
表4.5 第2次迭代节点计算结果
| 节点计算结果： |
| 节点 节点电压 节点相角（角度） 节点注入功率 |
| 1 242.000000 0.000000 201.788173 + j 151.690634 |
| 2 217.798165 5.202103 0.000000 + j 10.000000 |
| 3 214.356953 9.805041 180.000000 + j100.000000 |
| 4 208.562194 7.821131 50.000000 + j 30.000000 |
| 5 233.826862 2.400649 40.000000 + j 30.000000 |
表4.6 第2次迭代线路计算结果
| 线路计算结果： |
| 节点i 节点j 线路功率Sij 线路功率Sji 线路损耗△Sij |
| 1 2 181.461646 + j 158.835470 175.602645 + j127.554365 5.859001 + j 31.281105 |
| 3 2 174.680302 + j115.797851 175.445019 + j 137.783448 0.764716 + j 21.985597 |
| 3 4 5.081778 + j 15.992601 5.480116 + j 15.379773 0.398338 + j 0.612828 |
| 4 5 55.466037 + j 14.458867 60.375669 + j 22.012148 4.909632 + j 7.553280 |
| 5 1 20.300310 + j 8.106141 20.326527 + j 7.144836 0.026217 + j 0.961305 |
---------------------- 第3次迭代算，U最大误差为0.1839263364，角度最大误差为0.0002489062 -----------------------
表4.7 第3次迭代节点计算结果
| 节点计算结果： |
| 节点 节点电压 节点相角（角度） 节点注入功率 |
| 1 242.000000 0.000000 202.005528 + j 152.675770 |
| 2 217.698490 5.204339 0.000000 + j 10.000000 |
| 3 214.211734 9.819303 180.000000 + j100.000000 |
| 4 208.378268 7.830038 50.000000 + j 30.000000 |
| 5 233.711980 2.410727 40.000000 + j 30.000000 |
表4.8 第3次迭代线路计算结果
| 线路计算结果： |
| 节点i 节点j 线路功率Sij 线路功率Sji 线路损耗△Sij |
| 1 2 181.596428 + j 159.578749 175.708655 + j128.144025 5.887774 + j 31.434724 |
| 3 2 174.939910 + j116.053538 175.708279 + j 138.144162 0.768370 + j 22.090624 |
| 3 4 5.059613 + j 16.053722 5.460947 + j 15.436285 0.401335 + j 0.617438 |
| 4 5 55.460918 + j 14.563534 60.382919 + j 22.135842 4.922000 + j 7.572308 |
| 5 1 20.382884 + j 7.864213 20.409100 + j 6.902979 0.026215 + j 0.961234 |
采用牛顿-拉夫逊迭代法，总共迭代3次
根据给定的牛顿-拉夫逊迭代法潮流计算的结果，可以知道：
（1）节点电压误差：第1次迭代时，最大节点电压误差为18.4675047132，第2次迭代时为6.1929384291，第3次迭代时为0.1839263364。可以看出，随着迭代次数的增加，节点电压误差逐渐减小，说明迭代法逐渐趋向于收敛。
（2）节点相角误差：第1次迭代时，最大节点相角误差为0.1671410078，第2次迭代时为0.0048310350，第3次迭代时为0.0002489062。与节点电压误差类似，节点相角误差也逐渐减小，表明迭代法也在相角方向上逐渐收敛。
（3）节点注入功率：节点1注入功率逐渐增加，从第1次迭代的196.084682 + j 113.830126到第2次迭代的201.788173 + j 151.690634，最后到第3次迭代的202.005528 + j 152.675770。其他节点的注入功率基本保持稳定。
（4）线路功率：线路功率损耗也逐渐减小，第1次迭代时总线路功率损耗为10.089701 + j 55.312307，第2次迭代时为10.894090 + j 59.804637，第3次迭代时为11.005694 + j 60.680020。这表明迭代法能够准确计算线路功率分布。
综上所述，通过牛顿-拉夫逊迭代法进行潮流计算的结果显示，随着迭代次数的增加，节点电压和节点相角的误差逐渐减小，收敛性逐渐提高。同时，节点注入功率和线路功率的计算结果也逐渐趋于稳定。这说明迭代法能够有效解决潮流计算问题，得到准确的结果。
4.2PQ分解法潮流计算
4.2.1程序流程图
图4.2 PQ分解法潮流计算程序流程图
4.2.2结果分析
表4.9节点参数
| 节点 | 类型 | 发电机有功 | 发电机无功 | 负荷有功 | 负荷无功 | 电压幅值 | 电压相位 |
| 1 | 1 | 0 | 0 | 0 | 0 | 242 | 0 |
| 2 | 3 | 0 | 0 | 0 | 10 | 220 | 0 |
| 3 | 3 | 0 | 0 | -180 | -100 | 220 | 0 |
| 4 | 3 | 0 | 0 | -50 | -30 | 220 | 0 |
| 5 | 3 | 0 | 0 | 40 | 30 | 220 | 0 |
表4.10支路参数
| 入端 | 出端 | 电阻 | 电抗 | 电纳 | 有无调压变压器 | 变压器变比 |
| 1 | 2 | 5.9 | 31.5 | 0 | 0 | 0 |
| 3 | 2 | 0.8 | 23 | 0 | 0 | 0 |
| 3 | 4 | 65 | 100 | 0 | 0 | 0 |
| 4 | 5 | 65 | 100 | 0 | 0 | 0 |
| 5 | 1 | 3 | 110 | 0 | 0 | 0 |
表4.11 计算结果1
| 节点 | 类型 | 发电机有功 | 发电机无功 | 负荷有功 | 负荷无功 | 电压幅值 | 电压相位 |
| 1 | 1 | -176.7979027 | 254.6749786 | 0 | 0 | 242 | 0 |
| 2 | 3 | 0 | 0 | 0 | 10 | 220.7693827 | 0.120714478 |
| 3 | 3 | 0 | 0 | -180 | -100 | 220.9918527 | 0.202366266 |
| 4 | 3 | 0 | 0 | -50 | -30 | 221.0136374 | 0.18238969 |
| 5 | 3 | 0 | 0 | 40 | 30 | 220.2998793 | 0.025876774 |
表4.12 计算结果2
| 入端 | 出端 | 入端有功 | 入端无功 | 入段充电 | 出端有功 | 出端无功 | 出端充电 |
| 1 | 2 | -165.5730369 | 206.4597036 | 0 | -172.629 | 168.7871 | 0 |
| 3 | 2 | 173.1203522 | 3.183190009 | 0 | 172.6292 | -10.9363 | 0 |
| 3 | 4 | 6.881127909 | -4.423422941 | 0 | 6.792065 | -4.56044 | 0 |
| 4 | 5 | 56.79299342 | -29.3865728 | 0 | 51.3518 | -37.7576 | 0 |
| 5 | 1 | 11.35075676 | -43.60670897 | 0 | 11.22525 | -48.2087 | 0 |
通过PQ分解法潮流计算得到的结果，表4.5中节点1为发电机节点，其有功功率为负值，意味着该节点向系统中注入有功功率，而节点2-5均为负荷节点，需要从系统中吸收有功功率。此外，节点2-5均为无功负载节点，需要从系统中吸收无功功率。
从表4.5中可以看出，节点1的电压幅值为242V，高于220V的额定电压值，表明该节点存在电压过高的问题。而节点2-5的电压幅值基本达到了220V，符合电力系统的额定电压范围。节点1的电压相位为0，符号标志着电压相位为正序，而节点2-5的电压相位均为正值，说明电压相位处于正序相位。
表4.6中的结果可以用于计算电力系统的输电线路的潮流分布情况。从表4.6中可以看出，节点1和5之间的电力传输存在明显的潮流过载问题，节点1向节点2的潮流值为-165.57MW，节点5向节点4的潮流值为56.79MW。这说明电力系统需要进行调整，以增强输电线路的容量，以确保潮流过载并不会对电力系统的运行造成不利影响。
综上，通过PQ分解法进行电力系统潮流计算可以更好地了解电力系统的运行状况，发现问题并及时进行调整。但需要注意的是，在进行潮流计算时需要准确提供节点参数和电路参数，否则计算结果可能会出现偏差。
五、总结
本文基于牛顿-拉夫逊法和PQ分解法，对电力系统进行潮流计算，并利用MATLAB进行编程。通过对比分析两种计算方法的特点，我们得出了以下结论：
（1）牛顿-拉夫逊法适用于复杂的系统，其计算过程需要迭代多次，但能够保证计算结果的精度。
（2）PQ分解法适用于简单的系统，其计算过程直接，能够快速地计算出潮流结果，但其精度相对较低。
在陈珩教授的《电力系统稳态分析（第四板）》中的例4-4题目的解析中，我们根据给出的电气接线图和给定条件，设计了相应的矩阵计算方法，并使用MATLAB编程进行计算。通过分析计算结果，我们得出了以下结论：
（1）在牛顿-拉夫逊法和PQ分解法中，得出的潮流计算结果相对一致，但两种方法所需要的计算时间和精度略有不同。
（2）在给定的系统中，各支路的潮流分布相对均匀，系统总负荷得到了合理分配。
（3）在设计电力系统时，需要综合考虑系统的复杂程度、计算精度和计算效率，选择合适的计算方法进行潮流计算。
综上所述，本文基于牛顿-拉夫逊法和PQ分解法对电力系统进行了潮流计算，并对两种方法进行了比较分析。通过具体的题目解析和程序设计，我们得出了相应的结论和体会。在电力系统设计和维护中，潮流计算是非常重要的环节，需要细心认真地进行计算，以保证系统运行的稳定和可靠。同时，经过电力系统潮流计算上机实验课的培养训练后，我掌握了很多潮流计算的相关知识和论文的撰写方法，这其中包括：
对潮流计算理论到实际处理的进一步理解，加深了牛顿拉夫逊迭代法对解非齐次方程组的处理方式，理解了迭代收敛的前提条件；对功率流入流出关系的进一步认识，例将变压器励磁支路等效为一个无功功率等，加深理解了平衡节点的电压幅值固定、相角为参考相位的认识，平衡节点的存在简化了迭代步骤；编程时的简化，由于例题中只含有平衡节点和PQ节点，可以省略PV节点的寻找与运算，简化了编程；通过对电力系统相关学位论文的阅读，初步掌握了本专业相关论文的撰写方法，尤其是学会了使用AxMath进行公式的撰写，使得文本格式更具有可读性。
六、参考文献
[1]陈珩.电力系统稳态分析（第四板）[M].北京：中国电力出版社，2015.10
[2]李庚银.电力系统分析基础.北京：机械工业出版社，2011.8
[3] 何仰赞，温增银.电力系统分析.下（第四版）.武汉：华中科技大学出版社，2016.5
[4]曹井川.基于MATLAB的科研型潮流计算研究与设计[D].大连：大连海事大学，2017.
[5]蔡超豪.MATLAB在电力系统中的应用[M].北京：中国电力出版社，2022
[6]王增光. 科研目的使用的直角坐标牛顿法潮流计算研究与设计[D].大连海事大学,2019
[7]王震东,刘白杨,张随涵.基于Matlab的牛顿-拉夫逊法电力系统潮流计算[J].电工电气,2016(07):61-62.
[8]刘阳涵. 快速PQ分解法潮流计算方法研究[D].南昌大学,2016.
[9]李俊,肖宏,唐荣平,龙江.一种潮流计算的PQ改进算法研究[J].大众科技,2014,16(10):103-105.
[10]徐劲松,宁玉琳,杨永锋.基于Matlab的电力系统PQ分解法潮流计算研究[J].电气传动自动化,2011,33(02):10-18.
[11]罗杰.基于MATLAB的牛顿拉夫逊法电力潮流计算与实现[J].科技广场,2010(03):183-184.
七、附录
7.1牛顿拉夫逊法代码
%这段代码主要是实现了牛顿拉夫逊法潮流计算的过程。
% 首先从文件读取节点参数和支路参数，并进行一些初始化工作。
% 然后根据节点的类型统计了平衡节点、PV节点和PQ节点的个数，
% 并将这些节点分别存储在不同的数组中。接着重新排序节点，
% 并生成新旧节点对照表。对支路矩阵进行重新编号。
% 接下来，根据支路参数计算导纳矩阵Y，并根据支路类型进行不同的计算。
% 最后，对角线元素加入节点的功率注入信息。完成了迭代计算。
clc;clear;close all;%clc; % 清空命令行窗口 清除工作区的变量 % 关闭所有图形窗口
jj = sqrt(-1); %虚数单位
format long% 设置输出格式为长格式
Err = 1.0; % 设定误差精度1.0e-10
Max = 10; % 最大迭代次数
fileID = fopen('节点参数.txt','r');
formatSpec = '%d %f %f %f %f %d';% 定义读取格式
bus_ = (fscanf(fileID,formatSpec,[6 Inf]))';% 读取节点参数并进行转置操作
fileID = fopen('支路参数.txt','r');
formatSpec = '%d %d %f %f %f %f %f';
line_ = (fscanf(fileID,formatSpec,[7 Inf]))';% 读取支路参数并进行转置操作
delete '计算结果.txt'
Result_file = fopen('计算结果.txt','a');
% 设定最大迭代次数为Max，以便不收敛情况下及时跳出，显示每次迭代的结果
for T = 1:Max
bus=bus_;% 将初始节点参数赋值给当前节点参数
line=line_;% 将初始支路参数赋值给当前支路参数
[row_bus,column_bus]=size(bus); %bus的行数nb和列数mb,
[row_line,column_line]=size(line); %size的行数nl和列数ml
%统计三种节点的个数
nSW = 0; % nSW为平衡节点个数
nPV = 0; % nPV为PV节点个数
nPQ = 0; % nPQ为PQ节点个数
for i = 1:row_bus % row_bus为总节点数
type= bus(i,6);
if type == 3
nSW = nSW + 1;
elseif type == 2
nPV = nPV +1;
else
nPQ = nPQ + 1;
end
end
SW=zeros(nSW,6);% 创建空数组用于存储平衡节点
PV=zeros(nPV,6);% 创建空数组用于存储PV节点
PQ=zeros(nPQ,6);% 创建空数组用于存储PQ节点
oldnode=bus(:,1); % 保存旧的节点编号信息
%重新排序
nSW = 0;
nPV = 0;
nPQ = 0;
for i = 1:row_bus % row_bus为总节点数
type= bus(i,6);
if type == 3
nSW = nSW + 1;
SW(nSW,:)=bus(i,:); % 计算并储存平衡节点
elseif type == 2
nPV = nPV +1;
PV(nPV,:)=bus(i,:); % 计算并储存PV节点
else
nPQ = nPQ + 1;
PQ(nPQ,:)=bus(i,:);% 计算并储存PQ节点
end
end
bus=[PQ;PV;SW];% 将排序后的节点重新组合
newnode=bus(:,1);% 获取新的节点编号信息
nodenum=[oldnode newnode];% 生成新旧节点对照表
% 按排序以后的节点顺序对line矩阵重新编号
for i=1:row_line
for j=1:2
for k=1:row_bus
if line(i,j)==nodenum(k,2)
line(i,j)=nodenum(k,1);
break
end
end
end
end
Y=zeros(row_bus,row_bus); % 对导纳矩阵赋初值0
for t=1:row_line
i=line(t,1);
j=line(t,2);
Zt=line(t,3)+jj*line(t,4); % 线路阻抗参数
if Zt~=0
Yt=1/Zt; % 接地支路不计算Yt
end
Ym=line(t,5)+jj*line(t,6); % 线路导纳参数
k=line(t,7);
if (k==0)&&(j~=0) % 普通线路: K=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);
end
if (k==0)&&(j==0) % 对地支路: K=0,J=0,Z=0
Y(i,i)=Y(i,i)+Ym;
end
if k>0 % 变压器线路: Zt和Ym为折算到i侧的值,k在j侧
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);
end
if k<0 % 变压器线路: Zt和Ym为折算到i侧的值,k在i侧
Y(i,i)=Y(i,i)+Yt+Ym;
Y(j,j)=Y(j,j)+k*k*Yt;
Y(i,j)=Y(i,j)+k*Yt;
Y(j,i)=Y(i,j);
end
end
n = nPQ + nPV +1; % 总节点个数n = nPQ + nPV +1;
for times = 1:T % 开始迭代计算，迭代T次
% 对dP和dQ赋初值 PV节点不需计算dQ 平衡节点不参与计算
dP = bus(1:n-1,4);
dQ = bus(1:nPQ,5);
for i = 1:n-1
for j = 1:n
%△Pi=Pi-∑UiUj[Gijcos(∠i-∠j)+Bijsin(∠i-∠j)]，j从1到n
%△Qi=Qi-∑UiUj[GijSIN(∠i-∠j)-Bijcos(∠i-∠j)]，j从1到n
%其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
dP(i,1) = dP(i,1) - bus(i,2)*bus(j,2) * (real(Y(i,j))*cos(bus(i,3)-bus(j,3)) + imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
if i<=nPQ
dQ(i,1) = dQ(i,1) - bus(i,2)*bus(j,2) * (real(Y(i,j))*sin(bus(i,3)-bus(j,3)) - imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
end
end
end
% 将雅克比矩阵分块，即：J = [H N; K L]，并初始化其中，K对应课本的J
H = zeros(row_bus-1,row_bus-1);
N = zeros(row_bus-1,nPQ);
K = zeros(nPQ,row_bus-1);
L = zeros(nPQ,nPQ);
% 计算Qi和Pi
Qi = zeros(row_bus-1,1);
Pi = zeros(row_bus-1,1);
for i = 1:row_bus-1
for j = 1:row_bus
% Pi=∑UiUj[Gijcos(∠i-∠j)+Bijsin(∠i-∠j)]，j从1到n
% Qi=∑UiUj[GijSIN(∠i-∠j)-Bijcos(∠i-∠j)]，j从1到n
% 其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
Pi(i,1)=Pi(i,1) + bus(i,2)*bus(j,2) * (real(Y(i,j))*cos(bus(i,3)-bus(j,3)) + imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
Qi(i,1)=Qi(i,1) + bus(i,2)*bus(j,2) * (real(Y(i,j))*sin(bus(i,3)-bus(j,3)) - imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
end
end
for i = 1:row_bus-1
for j = 1:row_bus-1
% 分别计算H矩阵的非对角及对角元素
% Hij=UiUj[Gijsin(∠i-∠j)-Bijcos(∠i-∠j)]
% Hii=-(Qi+UiUiBii)
% 其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
if i~=j
H(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3)) - imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
else
H(i,j)=-(Qi(i,1) + imag(Y(i,j))*((bus(i,2))^2));
end
if j <= nPQ
% 分别计算N矩阵的非对角及对角元素
% Nij=UiUj[Gijcos(∠i-∠j)+Bijsin(∠i-∠j)]
% Nii=Pi+UiUiGii
% 其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
if i~=j
N(i,j)=bus(i,2)*bus(j,2) * (real(Y(i,j))*cos(bus(i,3)-bus(j,3)) + imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
else
N(i,j)=Pi(i,1) + real(Y(i,j))*((bus(i,2))^2);
end
end
if i <= nPQ
% 分别计算K矩阵的非对角及对角元素
% Kij=-UiUj[Gijcos(∠i-∠j)+Bijsin(∠i-∠j)]
% Nii=Pi-UiUiGii
% 其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
if i~=j
K(i,j)=-bus(i,2)*bus(j,2) * (real(Y(i,j))*cos(bus(i,3)-bus(j,3)) + imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));
else
K(i,j)=Pi(i,1) - real(Y(i,j))*((bus(i,2))^2);
end
if j <= nPQ
% 分别计算L矩阵的非对角及对角元素
% Lij=UiUj[Gijsin(∠i-∠j)-Bijcos(∠i-∠j)]
% Lii=Qi-UiUiBii
% 其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
if i~=j
L(i,j)=bus(i,2)*bus(j,2) * (real(Y(i,j))*sin(bus(i,3)-bus(j,3)) - imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
else
L(i,j)=Qi(i,1) - imag(Y(i,j))*((bus(i,2))^2);
end
end
end
end
end
j = [H N; K L]; % 生成雅克比矩阵
% 生成电压对角矩阵
UD = zeros(nPQ,nPQ);
for i = 1:nPQ
UD(i,i) = bus(i,2);
end
dAngU = j \ [dP;dQ];
dAng = dAngU(1:row_bus-1,1); % 计算相角修正量
dU = UD*(dAngU(row_bus:row_bus+nPQ-1,1)); % 计算电压修正量
bus(1:nPQ,2) = bus(1:nPQ,2) + dU; % 修正电压
bus(1:row_bus-1,3) = bus(1:row_bus-1,3) + dAng; % 修正相角
% 判断是否满足精度误差，如满足则跳出，否则返回继续迭代
if (max(abs(dU))<Err && max(abs(dAng))<Err)
break
end
end
% Pi=∑UiUj[Gijcos(∠i-∠j)+Bijsin(∠i-∠j)]，j从1到n 且已知UiUi[Gijcos(∠i-∠i)+Bijsin(∠i-∠i)]=bus(i,4)
% Qi=∑UiUj[GijSIN(∠i-∠j)-Bijcos(∠i-∠j)]，j从1到n 且已知UiUi[Gijsin(∠i-∠i)-Bijcos(∠i-∠i)]=bus(i,5)
% 其中Gij=real(Y(i,j)，Bij=imag(Y(i,j)
% 计算PV节点的无功注入
for i = nPQ+1:n-1
for j = 1:n
bus(i,5)=bus(i,5) + bus(i,2)*bus(j,2) * (real(Y(i,j))*sin(bus(i,3)-bus(j,3)) - imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));
end
end
% 计算平衡节点的有功注入及无功注入
for j =1:n
bus(n,4)=bus(n,4) + bus(n,2)*bus(j,2) * (real(Y(n,j))*cos(bus(n,3)-bus(j,3)) + imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));
bus(n,5)=bus(n,5) + bus(n,2)*bus(j,2) * (real(Y(n,j))*sin(bus(n,3)-bus(j,3)) - imag(Y(n,j))*cos(bus(n,3)-bus(j,3)));
end
% 按排序以后的节点顺序对line矩阵重新编号
bus_temp = zeros(row_bus,column_bus); % 矩阵用于临时存放bus矩阵的数据
k = 1;
for i = 1 :row_bus
for j = 1 : row_bus
if bus(j,1) == k
bus_temp(k,:) = bus(j,:);
k = k + 1;
end
end
end % 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中
bus = bus_temp; % 重新赋值回bus，完成bus矩阵的重新编号
for i=1:row_line
for j=1:2
for k=1:row_bus
if line(i,j)==nodenum(k,1)
line(i,j)=nodenum(k,2);
break
end
end
end
end
YtYm = zeros(row_line,5);
YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt，i侧的等效Ym，j侧的等效Ym
for time=1:row_line
i=line(time,1);
j=line(time,2);
Zt=line(time,3)+jj*line(time,4);
if Zt~=0
Yt=1/Zt;
end
Ym=line(time,5)+jj*line(time,6);
k=line(time,7);
if (k==0)&&(j~=0) % 普通线路: K=0
YtYm(time,3) = Yt;
YtYm(time,4) = Ym;
YtYm(time,5) = Ym;
end
if (k==0)&&(j==0) % 对地支路: K=0,J=0,R=X=0
YtYm(time,4) = Ym;
end
if k>0 % 变压器线路: Zt和Ym为折算到i侧的值,k在j侧
YtYm(time,3) = Yt/k;
YtYm(time,4) = Ym+Yt*(k-1)/k;
YtYm(time,5) = Yt*(1-k)/k/k;
end
if k<0 % 变压器线路: Zt和Ym为折算到K侧的值,k在i侧
YtYm(time,3) = -Yt*k;
YtYm(time,4) = Ym+Yt*(1+k);
YtYm(time,5) = Yt*(k^2+k);
end
end
% bus_result矩阵储存着节点计算结果
bus_result = zeros(row_bus,4);
bus_result(:,1:2) = bus(:,1:2);
bus_result(:,3) = bus(:,3) *180 / pi; % 相角采用角度制
bus_result(:,4) = bus(:,4) + jj*bus(:,5); % 注入功率
% line_res矩阵储存着线路潮流计算结果
line_result = zeros(row_line,5);
line_result(:,1:2) = line(:,1:2); % 前两列为节点编号，后三列为Sij，Sji，△Sij
for k = 1:row_line
i = line_result(k,1);
j = line_result(k,2);
if (j~=0)&&(i~=0)
%Sij=UiUi(Yij*+Yi*)-UiUj*Yij* *表示共轭
%Sji=UjUj(Yij*+Yj*)-UjUi*Yij*
%△Sij=Sij+Sji
line_result(k,3)=bus(i,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4))) - bus(i,2)*(cos(bus(i,3))+jj*sin(bus(i,3))) * bus(j,2)*(conj(cos(bus(j,3))+jj*sin(bus(j,3)))) * conj(YtYm(k,3));
line_result(k,4)=bus(j,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5))) - bus(j,2)*(cos(bus(j,3))+jj*sin(bus(j,3))) * bus(i,2)*(conj(cos(bus(i,3))+jj*sin(bus(i,3)))) * conj(YtYm(k,3));
line_result(k,5)=line_result(k,3) + line_result(k,4); % 利用公式计算非接地支路的潮流
elseif j==0
line_result(k,3)=bus(i,2)^2*conj(YtYm(k,4));
line_result(k,5)=line_result(k,3)+line_result(k,4);
else
line_result(k,4)=bus(j,2)^2*conj(YtYm(k,5));
line_result(k,5)=line_result(k,3)+line_result(k,4); % 利用公式计算接地支路的潮流
end
end
fprintf(Result_file,'\n\n---------------------- 第%d次迭代算，U最大误差为%12.10f，角度最大误差为%12.10f -----------------------\n',T,max(abs(dU)),max(abs(dAng)));
fprintf(Result_file,'\n节点计算结果：\n节点 节点电压 节点相角（角度） 节点注入功率\n');
for i = 1:row_bus
fprintf(Result_file,' %2.0f ',bus_result(i,1));
fprintf(Result_file,' %11.6f ',bus_result(i,2));
fprintf(Result_file,' %11.6f ',bus_result(i,3));
fprintf(Result_file,' %11.6f + j%11.6f\n',real(bus_result(i,4)),imag(bus_result(i,4)));
end
fprintf(Result_file,'\n线路计算结果：\n节点i 节点j 线路功率Sij 线路功率Sji 线路损耗△Sij\n');
for i = 1:row_line
fprintf(Result_file,' %2.0f ',line_result(i,1));
fprintf(Result_file,' %2.0f ',line_result(i,2));
fprintf(Result_file,' %11.6f + j%11.6f ',real(line_result(i,3)),imag(line_result(i,3)));
fprintf(Result_file,' %11.6f + j%11.6f ',real(line_result(i,4)),imag(line_result(i,4)));
fprintf(Result_file,' %11.6f + j%11.6f\n',real(line_result(i,5)),imag(line_result(i,5)));
end
% 判断是否满足精度误差，如满足则跳出，否则返回继续迭代
if (max(abs(dU))<Err && max(abs(dAng))<Err)
break
end
end
fprintf(Result_file,'\n采用牛顿-拉夫逊迭代法，总共迭代%d次\n',T);
fclose(Result_file);
7.2PQ分解法代码
%通过PQ分解法进行潮流计算。首先读取节点参数和支路参数，
% 然后根据节点类型对节点进行重新排序。
% 接下来，根据支路参数和节点类型，计算导纳矩阵的元素。
% 随后，根据节点的PQ不平衡量，进行迭代计算，
% 直到达到最大迭代次数或满足误差条件为止。
clear all
clc
%已清空工作区
n=5;%请输入节点数;
b=5;%请输入支路数;
ngnd=0;%请输入对地支路数
tree1=xlsread('支路参数');%读取支路参数
node1=xlsread('节点参数');%读取节点参数
node=node1;%支路参数备份（与节点重新排序有关）
tree=tree1;%节点参数备份（与节点重新排序有关）
% gnd=[9 0 0.19];%请输入对地支路（节点号 电导 电纳）
nitmax=1000;%最大迭代次数
err=1e-5;%最大允许误差
%对节点类型分类并排序
%1-平衡节点
%2-PV节点
%3-PQ节点
nPQ=0;nPV=0;nBAN=0;
%对节点重新排序
exch=zeros(1,n);%记录节点编号的变化
j=1;
for i=1:n
if node1(i,2)==3%找出PQ节点个数
nPQ=nPQ+1;
node(j,:)=node1(i,:);%i,j节点位置互换
exch(i)=j;
j=j+1;
end
end
for i=1:n
if node1(i,2)==2%找出PV节点个数
nPV=nPV+1;
node(j,:)=node1(i,:);%i,j节点位置互换
exch(i)=j;
j=j+1;
end
end
for i=1:n
if node1(i,2)==1%找出平衡节点
nBAN=i;
exch(i)=j;
node(n,:)=node1(i,:);%平衡节点放到最后的位置
end
end
%对应修改支路参数
for j=1:b
for m=1:2
tree(j,m)=exch(tree(j,m));
end
end
%修改对地支路参数
% for i=1:ngnd
% gnd(i,1)=exch(gnd(i,1));
% end
%初始化导纳矩阵
YR=zeros(n,n);%电导
YI=zeros(n,n);%电纳
Y=zeros(n,n);%导纳矩阵
%构造节点导纳矩阵，用于求解电力系统中的节点电压
for i=1:b%b：表示系统中支路（branch）的数量
if tree(i,6)==1%tree：一个 b 行 7 列的矩阵，表示每条支路的参数。
mo=tree(i,3)^2+tree(i,4)^2;
%自导纳，使用Π型等值（忽略励磁支路）
%k*在首端，Z*在末端
YR(tree(i,1),tree(i,1))=YR(tree(i,1),tree(i,1))+(tree(i,3)/mo)/tree(i,7)+(tree(i,3)/mo)*(1-tree(i,7))/tree(i,7)^2;
YI(tree(i,1),tree(i,1))=YI(tree(i,1),tree(i,1))+(-tree(i,4)/mo)/tree(i,7)+(-tree(i,4)/mo)*(1-tree(i,7))/tree(i,7)^2;
YR(tree(i,2),tree(i,2))=YR(tree(i,2),tree(i,2))+(tree(i,3)/mo)/tree(i,7)+(tree(i,3)/mo)*(tree(i,7)-1)/tree(i,7);
YI(tree(i,2),tree(i,2))=YI(tree(i,2),tree(i,2))+(-tree(i,4)/mo)/tree(i,7)+(-tree(i,4)/mo)*(tree(i,7)-1)/tree(i,7);
%互导纳
YR(tree(i,1),tree(i,2))=YR(tree(i,1),tree(i,2))-tree(i,3)/(tree(i,7)*mo);
YI(tree(i,1),tree(i,2))=YI(tree(i,1),tree(i,2))+tree(i,4)/(tree(i,7)*mo);
YR(tree(i,2),tree(i,1))=YR(tree(i,2),tree(i,1))-tree(i,3)/(tree(i,7)*mo);
YI(tree(i,2),tree(i,1))=YI(tree(i,2),tree(i,1))+tree(i,4)/(tree(i,7)*mo);
end
end
%线路的导纳元素，计算电力系统的节点导纳矩阵（Y矩阵）
for cc=1:b
if tree(cc,6)~=1%变量tree是一个大小为b×6的矩阵，表示系统中b个分支线路的信息
i=tree(cc,1);
j=tree(cc,2);
mo=tree(cc,3)^2+tree(cc,4)^2;%mo表示这条线路的阻抗模长的平方
%自导纳（不计对地电导）
YR(i,i)=YR(i,i)+tree(cc,3)/mo;
YI(i,i)=YI(i,i)-tree(cc,4)/mo+tree(cc,5);
YR(j,j)=YR(j,j)+tree(cc,3)/mo;
YI(j,j)=YI(j,j)-tree(cc,4)/mo+tree(cc,5);
%YR和YI分别表示实部和虚部，i和j分别表示当前线路连接的起点和终点节点编号。
%互导纳
YR(i,j)=YR(i,j)-tree(cc,3)/mo;
YI(i,j)=YI(i,j)+tree(cc,4)/mo;
YR(j,i)=YR(j,i)-tree(cc,3)/mo;
YI(j,i)=YI(j,i)+tree(cc,4)/mo;
end
end
%添加对地支路
% for cc=1:ngnd
% i=gnd(cc,1);
% YR(i,i)=YR(i,i)+gnd(cc,2);
% YI(i,i)=YI(i,i)+gnd(cc,3);
% end
Y=YR+1i*YI;
%节点导纳矩阵初始化结束
% 迭代初始化
%计算P-Q给定值
snet=zeros(1,n);%每个节点发电机与负荷的净注入功率
for i=1:n
snet(i)=node(i,3)+node(i,4)*1i-node(i,5)-node(i,6)*1i;
%功率的给定值
end
%（n-1）个P的给定值
%（mPQ）个Q的给定值
f=0;%迭代标志
nit=1;%当前迭代次数
%开始迭代
tic%计时开始
while(f==0&&nit<nitmax)%当f为0且nit小于nitmax时继续循环 nitmax：最大迭代次数
%计算P-Q不平衡量
str=zeros(1,n);%各节点净注入功率
for i=1:n% 遍历所有节点
psum=0;% 初始化节点i的有功注入
qsum=0;% 初始化节点i的无功注入
for j=1:n
%节点电压方程（极坐标形式）
psum=psum+node(i,7)*node(j,7)*YR(i,j)*cos(node(i,8)-node(j,8))+node(i,7)*node(j,7)*YI(i,j)*sin(node(i,8)-node(j,8));
%计算每个节点传输有功
qsum=qsum+node(i,7)*node(j,7)*YR(i,j)*sin(node(i,8)-node(j,8))-node(i,7)*node(j,7)*YI(i,j)*cos(node(i,8)-node(j,8));
%计算每个节点传输无功
end
str(i)=psum+qsum*1i;%str：各节点的净注入功率
end
DS=snet-str;%得到PQ不平衡量
%判断是否满足误差要求
DP=zeros(1,n-1);
DQ=zeros(1,nPQ);
%统计前n-1个节点的有功不平衡量绝对值
for i=1:n-1
DP(i)=abs(real(DS(i)));
end
%统计PQ节点的无功不平衡量绝对值
for i=1:nPQ
DQ(i)=abs(imag(DS(i)));
end
A=max(DP);
B=max(DQ);
MAX=max(A,B);
if MAX<err
f=1;%若误差满足要求，则标志位置1，表示停止迭代
end
%若不满足，继续迭代
DP=zeros(1,n-1);
DQ=zeros(1,nPQ);
%统计前n-1个节点的有功不平衡量
for i=1:n-1
DP(i)=real(DS(i));
end
%统计PQ节点的无功不平衡量
for i=1:nPQ
DQ(i)=imag(DS(i));
end
%形成雅可比矩阵
H=zeros(n-1,n-1);
L=zeros(nPQ,nPQ);
for i=1:n-1
for j=1:n-1
H(i,j)=node(i,7)*node(j,7)*YI(i,j);
end
end
for i=1:nPQ
for j=1:nPQ
L(i,j)=node(i,7)*node(j,7)*YI(i,j);
end
end
%求电压和相角的修正值
delt_p=zeros(1,n-1);
delt_a=zeros(1,nPQ);
Ud2=zeros(nPQ,nPQ);
for i=1:nPQ
for j=1:nPQ
if i==j
Ud2(i,j)=node(i,7);
end
end
end
delt_p=-H\DP';
delt_a=-Ud2*L\DQ';
%求得电压和相位的修正量
node(1:n-1,8)=node(1:n-1,8)+delt_p;
node(1:nPQ,7)=node(1:nPQ,7)+delt_a;
nit=nit+1;
end
toc%计时结束
%开始计算发电机功率
%负荷功率均为给定值
%str中记录了结果中每个节点注入的功率
for i=1:n
if node(i,2)==2%PV节点，需求解注入的无功
node(i,4)=imag(str(i))+node(i,6);
elseif node(i,2)==1%平衡节点，需求解注入的有功无功
node(i,3)=real(str(i))+node(i,5);
node(i,4)=imag(str(i))+node(i,6);
end
end
%给出节点电压的计算结果
node=sortrows(node)
xlswrite('计算结果1.xlsx',node,'A2:H6');
%计算线路中的功率分布
scc=zeros(b,8);%预存计算结果
%写入scc数组的前两列
for i=1:b
for m=1:2
scc(i,m)=tree1(i,m);
end
end
%给出节点电压的复数形式
Upha=zeros(1,n);
for i=1:n
Upha(i)=node(i,7)*cos(node(i,8))+1i*node(i,7)*sin(node(i,8));
end
%调整节点导纳矩阵
y1=zeros(n,n);
for i=1:n
for j=1:n
y1(i,j)=-Y(exch(i),exch(j));%节点之间线路的导纳（yij=-Yij）
end
end
%开始计算网络中每一条支路的入端和出端功率，利用一个 for 循环，对名为 scc 的矩阵进行赋值
for i=1:b%b 是一个变量，代表矩阵的行数
%入端功率
stemp=Upha(scc(i,1))*conj(Upha(scc(i,1))-Upha(scc(i,2)))*conj(y1(scc(i,1),scc(i,2)));
scc(i,3)=real(stemp);%Upha 是一个函数，输入参数为节点编号，返回值为该节点的电压幅值
scc(i,4)=imag(stemp);
scc(i,5)=Upha(scc(i,1))*Upha(scc(i,1))*tree1(i,5);
%出端功率
stemp=Upha(scc(i,2))*conj(Upha(scc(i,1))-Upha(scc(i,2)))*conj(y1(scc(i,1),scc(i,2)));
scc(i,6)=real(stemp);
scc(i,7)=imag(stemp);
scc(i,8)=Upha(scc(i,2))*Upha(scc(i,2))*tree1(i,5);
end
xlswrite('计算结果2.xlsx',scc,'A2:H7');%将变量scc写入第2行第1列到第7列（即A列到H列）中
895
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
