基于MATLAB电力系统潮流计算程序设计

来源：尚车旅游网

青岛某高校本科毕业设计（论文）

1 绪论

1．1 潮流计算

1.1.1 潮流计算概述

电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。所以潮流计算是研究电力系统的一种很重要也很基础的计算。

电力系统潮流计算也分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的随时监视及实时控制。

利用计算机进行电力系统潮流计算从 50 年代中期就已经开始。在这 20 年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。对潮流计算的要求可以归纳为下面几点：

（1）计算方法的可靠性或收敛性；（2）对计算机内存量的要求；（3）计算速度；

（4）计算的方便性和灵活性。电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题，其解

法都离不开迭代。因此，对潮流计算方法，首先要求它能可靠地收敛，并给出正确答案。由于电力系统结构及参数的一些特点，并且随着电力系统不断扩大，潮流计算方程式的阶数也越来越高，对这样的方程式并不是任何数学方法都能保证给出正确答案的。这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。

在用数字计算机解电力系统潮流问题的开始阶段，普遍采取以节点导纳矩阵为基础的逐次代入法。这个方法的原理比较简单，要求的数字计算机内存量比较低，适应 50 年代电子计算机制造水平和当时电力系统理论水平。但它的收

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

敛性较差，当系统规模变大时，迭代次数急剧上升，在计算中往往出现迭代不

收敛的情况。这就迫使电力系统计算人员转向以阻抗矩阵为基础的逐次代入法。 60 年代初，数字计算机已发展到第二代，计算机的内存和速度发生了很大的飞跃，从而为阻抗法的采用创造了条件。阻抗法要求的数字计算机储存表征系统接线和参数的阻抗矩阵需要较大的内存量。而且阻抗法每迭代一次都要求顺次取阻抗矩阵中的每一个元素进行运算，因此，每次迭代的运算量很大。这两种情况是过去电子管数字计算机无法适应的。阻抗法改善了系统潮流计算问题的收敛性，解决了导纳法无法求解的一些系统的潮流计算，在 60 年代获得了广泛的应用，曾为我国电力系统设计、运行和研究作出了很大的贡献。目前，我国电力工业中仍有一些单位采用阻抗法计算潮流。阻抗法的主要缺点是占用计算机内存大，每次迭代的计算量大。当系统不断扩大时，这些缺点就更加突出。一个内存 16K 的计算机在采用阻抗法时只能计算 100 个节点以下的系统，32K 内存的计算机也只能计算 150 个节点以下的系统。这样，我国很多电力系统为了采用阻抗法计算潮流就不得不先对系统进行相当的简化工作。

为了克服阻抗法在内存和速度方面的缺点，60 年代中期发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需要存储各个地区系统的阻抗矩阵及它们之间联络线的阻抗，这样不仅大幅度地节省了内存容量，同时也提高了计算速度。克服阻抗法缺点的另一途径是采用牛顿-拉夫逊法。这是数学中解决非线性方程式的典型方法，有较好的收敛性。在解决电力系统潮流计算问题时，是以导纳矩阵为基础的，因此，只要我们能在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿-拉夫逊法潮流程序的效率。自从 60 年代中期，在牛顿-拉夫逊法中利用了最佳顺序消去法以后，牛顿法在收敛性、内存要求、速度方面都超过了阻抗法，成为 60 年代末期以后广泛采用的优秀方法。与此同时，为了保证可靠的收敛，在我国还进行了网流法潮流计算的研究。随着电力系统的日益扩大和复杂化，特别是电力系统逐步实现自动控制的需要，对系统潮流计算在速度、内存以及收敛性方面都提出了更高的要求。70 年代以来，潮流计算方法通过不同的途径继续向前发展，其中比较成功的一个方法就是 P-Q 分解法。这个方法，根据电力系统的退热点，抓住主要矛盾，对纯数学的牛顿法进行了改进，从而在内存容量及计算速度方面都大大向前迈进内了一步。使一个 32K 内存容量的数字计算机可以计算 1000 个节点系统的潮流问题，此法计算速度已能用于在线计算，作系统静态安全监视。目前，我国很多电力系统都采用了 P-Q 分解法潮流

青岛某高校本科毕业设计（论文）

程序。

潮流计算灵活性和方便性的要求，对数字计算机的应用也是一个很关键的问题。过去在很长时间内，电力系统潮流计算是借助于交流台进行的。交流台模拟了电力系统，因此在交流计算台上计算潮流时，计算人员可以随时监视系统各部分运行状态是否满足要求，如发现某些部分运行不合理，则可以立即进行调整。这样，计算的过程就相当于运算人员对系统进行操作、调整的过程，非常直观，物理概念也很清楚。当利用数字计算机进行潮流计算时，就失去了这种直观性。为了弥补这个缺点，潮流程序的编制必须尽可能使计算人员在计算机计算的过程中加强对计算机过程的监视和控制，并便于作各种修改和调整。电力系统潮流计算问题并不是单纯的计算问题，把它当作一个运行方式的调整问题可能更为确切。为了得到一个合理的运行方式，往往需要不断根据计算结果，修改原始数据。在这个意义上，我们在编制潮流计算程序时，对使用的方便性和灵活性必须予以足够的重视。因此，除了要求计算方法尽可能适应各种修改、调整以外，还要注意输入和输出的方便性和灵活性，加强人机联系，以便使计算人员能及时监视计算过程并适当地控制计算的进行。

电力系统潮流计算是电力系统分析中的一种最基本的计算，是对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标是求取电力系统在给定运行状态的计算。即节点电压和功率分布，用以检查系统各元件是否过负荷、各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算的模型和方法有直接影响。实际电力系统的潮流技术主要采用牛顿-拉夫逊法。

