[基于牛顿—拉夫逊法的电力系统潮流计算毕业设计] 牛顿拉夫逊潮流计算
来源：司法考试 发布时间：2020-03-18 点击：
毕 业 设 计（论 文）
基于牛顿—拉夫逊法的电力系统潮流计算
专业年级 自动化
学
号
姓
名
指导教师
评 阅 人
20XX年X月
XX大 学
本科毕业设计（论文）任务书
（理 工 科 类）
Ⅰ、毕业设计（论文）题目：
基于牛顿拉夫逊法的电力系统潮流计算
Ⅱ、毕业设计（论文）工作内容（从综合运用知识、研究方案的设计、研究方法和手段的运用、应用文献资料、数据分析处理、图纸质量、技术或观点创新等方面详细说明）：
电力系统的基本任务是安全、可靠、经济地为用户提供电能。而对电力系统正常运行状态的分析和计算是其中十分重要的内容。
本课题拟采用牛顿拉夫逊法进行电力系统稳态，对复杂系统的潮流计算进行全面的分析，并得出结论。
开发工具采用matlab，对具体算例建立模型并进行仿真计算，结果图表显示。开发的软件具有一定的实用性。
参考文献来源于图书馆关于电力系统稳态计算的书籍和中国期刊网中文文献。外文期刊室查阅外文文献，并详细翻译1～2篇英文文献。
Ⅲ、进度安排：
第1～6周
资料翻阅，掌握计算的基本理论及matlab编程
第7～12周
找好实际系统，进行计算
第13～17周
编程和调试
第18～22周
分析及撰写论文
第23周
英文翻译
第24周
准备答辩
Ⅳ、主要参考资料：
1 电力工程基础
2 电力系统稳态分析
3 中国期刊网论文
4 英文文献
5 matlab编程
指导教师：
，
20XX
年
X
月
X 日
学生姓名：
，专业年级：
系负责人审核意见（从选题是否符合专业培养目标、是否结合科研或工程实际、综合训练程度、内容难度及工作量等方面加以审核）：
系负责人签字：
，
年
月
日
摘要
电力系统的潮流计算在电力系统稳态分析和电力系统设计中有很重要的作用，潮流计算也是电力系统暂态分析的基础。潮流计算是根据给定的系统运行条件来计算系统各个部分的运行状况，主要包括电压和功率的计算。到目前为止，利用电子计算机进行电力系统的潮流计算的算法已经出现了很多，其中应用最为广泛的是基于牛顿——拉夫逊法的潮流计算方法。
在利用计算机进行电力系统的潮流计算之前，需要对网络的节点进行划分和编号，建立电力网络的数学模型，即电力系统的网络方程式。本文主要介绍了节点导纳矩阵的形成方法。在形成节点导纳矩阵之前，需要将电力网络进行等值电路的变换，其中主要包括输电线路和变压器的等值电路的变换。
由于牛顿——拉夫逊潮流计算对于初值的给定有比较高的要求。因此在进行牛顿——拉夫逊迭代计算前，先采用高斯——赛德尔迭代法产生一组比较精确的初值。本文详细介绍了高斯——赛德尔法和牛顿——拉夫逊法迭代计算的过程。其中主要内容有迭代方程式的建立，雅克比矩阵的计算，功率和电压的计算，以及在迭代过程中PV节点转化为PQ节点时的处理方法。开发工具采用Matlab编程语言，采用读写Excel电子表格的方法进行数据的输入和输出。
本文采用一个5节点的网络进行实例分析，用Matlab开发的计算程序进行潮流计算，计算结果表明程序的算法具有良好的收敛性和实用性。
关键字：潮流计算，节点导纳矩阵，牛顿——拉夫逊，高斯——赛德尔，Matlab
Abstract
Power flow calculation has a very important role in power system steady-state analysis and power system design, and it is also the basis of transient analysis in power system. Flow calculation is based on given conditions of the power system and calculates the operational status of every part of the system, including voltage and power. So far, there are kinds of algorithm which use the electronic computer in power flow calculation, the most widely used algorithm is the Newton - Raphson power flow calculation method.
Before we the computer in power flow calculation, we need to need to have the nodes of the network classified and numbered and establish a mathematical model of power network, namely the power system network equations. This paper describes the formation of the node admittance matrix . In the formation of the node admittance matrix, we need to transform the power network to equivalent circuit, which includes transformation of transmission lines and transformers.
The Newton - Raphson power flow calculation has a relatively high demand for a given initial value. So before the Newton - Raphson iteration, we use Gauss - Seidel iterative method to produce a more precise initial value. This paper describes the process of Gauss - Seidel and Newton - Raphson iteration. The main contents are the establishment of iterative equation, the calculation of Jacobian matrix
and the calculation of power and voltage, as well as how to deal with the situation when a PV node transform to a PQ node iteration process.We use the Matlab programming language as development tools, the input and output of the data process in the Excel spreadsheets.
In this paper, we utilize a system contains 5 nodes to analyze , the result of
the calculation by the Matlab program shows that the algorithm is convergence and practice.
Key Words：Power flow calculation，node admittance matrix，Newton – Rap son，Gauss – Seidel，Matlab.
目 录
第一章 绪论 1
一、电力系统潮流计算的背景及意义 1
二、潮流计算的发展历史及现状 2
三、潮流计算的发展趋势 4
四、本文主要工作 5
第二章 电力网络的数学模型 6
一、节点电压方程 6
二、节点导纳矩阵的形成 7
（一）输电线路的等值电路 7
（二）变压器的等值电路 8
（三）节点导纳矩阵的计算
9
第三章 电力系统的潮流计算 11
一、迭代法简介 11
二、高斯——赛德尔潮流计算 11
（一）功率方程和变量、节点的分类 12
（二）高斯——赛德尔潮流计算 16
（三）算例分析 21
三、牛顿——拉夫逊潮流计算 24
（一）牛顿——拉夫逊法简介 24
（二）潮流计算时的修正方程 26
（三）算例分析 31
第四章 实例分析与程序设计 34
一、输入数据和输出数据 35
（一）输入数据 35
（二）输出数据 36
二、数学模型计算 36
（一）支路导纳矩阵的计算 36
（二）节点导纳矩阵的计算 38
三、潮流计算 38
（一）高斯——赛德尔潮流计算 38
（二）牛顿——拉夫逊潮流计算 40
四、程序设计 42
（一）主程序的设计 42
（二）子程序的设计 43
（三）数据的输入与输出 44
第五章 总结 45
参考文献 46
附录 48
附录1 源程序 48
1 高斯——赛德尔潮流计算源程序 48
2 牛顿——拉夫逊潮流计算程序 50
附录2 英文文献翻译 63
英文文献 63
中文翻译 73
第一章 绪论
电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件以及系统的界限情况确定整个电力系统各个部分的运行状态：各母线的电压。各元件中流过的功率，系统的功率损耗等等。电力系统的潮流计算的电力系统稳态分析、暂态分析和故障分析的基础。
一、电力系统潮流计算的背景及意义
在电力系统规划设计和现有的电力系统的运行方式的研究中，都需要用潮流计算来定量的分析比较供电方案或运行方式的合理性、可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。所以潮流计算是电力系统一种最重要最基本的运算。
电力系统的潮流计算也分为离线计算和在线计算两种，前者主要用于系统的规划设计和安排系统的运行方式，后者则用于正在运行系统的实时监视和控制。
在电网的设计规划阶段，通过潮流计算，合理的规划接入电源的容量和接入点，合理规划电网的结构，选择无功补偿方案，满足规划水平的大、小方式下的交流交换控制、调峰、调相、调压的要求。
在编年制运行方式时，在预计负荷增长及新设备投运基础上进行潮流计算，可以预计电网的运行情况，发现电网中的薄弱环节，供调度员日常调度控制参考，并为电网改造提供建议和依据。
正常检修以及特殊运行方式下的潮流计算，用于日常运行方式的编制，指导发电厂的开机方式，为有功、无功调整方案和负荷调整方案的制定提供依据，以满足电力网络正常运行的要求。
预测因事故或者电网负荷发生变化时，电网运行状态的变化，并以此来制定相应的处理方案。
二、潮流计算的发展历史及现状
在数字计算机出现之前，电力系统的潮流计算主要是借助于交流台通过人工计算完成，交流台模拟了电力系统，因此在交流计算台上计算潮流分布时，计算人员可以随时监视系统各个部分运行状态是否满足要求，如果发现某些部分不合理，则可以立即进行调整。这种方法直观，但是人工操作工作量大且易出错。
电力系统的潮流计算的计算量非常巨大，通过人来计算是非常困难的。随着电子计算机的产生和发展，人们开始探索利用计算机来进行潮流计算的方法。从50年代开始到现在，潮流计算曾采用了不同的方法，这些方法主要围绕着对潮流计算的一些基本要求进行的。对潮流计算的要求可以归纳为以下几点：
（1）计算方法的可靠性和收敛性
（2）对计算机内存量的要求
（3）计算速度
（4）计算的方便性和灵活性
其中第一个要求是最主要最基本的要求，即计算方法可行，计算的次数在计算过程中逐渐减少，而不是计算次数越来越多。电力系统的潮流计算在数学上是一组多元非线性方程式的求解问题，其方法都离不开迭代。因此，对潮流计算方法，首先要求它能可靠的收敛，并给出正确答案，由于电力系统结构及参数的一些特点，并且随着电力系统的不断扩大，潮流计算的方程式的阶数越来越多。(一般在几十阶甚至几百阶以上)，对这样的方程式不是任何是任何数学方法都能保证给出正确答案的。这种情况称为促使电力系统计算人员不断线的更可靠方法的重要因素。
在用数字计算机解电力系统潮流计算的开始阶段，普遍采用以节点导纳矩阵为基础的逐次带入法，即导纳法。这个方法的原理比较简单，要求数字计算机的内存比较小，适应50年代的电子计算机的制造水平和当时的电力系统理论水平。但它的收敛性比较差，当系统的规模增大时，迭代次数急剧上升，在计算中往往出现迭代不收敛的情况。这就迫使电力系统的计算人员转向以阻抗矩阵为基础的逐次代入法，即阻抗法。
60年代初期，数字计算机已发展到第二代，计算机的内存和速度发生了很大的飞跃，从而为阻抗法的采用创造了条件。阻抗法要求数字计算机贮存表征系统接线和参数的阻抗矩阵，这就需要大量的内存。而且阻抗法每迭代一次都要求顺次取阻抗矩阵中的每一个元素进行计算，因此，每次迭代的运算量很大。这两种情况都是过去电子管计算机无法适应的。
阻抗法改善了系统潮流计算的收敛性问题，解决了导纳法无法求解的一些系统的潮流计算，在60年代获得了广泛的应用，曾为我国的电力系统的设计、运行和研究做出了很大的贡献。
阻抗法的缺点是占用的计算机 的内存比较大，每次迭代的计算量大。当系统不断扩大时，这些缺点就更加突出。一个内存16K的计算机在采用阻抗法时只能计算100个节点以下的系统。这样，我国很多电力系统为了采用阻抗法潮流计算就不得不对系统进行相当的简化工作。
为了克服阻抗法在内存和速度上的缺点，60年代中期发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需要存储各个地区系统的阻抗矩阵及他们之间的联络线的阻抗，这样大幅度的节省了内存容量，提高了计算速度。
克服阻抗法的缺点的另一个方法是采用牛顿—拉夫逊法。牛顿—拉夫逊法是数学中解决非线性方程式的典型方法，有较好的收敛性。在解决电力系统的潮流计算问题时，是以导纳矩阵为基础的，因此只要我们能在迭代过程中尽可能保持方程式的系数矩阵的稀疏性，就可以大大提高牛顿—拉夫逊法潮流程序的效率。自从60年代中期，在牛顿—拉夫逊法中利用了最佳顺序消去法以后，牛顿法在收敛性、内存要求和速度方面都超过了阻抗法，成为60年代末期以后广泛采用的优秀方法。
与此同时，为了保证可靠的收敛，在我国还进行了利用非线性规划法计算潮流计算的研究。
随着电力系统的日益扩大和复杂化，特别是电力系统逐步实现自动控制的需要，对系统潮流计算在速度、内存以及收敛性的方面都提出了更高的要求。
70年代以来，潮流计算方法通过不同的途径继续向前发展，其中比较成功的就是P—Q分解法。这个方法，根据电力系统的特点，抓住主要矛盾，对出数学的牛顿—拉夫逊法进行了改进，从而在内存容量以及计算速度方面都大大向前
迈进了一步。使一个32K内存容量的数字计算机可以计算1000个节点的潮流计算问题，此方法计算速度以能用于在线计算，做系统静态安全监测。目前，我国很多电力系统都采用了P—Q分解法潮流程序。
近20多年来，潮流算法的研究仍然非常活跃，但是大多数研究都是围绕改进牛顿法和P-Q分解法进行的。此外，随着人工智能理论的发展，遗传算法，人工神经网络和模糊算法也被逐渐引入到潮流计算当中，但是到目前为止，这些新的方法还不能取代牛顿——拉夫逊法和P-Q分解法的地位。由于电力系统规模的不断扩大，对计算速度的要求不断提高，计算机的并行计算技术也将在潮流计算中得到广泛的应用，成为重要的研究领域。
三、潮流计算的发展趋势
现在应用最为广泛的牛顿——拉夫逊法是将非线性的潮流方程逐次线性化，为了进一步提高算法的收敛性和计算速度，人们考虑采用将泰勒级数的高阶项或非线性项也考虑进来，于是产生了二阶潮流算法。后来又提出了根据直角坐标形式的潮流方程是一个二次代数方程的特点，提出了采用直角坐标的保留非线性快速潮流算法。
对于一些病态系统，应用非线性潮流计算方法往往会造成计算过程的振荡或者不收敛，从数学上讲，非线性的潮流计算方程组本来就是无解的。这样，人们提出来了将潮流方程构造成一个函数，求此函数的最小值问题，称之为非线性规划最优潮流的计算方法。优点是原理上保证了计算过程永远不会发散。如果将数学规划原理和牛顿潮流算法有机结合一起就是最优乘子法。另外，为了优化系统的运行，从所有以上的可行潮流解中挑选出满足一定指标要求的一个最佳方案就是最优潮流问题。最优潮流是一种同时考虑经济性和安全性的电力网络分析优化问题。OPF 在电力系统的安全运行、经济调度、可靠性分析、能量管理以及电力定价等方面得到了广泛的应用。
另外随着直流输电技术的研究和发展，直流输电网络和交流混合电力系统的潮流计算也有了一定发展，随着直流输电技术的不断应用，混合电力系统的潮流计算
必将获得一个广阔的发展空间。
四、本文主要工作
本文主要工作主要是详细介绍牛顿——拉夫逊法的原理、算法设计和Matlab程序的编写。以及牛顿拉夫逊法的优缺点以及对于其缺点的改进方法。
牛顿—拉夫逊法，把非线性方程式的求解过程变成反复对相对应的线性方程式的求解过程，通常称为逐次线性化过程，这是牛顿—拉夫逊法的核心。每一次的迭代都要先解修正方程，然后用解得的各节点电压变量（修正量）求个节点的新值（修正后值）。步骤是：
（1）设一组结点电压；
（2）求功率和电压的不平衡量；
（3）求雅克比矩阵的各个元素；
（4）解修正方程式。
牛顿——拉夫逊法最重要的一步是计算雅克比矩阵。这是将非线性潮流方程线性化的关键步骤。牛顿——拉夫逊法具有很好的收敛性，计算速度快，计算结果准确。但是牛顿——拉夫逊法对于初值有比较高的要求，当给定的初值与精确值相差较大时，计算结果会产生很大的误差，甚至不能收敛。为了解决这个问题，通常先利用高斯——赛德尔法进行计算，将计算得到的结果作为牛顿——拉夫逊法的初值进行计算。
本文还介绍了电力网络数学模型的建立，这是利用数字计算机进行潮流计算的基础，主要包括输电线路和变压器数学模型的建立和节点导纳矩阵的计算。
第二章 电力网络的数学模型
电力网络的数学模型指的是将网络的有关参数和变量及其相互关系归纳起来所组成的，可以反映网络性能的数学方程式组。在利用计算机的复杂电力系统的潮流计算中用的最多的是节点电压方程。节点电压方程用节点电压和来表示支路电流，根据基尔霍夫电流定律列出方程组。
一、节点电压方程
如图2.1所示是一个电力系统的等值网络图，它共有三个节点。
图2.1 电力系统的等值网络图
根据基尔霍夫电流定律，对该电路图列写节点电压方程得：
（2.1）
式（2.1）就是图2.1所示电力网络等值电路的数学模型。将其用矩阵形式表式为：
(2.2)
式（2.2）可以展开为
(2.3)
在式（2.2）中是节点注入电流的列向量。在电力系统计算中，节点注入电流可以理解为与该节点相连的正电流源与负电流源之和，其中规定注入该节点的电流为正，流
出该节点的电流为负。有的节点不与电流源相连，则其注入电流为零，如图2.1中的节点3。是节点电压列向量，节点电压是该节点对参考地的电压。是一个阶的节点导纳矩阵，其就等于网络中出参考地之外的节点数。
二、节点导纳矩阵的形成
节点导纳矩阵在利用计算机进行潮流计算中具有十分重要的地位，它是电力网络的数学表示形式。电力网络拓扑结构经过一系列的等效之后，形成类似于图2.1所示的等值电路。在等值电路的求解过程中，输电线路和变压器的等值电路的求取是主要的工作。
（一）输电线路的等值电路
输电线路等值电路一般以图2.2所示的形式表示。
图2.2 输电线路的等值电路
、和都是导纳，单位是西门子。在计算输电线路的等值电路时，一般要先知道线路的长度，单位长度线路的阻抗和单位长度线路的导纳。对于小于100公里的输电线路一般忽略其导纳，即在图2.2中的和都等于零，等于线路阻抗的倒数。
当线路长度在100公里到300公里之间时，要考虑线路的导纳，导纳用集中参数表示。此时、和的计算方法为：
(2.4)
当线路长度大于300公里时，则要考虑线路的分布参数。此时线路的等值电路的计算方法为：
(2.5)
其中称为线路特性阻抗，成为线路传播常数。
(2.6)
（二）变压器的等值电路
在变压器的等值电路的计算中，我们将变压器看做是一个理想变压和一个阻抗的串联，并且忽略了变压器的漏抗。理想变压器的变比是实际的变压器变比，而阻抗一般是折算到低压侧的变压器短路阻抗，串联在理想变压器的低压侧。如图2.3所示。
图2.3
变压器的等值电路
图2.3所示的等值电路仍然不能直接形成节点导纳矩阵，需要进一步化为图2.2所示的等值电路，此时的计算方法为：
(2.7)
而对于三绕组的变压器，则可以将其中一个绕组看成是低压绕组，另外两个绕组看成为高压绕组，把绕组变压器化为两个双绕组的变压器。如图2.4所示。
图2.4
三绕组变压器的等值电路
将三绕组变压器化为两个双绕组变压器后就可以按照双绕组变压器等值电路的求法进行进一步的化简。
（三）节点导纳矩阵的计算
将电力网络的拓扑化为如图2.1所示的等值电路后，就可以进行节点导纳矩阵的计算了。节点导纳矩阵求取时，注意一下几点。
（1）节点导纳矩阵是方阵，其阶数就等于网络中除参考地之外的节点数。
节点导纳矩阵是稀疏矩阵，其各非零非对角元数就等于与该行相对应节点所连接的不接地支路数。如图2.1所示，与节点2对应的第二行非零非对角元数为2.
（2）节点导纳矩阵的对角元就等于与该节点所连接的导纳的总和。如图2.1中，与节点2对应的对角元。
（3）节点导纳矩阵的非对角元就等于连接节点，支路导纳的负值。如图2.1所示，
。
（4）节点导纳矩阵一般是对称矩阵。
第三章 电力系统的潮流计算
一、迭代法简介
在电力系统的计算中，迭代法主要用于求解非线性方程。但在某些情况下，迭代法也适用于线性方程，下面以一个线性方程为例介绍迭代法。
【例 3.1】解方程组
【解】用消去法求得它的解是，。下面用迭代法求解。把该方程组变形为：
取，的初值，代入上式的右边，即可从左边解得 ，。然后再将，代入上式，解得又一组解，经过三次迭代后可得
，，误差已经减小到了千分之一。
二、高斯——赛德尔潮流计算
建立了节点导纳矩阵，就可以进行潮流计算。潮流计算的额主要思想是迭代法。其关键是要建立一个方程组。但是由于在工程实践中，通常我们已知的既不是节点电压，也不是节点电流，而是节点功率，实际计算中几乎无一例外的要迭代解非线性的节点电压方程（*表示复数的共轭）。
（一）功率方程和变量、节点的分类
设有简单系统如图3.1所示。图中，分别为母线1、2的等值电源功率；，分别为母线1、2的负荷功率；它们的合成，分别为母线1、2的注入功率，与之对应的电流，分别为母线1、2的注入电流。于是，
(3.1)
图3.1 简单系统及其等值网络
（a）简单系统；（b）简单系统的等值网络；（c）注入功率和注入电流
(3.2)
如令
(3.3)
并将式（3.3）代入式（3.2）展开，将有功功率和无功功率分列，可得
(3.4)
式（3.4）就是这个简单系统的功率方程。
由式（3.4）可见，在功率方程中，母线电压的相位角是以差=的形式出现的，亦即决定功率大小的是相对相位角或相对功率角，而不是绝对相位角或绝对功率角。
由式（3.4）可得
(3.5)
它们都是母线电压，和相位角，相对相位角或的函数。
由式（3.4）还可见，在这四个一组的功率方程组中，除网络参数、、、外，共有十二个变量，它们是：
负荷消耗的有功、无功功率——、、、；
电源发出的有功、无功功率——、、、；
母线或节点的大小和相位角——、、、；
因此，除非已知其中的八个变量，否则将无法求解。
在这十二个变量中，负荷消耗的有功、无功功率无法控制，因它们取决于用户。它们就称为不可控变量或扰动变量。之所以称扰动变量是由于这些变量出现事先没有预计的变动时，系统将偏离它们的原始运行状况。不可控变量或扰动变量用列向量表示。
余下的八个变量中，电源发出的有功、无功功率是可以控制的自变量，因而它们就称为控制变量。控制变量用列向量表示。
最后余下的四个变量——母线电压或节点电压的大小和相位角——是受控制变量控制的因变量。其中、主要受、的控制；、主要受、的控制。这四个变量就是系统的状态变量。状态变量一般用列向量表示。
变量的这种分类也适用于个节点的复杂系统。只是对于这种复杂的系统，变量数将增加到个，其中扰动变量、控制变量、状态变量各个。换言之，扰动向量、控制向量、状态向量，都是阶的列向量。
看来似乎将变量做如上分类后，只要已知或给定扰动变量和控制变量就可以运用功率方程式（3.4）解出状态变量。其实不然，因已如上述，功率方程中，母线或节点电压的相位角是以相对值出现的，以致式（3.4）中的、变化同样大小时，功率的数值不变。从而也不可能运用它们求取绝对相位角，也如上所述，系统的功率损耗本身就是状态变量的函数，在解得状态变量前，不可能确定这些功率损耗。
为克服上述困难，可对变量的给定稍作调整：
在一个有个节点的系统中，只给定对控制变量、，余下的一对控制变量、待定。这一对控制变量、将使系统的功率保持平衡。
在这系统中，给定一对状态变量、，只要求确定对状态变量、。给定的通常为零。这实际上就是相当于取节点s的电压相量为参考轴。
这样，原则上可以从个方程中解出个未知量。但是实际上，这个解还应满足以下的约束条件，这些约束条件是保证系统正常运行所不可缺少的。其中对控制变量的约束条件是
；
对没有电源的节点则为
；
这些限制条件取决于一系列技术经济因素，应根据实际情况而定。
对于状态变量的约束条件则是：
这条件表示，系统中各节点的电压大小不得越出一定的范围，因系统运行的基本要求之一就是要保证良好的电压质量。
对于某些状态变量还有如下的约束条件
这条件主要是保证系统运行的稳定性所要求的。
对于扰动变量、不可控，对它们没有约束。
考虑到这些约束条件后，对于某些节点，不是给定控制变量、而留下状态变量、待求，而是给定这些节点的和而留下和待求。这其实意味着让这些电源调节它们发出的无功功率以保证与之连接的节点电压为定值。
这样，系统的节点就因给定变量的不同而分为三类。
第一类称PQ节点。对这类节点，等值负荷功率、和等值电源功率、是给定的，从而注入功率、是给定的。待求的则是节点电压的大小和相位角。属于这一类的节点又按给定有功、无功功率的发电厂母线和没有其它电源的变电所母线。
第二类节点称PV节点。对这类节点，等支负荷和等值电源的有功功率、是给定的，从而注入有功功率是给定的。等值负荷的无功功率和节点电压的大小也是给定的。待求的则是等值电源的无功功率，也就是要求注入无功功率和节点电压的相位角。有一定无功功率储备的发电厂和有一定无功功率电源的变电所母线都可选作PV节点。
第三类节点称平衡节点。潮流计算时，一般都只设一个平衡节点。对这个节点，等值负荷功率、是给定的，节点电压的大小和相位角是给定的。待求的则是等值电源功率、，从而要求注入功率、。担负调整系统频率任务的发电厂母线往往被选作平衡节点。
进行计算时，平衡节点是不可少的；PQ节点是大量的，PV节点是较少的，甚至可以没有。
（二）高斯——赛德尔潮流计算
自1956年成功地运用数字计算机计算潮流分布以来，曾先后出现过许多种结算方法。目前常用的有运用节点导纳矩阵的牛顿——拉夫逊法和由该法派生的P-Q分解法。但是由于牛顿——拉夫逊法对初值的选取要求严格，某些程序的第一、二次迭代又往往采用高斯——赛尔德法估计初值。以下是对高斯——赛尔德的法的具体介绍。
高斯——赛尔德法比较简单，是由于它可以直接迭代解节点电压方程。因将节点电压方程展开，可得
(3.6)
移项后，又可得
(3.7)
将式（3.7）进一步展开后，就可以用高斯——赛尔德法迭代求解。例如，对节点1为平衡节点，其余都是PQ节点的网络，上式可展开如式（3.8）。
(3.8)
式中 ——给定的各节点注入功率的共轭值
——给定的平衡节点电压
——迭代次数
上式是按高斯——塞尔德法解方程式组时的标准模式书写的，按这种模式，式中等号右侧的采用经次迭代后的值；等号右侧的，当时，采用次迭代后的值；当，采用次迭代后的值。
迭代解式（3.8）的步骤是：先假设一组，一般可设，将它们代入第一式，可解得。然后将、、…代入第二式，可解得。再将、、…代入第三式，又可解得。依次类推，直至解得。这就是第一次迭代。第一次迭代结束时，解得了所有的。
再将解得的这组再一次代入式（3.8）进行第二次迭代。先代入第一式，解得。然后将、、…代入第二式，可解得
。再将、、…代入第三式，又可解得。依次类推，直至解得。这就是第二次迭代。第二次迭代结束时，解得了所有的。
如此不断迭代，直至某一次迭代后的解得的与前一次迭代后的值相差小于事先给定的允许误差，即而停止。因这个条件满足就是迭代结束的标志。
网络中往往还会有PV节点，而且，这中PV节点的注入无功功率由于受到电源供应的无功功率的限制，对于这些节点，如上的计算步骤应修改如下。
设节点是PV节点，则由于已经给定，在每次迭代求得后，应首先将求得的修正为，即将求得的电压大小由改为，而求得的相位则不改动。然后，将它和其它节点一起代入式（3.9）以求取,亦即按式（3.10）计算
(3.9)
(3.10)
求得后，再将其代入下式以求取
(3.11)
显然，上两式中的都应为修正后的值，而列出上两式时设PV节点的编号。
迭代过程中往往会出现越限，即按式（3.10）求得的不能满足的情况。考虑到实践中对节点电压的限制不如对节点功率的限制严格，出现这种情况时，只能用或代入式（3.11）以求取，而在求得后，不在像对那样修正它的值。换言之，这时只能满足约束条件而不能满足=给定值。而事实上，这时PV节点以变为了PQ节点。
迭代收敛后，就可以计算平衡节点的功率
(3.12)
并计算各线路上流动的功率
(3.13a)
(3.13b)
以及各线路上的功率损耗
(3.14)
这样由式（3.11）求得了所有PQ节点的电压大小和相位；由式（3.10）、（3.11）求得了PV节点的无功功率和电压的相位角；由式（3.12）求得了平衡节点的视在功率，由式（3.13）、（3.14）求得了所有线路上流动的功率；换言之，网络中所有支路的功率和功率损耗都已经确定，潮流分布的计算已经完成。
图3.2 高斯——赛德尔程序框图
高斯——赛德尔潮流计算的程序流程图如图3.2所示。对于图3.2，作如下说明。
（1）设节点1、2、3…m中，除去一个平衡节点s外，其余的为PQ节点。节点m+1、m+2…n则全为PQ节点。在输入数据时，节点的编号要按照上述要求进行编号。
（2）框1和框2的步骤在调用高斯——塞尔德子程序之前完成，在调用高斯塞尔德程序时输入的数据有：
平衡节点的电压；
PQ节点的视在功率；
PV节点的无功功率和节点电压的大小，；
PV节点的无功功率限额，；
（3）对于框4，迭代次数k的设置，是为了防止迭代不收敛时来终止迭代的。变量f是标志变量，当某一个PV节点在计算过程中因没有满足限制条件而转化为PQ节点时，f中的对应位变为1，否则为0；
（4）对于框8，在修正时，应保持其大小始终等于给定的PV节点的电压大小，只对其相位角进行修正。
（三）算例分析
【例3.2】三节点等值网络如图3.3所示。图中，节点1为平衡节点，给定=1.0+j0，节点2为PQ节点，给定；节点3为PV节点，给定，。各支路导纳以示于图3.3。
图3.3 三节点等值网络
【解】图3.3是一个已经经过转换的等值网络，我们先建立一个矩阵，是一个阶的矩阵，其中第1、2、3节点分别对应的第1、2、3行（列），参考地对应第4行（列），的每一个元素表示连接第个节点和第个节点的导纳，或者表示连接第个节点和地之间的导纳的值。我们称矩阵为支路导纳矩阵，以区别于节点导纳矩阵，于是，由图3.3我们可得矩阵为
然后根据前文1.3.3节中介绍的节点导纳矩阵的形成方法，得到节点导纳矩阵，
我们设定精确度为0.00001，=10，=0.02，将数据代入根据图3.2的流程图编写的程序中进行计算，以下是计算结果。
表3.1
第一次迭代后的节点电压
节点编号
1
2
3
节点电压
1.0000+j0.0000
0.9680-j0.0260
1.0988+j0.0518
最大电压变化量dumax为：dumax=0.0519;
表3.2
第二次迭代后的节点电压
节点编号
1
2
3
节点电压
1.0000+j0.0000
0.9662-j0.0260
1.0987+j0.0534
最大电压变化量dumax为：dumax=0.0018;
表3.3
第三次迭代后的节点电压
节点编号
1
2
3
节点电压
1.0000+j0.0000
0.9661-j0.0260
1.0987+j0.0534
最大电压变化量dumax为：dumax=;
表3.4
第四次迭代后的节点电压
节点编号
1
2
3
节点电压
1.0000+j0.0000
0.9661-j0.0260
1.0987+j0.0534
最大电压变化量dumax为：dumax=;
由以上计算结果可以看出，经过四次迭代后，计算结果的精确度就达到了，从dumax可以看出，迭代计算能很快的收敛。
迭代计算结束后，再根据式（3.10）和式（3.12）计算PV节点的无功功率和平衡节点的视在功率，并根据式组（3.13）计算线路上的视在功率，根据式（3.14）计算线路上的功率损耗。计算结果如下：
表3.5
节点注入功率
节点编号
1
2
3
节点注入功率
0.4437+j0.2405
-0.8000-j0.6000
0.4000+j0.0624
表3.6
各线路上的流动功率
i
j
1
2
3
1
0
0.8107+j0.6429
-0.3670-j4024
2
-0.8000-j.6000
0
0
3
0.3818+j0.0624
0-j0.3993
0
表3.7
各线路上的损耗功率
i
j
1
2
3
1
0
0.0107+j0.0429
0.0147-j0.3400
2
0.0107+j0.0429
0
0-j0.3993
3
0.0147-j0.3400
0-j0.3993
0
三、牛顿——拉夫逊潮流计算
牛顿——拉夫逊是广泛采用的解非线性方程式组的方法，也是当前广泛采用的计算潮流分布的方法。这里的非线性方程式组就是非线性功率方程。
（一）牛顿——拉夫逊法简介
设有非线性方程式组如下
(3.15)
其近似解为、、…。设近似解与精确解分别相差、、…
则如下的关系式应该成立
(3.16)
上式中任何一式都可按泰勒级数展开。以第一式为例
式中、…分别表示以、、…代入这些偏导数的表示式时计算所得；则是一包含、、…的高次方与的高阶偏导数乘积的函数，如果初值与精确解相差不大，则可以略去。于是式（3.16）可以写为：
(3.17)
这是一组线性方程式组或线性化了的方程式组，成为修正方程式组。将它改写为矩阵形式为：
(3.18)
或简写为
(3.19)
式中，称为雅克比矩阵。将代入，可得、中的各元素。然后运用解线性矩阵方程的解法，可求的，从而可以求得经第一次迭代后的新值=
。再将代入，从而可以求得。如此不断迭代，直至满足对精确度的要求。
运用这种方法计算时，的初值要选择得比较接近它们的精确解，否则迭代的过程可能不收敛。正因如此，在运用牛顿——拉夫逊法的潮流计算中，第一、二次的迭代先采用高斯——赛德尔法得到一个比较精确的初值，然后在用牛顿——拉夫逊法进行迭代。
（二）潮流计算时的修正方程
运用牛顿——拉夫逊法计算潮流分布时，节点导纳矩阵的形成，平衡节点和线路功率的计算都和运用高斯——赛德尔法时相同，不同的知识迭代过程。迭代过程中，两种方法应用的基本方程式又都是，只是高斯——赛德尔将其展开为电压方程式，用牛顿——拉夫逊法时，将其展开为功率方程式。
(3.20)
式中，第一部分为给定的节点注入功率，第二部分为由节点电压求得的节点注入功率，它们二者的差就是节点功率的不平衡量。有待解决的问题就是各节点功率不平衡量都趋近于零时，各节点电压应具有何值。将式（3.20）与式（3.15）对照，式（3.15）中的就对应这里的节点功率不平衡量，而式（3.15）中的、…，则对应于这里的节点电压。
建立了这种对应关系，就可以仿照式（3.18）列出修正方程式，并迭代求解。当节点电压以直角坐标系的形式表示时，其修正方程的建立如下所述。
节点电压以直角坐标表示时，令，，并将式（3.20）改写为
(3.21)
将实部和虚部分开写，可得
(3.22a)
(3.22b)
鉴于系统中还有电压大小给定的PV节点，还应由列出
(3.22c)
式中，第一部分为给定的PV节点的节点电压的平方，第二部分为求得的节点电压的平方，它们二者的差可看做为节点电压大小的不平衡量。
对于一个具有个节点的网络，式组（3.22）共有个方程式。如仍按前述的节点编号划分法，式（3.22a）共有个，包括平衡节点外所有节点有功功率不平衡量的表示式，即；式（3.22b）有，包括所有PQ节点无功功率不平衡量的表示式，即；式（3.22c）类型的有个，包括所有PV节点电压大小的不平衡量的表示式，即。平衡节点的功率和电压之所以不包括在这方程式组内是由于平衡节点的注入功率不可能事先给定，而平衡节点的电压已知，不必列节点电压不平衡量的方程。
这样就可建立类似式（3.18）的修正方程式如下
(3.23)
式中的、、分别就是如式组（3.22）所示。而式中的雅克比矩阵的各个元素分别为
(3.24)
当时：
(3.25a)
当时，为使这些偏导数的表示更加简洁，先引入节点注入电流的表示式如下
(3.26)
然后可得
(3.25b)
以上是牛顿——拉夫逊潮流计算的方法，牛顿——拉夫逊潮流计算的程序框图如图3.4所示
对于图3.4作以下几点说明
（1）原始数据的输入和节点导纳矩阵的形成是在进入牛顿——拉夫逊潮流计算子程序之前完成的。
（2）框2中为了得到比较精确的节点电压初值，要先调用高斯——赛德尔潮流计算子程序，以得到比较精确的节点电压初值。然后在进行牛顿——拉夫逊潮流计算。
（3）迭代次数k设为1000，每次迭代后减1，直至为零，如果迭代不收敛，则当k减为零时，程序自动终止。
（4）雅克比矩阵计算过程中，注意PV节点转化为PQ节点的情况。处理方法与高斯——赛德尔中的相同。
图3.4 牛顿——拉夫逊潮流计算流程图
（三）算例分析
【例3.3】五节点等值网络如图3.5所示。图中，节点1为平衡节点，给定；其它节点都是PQ节点，给定，，，。节点导纳矩阵已知，在直角坐标系下表示各节点电压，进行牛顿——拉夫逊潮流计算。
图3.5 五节点等值网络
【解】根据已知的条件和限制，将输入数据输入到matlab的工作空间。然后调用牛顿——拉夫逊潮流计算程序进行计算。
在进行牛顿——拉夫逊潮流计算之前，我们先调用高斯——赛德尔潮流计算程序计算初值。我们将高斯——赛德尔潮流计算的精确度定为，将牛顿——拉夫逊潮流计算的精确度定为，进过计算后得出以下结果。
高斯——赛德尔潮流计算结果：
迭代次数23次，精确度；
表3.8
节点注入功率：
节点编号i
1
2
3
4
5
节点功率
1.2902
-j0.0745
0.2000
+j0.2000
-0.4500
-j0.1500
-0.4000
-j0.0500
-0.6000
-j0.1000
表3.9
节点电压：
节点编号i
1
2
3
4
5
节点功率
1.0600
1.0463
-j0.0511
1.0204
-j0.0889
1.0193
-j0.0947
1.0122
-j0.1088i
牛顿——拉夫逊潮流计算结果：
迭代次数2次，精确度；
表3.10
节点注入功率：
节点编号i
1
2
3
4
5
节点功率
1.2947- -j0.0743
0.2000
+j0.2000
-0.4500
-j0.1500
-0.4000
-j0.0500
-0.6000
-j0.1000
表3.11
节点电压：
节点编号i
1
2
3
4
5
节点功率
1.0600
1.0462
-j0.0512
1.0203
-j0.0892
1.0192
-j0.0950
1.0121
-j0.1090
从以上计算结果可以看出，牛顿——拉夫逊潮流计算的计算精确度和收敛速度远大于高斯——赛德尔潮流，但是牛顿——拉夫逊潮流计算对初值的要求比较高，当高斯——赛德尔潮流计算输出结果的精度为0.01时，牛顿——拉夫逊潮流计算在进行了1000次迭代后仍然不能收敛，所以在应用牛顿——拉夫逊潮流计算时，为了保证其能够收敛，必须要提供足够精确的初值，这可以通过提高高斯——赛德尔潮流计算的输出精确度来实现。通过将高斯——赛德尔潮流计算和牛顿——拉夫逊潮流计算进行配合使用，既可以减小迭代的次数，保证计算的收敛，又可以提高其计算的精确度。
根据式（3.10）和（3.12）计算PV节点的无功功率和平衡节点的视在功率，根据（3.13）计算线路上的视在功率，根据式（3.14）计算线路上的功率损耗。
表3.12
各线路上的流动功率
i
j
1
2
3
4
5
1
0
0.8877
-j0.1140
0.4070
-j0.0221
0 - j0.0618
0 - j0.0618
2
-0.8747
+j0.0011
0
0.2459
-j0.0359i
0.2783
-j0.0417
0.5472
-j0.0034
3
-0.3951
-j0.0616
-0.2435
-j0.1045
0
0.1886
-j0.0992
0 - j0.0577
4
-0.0000
-j0.0576
-0.2750
-j0.0959
-0.1883
-j0.0150
0
0.0633
-j0.0543
5
-0.0000
-j0.0414
-0.5370
-j0.0976
-0.0000
-j0.0414
-0.0630
-j0.0439
0
表3.13
各线路上的流动功率
i
j
1
2
3
4
5
1
0
0.0130
-j0.1128
0.0119
-j0.0838
-0.0000
-j0.1194
-0.0000
-j0.1032
2
0.0130
-j0.1128
0
0.0024
-j0.1404
0.0033
-j0.1376
0.0102
-j0.1009
3
0.0119
-j0.0838
0.0024
-j0.1404
0
0.0004
-j0.1143
-0.0000
-j0.0991
4
-0.0000
-j0.1194
0.0033
-j0.1376
0.0004
-j0.1143
0
0.0003
-j0.0982
5
-0.0000
-j0.1032
0.0102
-j0.1009
-0.0000
-j0.0991
0.0003
-j0.0982
0
第四章 实例分析与程序设计
前面两章已经介绍了电力系统数学模型的建立和潮流计算的理论和步骤，本章将通过一个实例，利用上述方法对电力网络进行实例分析。
【例4.1】如图4.1所示，是一个电力网络的拓扑结构。其中节点1是平衡节点，节点2、3、4是PQ节点，节点5是PV节点。各节点之间输电线路的长度，输电线路单位长度（1千米）的阻抗、导纳以及变压器的变比和阻抗（折算到低压侧）的值如表4.1所示。
图4.1 电力网络实例
表4.1
电力网络参数
节点编号i
节点编号j
阻抗()
导纳
长度
变比
3
1
0+j0.03
0
0
1.05
2
3
0.01+j0.035
0
10
0
2
4
0.0002+j0.00125
0+j0.0025
200
0
3
4
0.0004+j0.0015
0+j0.0025
200
0
4
5
0+j0.0015
0
0
1.05
表4.1中如果节点i和j之间是变压器，则阻抗为变压器折算到低压侧的阻抗，并且忽略了变压器的漏抗，这时其线路长度为0，变比不为0。电力网络中的同一支路不能再表中重复出现。
例图4.1中，平衡节点1的电压为1.05+j0；PV节点5的有功功率为5，电压幅值为1.05；PQ节点2、3、4的视在功率分别为1.6+j0.8、3.7+j1.3、2+j1;
一、输入数据和输出数据
输入数据和输出数据都以数组的形式在Excel文档中读写，组织好输入数据和输出数据是编写程序的关键一步。
（一）输入数据
潮流计算的输入数据有：
电力网络的参数矩阵data，数据组织形式和数据含义如表4.1所示；
PQ节点的视在功率向量S_PQ，其中包括平衡节点的视在功率（暂定为0）；
PV节点的有功功率向量P_PV；
PV节点的电压幅值向量U_PV；
电网节点数n；
PQ节点数m（包括平衡节点）；
平衡节点编号s；
平衡节点电压us；
PV节点无功功率最大值Qmax；
PV节点无功功率最小值Qmin；
精确度epxl。
根据例4.1，我们组织输入数据如下：
；
；
；
；
n=5；
m=4；
；
Qmax=10；
Qmin=0.1；
epxl=0.00001；
（二）输出数据
输出数据主要有:
节点注入功率s_node，向量，表示各个节点的注入功率
节点电压u_node，向量，表示各个节点的节点电压；
线路流动功率Sij，矩阵，表示有节点i流向节点j的视在功率；
线路损耗功率dSij，矩阵，表示节点i和节点j之间线路上的功率损耗；
二、数学模型计算
建立电力网络数学模型的目的是得到电力网络的节点导纳矩阵，其中主要的工作包括输电线路等值电压和变压器等值电路的计算，支路导纳矩阵的计算，节点导纳矩阵的计算。其中的支路导纳矩阵是一个(n+1)阶的方阵，其前n行（列）对应着网络的n个节点，第(n+1)行（列）对应着网络的参考地，的每一个元素表示连接第个节点和第个节点的导纳，或者表示连接第个节点和地之间的导纳的值。
（一）支路导纳矩阵的计算
按照4.1.1节组织的输入数据输入到潮流计算的程序之后，第一步数要计算支路导纳矩阵。计算支路导纳矩阵的流程图如图4.2所示。
图4.2 支路导纳矩阵计算流程图
将4.1.1节的数据输入到潮流计算的程序后，得到的支路导纳矩阵为：
（二）节点导纳矩阵的计算
得到支路导纳矩阵后，节点导纳矩阵Y就能很容易的求出来了，按照2.2.3节介绍的节点导纳矩阵的求法，可以求得例4.1的节点导纳矩阵为
三、潮流计算
计算出电力网络的节点导纳矩阵后，就可以调用潮流计算子程序进行潮流计算，计算电力网络的节点注入功率，节点电压，线路流动功率和线路损耗功率。
（一）高斯——赛德尔潮流计算
设定高斯——赛德尔潮流计算的精确度epxl为，输入数据为Y、S_PQ、P_PV、U_PV、n、m、s、us、Qmax、Qmin、epxl；输出数据为s_node，u_node、Sij、dSij，经过计算后得出的结果为：
迭代次数40，精确度
表4.2
节点注入功率s_node
节点编号
1
2
3
4
5
节点注入功率s_node
3.3133
+ j1.9998
-1.6000
- j0.8000
-3.7000
- j1.3000
-2.0000
-j 1.0000
5.0000
+ j1.5104
表4.3
节点电压u_node
节点编号
1
2
3
4
5
节点电压u_node
1.0500
0.8670
-j 0.1726
1.0425
- j0.0994
1.0736
+j 0.1317
1.0331
+ j0.1876
i
表4.4 线路流动功率Sij
Sij
j
1
2
3
4
5
1
0
0 + j1.7500
3.3133
+ j1.9998
0 + j1.7500
0 +j 1.7500
2
-0.0000
- j0.1954
0
-0.3531
- j0.4932
-1.2472
- j0.5021
-0.0000
- j0.1954
3
-3.3133
- j1.8664
0.3804
- j1.5387
0
-0.7673
- j1.7589
0.0000
- j1.9321
4
-0.0000
- j4.1222
1.3317
- j3.2877
0.8124
- j4.1261
0
-4.1498
- j1.8304
5
0 + j3.5000
0 + j3.5000
0 + j3.5000
4.1498
+ j1.5108
0
表4.5
线路损耗功率dSij
i
dSij
j
1
2
3
4
5
1
0
-0.0000
+ j1.5546
0 +j 0.1334
-0.0000
- j2.3722
0 + j5.2500
2
-0.0000
+ j1.5546
0
0.0273
- j2.0319
0.0844
- j3.7898
-0.0000
+ j3.3046
3
0 +j 0.1334
0.0273
- j2.0319
0
0.0451
- j5.8850
0.0000
+ j1.5679
4
-0.0000
- j2.3722
0.0844
- j3.7898
0.0451
- j5.8850
0
0 - j0.3196
5
0 + j5.2500
-0.0000
+ j3.3046
0.0000
+ j1.5679
0 - j0.3196
0
（二）牛顿——拉夫逊潮流计算
设定牛顿——拉夫逊潮流计算的精确度epxl为，输入数据为Y、S_PQ、P_PV、U_PV、n、m、s、us、Qmax、Qmin、epxl；其电压初值为4.3.1节中计算得出的节点电压u_node；输出数据为s_node，u_node、Sij、dSij，经过计算后得出的结果为：
迭代次数1，精确度
表4.5
节点注入功率s_node
节点编号
1
2
3
4
5
节点注入功率s_node
3.3133
+ j1.9998
-1.6000
- j0.8000
-3.7000
- j1.3000
-2.0000
-j 1.0000
5.0000
+ j 0.2000
表4.6
节点电压u_node
节点编号
1
2
3
4
5
节点注入功率s_node
1.0500
0.8670
-j 0.1726
1.0425
- j0.0994
1.0736
+j 0.1317
1.0331
+ j0.1876
i
表4.7
线路流动功率Sij
Sij
j
1
2
3
4
5
1
0
0 + j1.7500
3.3133
+ j1.9998
0 + j1.7500
0 +j 1.7500
2
-0.0000
- j0.1954
0
-0.3531
- j0.4932
-1.2472
- j0.5021
-0.0000
- j0.1954
3
-3.3133
- j1.8664
0.3804
- j1.5387
0
-0.7673
- j1.7589
0.0000
- j1.9321
4
-0.0000
- j4.1222
1.3317
- j3.2877
0.8124
- j4.1261
0
-4.1498
- j1.8304
5
0 + j3.5000
0 + j3.5000
0 + j3.5000
4.1498
+ j1.5108
0
表4.8
线路损耗功率dSij
i
dSij
j
1
2
3
4
5
1
0
-0.0000
+ j1.5546
0 +j 0.1334
-0.0000
- j2.3722
0 + j5.2500
2
-0.0000
+ j1.5546
0
0.0273
- j2.0319
0.0844
- j3.7898
-0.0000
+ j3.3046
3
0 +j 0.1334
0.0273
- j2.0319
0
0.0451
- j5.8850
0.0000
+ j1.5679
4
-0.0000
- j2.3722
0.0844
- j3.7898
0.0451
- j5.8850
0
0 - j0.3196
5
0 + j5.2500
-0.0000
+ j3.3046
0.0000
+ j1.5679
0 - j0.3196
0
四、程序设计
程序编写用的是Matlab语言，程序编写的主要工作是主程序的设计与编写，子程序的设计与编写，数据的输入与输出，计算结果的保存。
（一）主程序的设计
主程序主要完成数据的输入与输出，调用各个子程序，协调个子程序的运行。主程序的流程图如图4.2所示。
图4.2 主程序流程图
在步骤1中，需要调用输电线路和变压器数学模型的计算子程序，在步骤3中需要调用高斯——斯德尔潮流计算子程序。输入数据的初步处理包括根据数据获得电网的基本参数，如支路数，将电网参数转化为潮流计算要求的形式。
（二）子程序的设计
子程序主要包括：
（1） 输电线路数学模型计算子程序：line_para；
（2） 变压器数学模型计算子程序：transfer_para；
（3） 支路导纳矩阵计算子程序：Jmatrix；
（4） 节点导纳矩阵计算子程序：Ymatrix；
（5） 牛顿——拉夫逊潮流计算子程序：DYnewton；
（6） 高斯——赛德尔潮流计算子程序：DYgaosai；
（7） 写文件子程序：write；
其中输电线路数学模型计算子程序的算法在2.2.1节中作了详细介绍，变压器数学模型建立子程序的算法在2.2.2节中作了详细介绍。支路导纳矩阵计算子程序的流程图如图4.2所示，牛顿——拉夫逊潮流计算子程序流程图如图3.4所示，高斯——赛德尔潮流计算子程序流程图如图3.2所示。写文件子程序，主要完成计算结果的保存，计算结果保存在excel文件中。在计算完支路导纳矩阵后，节点导纳矩阵Y的计算就很容易了，按照2.2.3节介绍的节点导纳矩阵的求法，可以编写Ymatrix子程序如下：
function Y=Ymatrix(J,n)
%n:网络节点数
%J：支路导纳矩阵
fprintf('节点导纳矩阵计算开始\n')
Y=zeros(n,n);
for i=1:n
for j=1:n
if(i==j)
Y(i,j)=sum(J(i,:));
else
Y(i,j)=-J(i,j);
end
end
end
fprintf('节点导纳矩阵计算结束\n')
以上是所有子程序的设计方法，每一个子程序分别完成了不同的工作，从而将整个潮流计算完成。
（三）数据的输入与输出
输入数据首先卸载excel文件中，其组织形式如4.1.1节所述，计算结果保存在excel文件中，主要内容如4.1.2节所述。在excel中加载Matlab宏，即可在excel中调用Matlab程序，完成计算。
第五章 总结
本文详细的介绍了牛顿——拉夫逊潮流计算方法，并对其中存在的缺点进行了完善。从3.3.3节的算例分析可以看出，牛顿——拉夫逊法具有很好的收敛性，而且计算速度快。特别将牛顿——拉夫逊法和高斯——赛德尔法进行综合运用时，计算的速度和精度能满足更高的要求。
在进行电力系统的其他分析，比如输电线路的损耗，电力网络的无功优化，电力系统的故障分析，在这些分析中，潮流计算也可以得到广泛应用。在暂态类得分析中，电力网络的拓扑结构发生变化，其节点导纳矩阵需要进行重新计算，输入数据也需要进行更新，然后在进行潮流计算，分析计算结果，得出结论。为了加快计算速度，可以对节点导纳矩阵进行部分修改，而不必进行重新计算。
牛顿——拉夫逊法在迭代过程中，如果雅克比矩阵行列式的值比较小，会使得其精确度下降，这时，可以只用高斯——赛德尔，而不用牛顿——拉夫逊法，总之，要充分利用两者的优点，使得计算更迅速准确。也可以采用优化的牛顿——拉夫逊法，来弥补牛顿——拉夫逊法存在的缺点。
参考文献
[1] 陈珩. 电力系统稳态分析 [M]. 南京：南京工学院，1985-05-1.
[2] 王锡凡. 现代电力系统分析 [M]. 北京：科学出版社，2003-3.
[3] 西安交通大学. 电子数字计算机的应用：电力系统计算 [M]. 西安：西安交通大学，1987-10-1.
[4] 南京工学院. 电力系统 [M]. 南京：南京工学院，1980-08-1.
[5] 吴天明，彭彬.
MATLAB电力系统设计与分析 [M]. 国防工业出版社，2004-01.
[6] 林晖. VSC-MTDC和交流混合的电力系统潮流计算 [D]. 湖南大学电气与信息工程学院，2007.4.20.
[7] 梁海平. 智能潮流分析系统的研究与开发 [D]. 华北电力大学电力系统及其自动化专业，2004.12.18.
[8] 许建明. 经典法和遗传算法在电力系统有功优化中的研究 [D]. 南昌大学信息工程学院，2006.6.
[9] 赵芳 陈士方 张圣集. 保留非线性的电力系统概率潮流计算 [J]. 电力情报，2000年03期.
[10] 成海滨 沈茂亚. 电力系统最优潮流算法研究综述 [J]. 电力电气，2006年第25卷第11期.
[11] 罗玮. 电力市场环境下基于最优潮流的输电容量充裕度研究 [D]. 上海交通大学电力系统及其自动化专业，2008.1.1
[12] 郭金伟. 电力市场环境下基于最优潮流的节点实时电价和购电份额研究 [D]. 北京交通大学，2008.
[13] 刘方. 电力系统动态最优潮流的模型与算法研究[D]. 重庆大学电气工程学院，2007.4.1
[14] 郗忠梅 李有安 赵法起 张博. 基于Matlab的电力系统潮流计算[J]. 山东农业大学学报，2010年41卷第2期.
[15] 樊宇璐 李世作 张志斌. 基于拟牛顿法的电力系统潮流计算[J]. 电气开关，2009年第2期.
[16] Yan Sun, Nai-Shan Han.gPower Flow Calculation Method for Islanded Power Network [J]. Department of Electrical Engineering Guangxi University ,Nanning, China,2007.10.3.
[17] Andrey Pazderin, Sergey Yuferev.Power Flow Calculation by Combination of
Newton-Raphson Method and Newton’s Method in Optimization.[J].URAL STATE TECHNICAL UNIVERSITY – UPI,2008.11.3
[18]D. A. Arzamastsev, P. I. Bartolomej, A. M. Holjan, Automatic control system and optimization of modes of electric power systems, Moscow:Vysshaya Shkola, 1983.
[19]V. I. Tarasov, Nonlinear minimization methods for calculating steady states of electric power systems, – Novosibirsk: Nauka, 2001, p. 214.
[20]A. Z. Gamm, Statistical methods of power system state estimation, Moscow: Nauka, 1976, p. 220.
[21]V. M. Gornshtejn, B. P. Miroshnichenko, A.V. Ponomarev, Methods of power system mode optimization, Moscow: Energiya, 1981, p. 336.
附录
附录1 源程序
1 高斯——赛德尔潮流计算源程序
%高斯——赛尔德潮流计算
%[s_node,u_node]=gaosai(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,expl)
%返回值：s_node:
节点功率
%
u_node：
节点电压
%输入值： Y：
导纳矩阵
%
m：
PQ节点数,包含平衡节点
%
n：
总结点数
%
s：
平衡节点编号，一般定为1号编号
%
S_PQ: PQ节点的视在功率[1:m],其中平衡节点的视在功率暂定为0；
%
P_PV：PV节点的有功功率[m+1:n]
%
U_PV: PV节点的电压给定值,只有幅值，没有相位角，是实数[m+1:n]
%
us:
平衡节点的电压给定值，复数。
%
Qmax；PV节点无功功率最大值
%
Qmin：PV节点无功功率最小值
%
epxl: 最大电压变量dumax的最大值，即规定的潮流计算的精度，实数。
function [u_node]=DYgaosai(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,epxl)
%**************************************************************************
%以下为高斯塞尔德潮流计算程序
fprintf('高斯——塞尔德潮流计算开始\n')
U=ones(1,n);%节点电压的初始值；
U(s)=us;%将电压初始值的平衡节点的初始电压设为给定值
k=1000;%设置迭代次数
dumax=0;%最大电压变化量，实数
flag=zeros(1,n-m);%PV节点在计算过程中转化为PQ节点的标志，PV转化为PQ节点时，相应的标志位置1
Q_PV=ones(1,n-m);%PV节点的无功功率；
Q_PV=Q_PV.*0.2;
Ss=0;%平衡节点的视在功率
for a=m+1:n
U(a)=U_PV(a-m);%第a个节点是PV节点，暂时将PV节点电压给定值赋给节点电压初始%值,其相角暂时为0
end
while(k>=1)
for i=1:n
u=U(i);%保存U(i)的初值，为求du和dumax做准备；
if(i==s)%第i个节点是否为平衡节点？
% i=i+1;%d第i个节点是平衡节点，跳过。
else
%第i节点不是平衡节点，进行以下处理
A=0;
for a=1:n
if(a~=i)
A=A+Y(i,a)*U(a);
end
end
a=1;
if(i>=m+1 && flag(i-m)==0) %第i各节点是否为PV节点？
Q_PV(i-m)=-imag(U(i)'*Y(i,i)*U(i)+U(i)'*A);
if(Q_PV(i-m)<=Qmax && Q_PV(i-m)>=Qmin)
Upv=((P_PV(i-m)-Q_PV(i-m)*j)/U(i)-A)/Y(i,i);
ang=angle(Upv);
U(i)=abs(U(i))*(cos(ang)+sin(ang)*j);
else
if(Q_PV(i-m)>Qmax)
Q_PV(i-m)=Qmax;
flag(i-m)=1;
fprintf('PV节点%gQ大于Qmax转化为PQ节点\n',i);
elseif(Q_PV(i-m)<Qmin)
Q_PV(i-m)=Qmin;
flag(i-m)=1;
fprintf('PV节点%gQ小于Qmin转化为PQ节点\n',i);
end
U(i)=((P_PV(i-m)-Q_PV(i-m)*j)/U(i)'-A)/Y(i,i);
end
else
if(i<=m)%第i个节点是PQ节点
U(i)=(S_PQ(i)'/U(i)'-A)/Y(i,i);
elseif(i>=m+1 && flag(i-m)==1)
U(i)=((P_PV(i-m)-Q_PV(i-m)*j)/U(i)'-A)/Y(i,i);
end
end
end
du=abs(U(i)-u);
if(du>dumax)
dumax=du;
end
end
i=1;
if(dumax<epxl)
fprintf('迭代次数%g\n',1000-k+1)
break;
else
k=k-1;
dumax=0;
end
end
if(k==0)
fprintf('高斯——塞尔德潮流计算不能收敛！！！\n')
end
for i=1:n %计算平衡节点的视在功率
Ss=Ss+Y(s,i)'*U(i)';
end
Ss=Ss*U(s);%平衡节点的视在功率
for i=1:n%计算各个节点的视在功率
if(i==s)
s_node(i)=Ss;
else
if(i<=m)
s_node(i)=S_PQ(i);
elseif(i>=m+1)
s_node(i)=P_PV(i-m)+Q_PV(i-m)*j;
end
end
end
s_node
u_node=U
2 牛顿——拉夫逊潮流计算程序
%牛顿——拉夫逊潮流计算。
%直角坐标系
%输出值：s_node：节点视在功率
%
u_node：节点电压
%输入值：Y：
导纳矩阵
%
m：
PQ节点数,包含平衡节点
%
n：
总结点数
%
s：
平衡节点编号，一般定为1号编号
%
S_PQ: PQ节点的视在功率[1:m],其中平衡节点的视在功率暂定为0；
%
P_PV：PV节点的有功功率[m+1:n]
%
U_PV: PV节点的电压给定值,只有幅值，没有相位角，是实数[m+1:n]
%
us:
平衡节点的电压给定值，复数。
%
Qmax；PV节点无功功率最大值
%
Qmin：PV节点无功功率最小值
%
epxl: 最大电压变量dumax的最大值，即规定的潮流计算的精度，实数。
function [s_node,u_node]=DYnewton(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,epxl)
%****************************************************************
%以下为牛顿——拉夫逊法的程序
%以下为牛顿——拉夫逊法的程序
U=DYgaosai(Y,m,n,s,S_PQ,P_PV,U_PV,us,Qmax,Qmin,0.000001);%利用高斯-塞尔德法得到比较准确的节点电压的初值。
fprintf('牛顿——拉夫逊潮流计算开始\n')
dPQV2=zeros(2*(n-1),1);%不平衡量向量；
Ykb=zeros(2*(n-1),2*(n-1));%雅克比矩阵；
I=zeros(1,n);
dU=zeros(2*(n-1),1);%电压不平衡量
dumax=0;
k=1000;%设置迭代次数
Ss=0;%平衡节点的视在功率
S=zeros(1,n);
Q_PV=ones(1,n-m);%PV节点的无功功率；
Q_PV=Q_PV.*0.2;
i=1;
flag=zeros(1,n-m);%PV节点在计算过程中转化为PQ节点的标志，PV转化为PQ节点时，相应的标志位置1
while(k>=1)%迭代次数控制
%计算节点的不平衡量的dPQV2
for i=1:n %i:节点编号；
%*******************************
A=0;
for a=1:n
A=A+(Y(i,a)'*U(a)');
end
a=1;
A=U(i)*A;
%***************************************
if(i<=s-1)%如果i号节点是平衡节点之前的PQ节点
dPQV2(2*i-1)=real(S_PQ(i))-real(A);%PQ节点有功功率不平衡量
dPQV2(2*i)=imag(S_PQ(i))-imag(A);%PQ节点无功功率不平衡量
else
if(i>=s+1 && i<=m)%如果i号节点是平衡节点之后的PQ节点
dPQV2(2*(i-1)-1)=real(S_PQ(i))-real(A);%PQ节点有功功率不平衡量
dPQV2(2*(i-1))=imag(S_PQ(i))-imag(A);%PQ节点无功功率不平衡量
elseif(i>=m+1 && i<=n)%如果i号节点是PV节点
dPQV2(2*(i-1)-1)=P_PV(i-m)-real(A);%PV节点有功功率不平衡量
dPQV2(2*(i-1))=U_PV(i-m)^2-abs(U(i))^2;%PV节点电压平方不平衡量
end
end
end
i=1;
%节点的不平衡量的dPQV2计算结束
for i=1:n%计算节点注入电流I；
if(i<=s-1)
for a=1:n
I(i)=I(i)+Y(i,a)*U(a);
end
a=1;
elseif(i>=s+1)
for a=1:n
I(i)=I(i)+Y(i,a)*U(a);
end
a=1;
end
end
i=1;
for i=1:n
%计算雅克比矩阵
if(i<=s-1)%如果i号节点是平衡节点之前的PQ节点
for a=1:n
if(a~=i)
if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点
Ykb(2*i-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));
%PQ节点Hia参数
Ykb(2*i-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i)); %PQ节点Nia参数
Ykb(2*i,2*a-1)=-Ykb(2*i-1,2*a);%PQ节点的Jia参数=-Nia
Ykb(2*i,2*a)=Ykb(2*i-1,2*a-1);%PQ节点的Lia参数=Hia
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点,Ykb中的a->a-1;i不变
Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数
Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数
Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1));%PQ节点的Jia参数=-Nia
Ykb(2*i,2*(a-1))=Ykb(2*i-1,2*(a-1)-1);%PQ节点的Lia参数=Hia
elseif(a>=m+1 && a<=n)%Ykb中的a->a-1;i不变
Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数
Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数
Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1));%PV节点Jia参数=-Nia
Ykb(2*i-1,2*(a-1))=Ykb(2*i-1,2*(a-1)-1);%PV节点Lia参数=Hia
end
end
elseif(a==i)
if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点
Ykb(2*i-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*i-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*i,2*a-1)=-Ykb(2*i-1,2*a)-2*real(I(i));%PQ节点的Jia参数
Ykb(2*i,2*a)=Ykb(2*i-1,2*a-1)+2*imag(I(i));%PQ节点的Lia参数
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点
Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1))-2*real(I(i));%PQ节点的Jia参
数=-Nia-2real（I(i)）；
Ykb(2*i,2*(a-1))=Ykb(2*i-1,2*(a-1)-1)+2*imag(I(i));%PQ节点的Lia参数=Hia+2imag(I(i));
elseif(a>=m+1 && a<=n)
Ykb(2*i-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*i-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*i,2*(a-1)-1)=-Ykb(2*i-1,2*(a-1))-2*real(I(i));%PQ节点的Jia参数=-Nia-2real（I(i)）；
Ykb(2*i-1,2*(a-1))=Ykb(2*i-1,2*(a-1)-1)+2*imag(I(i));%PQ节点的Lia参数=Hia+2imag(I(i));
end
end
end
end
a=1;
else
if(i>=s+1 && i<=m)%如果i号节点是平衡节点之后的PQ节点;Ybk中的i->i-1;
for a=1:n
if(a~=i)
if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点;Ybk中的i->i-1;a不变
Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数
Ykb(2*(i-1),2*a-1)=-Ykb(2*(i-1)-1,2*a);%PQ节点的Jia参数=-Nia
Ykb(2*(i-1),2*a)=Ykb(2*(i-1)-1,2*a-1);%PQ节点的Lia参数Hia
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点;Ybk中的i->i-1;a->a-1;
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1));%PQ节点的Jia参数=-Nia
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1);%PQ节点的Lia参数=Hia
elseif(a>=m+1 && a<=n)
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1));%PQ节点的Jia参数=-Nia
Ykb(2*(i-1)-1,2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1);%PQ节点的Lia参数
=Hia
end
end
elseif(a==i)
if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点;Ybk中的i->i-1;a不变
Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*a-1)=-Ykb(2*(i-1)-1,2*a);%PQ节点的Jia参数
Ykb(2*(i-1),2*a)=Ykb(2*(i-1)-1,2*a-1)+2*imag(I(i));%PQ节点的Lia参数
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1))-2*real(I(i));%PQ节点的Jia参数=-Nia-2real(I(i))
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1)+2*imag(I(i));%PQ节点的Lia参数=Hia+2*imag(I(i))
elseif(a>=m+1 && a<=n)
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1))-2*real(I(i));%PQ节点的Jia参数=-Nia-2real(I(i))
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1)+2*imag(I(i));%PQ节点的Lia参数=Hia+2*imag(I(i))
end
end
end
end
a=1;
%**************************************************************
else
if(i>=m+1 && flag(i-m)==0)%如果i号节点是PV节点
for a=1:n
if(a~=i)
if(a<=s-1)%如果a号节点在平衡节点之前;Ykb中的i->i-1,a不变
Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PV节点Hia参数
Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PV节点Nia参数
Ykb(2*(i-1),2*a-1)=0;%PV节点的Ria参数
Ykb(2*(i-1),2*a)=0;%PV节点的Sia参数
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点；Ykb中的i->i-1,a->a-1
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PV节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PV节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=0;%PV节点的Ria参数
Ykb(2*(i-1),2*(a-1))=0;%PV节点的Sia参数
elseif(a>=m+1 && a<=n)
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PV节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PV节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=0;%PV节点的Ria参数
Ykb(2*(i-1),2*(a-1))=0;%PV节点的Sia参数
end
end
elseif(a==i)
if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点
Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PV节点Hia参数
Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PV节点Nia参数
Ykb(2*(i-1),2*a-1)=-2*imag(U(i));%PV节点的Ria参数=-2*imag(U(i))
Ykb(2*(i-1),2*a)=-2*real(U(i));%PV节点的Sia参数=-2*real(U(i))
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-2*imag(U(i));%PV节点的Ria参数=-2*imag(U(i))
Ykb(2*(i-1),2*(a-1))=-2*real(U(i));%PV节点的Sia参数=-2*real(U(i))
elseif(a>=m+1 && a<=n)
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-2*imag(U(i));%PV节点的Ria参数=-2*imag(U(i))
Ykb(2*(i-1),2*(a-1))=-2*real(U(i));%PV节点的Sia参数=-2*real(U(i))
end
end
end
end
a=1;
elseif(i>=m+1 && flag(i-m)==1)%如果该PV节点已经转化为了PQ节点。将按照PQ节点的计算方法计算雅克比矩阵
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
for a=1:n
if(a~=i)
if(a<=s-1)%如果a号节点在平衡节点之前;Ykb中的i->i-1,a不变
Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PV-PQ节点Hia参数
Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PV-PQ节点Nia参数
Ykb(2*(i-1),2*a-1)=-Ykb(2*(i-1)-1,2*a);%PV-PQ节点的Jia参数=-Nia
Ykb(2*(i-1),2*a)=Ykb(2*(i-1)-1,2*a-1);%PV节点的Lia参数=Hia
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点；Ykb中的i->i-1,a->a-1
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PV-PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PV-PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1));%PV-PQ节点的Jia参数=-Nia
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1);%PV-PQ节点的Lia参数=Hia
elseif(a>=m+1 && a<=n)
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i));%PV-PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i));%PV-PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1));%PV-PQ节点的Jia参数=-Nia
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1);%PV-PQ节点的Lia参数=Hia
end
end
elseif(a==i)
if(a<=s-1)%如果a号节点是平衡节点之前的PQ节点
Ykb(2*(i-1)-1,2*a-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PV-PQ节点Hia参数
Ykb(2*(i-1)-1,2*a)=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PV-PQ节点Nia参数
Ykb(2*(i-1),2*a-1)=-Ykb(2*(i-1)-1,2*a)-2*real(I(i));%PV-PQ节点的Jia参数=-Nia-2*real(I(i))
Ykb(2*(i-1),2*a)=Ykb(2*(i-1)-1,2*a-1)+2*imag(I(i));%PV-PQ节点的Lia参数=Hia+2*imag(I(i))
else
if(a>=s+1 && a<=m)%如果a号节点是平衡节点之后的PQ节点
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1))-2*real(I(i));%PV-PQ节点的Jia参数=-Nia-2*real(I(i))
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1)+2*imag(I(i));%PV-PQ节点的Lia参数=Hia+2*imag(I(i))
elseif(a>=m+1 && a<=n)
Ykb(2*(i-1)-1,2*(a-1)-1)=imag(Y(i,a))*real(U(i))-real(Y(i,a))*imag(U(i))-imag(I(i));%PQ节点Hia参数
Ykb(2*(i-1)-1,2*(a-1))=-real(Y(i,a))*real(U(i))-imag(Y(i,a))*imag(U(i))-real(I(i));%PQ节点Nia参数
Ykb(2*(i-1),2*(a-1)-1)=-Ykb(2*(i-1)-1,2*(a-1))-2*real(I(i));%PV-PQ节点的Jia参数=-Nia-2*real(I(i))
Ykb(2*(i-1),2*(a-1))=Ykb(2*(i-1)-1,2*(a-1)-1)+2*imag(I(i));%PV-PQ节点的Lia参数=Hia+2*imag(I(i))
end
end
end
end
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%PV节点已经转化为了PQ节点。将按照PQ节点的计算方法计算雅克比矩阵&&&结束！！！
end
%**************************************************************
end
end
end
i=1;
%雅克比矩阵计算结束
dU=Ykb\dPQV2;%计算电压变化量
for i=1:2*(n-1)
if(abs(dU(i))>dumax)
dumax=abs(dU(i));
end
end
i=1;
if(dumax<=epxl)
fprintf('迭代次数%g\n',1000-k+1)
break;
else
for i=1:n%计算新的电压U
if(i<=s-1)
U(i)=U(i)-(dU(2*i)+dU(2*i-1)*j);
elseif(i>=s+1 && i<=n)
U(i)=U(i)-(dU(2*(i-1))+dU(2*(i-1)-1)*j);
end
end
%检验PV节点的无功功率是否超出规定范围，即PV是否要转化成为PQ
for i=m+1:n
A=0;
for a=1:n
if(a~=i)
A=A+Y(i,a)*U(a);
end
end
a=1;
if(flag(i-m)==0) %第i各节点是否为PV节点？
Q_PV(i-m)=-imag(U(i)'*Y(i,i)*U(i)+U(i)'*A);%计算PV节点的Q
if(Q_PV(i-m)<=Qmax && Q_PV(i-m)>=Qmin)%经过计算PV节点的Q在规定范围内
Upv=((P_PV(i-m)-Q_PV(i-m)*j)/U(i)-A)/Y(i,i);%修正PV节点电压的相角，幅值不变
ang=angle(Upv);%修正PV节点电压的相角，幅值不变
U(i)=abs(U(i))*(cos(ang)+sin(ang)*j);%修正PV节点电压的相角，幅值不变
%
Q_PV(i-m)=-imag(U(i)'*Y(i,i)*U(i)+U(i)'*A);
else%经过计算PV节点的Q不在规定范围内，PV节点变为PQ节点
if(Q_PV(i-m)>Qmax)%PV的Q用规定的Qmax代替
Q_PV(i-m)=Qmax;
flag(i-m)=1;%PV转化为PQ
fprintf('PV节点%gQ大于Qmax转化为PQ节点\n',i);
elseif(abs(Q_PV(i-m))<Qmin)%PV的Q用规定的Qmin代替
Q_PV(i-m)=Qmin;
flag(i-m)=1;%PV转化为PQ
fprintf('PV节点%gQ小于Qmin转化为PQ节点\n',i);
end
end
end
end
i=1;
k=k-1;
dumax=0;
end
end
if(k==0)
fprintf('牛顿——拉夫逊潮流计算不能收敛！！！\n')
end
for i=1:n %计算平衡节点的视在功率
Ss=Ss+Y(s,i)'*U(i)';
end
Ss=Ss*U(s);%平衡节点的视在功率
for i=1:n%计算各个节点的视在功率
if(i==s)
S(i)=Ss;
else
if(i<=m)
S(i)=S_PQ(i);
elseif(i>=m+1)
S(i)=P_PV(i-m)+Q_PV(i-m)*j;
end
end
end
i=1;
u_node=U;
s_node=S;
fprintf('牛顿——拉夫逊潮流计算结束，计算结果如下\n')
附录2 英文文献翻译
英文文献
Power Flow Calculation by Combination of Newton-Raphson Method and Newton’s Method in Optimization.
Andrey Pazderin, Sergey Yuferev
URAL STATE TECHNICAL UNIVERSITY – UPI
E-mail: pav@daes.ustu.ru, usv@daes.ustu.ru
Abstract--In this paper, the application of the Newton’s method in optimization for power flow calculation is considered.
Convergence conditions of the suggested method using an example of a three-machine system are investigated. It is shown, that the method allows to calculate non-existent state points and automatically pulls them onto the boundary of power flow existence domain. A combined method which is composed of Newton-Raphson method and Newton’s method in optimization is presented in the paper.
Index Terms—Newton method, Hessian matrix, convergence of numerical methods, steady state stability
Ⅰ.INTRODUCTION
The solution of the power flow problem is the basis on which other problems of managing the operation and development of electrical power systems (EPS) are solved. The complexity of the problem of power flow calculation is attributed to nonlinearity of steady-state equations system and its high dimensionality, which involves iterative methods. The basic problem of the power flow calculation is that of the solution feasibility and iterative process convergence [1].
The desire to find a solution which would be on the boundary of the existence domain when the given nodal capacities are outside the existence domain of the solution, and it is required to pull the state point back onto the feasibility boundary, motivates to develop methods and algorithms for power flow calculation, providing reliable convergence to the solution.
The algorithm for the power flow calculation based on the Newton's method in optimization allows to find a solution for the situation when initial data are outside the existence domain and to pull the operation point onto the feasibility boundary by an optimal path. Also it is possible to estimate a static stability margin by utilizing Newton's method in optimization.
As the algorithm based on the Newton’s method in optimization has considerable computational cost and power control cannot be realized in all nodes, the algorithm based on the combination of the Newton-Raphson methods and the Newton’s method in optimization is offered to be utilized for calculating speed, enhancing the power flow calculation.
II. THEORETICAL BACKGROUND
A．Steady-state equations
The system of steady-state equations, in general, can be expressed as follows:
(1)
where is the vector of parameters given for power flow calculation. In power flow calculation, real and reactive powers are set in each bus except for the slack bus. Ingeneration buses, the modulus of voltage can be fixed. W(X,Y) is the nonlinear vector function of steady-state equations. Variables Y define the quasi-constant parameters associated with an equivalent circuit of an electrical network. X is a required state vector, it defines steady state of EPS. The dimension of the state vector coincides with the number of nonlinear equations of the system (1). There are various known forms of notation of the steady-state equations. Normally, they are nodal-voltage equations in the form of power balance or in the form of current balance. Complex quantities in these equations can be presented in polar or rectangular coordinates, which leads to a sufficiently large variety forms of the steady-state equations notation. There are variable methods of a nonlinear system of steady-state equations solution. They are united by the incremental vector of independent variables
ΔX
being searched and the condition of convergence being assessed at each iteration.
B. The Newton's method in optimization
Another way of solving the problem of power flow calculation is related to defining a zero minimum of objective function of squares sum of discrepancies of steady-state
equations:
(2)
The function minimum (2) is reached at the point where derivatives on all required variables are equal to zero:
(
(3)
It is necessary to solve a nonlinear set of equations (3) to find the solution for the problem
. Calculating the power flow, which is made by the system of the linear equations with a Hessian matrix at each iteration, is referred to as the Newton's
method in optimization [4]:
(
(4)
The Hessian matrix contains two items:
(5)
During the power flow calculation, the determinant of Hessian matrix is positive round zero and negative value of a determinant of Jacobian .This allows to find the state point during the power flow calculation, when initial point has been outside of the existence domain.
The convergence domain of the solution of the Newton's optimization method is limited by a positive value of the Hessian matrix determinant. The iterative process even for a solvable operating point can converge to an incorrect
solution if initial approximation has been outside convergence domain. This allows to estimate a static stability margin of the state and to find the most perilous path of its weighting.
III. INVESTIGATIONS ON THE TEST SCHEME
Convergence of the Newton's method in optimization with a full Hessian matrix has been investigated. Calculations were made based on program MathCAD for a network comprising three buses the parameters of which are presented in Figure 1.Dependant variables were angles of vectors of bus voltage 1 and 2 ，independent variables were capacities in nodes 1 and 2, and absolute values of voltages of nodes 1, 2 and 3 were fixed.
Fig. 1 – The Test scheme
In Figure 2, the boundary of existence domain for a solution of the steady-state is presented in angular coordinates δ1-δ2. This boundary conforms to a positive value of the Jacobian determinant:
(6)
As a result of the power flow calculation based on the Newton method in optimization, the angle values have been received, these values corresponding to the given capacities in Fig.2 (generation is positive and loading is negative).
For the state points which are inside the existence domain, the objective function (2) has been reduced to zero. For the state points which are on the boundary of the existence domain, objective function (2) has not been reduced to zero and the calculated values of capacities differed from the given capacities.
Fig. 2 – Domain of Existence for a Solution
Fig.3 - Boundary of existence domain
In Fig.3, the boundary of the existence domain is presented in coordinates of capacities P1-P2. State points occurring on the boundary of the existence domain (6) have been set by the capacities which were outside the existence domain. As a
result of power flow calculation by minimization (2) based on the Newton's method in optimization, the iterative process converges to the nearest boundary point. It is due to t
he fact that surfaces of the equal level of objective function (2) in coordinates of nodal capacities are proper circles (for threemachine system) having the centre on the point defined by given values of nodal capacities . The graphic interpretation of surfaces of the equal level of objective function for operating point state with 13000 MW loading bus 1 and 15000 MW generating bus 2 is presented in Fig.3.
Hessian matrix is remarkable in its being not singular on the boundary of existence domain. The determinant of a Hessian matrix (5) is positive around zero and a negative value of the Jacobian matrix determinant. This fact allows the power flow to be calculated even for the unstable points which are outside existence domain. The iterative process based on the system of the linear equations (4) solution has converged to the critical stability point within 3-5 iteration. Naturally, the iterative process based on Newton-Rapson method is divergent for such unsolvable operating points.
The convergence domain of the method under consideration has been investigated. What is meant is that not all unsolvable operating points will be pulled onto the
boundary of existence domain. A certain threshold having been exceeded the iterative process has begun to converge to the imaginary solution with angles exceeding 360°.
It is necessary to note that to receive a critical stability operating point in case when initial nodal capacities are set outside the boundary of the existence domain, there is no necessity to make any additional terms as the iterative process converges naturally to the nearest boundary point.
Pulling the operation point onto feasibility boundary is not always possible by the shortest and optimal path. There are a number of constraints, such as impossibility of load (consumption) increase at buses, constraints of generation shedding/gaining at stations. Load following capability of generator units is various, consequently for faster pulling the operation point onto the feasibility boundary it is necessary to
carry out this pulling probably by longer, but faster path.
The algorithm provides possibility of path correction of pulling. It is carried out by using of the weighting coefficients, which define degree of participation of each
node in total control action. For this purpose diagonal matrix A of the weighting coefficients for each node is included into the objective function (2):
(7)
All diagonal elements of the weighting coefficient matrix A should be greater-than zero:
When initial approximation lies into the feasibility domain, coefficients are not influence on the computational process and on the result.
In the figure 4 different paths of the pulling the same operation point onto feasibility boundary depending on the weighting coefficients are presented. Paths are presented for two different operating points.
In tables I and II effect of weighting coefficients on the output computation is presented. In tables I and II k1 and k2 are weighting coefficient for buses 1 and 2, respectively.
TABLE I
WEIGHTING COEFFICIENT EFFECT ON OUTPUT COMPUTATION FOR INITIAL SET CAPACITIES P1= -13000 MW AND P2= 15000 MW
Coefficients
,MW
,MW
,deg
,deg
=1,=1
-7800
9410
-45
55
=5,=1
-8600
8080
-69
25
=0.005,=1
-5700
10140
-1
93
TABLE II
WEIGHTING COEFFICIENT EFFECT ON OUTPUT COMPUTATION FOR INITIAL SET CAPACITIES P1= -8000 MW AND P2= -5000 MW
Coefficients
,MW
,MW
,deg
,deg
=1,=1
-4360
-1680
-92
-80
=0.01,=1
-1050
-4920
-76
-94
=1,=0.35
5800
0
-99
-71
Fig.4 - Paths of pulling the operation point onto the feasibility boundary
IV. COMBINATION OF METHODS
If to compare the Newton’s method in optimization for power flow calculation with newton-Raphson using a Jacobian matrix, the method computational costs on each
iteration will be several times greater as the property of Hessian matrix being filled up by nonzero elements 2.5-3 times greater than with Jacobian one. Each row of Jacobian matrix corresponding to any bus contains nonzero elements corresponding to all incident buses of the scheme. Each row of Hessian matrix contains nonzero elements in the matrix corresponding not only to the neighboring buses, but also their neighbors. However, it is possible to compensate this disadvantage through the combination Newton-
Rap son method with Newton’s method in optimization. It means that the part of nodes can be calculated by conventional Newton method, and the remaining buses will be computed by Newton’s method in optimization. The first group of passive nodes consists of buses in which it is not possible to change
nodal capacity or it is not expedient. Hence, emergency control actions are possible only in a small group of buses supplying with telecontrol. Most of the nodes including
purely transit buses are passive. Active nodes are generating buses in which operating actions are provided. Such approach allows to fix nodal capacity for all passive buses of the scheme which have been calculated by Newton-Rap son method. In active buses which have been calculated by Newton’s method in optimization, deviations from set values of nodal capacity are possible. These deviations can be considered as control action. The power flow calculation algorithm based on combination Newton – Raphson method and Newton’s method in optimization can be presented as follows:
1.The linear equation system with Jacobian matrix is generated for all buses of the scheme.
2. The solution process of the linear equation system with Jacobian is started by utilizing the Gauss method for all passive buses. Factorization of the linear equations system is terminated when all passive buses are eliminated. Factorized
equations are kept.
3.The nodal admittance matrix is generated from not factorized the part of Jacobian matrix corresponding to active buses. This admittance matrix contains parameters of the equivalent network which contains only active buses.
4.The linear equation system with Hessian matrix (4) is generated for the obtained equivalent by Newton’s method in optimization.
5.The linear equation system with Hessian matrix is calculated and changes of independent variables
are defined for active buses.
6.Factorized equations of passive buses are calculated, and changes of independent variables
are defined for passive buses.
7.The vector of independent variables is updated using the changes of independent variables for all buses.
8. New nodal capacities in all buses of the network are defined; constraints are checked; if it necessary, the list of active buses will be corrected.
9. Convergence of the iterative process is checked. If changes of variables are significant, it is necessary to return to item 1.
Taking into account the number of active buses in the network aren’t large, computational costs of such algorithm slightly exceed computational costs of the Newton-Rapson method.
V. CONCLUSION
1. The power flow calculation of an electric network by minimizing the square sum of discrepancies of nodal capacities based on the Newton's method in optimization
materially increases the productivity of deriving a solution for heavy in terms of conditions of stability states and the unstable states outside the existence domain of the solution.
2. During the power flow calculation, the determinant of Hessian matrix is positive around zero and negative value of the Jacobian matrix determinant. The iterative process naturally converges to the nearest marginal state point during the power flow calculation, when the initial operating point has been outside of the existence domain.
3. There is a possibility of control action correction for the pulling operation point onto feasibility boundary by using matrix of weighting coefficients.
4. Utilization of the combined method for power flow calculation allows to use all advantages of Newton’s method in optimization and to provide high calculating speed.
5. In case when the setting nodal powers are outside the existence domain, there are discrepancies in the active buses, which can be considered as control actions for pulling the state point onto the feasibility boundary. When the initial state point is inside the existence domain, the iterative process converges with zero discrepancies for both active and passive buses.
中文翻译
基于优化的牛顿——拉夫逊法和牛顿法的潮流计算
摘要——在本文中，考虑到了优化的牛顿法在潮流计算中的应用。我们将在一个三节点系统的例子中探究这个方法的收敛条件。这个例子表明，这个方法可以计算不存在的状态点并自动将这个点放到功率潮流域的边界上。本文将介绍一个由经过优化的牛顿——拉夫逊法和牛顿法相结合的方法。
关键词——牛顿法，Hessian矩阵，数值方法，静态稳定
1.简介
解决潮流计算问题是解决电力系统运行和发展问题的基础。由于静态方程的非线性特征，以及静态方程的阶数较高，潮流计算往往比较复杂，于是便引进了迭代法。潮流计算的基本问题在于方案的可行性和迭代过程的收敛性。
当给定节点的容量在存在域之外时，我们总是希望能找到一个方法将这个状态点拖回到可行性的边界，研究方法和算法的动机为我们提供了具有可靠收敛性的方案。
当给定的初始值不在存在域中的时候，基于优化牛顿法的潮流计算算法能够找到一条最优路径将工作点移到可行性的边界。它还可以利用优化的牛顿法来估计稳定裕度。
由于优化的牛顿法已经考虑了计算成本和那些不能实现功率控制的节点，这种算法是以优化的牛顿——拉夫逊法和牛顿法相结合的方法为基础的，它可以用来提高计算速度，优化潮流计算。
2. 理论基础
2.1 静态方程
系统的静态方程一般情况下可以如下表示：
（1）
是给定的用来进行潮流计算的向量参数。在潮流计算中，除了平衡节点外，其它节点的有功功率和无功功率都是给定的，PV节点的电压幅值是给定的。W(X,Y)是静态方程组的非线性向量函数。变量Y是一个与电力网络的等值电路相关的准常量参数，
X 是要求的状态向量，它表示电力系统的稳定运行时的状态。这个状态向量的维数与系统非线性方程式（1）的个数一致。现在已知的表示静态方程的形式有多种，一般情况下是以功率平衡或者电流平衡为依据的节点电压方程。这些方程中的复数可以在极坐标或者直角坐标系下表示，这使得静态方程的表示形式多种多样。解非线性系统静态方程的方法有很多,它们都会在每次迭代结束后计算增量向量ΔX，并评估其收敛性。
2.2 优化牛顿法
另外一个解决潮流计算的办法是定义一个静态方程平方和之差零最小的目标函数：
（2）
当所有的有关变量都为零时，函数取得的最小值。
（3）
为了解决问题，我们必须对式（3）的一系列非线性方程进行求解。计算潮流分布涉及到优化牛顿法，每一次迭代计算都是解一组由Hessian矩阵组成的线性方程。
（4）
Hessian矩阵包括两个部分：
（5）
在潮流计算过程中，Hessian矩阵的行列式是一个接近于零的正数或者是雅克比矩阵的行列式的相反数。这样当初始点在存在域之外时可以在潮流计算时找到状态点。
优化牛顿法的收敛于被Hessian矩阵行列式的正值所限制。在迭代过程中，如果初试近似值在收敛域外，为了能找到一个有解的工作点，迭代可能会收敛到一个不正确的解。这能使我们估计静态稳定裕度并找出权重最危险的路径。
3. 试验方案的调查研究
我们已经研究了带有满Hessian矩阵的优化牛顿法的收敛性。如图1所示，我们以一个三个节点的网络为例，以MathCAD程序为工具，进行了计算。不独立的变量是节点1和2的节点电压的相
位角向量，独立变量是节点1和2的容量，节点1、2的电压幅值是给定的。
图1 试验网络
在图2中，稳态解得存在域边界在角坐标系δ1-δ2中表示出来。这个边界与雅克比矩阵行列式的正值相符合：
（6）
基于优化就顿法的潮流计算的结果是，相位角的值被考虑进去，这些值与图2（供电节点是正的，负载节点是负的）中所示的容量相对应。
对于状态点都在存在域内的，目标函数（2）已经减少到零。对于那些状态点在不在存在域内的，目标函数（2）还没有减少到零，计算出来的容量也与给定的容量不同。
图2 一个解的存在域
图3 存在域边界
在图3中，存在域的边界在坐标系P1-P2中表示出来。存在域（6）边界上的状态点受到存在域外容量的限制。以优化牛顿法为基础对最小化函数（2）进行潮流计算导致了迭代过程收敛到最近的边界点上，这是因为在同一水平上的目标函数（2）的边界是一个圆（对于三节点的系统），这个圆的圆心与给定节点容量
的点重合。在图3中，我们对拥有13000MW负载的节点1和15000MW注入功率的节点2工作点的同水平目标函数边界作了图形解释。当存在域边界奇异时，Hessian是值得注意的。Hessian矩阵（5）的行列式是一个接近于0的正数或者是雅克比矩阵行列式的负值，这使得对于那些在存在域之外的不稳定点也能计算其潮流分布。以线性方程组（4）为基础进行迭代计算，经过3到5次迭代后就可收敛到临界稳定点。当然，对于那些没有解的工作点，以牛顿——拉夫逊为基础的迭代过程是发散的。
我们已经研究了这个方法的收敛域。那意味着并不是所有的无解工作点都可以移到存在域的边界上。在迭代过程中如果某一阈值被超过，则迭代过程就会收敛到假象的解决方案，同时角度会超过360°。我们必须注意到，当初始节点容量在存在域边界外时，为了得到一个稳定的工作点，我们没有必要添加任何的额外条件，因为迭代过程会向最近的边界点收敛。
按照最短最优路径将工作点移到可行性边界并不总是可行的，这个过程有许多限制因素。比如节点负载增长的不可能性，发电厂发电量的限制。发电机组的负载是各种各样的，因此为了能够更快地将工作点移动到可行性边界上，必须要通过更长但是更快地路径。
这个算法可以修改移动路径，通过计算权重系数，权重系数表示每个节点在总量控制中参与的程度。为了实现这个目标，表示每个节点权重系数的对角矩阵被包括到了目标函数（2）中。所有权重矩阵A的对角线元素都必须大于零。
当初始近似值在可行域内时，在计算过程中权重矩阵不会对计算过程和计算结果产生影响。
在图4中，我们可以看到在不同的权重系数下将同一个工作点按照不同的路径移到可行性边界。图中表示了两个工作点的移动路径。
在表1和表2中，可以看出权重系数对计算结果的影响。表1和表2中的k1和k2是节点1和2的权重系数。
表1：权重系数对计算输出的影响
初始容量：P1=-13000MW，P2=15000MW
系数
,MW
,MW
,deg
,deg
=1,=1
-7800
9410
-45
55
=5,=1
-8600
8080
-69
25
=0.005,=1
-5700
10140
-1
93
表2：权重系数对计算输出的影响
初始容量：P1=-8000MW，P2=-5000MW
系数
,MW
,MW
,deg
,deg
=1,=1
-4360
-1680
-92
-80
=0.01,=1
-1050
-4920
-76
-94
=1,=0.35
5800
0
-99
-71
图4 工作点移到到可行性边界的路径
4. 组合方法
在潮流计算中如果将优化的牛顿法和利用雅克比矩阵的牛顿——拉夫逊法进行比较，我们会发现每一种方法的计算成本都是利用Hessian矩阵金子那个非零元素填充的方法的2.5到3倍。雅克比矩阵的每一行都与每一个节点对应，它的每一个非零元素都和某一个节点有关系。Hessian矩阵的每一行的非零元素不仅与邻近节点相对应，而且与邻近节点的礼金节点对应。然而，我们可以通过综合应用牛顿——拉夫逊法和优化牛顿法来补偿这些不足。这意味着一部分节点可以用传统的牛顿法来计算，其它的节点可以用优化牛顿法来计算。第一组的节点包括那些节点容量一般不会改变的母线或者不是用来临时改变状态的节点。因此，紧急情况下用来控制电网的那些节点是少数的。大多数的节点包括纯粹的转换母线，容量都是正的。活动节点是发电机母线用来提供对电网的控制。这种方法可以使牛顿法计算过的电网的所有正节点的容量保持不变，而由优化牛顿法计算的活动节点，它们的容量可能与给定值有偏差。这些偏差可以看作是操作输入量。优化牛顿法和牛顿——拉夫逊法组合方法的步骤如下所述：
1. 线性方程的雅克比矩阵由网络的所有节点生成。
2. 在求解带有雅克比矩阵线性方程的过程中必须从对所有正节点利用高斯法开始，当所有的正节点都被删除后，线性方程的因式分解才会停止。因式分解的等式则被保留。
3. 节点导纳矩阵是由与活动节点相对应的未被分解的雅克比矩阵的那部分生成。这个节点导纳矩阵包含只含有活动节点的等值电路的参数。
4. 生成带有Hessian矩阵的的线性方程是为了能从优化牛顿法中得到等式。
5. Hessian矩阵线性方程是为活动节点定义耳朵，它随着独立变量的变化而变化，并被重新计算。
6. 因式分解方程是为正节点定义的，它随着独立变量的变化而变化，并被重新计算。
7. 独立变量向量随着所有节点的变化而更新。
8. 所有节点的容量将被重新计算，限制条件被检验，活动节点列表被修正。
9. 检验迭代过程的收敛性，如果变量变化明显，要回到第一步重新计算。
活动节点的数目不会很大，优化后的算法仅仅比优化前计算成本大一点点。
5. 总结
1.基于优化牛顿法的电力系统潮流计算以最小平方和之差为目标函数，它实质上是拓宽了解决在超出存在域范围的稳态和非稳态条件下解决问题的途径。
2.在潮流计算中，Hessian矩阵行列式是一个接近于零的正数或者是雅克比矩阵的负值。当初始值在存在域外部时，迭代过程自动向最近的区域状态点靠近。
3.通过利用权重系数矩阵可以控制工作点向稳定性边界移动的路径。
4.利用综合方法进行潮流计算可以让我们充分利用优化牛顿法的优点并且提高计算速度。
5.当设定的节点容量在存在域之外时，活动点会存在差别，这个差别可以看作是对移动工作点到可行性边界路径的操作。当初始状态在存在域之内时，所有的活动节点和正节点都会在迭代过程中向零差异收敛。
推荐访问:毕业设计 电力系统 拉夫