在运行方式管理中，潮流是确定电网运行方式的基本出发点；在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流保证了在预想操作情况下电网的潮流分布以及校验运行可靠性。在电力系统调度运行的多个领域都涉及到电网潮流计算。潮流是确定电力网络运行状态的基本因素，潮流问题是研究电力系统稳态问题的基础和前提。

牛顿-拉夫逊法早在 50 年代末就已应用于求解电力系统潮流问题，但作为一种实用的，有竞争力的电力系统潮流计算方法，则是在应用了稀疏矩阵技巧和高斯消去法求修正方程后。牛顿-拉夫逊法是求解非线性代数方程有效的迭代计算。P-Q 分解法进行电力系统分析的潮流计算程序的编制与调试，获得电力

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

系统中各节点电压，为进一步进行电力系统分析作准备。通过本文加深对电力

系统潮流计算原理的理解和计算，初步学会运用计算机知识解决电力系统的问题，掌握潮流计算的过程及其特点。熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。

1.1.2 潮流计算的意义

电力系统潮流计算是电力系统最基本的计算，也是最重要的计算。所谓潮流计算，就是已知电网的接线方式与参数及运行条件，计算电力系统稳态运行各母线电压、各支路电流与功率及网损。对于正在运行的电力系统，通过潮流计算可以判断电网母线电压、支路电流和功率是否越限，如果有越限，就应采取措施，调整运行方式。对于正在规划的电力系统，通过潮流计算，可以为选择电网供电方案和电气设备提供依据。潮流计算还可以为继电保护和自动装置定整计算、电力系统故障计算和稳定计算等提供原始数据。

在运行方式管理中，潮流是确定电网运行方式的基本出发点；在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了多个在预想操作情况下电网的潮流分布以校验运行可靠性。在电力系统调度运行的多个领域都涉及到电网潮流计算。潮流是确定电力网络运行状态的基本因素，潮流问题是研究电力系统稳态问题的基础和前提。

1.2 MATLAB

1.2.1 MATLAB 概述

目前电子计算机已广泛应用于电力系统的分析计算，潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流计算方法有五方面的要求：（1）计算速度快（2）内存需要少（3）计算结果有良好的可靠性和可信性（4）适应性好，亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强（5）简单。

MATLAB 是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB 程序设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率

青岛某高校本科毕业设计（论文）

地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB 与 C 语言和 FORTRAN 语言相比更容易被掌握。通过 M 语言，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。另外，MATLAB 提供了一种特殊的工具：工具箱（TOOLBOXES）.这些工具箱主要包括：信号处理（SIGNAL PROCESSING）、控制系统（CONTROL SYSTEMS）、神经网络（NEURAL NETWORKS）、模糊逻辑(FUZZY LOGIC)、小波(WAVELETS)和模拟（SIMULATION）等等。不同领域、不同层次的用户通过相应工具的学习和应用，可以方便地进行计算、分析及设计工作。MATLAB 设计中，原始数据的填写格式是很关键的一个环节，它与程序使用的方便性和灵活性有着直接的关系。原始数据输入格式的设计，主要应从使用的角度出发，原则是简单明了，便于修改。

1.2.2 MATLAB 的优缺点

MATLAB 语言最大的特点是简单和直接，它主要有如下特点：编程效率

高

MATLAB 是一种面向科学与工程计算的高级语言，允许使用数学形式的语言编写程序，且比 BASIC、FORTRAN 和 C 等语言更加接近我们书写计算公式的思维方式，用 MATLAB 编写程序犹如在演算纸上排列出公式与求解问题。因此，MATLAB 语言也可通俗地称为演算纸式科学算法语言。由于它编写简单，所以编程效率高，易学易懂。

用户使用方便

MATLAB 语言是一种解释执行的语言（在没被专门的工具编译之前），它灵活、方便，其调试程序手段丰富，调试速度快，需要学习时间少。人们用任何一种语言编写程序一般都要经过四个步骤：编辑、编译、链接，以及执行和调试。各个步骤之间是顺序关系，编程的过程就是在它们之间做瀑布型的循环。 MATLAB 语言与其他语言相比，较好的解决了上述问题，把编辑、编译、链接和执行融为一体。它能在同一画面上进行灵活操作，快速排除输入程序中的书写错误、语法错误以至语义错误，从而加快了用户编写、修改和调试程序的速度，可以说在编程和调试过程中它是一种比 VB 还要简单的语言。

具体的说，MATLAB 运行时，如直接在命令行输入 MATLAB 语句（命令），包括调 M 文件的语句，每输入一条语句，就立即对其进行处理，完成编译、链接和运行的全过程。又如，将 MATLAB 源程序编辑为 M 文件，由于 MATLAB

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

磁盘文件也是 M 文件，所以编辑后的源文件就可以直接运行，而不需要进行编

译和链接。在运行 M 文件时，如果有错，计算机屏幕上会给出详细的出错信息，用户经修改后再执行，直到正确为止。所以可以说，MATLAB 语言不仅是一种语言，广义上讲是一种该语言的开发系统，即语言调试系统。

扩充能力强，交互性好高版本的 MATLAB 语言有丰富的库函数，在进行复杂的数序运算时可以直接调用，而且 MATLAB 的库函数同用户文件在形成上一样，所以用户文件也可作为 MATLAB 的库函数来调用。因而，用户可以根据自己的需要方便地建立和扩充新的库函数，以便提高 MATLAB 的使用效率和扩充它的功能。另外，为了充分利用 FORTRAN、C 等语言的资源，包括用户已编好的 FORTRAN、C 语言程序，通过建立 Me 文件的形式，混合编程，方便地调用有关的 FORTRAN、C 语言的子程序，还可以在 C 语言和 FORTRAN 语言中方便地使用 MATLAB 的数值计算功能。这样良好的交互性使程序员可以使用以前编写过的程序，减少重复性工作，也使现在编写的程序具有重复利用的价值。

移植性和开放性很好 MATLAB 是用 C 语言编写的，而 C 语言的可移植性很好。于是 MATLAB 可以很方便地移植到能运行 C 语言的操作平台上。MATLAB 合适的工作平台有：Windows 系列、UNIX、Linux、VMS6.1 和 PowerMac。除了内部函数外， MATLAB 所有的核心文件和工具箱文件都是公开的，都是可读可写的源文件，用户可以通过对源文件的修改和自己编程构成新的工具箱。

语句简单，内涵丰富 MATLAB 语言中最基本最重要的成分是函数，其一般形式为 [a,b,c，...]=fun(d,e,f,...),即一个函数由函数名，输入变量 d,e,f,...和输出变量 a,b,c,... 组成，同一函数名 F，不同数目的输入变量（包括无输入变量）及不同数目的输出变量，代表着不同的含义（有点像面向对象中的多态性）。这不仅使 MATLAB 的库函数功能更丰富，而且大大减少了需要的磁盘空间看，使得 MATLAB 编写的 M 文件简单、短小而高效。

高效方便的矩阵和数组运算 MATLAB 语言像 BASIC、FORTRAN 和 C 语言一样规定了矩阵的算术运算符、关系运算符、罗技运算符、条件运算符，而且这些运算符大部分可以毫无改变地照搬到数组间的运算，有些如算术运算符只要增加“.”就可用于数组间的运算。另外，它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库

青岛某高校本科毕业设计（论文）

函数，使之在求解诸如信号处理、建模、系统识别、控制、优化等领域的问题时，显得大为简捷、高效、方便，这是其他高级语言所不能比拟的。在此基础上，高版本的 MATLAB 一定能名副其实地称为“万能演算纸”式的科学算法语言。

方便的绘图功能

MATLAB 的绘图是十分方便的，它有一系列绘图函数（命令），例如线性坐标、对数坐标、半对数坐标及极坐标，均只需调用不同的绘图函数（命令），在图上标出图题、XY 轴标注，格（栅）绘制也只需调用相应的命令，简单易行。另外，在调用绘图函数时调整自变量可绘出不变颜色的点、线、复线或多重线。这种为科学研究着想的设计是通用的编程语言所不能及的。

MATLAB 的缺点是，它和其他高级程序相比，程序的执行速度较慢。由于 MATLAB 的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢，并且 MA TLAB 的界面功能比较弱, 不能实现交互界面、数据采集和端口操作等功能。

1.3 总结

本文介绍了图形化潮流计算软件的开发设计思想和总体结构，阐述了该软件所具备的功能和特点。结合电力系统的特点，软件采用 MATLAB 语言运行于 WINDOWS 操作系统的图形化潮流计算软件。本系统的主要特点是操作简单，图形界面直观，运行稳定，计算准确。计算中，算法做了一些改进，提高了计算速度，各类的有效封装又使程序具有很好的模块性、可维护性和重用性。

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

2 电力系统的潮流计算

2.1 概述

潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种。自从 20 世纪 50 年代计算机应用于电力系统以来，当时求解潮流的方法是以节点导纳矩阵为基础的逐次代入法（导纳法），后来为解决导纳法的收敛性较差的问题，出现了以阻抗矩阵为基础的逐次代入法（阻抗法），到 20 世纪 60 年代，针对阻抗法占用计算机内存大的问题又出现了分块阻抗法及牛顿 —拉夫逊法。牛顿-拉夫逊法是数学上解非线性方程式的有效方法，有较好的收敛性，将 N—R 法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧，使 N—R 法在收敛性、占用内存、计算机速度方面的优点都超过了阻抗法，成为 20 世纪 60 年代末期以后普遍采用的方法，同时国内外广泛研究了诸如非线性规划法、直流法、交流法等各种不同的潮流计算方法。20 世纪 70 年代以来，又涌现出了更新的潮流计算方法，其中有 1974 年由 B.Stoot 和 O.Alsac 提出的快速分解法以及 1978 年由岩本伸一等提出的保留非线性的高速潮流计算法，其中快速分解法从 1975 年开始已在国内使用，并习惯称之为 PQ 分解法。由于 PQ 分解法在计算速度上大大超过了 N—R 法，不但能应用于离线潮流计算，而且也能应用于在线潮流计算。

2.2 潮流计算的基本方程

2.2.1 节点的分类

用一般的电路理论求解网络方程，目的是给出电压源（或电流源）研究网络内的电流（或电压）分布，作为基础的方程式，一般用线性代数方程式表示。然而在电力系统中，给出发电机或负荷连接母线上电压或电流（都是向量）的情况是很少的，一般是给出发电机母线上发电机的有功功率（P）和母线电压的幅值（U），给出负荷母线上负荷消耗的有功功率（ P）和无功功率（Q）。主要目的是由这些已知量去求电力系统内的各种电气量。所以，根据电力系统中各节点性质的不同，很自然的把节点分成三种类型。

PQ 节点对这一类节点，事先给定的是节点功率（P、Q），待求的未知量是节点电

青岛某高校本科毕业设计（论文）

压向量（U、θ），所以叫“PQ 节点”。通常变电所母线都是 PQ 节点，当某些发电机的输出功率 P、Q 给定时，也作为 PQ 节点。PQ 节点上的发电机称之为 PQ 机（或 PQ 给定型发电机）。在潮流计算中，系统大部分节点属于 PQ 节点。

PU 节点

这类节点给出的参数是该节点的有功功率 P 及电压幅值 U，待求量为该节点的无功功率 Q 及电压向量的相角θ。这类节点在运行中往往要有一定可调节的无功电源，用于维持给定的电压值。通常选择有一定无功功率储备的发电机母线或者变电所有无功补偿设备的母线作 PU 节点处理。PU 节点上的发电机称之为 PU 机。

平衡节点在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在

计算中取该节点电压向量的方向作为参考轴，相当于给定该点电压向量的角度

为零。也就是说，对平衡节点给定的运行参数是 U 和θ，因此又称为 Uθ节点，而待求量是该节点的 P，Q，整个系统的功率平衡由这一节点承担。

关于平衡节点的选择，一般选择系统中担任调频调压的某一发电厂（或发电机），有时也可能按其他原则选择，例如，为提高计算的收敛性，可以选择出线数多或者靠近电网中心的发电厂母线作平衡节点。

以上三类节点 4 个运行参数 P、Q、U、θ中，已知量都是两个，待求量也是两个，只是类型不同而已。

2.2.2 基本方程式

在潮流计算中任何复杂的电力系统都可以归结为以下元件（参数）组成：发电机（注入电流或功率）；

负荷（注入负的电流或功率）；输电线支路（电阻、电抗）；变压器支路（电阻、电抗、变比）；母线上对地支路（阻抗和导纳）；线路上的对地支路（一般为线路充电电容导纳）。集中了以上各种类型的元件的简单网络如图 2.1（a）所示。必须指出，如果仅研究稳态情况下的潮流而不涉及暂态过程的计算，则不需要发电机和负荷的阻抗参数，只需要给出发电机和负荷的注入功率和电流，

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

并且规定发电机和负荷的注入功率和电流取正，而负荷取负。

(a)

潮流计算用的电网结构图

Pg1 jQg1 Pg 2 jQg 2 0 j0 (a) The flow calculation grid chart

I 1 I2 I3

无源线性网络（可用导纳矩阵或阻抗矩阵表示）

 I PR 4 jQR 4 4 

I5 P  jQ

R5R5

I6 PR6 jQR6

(b)

潮流计算等值网络

(b) the trend to calculate the equivalent network

图 2.1 潮流计算用的电网结构图和等值网络 equivalent network

Figure 2.1 The flow calculation grid structure diagram and the 将图 2.1（a）中的发电机和负荷节点用无阻抗线从网络中抽出（为不失一

般性，将既非发电机又非负荷的浮动节点当作零注入功率的母线抽出网络之外），剩下的部分既由接地和不接地支路组成一个无源线性网络如图 2.1（b）所示。

对于这个无源线性网络可用相应的导纳矩阵（或阻抗矩阵）来描述，采用导纳矩阵时，节点注入电流和节点电压构成以下线性方程组

I YU

青岛某高校本科毕业设计（论文）

I1 I2 其中 I In 

可展开为如下形式：

U1 UU 2 

Un 

1,IY U2,..., ) n i  ij j (i

（2.1）

j 1

若 U ZI

可展开为如下形式：

 (iUZ I1,2,..., ) n i  ij j

（2.2）

j 1

式中 n 为网络节点数。由于实际电网中测量的节点注入量一般不是电流而

是功率，因此必须将式（2.1）中的注入电流用节点注入功率来表示。根据电工

理论，节点功率与节点电流之间的关系为

P jQ USI i i i i i （2.3）式中，

Pi PGi PLDi , Qi QGi QLDi

因此用导纳矩阵式（2.1）时，PQ 节点可以表示为

Si PjQi Ii i Ui Ui

把这个关系式代入式（2.1）中，得

Pi jQi Y U (i1,2,..., ) n ij j Ui j 1

（2.4）

（2.5）

比较式（2.1）和（2.5），由于功率代替电流的结果，使式（2.1）电流电压

的线性方程组变化为功率和电压的非线性方程组，这个非线性方程组就是潮流计算的基本方程。

式（2.5）是一组共有 n 个非线性方程组组成的复数方程组，如果把实部和虚部分开便得到 2n 个实数方程，因此由该方程组可解出 2n 个运行参数。但是已经知道每一个节点都有 4 个运行变量，即节点的功率 P 、 Qi ，以及节点电压

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

的幅值和相位（或对应于某一选定参考直角坐标的实部和虚部），记作（ Pi 、Qi 、 Ui 、i ）或（ Pi 、 Qi 、 ei 、 fi ），当节点数为 n 时，则共有 4n 个运行参数。由 2n 个方程式要求出 4n 个运行参数是不可能的，只能求出 2n 个运行参数，而其余 2n 个参数应作为原始数据事先给定。这就得根据节点的分类，将每一个节点的 4 个运行参数中的两个作为原始数据，而将另外两个作为待求量。式（2.5）含有 n 个非线性复数方程，是潮流计算问题的基本方程，对这个方程不同的应用和处理就形成了不同的潮流程序。

若用极坐标形式表示则有：

Pi jQi UiGij jBij U j j Ui Gij jBij U j cosij j sinij （2.6）

ji ji

其中 Pi ， Qi 可分别表示为：

Pi Ui U j Gij cosij jBij sinij 

 ji

i 1, 2n（2.7） 

Qi Ui U j Gij sinij jBij cosij 

ji 

2.2.3 牛顿-拉夫逊

牛顿-拉夫逊算法产生于 50 年代末期，是一种实用且有竞争力的电力系统潮流计算方法。在稀疏矩阵技巧和高斯消去法被应用以后，其真正的价值才体现出来。牛顿—拉夫逊法是求解非线性代数方程有效的迭代计算方法，已经成为求解电力系统潮流问题应用最为广泛的一种方法。方程式（2.7）为用极坐标形式表示的牛顿-拉夫逊潮流方程，将这两个方程改写成残差的形式，即：

Pi =Pis Ui Ui Gij cosij Bij sinij  ji

i 1, 2n（2.8） 

Qi Qis Ui U j Gij sinij Bij cosij 

其中：ij 为节点 i 和 j 之间的电压相角差；Gij 和 Bij 分别为支路电导和电纳； Ui 和U j 分别为节点 i 和 j 的电压向量。对式（2.8）进行泰勒级数展开，取一次项近似，即可得到牛顿法潮流计算的修正方程式，即：

青岛某高校本科毕业设计（论文）

P  P / P / U   H N     

J       （2.9） Q Q / Q / U U / U  M L U / U  U / U 

其中： Pi 和 Qi 为潮流方程的有功功率和无功功率残差向量，共（2n－2）维； U 和 △θ 为母线电压修正量，共（2n－2）维；系数 J 为雅可比矩阵

对方程式（2.9）进行变换即可得到变量 U 和 的求解公式，即：

P 

 J 1 （2.10） U / U  Q   

雅可比矩阵各元素可表达为：

UiU j Gij sinij Bij cosij i j Pi （2.11） H U 2 B Q i j i ii i UiU j Gij cosij Bij sinij i j Pi ij （2.12） N U U 2G  P i j i  ii i 

Q UiU j Gij cosij Bij sinij i j i jM ij   U （2.13） U 2G  P i j 

j ii i 

UiU j Gij sinij Bij cosij i j Pi （2.14） 2

L Uj  Q i j ij U B  i ii i 

根据公式（2.8）至公式（2.13），可得牛顿法潮流计算的具体步骤。 ① 输入原始数据计算节点导纳矩阵；

② 给出各节点电压初值U 0；

③ 将电压初值代入式（2.8），求出 P0Q0。判断是否满足收敛条件，

如果满足，则停止计算。否则，继续进行下面的步骤；

④ 将电压初值代入式（2.11）至式（2.14）中求出雅可比矩阵 J；⑤ 解式

（2.10）中潮流残差方程，求出节点电压的修正量 U 0；

100

U U U ⑥ 修正节点电压向量；

⑦ 返回步骤③继续迭代；

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

⑧ 判断是否满足收敛条件，如果满足，则停止计算，否则，再以 U 0为初值，返回第③步进行下一次迭代。

使用牛顿-拉夫逊法有以下优点:

① 收敛速度快，具有平方收敛特性，迭代次数与系统规模基本无关，若初值选择得较好，一般迭代几次就能收敛；

② 对于有些病态条件的问题，也能利用该方法求解；

③ 应用了稀疏矩阵技巧，所需计算机内存适中。牛顿-拉夫逊法虽是一种广泛使用的方法，但也存在以下缺点。

① 编程比较复杂，且收敛速度的快慢和迭代次数与初始值的好坏有很大的关系，如果初始值好，可以大大减小迭代次数与收敛速度，如果选择不合适有可能永远不收敛；

② 因非对称的雅可比矩阵不是固定的，每次迭代都需要重新计算，大量的求导运算，计算量很大，降低了计算速度。

针对牛顿-拉夫逊法计算速度方面存在的不足和电力系统实现在线控制的要求，在改进牛顿-拉夫逊法的基础上，提出了快速解耦算法。快速解耦算法派生于牛顿-拉夫逊法的极坐标形式，又称为 PQ 分解法。其基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，把有功功率误差作为修正电压向量角度的依据，把无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开进行。它密切地结合了电力系统的固有特点，无论是内存占用量还是计算速度方面都比牛顿-拉夫逊法有了较大的改进。

青岛某高校本科毕业设计（论文）

3 PQ 分解法的潮流计算

3.1 PQ 分解法的基本方程式

20 世纪 60 年代以来 N—R 法曾经是潮流计算中应用比较普遍的方法，但随着网络规模的扩大以及计算机从离线计算向在线计算的发展，N—R 法是在内存需要量及计算速度方面越来越不适应要求。20 世纪 70 年代中期出现的快速分解法比较成功的解决了上述问题，使潮流计算在 N—R 法的基础上向前迈进了一大步，成为取代 N—R 法的算法之一。

快速分解法（又称 PQ 分解法）是从简化牛顿法极坐标形式计算潮流程序的基础上提出来的。它的基本思想是根据电力系统实际运行特点：通常网络上

U 对母线有功功率的明的电抗远大于电阻值，则系统母线电压幅值的微小变化

Q ，因此，节点功率方程在用极坐标形式表示时，它的修正方程式可简显改变化为:

P  H 0   （3.1）    

Q 0 LU /U 

这就是把 2(n-1)阶的线性方程组变成了两个 n-1 阶的线性方程组，将 P 和 Q 分开来进行迭代计算，因而大大地减少了计算工作量。但是 H，L 在迭代过程中仍然在不断的变化，而且又都是不对称的矩阵。对牛顿法的进一步简化，即把式（3.1）中的系数矩阵简化为在迭代过程中不变的对称矩阵。

在一般情况下，线路两端电压的相角ij 是不大的，因此，可以认为



Gij sinij Bij 

部，即

cosij 1

（3.2）

此外，与系统各节点无功功率相应的导纳 BLDj 远远小于该节点自导纳的虚

BLDj i2 Bii

U i Q U 2 B

因而（3.3）

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

考虑到以上关系，式（3.1）的系数矩阵中的各元素可表示为

Hij UiU j Bij (i, j 1, 2,..., n 1) （3.4）

Li, j 1, 2,..., m)

ij UiU j Bij

(而系数矩阵 H 和 L 则可以分别写成：



 U1B11U1 U1B12U2 ... U1B1,n1Un1 

H U2 B21U1 U2 B22U2 ... U2 B2,n1Un1 



Un1Bn1,1U1

Un1Bn1,2U

2 ... Un1Bn1,n1Un1 

U 1

 B11 B

12 ... B1,n1 U1 B21 B 22 ... B2,n1  U2

 U2

   U  n1 Bn1,1 Bn1,2 ... Bn1,n1





U D1B 'U D1

（3.6）

 U1B11U1 U1B12U2 ... U1B1mUm 

L U2 B21U1 U2 B22U2 ... U2 B2mUm  UUm Bm1Um Um Bm2U2 ... m BmmUm 

U 1  B11 B12 ... B1m U1  U2

B21

B 

22 ... B2m U2

  

 

U  m Bm1 Bm2 ... Bmm Un1 

U D 2 B ''U D 2

（3.7）将式（3.6）和（3.7）代入式（3.1）中，得到

PU D1 B 'U D1 

QU D 2 B ''U 

用 U 1 D

和 U 1 1

D 2 分别左乘以上两式，便得

（3.5）



Un1 

青岛某高校本科毕业设计（论文）

U D1 PB 'U D1 U D 2 1 1 （3.8）（3.9）

QB ''U 

这就是简化了的修正方程式，它们也可展开写成

P U1 

1 B B ... B U  

11 12 1,n1 1 P 2  B  U 1 B B ... U 21

2 2 22 2,n1 （3.10）

2          Bn1,1 Bn1,2 ... Bn1,n1 Un1n1 P  n 1  Un1 Q1   U1 B B BU  ... 11 12 1m 1 Q2  B B ... B U U 21 （3.11） 22 2m  2 

2           BBBUQ  m1 m2 ... mm m  m Um 

在这两个修正方程式中系数矩阵元素就是系统导纳矩阵的虚部，因而系数矩阵是对称矩阵，且在迭代过程中保持不变，这就大大减少了计算工作量。用

极坐标表示的节点功率增量为:

ij sBin)0ij  ij Pi Pis Ui U j (Gcoijs

j 1 （3.12）

Qi Qis Ui U j (Gsinij

j 1

ij cBo0ij ijs )

式（3.10）（3.11）（3.12）构成了 PQ 分解法迭代过程的基本方程式。

3.2 计算步骤

)

（1）给定各节点电压的初始值(0)(,0U

i i

Pi （2）代入式（3.12）计算各节点有功功率 Pi ，并求出 Ui

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

（3）解修正方程式（3.10），得出各节点电压相角修正量 i

kk)()() （4）修正各节点电压的相角i : (1

i k i i

Qi

（5）根据式（3.12）求得各节点无功功率误差 Qi ，并求出 Ui （6）求解修正方程式（3.11），得出各节点电压幅值的修正量 Ui

k)()() （7）修正各节点电压的幅值Ui : U (1

U k U k

（8）返回（2）进行迭代，直到各节点功率误差 Pi 及 Qi 都满足收敛条件。

3.3 PQ 分解法程序框图

开始

输入原始数据

形成矩阵 B ' 和 B，并进行三角分解

设 PQ 节点电压初值，各节点电压相角初值

置迭代记数 k=0

k p 1, kQ 1 用式（3.12）计算不平衡功率 P, 计算 Pi ()k ()k()

i / U k i

青岛某高校本科毕业设计（论文）

()k ? max Pi P 置 kP 0

解修正方程式（3.10），求 i ()k kQ 0 ?

k(1k)()() i k i i

置 kQ 1

()k() k 用式（3.12）计算不平衡功率 Q 计算 Q /U i i i ()k

max Q ()k

i ? Q ()k 置 kQ 0

解修正方程式（3.11），求 U i kP 0 ?

k)()() kU i (1Ui k U i

置 kP 1 计算平衡节点功率及全部线路功率

k+1→k

结束

图 3.1 流程图

Figure 3.1 flowchart

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

4 MATLAB 编程及其实例仿真

基于 MATLAB 的 PQ 分解法潮流计算程序 MATLAB 已广泛应用于自动控制、数学运算、信号分析、计算机技术、图像信号处理、财务分析、航天工业和生物医学工程等领域。由于 MATLAB 语言功能强大、人际界面友好、编程效率高、强大而智能化的作图功能，且具有编程语句简洁、灵活、表达和运算能力强等显著特点。

4.1 程序清单

%本程序的功能是用 PQ 分解法进行潮流计算 n=input(‘请输入节点数：n=’); nl==input(‘请输入支路数：nl==’); isb==input(‘请输入平衡母线节点号：isb=’); pr==input(‘请输入误差精度:pr=’); B1==input(‘请输入由支路参数形成的矩阵:B1=’); B2==input(‘请输入由节点参数形成的矩阵:B2=’); X==input(‘请输入由节点号及其对地阻抗形成的矩阵:X=’); na==input(‘请输入 PQ 节点号:na=’);

Y==zeros(n);YI= zeros(n);e= zeros(1,n);f= zeros(1,n);V== zeros(1,n);O= zeros(1,n);

for i=1:n if X(i,2)~=0; p=x(i,1); Y(p,p)=1./X(i,2); end end for i=1:nl if B1(i,6)= =0 p=B1(i,1);q=B1(i,2); else p= B1(i,2);q=B1(i,1); end

青岛某高校本科毕业设计（论文）

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); YI(p,q)= YI(p,q) -1./B1(i,3); Y(q,p) =Y(p,q); YI(q,p) = YI(p,q);

Y(q,q)= Y(q,q)+1./(B1(i,3)* B1(i,5)~2)+B1(i,4)./2; YI(q,q)= YI(q,q) +1./B1(i,3);

Y(p,p)= Y(p,p) +1./B1(i,3)+ B1(i,4)./2; YI(p,p)= YI(p,p) +1./B1(i,3); end %求导纳矩阵 G=real(Y);B=imag(YI);BI=imag(Y); for i=1:n

S(i)=B2(i,1)-B2(i,2); BI(i,i)= BI(i,i)+ B2(i,5); end

P=real(S);Q=imag(S); for i=1:n e(i)=real(B2(i,3)); f(i)=imag(B2(i,3)); V(i)=B2(i,4); end for i=1:n if B2(i,6)= =2

V(i)=sqrt(e(i)~2+f(i)~2); O(i)=atan(f(i)./e(i)); end end for i=2:n if i= =n B(i,i)=1./B(i,i);

else IC1=i+1; for j1=IC1:n

B(i,j1)= B(i,j1)./B(i,i);

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

end

B(i,i)=1./B(i,i); for k=i+1:n

for j1=i+1:n

B(k,j1)=B(k,j1)-B(k,i)*B(i,j1); end end end end p=0;q=0; for i=1:n if B2(i,6)= =2 p=p+1;k=0; for j1=1:n if B2(i1,6)= =2

k=k+1; A(p,k)=BI(i,j1); end end end end for i=1:na if i= =na

A(i,i)=1./A(i,i); else k=i+1;

for j1=k:na

A(i,j1)= A(i,j1)./A(i,i); end

A(i,i)=1. /A(i,i); for k=i+1:na for j1=i+1:na

A(k,j1)= A(k,j1)-A(k,i)*A(i,j1);

青岛某高校本科毕业设计（论文）

end end end

end

ICT2=1;ICT1=0;kp=1,kq=1;K=1;DET=0;ICT3=1; while ICT2~=0|ICT3~=0 ICT2=0;ICT3=0; for i=1:n if i~=isb C(i)=0; for k=1:n

C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k))); end

DP1(i)=P(i)-V(i)*C(i); DP(i)=DP1(i)./V(i); DET=abs(DP1(i)); if DET>=pr ICT2=ICT2+1; end end end

Np(K)=ICT2;

if ICT2~=0 for i=2:n

DP(i)=B(i,i)*DP(i); if i~ =n IC1=i+1; for k=ic1:n

DP(k)=DP(k)-B(k,i)*DP(i); end else for LZ=3:i

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

L=i+3-LZ; IC4=L-1;

for MZ=2:IC4 I=IC4+2-MZ;

DP(I)=DP(I)-B(I,L)*DP(L);

end end end

end for i=2:n

O(i)=O(i)-DP(i); end kq=1;L=0; for i=1:n if B2(i,6)= =2 C(i)=0;L=L+1; for k=1:n

C(i)= C(i) +V(k)*(G(i,k)*sin(O(i)-O(k))+BI(i,k)*cos(O(i)-O(k))); end

DQ1(i)=Q(i)-V(i)*C(i); DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i)); if DET>=pr ICT3=ICT3+1; end end end end else kp=0; if kq~=0; L=0; for i=1:n

青岛某高校本科毕业设计（论文）

if B2(i,6)= =2 C(i)=0;L=L+1; for k=1:n

C(i)= C(i) +V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k))); end

DQ1(i)=Q(i)-V(i)*C(i); DQ(L)= DQ1(i)./V(i); DET=abs(DQ1(i)); end end end end

Nq(K)=ICT3; if ICT3~=0 L=0; for i=1:na DQ(i)=A(i,i)*DQ(i); if i= =na for LZ=2:i

L=i+2-LZ; IC4=L-1;

for MZ=1:IC4 I=IC4+1-MZ;

DQ(I)=DQ(I)-A(I,L)*DQ(L); end end else IC1=i+1; for k=IC1:na

DQ(k)=DQ(k)-A(k,i)*DQ(i); end end

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

end L=0; for i=1:n

if B2(i,6)= =2 L=L+1; V(i)=V(i)-DQ(L); end end kp=1; K=K+1; else

kq=0; if kp~=0 K=K+1; end end for i=1:n Dy(K-1,i)=V(i); end end

disp(‘迭代次数‘)；

disp（K）； disp（‘每次没有达到精度要求的有功功率个数为‘）； disp（NP）； disp（‘每次没有达到精度要求的无功功率个数为‘）； disp（Nq）;

for k=1:n

E(k)=V(k)*cos(O(K))+V(k)*sin(O(k))*j; O(k)=O(k)*180./pi;

end

disp(‘各节点的电压标么值 E 为（节点号从小到大排）：‘); disp(E);

青岛某高校本科毕业设计（论文）

disp(‘各节点的电压 U 大小（节点号从小到大排）：‘)为；

disp(U);

disp(‘各节点的电压相角 O（节点号从小到大排）：‘)为； disp(O); for p=1:n C(p)=0;

for q=1:n

C(p)= C(p)+conj(Y(p,q))*conj(E(q)); end

S(p)=E(p)*C(p);

end

disp(‘各节点的功率 S（节点号从小到大排）：‘)为； disp（S）； disp(‘各条支路的首端功率 Si 为（顺序同您输入 B1 时的一样）：‘）为；

for i=1:nl if B1(i,6)= =0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*con j(1./(B1(i,3)*B1(i,5))));

disp(Si(p,q)); end

disp(‘各条支路的末端功率 Sj（顺序同您输入 B1 时的一样）：‘）为； for i=1:nl if B1(i,6)= =0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))-conj(E(p)))*con j(1./(B1(i,3)*B1(i,

5))));

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

disp(Sj(q,p)); end

disp(‘各条支路的功率损耗 DS（顺序同您输入 B1 时的一样）：‘）为； for i=1:nl if B1(i,6)= =0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end

DS(i)=Si(p,q)+Sj(q,p); disp(DS(i)); end for i=1:K Cs(i)=i; for j=1:n Dy(K,j)=Dy(K-1,j); end end

disp(‘以下是每次迭代后各节点的电压值（如图所示）‘)；

plot(Cs,Dy),xlable(‘迭代次数‘)，ylable(‘电压‘),title(‘电压迭代次数曲线’)；

4.2 设计验证

现以一个例子来说明该软件得使用方法: 电力系统接线，如 4.1 所示,试求潮流分布

（1）n 为节点数、nl 为支路数、isb 为平衡母线节点号（固定为 1）、pr 为误差精度。

（2）输入由支路参数形成的矩阵 B1 矩阵 B1 的每行是由下列参数构成的： ① 某支路的首端号 P； ② 末端号 Q，且 P 28

青岛某高校本科毕业设计（论文）

⑥ 折算到哪一侧的标志（如果支路的首端 P 处于高压侧则请输入“1”，否则请输入“0”）。

P5 5

1:1.05 j0.015

0.08 j0.30 1.05 :1

U5 1.05 j0.25 j0.25 2 j1 0.04 j0.25 j0.25 j0.03 U1 1.05 1 0

0.1j0.35

j0.25 3.7 j1.3

1.6 j0.8 图 4.1 电路图 Figure 4.1 schematic

（3）请输入各节点参数形成的矩阵 B2 矩阵 B2 的每行是由下列参数构成的： ① 节点所接发电机的功率 SG； ② 节点负荷的功率 SL； ③ 节点电压的初始值；

④ PU 节点电压 U 的给定值；

⑤ 节点所接的无功补偿设备的容量； ⑥ 节点分类标号 igl。输入数据：请输入节点数：n=5 请输入支路数：nl=5 请输入平衡母线节点号：isb=1 请输入误差精度：pr=0.00001

请输入由支路参数形成的矩阵：B1=[1 2 0.03i 0 1.050;2 3 0.08+0.3i 0.5i 1 0;2

4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1]

请输入由节点参数形成的矩阵：B2=[0 0 1.05 1.05 0 1;0 3.7+1.3i 1.05 0 0 2;0

2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0 0 2;5 0 1.05 1.05 0 3]

请输入由节点号及其对地阻抗形成的矩阵：X=[1 0;2 0;3 0;4 0;5 0] 请输入 PQ 节点数 na=3

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

输入如图 4.2 所示:

4.2 输入框图

4.2 Input Block Diagram 在

输入需要的数据后,点击运算,出现结果如下：结果：迭代次数

每次没有达到精度要求的有功功率个数为

4 4 4 4 4 4 4 4 每次没有达到精度要求的无功功率个数为

3 0

3 3 3 3 3 3 3 3 1 0 各节点的电压标么值 E（节点号从小到大排）为：

1.0500 1.0335-0.0774i 1.0260+0.3305i 0.8592-0.0718i

0.9746+0.3907i 各节点的电压 U 大小（节点号从小到大排）为： 1.0500

1.0364

1.0779

0.8662

1.0500

青岛某高校本科毕业设计（论文）

各节点的电压相角 O（节点号从小到大排）为：

0 -3.519 17.8535 -4.7785 21.8433 各节点的功率 S（节点号从小到大排）为：

2.5794+2.2994i -3.7000-1.3000i -2.0000-1.0000i -1.6000-0.8000i

5.0000+1.8131i

各条支路的首端功率 Si（顺序同您输入 B1 时的一样）为：

2.5794+2.2994i -1.2774+0.2032i 0.1568+0.4713i 1.5845+0.6725i

5.0000+1.8131i 各条支路的末端功率 Sj（顺序同您输入 B1 时的一样）为：

-2.8794-1.9745i 1.4154-0.2443i

-0.1338-0.3909i -1.4662-0.4091i

-5.0000-1.4282i 各条支路的功率损耗 DS（顺序同您输入 B1 时的一样）为：

0.0000+0.3249i 0.1381-0.0412i 0.0230+0.0804i 0.1184+0.2635i

-0.0000+0.3849i

点击计算后软件运行情况如图 4.3 所示：

基于 MATLAB 的电力系统 P-Q 潮流计算程序设计

图 4.3 仿真图

Figure 4.3 Simulation diagram

经过手动计算，得 1 至 5 号节点功率分别为：

2.5794+2.2994i -3.7000-1.3000i -2.0000-1.0000i -1.6000-0.8000i 5.0000+1.8131i

举例验算结论：通过手动运算，发现运算结果与手动结果运算结果一样，

故软件计算正确。

青岛某高校本科毕业设计（论文）

5 总结

通过学习潮流计算原理参阅一些参考资料及文献，对传统潮流死算有所理解，对一些新的潮流计算技术有了整体的认识。

潮流计算是电力系统的一项基本计算，它根据给定的运行条件及系统接线情况确定系统的运行状态、母线的电压，各元件中流过的功率，系统的功率损耗等，进而比较运行方式的合理性、可靠性和经济性。

本文分析和比较了电力系统潮流计算的各种方法，其中 PQ 分解法具有计算速度快、所占内存少等优势被广泛采用。本设计首先进行了 PQ 分解法的理论研究，然后应用 MATLAB 软件编制 PQ 分解法潮流程序，并通过一个电力系统实例对所编程序进行验证，得出结论。

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文