牛顿拉斐逊法潮流计算

1、攀枝花学院本科专业综合实训电气工程与自动化专业综合实训一学生姓名： 王 佳 唐 倩 学生学号：201210502066 201210502058院（系）： 电气信息工程学院 年级专业： 2012级电气工程与自动化 指导教师： 方重秋 副教授 二一五年五月摘 要潮流计算是电力系统非常重要的分析计算，用以研究系统规划和运行中提出的各种问题。对规划中的电力系统，通过潮流计算可以检验所提出的电力系统规划方案能否满足各种运行方式的要求；对运行中的电力系统，通过潮流计算可以预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以内，系统中各种元件(线路、变压器等)是否会

2、出现过负荷，以及可能出现过负荷时应事先采取哪些预防措施等。潮流计算是电力系统分析最基本的计算。除它自身的重要作用之外，在电力系统分析综合程序(PSASP)中，潮流计算还是网损计算、静态安全分析、暂态稳定计算、小干扰静态稳定计算、短路计算、静态和动态等值计算的基础。传统的潮流计算程序缺乏图形用户界面，结果显示不直接难与其他分析功能集成。网络原始数据输入工作大量且易于出错。本文采用MATLAB语言运行WINDOWS操作系统的潮流计算软件。而采用MATLAB界面直观，运行稳定，计算准确。关键词：电力系统潮流计算；牛顿拉夫逊法潮流计算；MATLAB目 录1.概述1.1 设计目的与要求.31.1.1 设

3、计目的.31.1.2 设计要求.31.2 设计题目.31.3 设计内容.32.电力系统潮流计算概述2.1 电力系统简介.32.2 潮流计算简介.2.3 潮流计算的意义及其发展. .3.潮流计算设计题目3.1 潮流计算题目.3.2 对课题的分析及求解思路.4.潮流计算算法及手工计算4.1 变压器的型等值电路.74.2 节点电压方程.94.3 节点导纳矩阵.114.4 导纳矩阵在潮流计算中的应用.124.5 潮流计算的手工计算.135.Matlab概述5.1 Matlab简介.175.2 Matlab的应用.175.3 矩阵的运算.185.3.1 与常数的运算.185.3.2 基本数学运算.185

4、.3.3 逻辑关系运算.185.4 Matlab中的一些命令.186.潮流计算流程图及源程序6.1 潮流计算流程图.226.2 潮流计算源程序图.236.3 运行计算结果.30总结参考文献1.1 设计目的与要求1.1.1设计目的1. 掌握电力系统潮流计算的基本原理；2. 掌握并能熟练运用一门计算机语言（MATLAB语言或C语言或C+语言）；3. 采用计算机语言对潮流计算进行计算机编程。1.1.2 设计要求1. 程序源代码；2. 给定题目的输入，输出文件；3. 程序说明；4. 给定系统的程序计算过程；5. 给定系统的手算过程（至少迭代2次）。1.2设计题目电力系统潮流计算（牛顿-拉夫逊法、P-Q

5、分解法）1.3设计内容1. 根据电力系统网络推导电力网络数学模型，写出节点导纳矩阵；2. 赋予各节点电压变量（直角坐标系形式）初值后，求解不平衡量；3. 形成雅可比矩阵；4. 求解修正量后，重新修改初值，从2开始重新循环计算；5. 求解的电压变量达到所要求的精度时，再计算各支路功率分布、功率损耗和平衡节点功率；6. 上机编程调试；7. 计算分析给定系统潮流分析并与手工计算结果做比较分析；8. 书写课程设计说明书，准备答辩。2.1电力系统叙述 电力工业发展初期，电能是直接在用户附近的发电站（或称发电厂）中生产的，各发电站孤立运行。随着工农业生产和城市的发展，电能的需要量迅速增加，而热能资源和水能

6、资源丰富的地区又往往远离用电比较集中的城市和工矿区，为了解决这个矛盾，就需要在动力资源丰富的地区建立大型发电站，然后将电能远距离输送给电力用户。同时，为了提高供电的可靠性以及资源利用的综合经济性，又把许多分散的各种形式的发电站，通过送电线路和变电所联系起来。这种由发电机、升压和降压变电所，送电线路以及用电设备有机连接起来的整体，即称为电力系统。现代电力系统提出了“灵活交流输电和新型直流输电”的概念。灵活交流输电技术是指运用固态电子器件与现代自动控制技术对交流电网的电压、相位角、阻抗、功率以及电路的通断进行实时闭环控制，从而提高高压输电线路的诉讼能力和电力系统的稳态水平。新型直流输电技术是指应用

7、现电力电子技术的最新成果，改善和简化变流站的造价等。运营方式管理中，潮流是确定电网运行方式的基本出发点：在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了电网在预想操作预想下的电网的潮流分布以及校验运行的可靠性。在电力系统调度运行的多个领域都涉及到电网潮流计算。潮流是确定电力网咯运行状态的基本因素，潮流问题是研究电力系统稳态问题的基础和前提。1.2潮流计算简介电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压。各元件中流过的功率，系统的功率损耗等等。在电力系统规划的设计和现有电力

8、系统运行方式的研究中，都需要利用潮流计算来定量的分析比较供电方案或运行方式的合理性。可靠性和经济性。此外，电力系统的潮流计算也是计算机系统动态稳定和静态稳定的基础，所以潮流计算是研究电力系统的一种和重要和基础的计算。电力系统潮流计算也分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。利用电子数字计算机进行潮流计算从50年代中期就已经开始了。在这20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的，对潮流计算的要求可以归纳为以下几点：（1） 计算方法的可靠性或收敛性；（2） 对计算机内

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

10、。运行方式管理中，潮流是确定电网运行方式的基本出发点；在规划领域，需要进行潮流分析验证规划方案的合理性；在实时运行环境，调度员潮流提供了多个在预想操作情况下电网的潮流分布以及校验运行可靠性。在电力系统调度运行的多个领域问题是研究电力系统稳态问题的基础和前提。在用数字见算计算机解电力系统潮流问题的开始阶段，普遍采取以节点导纳矩阵为基础的逐次代入法。这个方法的原理比较简单，要求的数字计算机内存量比较差下，适应50年代电子计算机制造水平和当时电力系统理论水平，但它的收敛性较差，当系统规模变大时，迭代次数急剧上升，在计算中往往出现迭代不收敛的情况。这就迫使电力系统的计算人员转向以阻抗矩阵为基础的逐次代

11、入法。阻抗法改善了系统潮流计算问题的收敛性，解决了导纳无法求解的一些系统的潮流计算，在60年代获得了广泛的应用，阻抗法德主要缺点是占用计算机内存大，每次迭代的计算量大。当系统不断扩大时，这些缺点就更加突出，为了克服这些缺点，60年代中期发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需要存储各个地区系统的阻抗矩阵及它们之间联络的阻抗，这样不仅大幅度的节省了内存容量，同时也提高了计算速度。克服阻抗法缺点是另一个途径是采用牛顿-拉夫逊法。这是数学中解决非线性方程式的典型方法，有较好的收敛性。在解决电力系统潮流计算问题时，是以导纳矩阵为基础的，因此，只要我

12、们能在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿法潮流程序的效率。自从60年代中期，牛顿法中利用了最佳顺序消去法以后，牛顿法在收敛性。内存要求。速度方面都超过了阻抗法，成为了60年代末期以后广泛采用的优秀方法。3.1 潮流计算题目牛顿拉斐逊法计算潮流参数如下：G1、G2：SN=35MVA，VN=10.5kV，X=0.35；T1: SN=31.5MVA，Vs%=10， k=10.5/121kV,Ps=180kW, Po=33kW,Io%=0.8；YN/d-11T2: SN=31.5MVA，Vs%=10， k=10.5/121kV,Ps=170kW, Po=33kW,Io%=0

13、.8；YN/d-11L1:线路长60km，电阻0.25/km，电抗0.41/km，对地容纳3.00×10-6S/km；L2:线路长100km，电阻0.15/km，电抗0.38/km，对地容纳2.78×10-6S/km；L3:线路长80km，电阻0.18/km，电抗0.4/km，对地容纳2.78×10-6S/km；负荷：S3=50MVA，功率因数均为0.85。G1为平衡节点；G2按额定功率工作， PN、UN(PV节点)3.2 对课题的分析及求解思路此电力系统是一个5节点，5支路的电力网络。采用牛顿拉夫逊法的极坐标法进行求解。4.1变压器的型等值电路在电力系统潮流计算

14、中，往往要计算节点导纳矩阵，而我们计算节点导纳矩阵采用节点电压法来实现，如在变压器构成的电力系统中，需要将变压器模型转变成变压器型等值电路（见图1-1），在利用电路知识列节点电压方程，从而导出所需的导纳矩阵。图3-1双绕组变压器的型等值电路（i，j为节点）而在电力系统潮流计算中一般采用标幺值进行计算，标幺值公式如下：如果采用标么值计算，元件参数都应归算到同一基准值时得标么值，才能在同一个等值电路上分析和计算。所以，变压器转变成型等值电路时，我们采用标幺值计算，使所求参数为变压器变比k的函数。而在一个已经归算好的电力系统网中，若改变变压器的分接头来进行调压，这时变压器的等值电路参数也会相应得改变

15、，此时采用型等值电路进行折算就显得较为方便。下面是变压器的型等值电路分析过程： 如不计励磁支路的影响，双绕组变压器可用其阻抗与一个理想变压器串联的电路表示，如图所示。理想变压器只有一个参数，那就是变比k=。现以变压器阻抗按实际变比归算到低压侧的情况为例，推导出双绕组变压器的型等值电路。流入和流出理想变压器的功率相等： (、分别为变压器高、低绕组的实际电压） （3-1） （3-2）联立（3-1）、（3-2）两个公式解得： （3-3） （3-4）根据电路原理节点1、2的节点电流方程具有如下形式： （3-5）将式（3-3）、（3-4）与式（3-5）比较得(3-6): (3-6)因此可以的得到各支路导

16、纳为 （3-7） 4.2 节点电压方程在电路中我们学过利用节点电压方程求解某几条支路的电流，现以下图3-2-1与图3-2-2为例推导节点电压方程组。 图3-2-1节点电压法为例 图3-2-2用电流源代替电压源为例 图3-2-1表示了一个具有两个电源和你一个等值负荷的系统。、为电源电势，、为电源的内部导纳，为负荷的等值导纳，、为各支路的导纳。如果以地为电压参考点，设节点1、2、3的电压为，根据基尔霍夫电流KCL法对节点1、2、3列节点电流方程得式（3-8）： （3-8）上式中左端为节点1、2、3流出的电流，右端为注入个节点的电流。由上式可以得到一个等效的等值电路图3-2-2。图3-2-2中利用了

17、电流源代替的电压源。在图1-2-2中可知的式（3-9）： （3-9） 为等值电流源向网络注入的电流。将与式（3-8）联立得式（3-10）： (3-10）上式中称为节点1、2、3的自导纳，称为相应节点之间的互导纳。 因此，在一般情况下，在电力网络中有n个节点，则可以按式（3-10）的形式列出n个节点方程式，也可用矩阵的形式表示。其中 分别为节点注入电流列向量及节点电压列向量；为节点导纳矩阵，其中对角元素为节点i的自导纳，非对角线为节点i与节点j之间的互导纳。 4.3 节点导纳矩阵节点导纳矩阵既可根据自导纳和互导纳的定义直接求取，也可根据电路知识中找出改网络的关联矩阵，在节点电压方程的矩阵形式进行

18、求解。本章节我们主要讨论的是直接求解导纳矩阵。根据节点电压方程章节我们知道，在利用电子数字计算机计算电力系统运行情况是，多采用形式的节点方程式。其中阶数等于电力网络的节点数。从而可以得到n个节点时的节点导纳矩阵方程组（3-11）如下： (3-11） 由此可以得到n个节点导纳矩阵：它反映了网络的参数及接线情况，因此导纳矩阵可以看成是对电力网络电气特性的一种数学抽象。由导纳短阵所联系的节点方程式是电力网络广泛应用的一种数学模型。 通过上面的讨论，可以看出节点导纳矩阵的有以下特点： （1）导纳矩阵的元素很容易根据网络接线图和支路参数直观地求得，形成节点导纳矩阵的程序比较简单。 （2）导纳矩阵为对称矩

19、阵。由网络的互易特性易知。（3）导纳矩阵是稀疏矩阵。它的对角线元素一般不为零，但在非对角线元素中则存在不少零元素。在电力系统的接线图中，一般每个节点与平均不超过34个其他节点有直接的支路连接。因此，在导纳矩阵的非对角线元素中每行仅有34个非零元素，其余的都是零元素，而且网络的规模越大，这种现象越显著。导纳矩阵的对称性和稀疏性对于应用计算机求解电力系统问题有很大的影响。如果能充分地利用这两个特点，如在程序设计中储存导纳矩阵的对角元素和上三角元素（或下三角元素），排除零元素的储存和运算，就可以大大地节省储存单元和提高计算速度。节点导纳矩阵的形式可归纳如下：（1）导纳矩阵的阶数等于电力网络的节点数。

20、（2）导纳矩阵各行非对角元素中非零元素的个数等于对应节点所连得不接地支路数。（3）导纳矩阵各对角元素，即节点的自导纳等于相应节点之间的支路导纳之和。（4）导纳矩阵非对角元素，即节点之间的互导纳等于相应节点之间的支路导纳的负值。而在电力系统中进行潮流计算时，往往要计算不同接线下的运行状况，例如，改变变压器主抽头时，潮流分布也随之变化，以及改变其他设备参数进行计算潮流分布，此时就需要导出变化时的导纳矩阵就需要对所设计的程序进行参数设定，而不需要重复上述步骤去导出所求的导纳矩阵。4.4 导纳矩阵在潮流计算中的应用导纳矩阵在潮流计算中的应用起到重要的作用，前面我们介绍了根据系统网络的接线盒参数形成节点

21、导纳矩阵的方法。尽管形成节点导纳矩阵的原理是简单的，但如果采用手算的方法，即使节点数不多的系统也仍然有相当大的工作量。因此只有应用计算机才能快速而准确地完成这些计算任务。本章节我们介绍形成系统节点导纳矩阵的实用程序。为了形成节点导纳矩阵，必须知道电力系统的接线图。从前面的讨论知道，网络接线由节点及连接两个节点的支路确定。实际上，只有输入了各支路两端的节点号，就相当输入了系统的接线图。在计算潮流分布时，我们必须先导出该网络的导纳矩阵，而进行潮流计算时解非线性的节点电压方程的有关方法中，高斯-塞德尔迭代法和牛顿-拉夫逊迭代法是计算机潮流计算中常用的基本方法。这两种方法既可用以解线性方程组，也课用以

22、解非线性方程组。高斯-塞德尔迭代法由于其简单而在早期的潮流计算程序中得以采用。但嗣后就逐渐被牛顿型算法所取代。目前这种方法多半与牛顿型算法配合使用以弥补后者的不足。牛顿-拉夫逊法的收敛性较好，但对初值的要求比较严格，是当前广泛采用的计算机潮流算法。运用计算机进行潮流计算，一般要完成以下几个步骤：建立数学模型、确立计算方法、制定计算机流程并编制程序、上机计算及对计算结果进行分析。因此我们可以知道导纳矩阵在潮流计算中是很重要的。本节只是对导纳矩阵在潮流分布的计算机算法一些简单的描述，我们将在下一章对其进行详细讲解。4.5 潮流计算的手工计算在计算电力系统网络的潮流分布时，我们需要把变压器转化成变压

23、器的型等值电路来进行计算器等效导纳，以下是变压器的型等值电路（见图3-1），变压器的型等值电路（见图3-2），其等效导纳的计算过程如下：变压器的型等值电路0.95:1 图3-1所求变压器参数数据如下：=Y=-=j3.3158y=变压器的型等值电路1:0.95 图3-2所求变压器参数数据如下：=Y=-=j3.3158Y=所求电力系统网络的各个节点的相关参数如下：潮流计算极坐标形式计算如下：（1）各个节点电压及功率的标幺值为： (2)各节点的功率方程： 带入各数据得： （3）按泰勒级数展开，可得到修正方程式为： 其中由 得= = =由以上各向量解得= =由已知进行第一次迭代可得： 已知PV节点为节

24、点5，求节点5得无功：由公式 式中P为PV节点，节点1为平很节点带入已知的数据 =-2.83845.1 Matlab简介MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB和Mathematica、Maple并称为三

25、大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 5.2 Matlab的应用MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且mathwork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C+ ，JAVA的支持。可以

26、直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。MATLAB 产品族可以用来进行以下各种工作： 数值分析 数值和符号计算 工程与科学绘图 控制系统的设计与仿真 数字图像处理技术 数字信号处理技术 通讯系统设计与仿真 财务与金融工程 MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了 MATLAB 环境，以解决这些应用领域内特定类型的问题。 5.

27、3.1 与常数的运算 常数与矩阵的运算即是同该矩阵的每一元素进行计算。但需注意进行数除时，常数通常只能做除量。 基本函数运算中，矩阵的函数运算是矩阵预算中最实用的部分，常用的主要有以下几个： det(a) 求矩阵a的行列式 eig(a) 求矩阵a的特征值 inv(a)或a(-1) 求矩阵a的逆矩阵 rank(a) 求矩阵a的秩 trace(a) 求矩阵a的迹（对角线元素之和） 我们进行工程计算时常常遇到矩阵对应元素之间的运算。这种运算不同于前面讲的数学运算，为有所区别，我们称之为数组运算。5.3.2 基本数学运算 数组的加、减与矩阵的加、减运算完全相同。而乘除法运算有相当大的区别，数组的乘除法

28、是指两同维数组对应元素之间的乘除法，它们的运算符为“.*”和“./”或“.。”前面讲过常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定，数组与常数之间的除法运算没有任何限制。 另外，矩阵的数组运算中还有幂运算（运算符.）、指数运算（exp）、对数运算(log)、和开方运算（sqrt）、等，有了“对应元素”的规定，数组的运算实质上就是针对数组内部的每个元素进行的。矩阵的幂运算与数组的幂运算有很大的区别。5.3.4. 逻辑关系运算 逻辑运算是MATLAB中数组运算所特有的一种运算形式，也是几乎所有的高级语言普遍适用的一种运算。5.4. Matlab中的一些命令1）一般MAT

29、LAB命令格式为 输出参数1，输出参数2，=（命令名）（输入参数1，输入参数2，） 输出参数用方括号，输入参数用圆括号如果输出参数只有一个可不使用括号。 2）可用、键来重现已输入的数据或命令。用、键来移动光标进行修改。3）所有MATLAB命令都用小写字母。大写字母和小写字母分别表示不同的变量。4）常用数有特定的名字，如pi（=3.141596） 、Inf（=）、NaN则表示不定型求得的结果（如0/0）。5）矩阵的输入要一行一行的进行，每行各元素用空格或（，）分开，每行用（；）分开。 6）MATLAB书写格式为A=1 2 3 ；4 5 6 ；7 8 9 在MATLAB中运行如下程序可得到A矩阵

30、a=1 2 3;4 5 6;7 8 9 a = 1 2 3 4 5 6 7 8 97）需要显示命令的计算结果时，则语句后面不加“；”号，否则要加“；”号。 运行下面两种格式可以看出他们的区别 a=1 2 3;4 5 6;7 8 9 a=1 2 3;4 5 6;7 8 9; a = （不显示计算结果） 1 2 3 4 5 68）当输入语句过长需要换行时，应加上“”后再回车，则可连续输入。9）diary命令使用该命令可以在窗口中以ASCII码形式记录所有的输入和输出。但这个命令不是存储数据，而是存储输入与屏幕上输出的内容。它可以记录下工作的过程。在每个工作过程之前使用该命令，工作结束后使用diar

31、y off 则能将整个工作过程记录下来。 格式diary(文件名)（扩展名）. diary off一般来说扩展名可取，m这样就可在MATLABBIN目录下存入该文件。 10）save 命令该命令存储定义的变量或演算结果，也可以用来存储指定的变量。 命令格式为 save 文件名 . 扩展名 11）what 命令 该命令可以在当目录下显示MATLAB文件和MAT数据文件12）dir 命令 显示当前目录下的所有文件.13）clear命令 14）d1,d2,d3,.=size(a) 求矩阵的大小，对m*n二维矩阵，第一个为行数m,第二个为列数n。如果输入 calear a b c ，则表示清除工作空间

32、中指定变量a,b,c；如果仅仅输入calear命令，则清除整个工作空间。 与此同时，MATLAB具有强大的矩阵运算功能, 但由于我们在求节点导纳矩阵时用的不多, 因此这里我们只作简单介绍。1) 在MATLAB中表示一个矢量要用方括号， 而列矢量的输入只需在行矢量输入格式基础上加转置符()即可。 如 x=1 2 3;4 5 6 x = 1 2 3 4 5 6 而 x=1 2 3;4 5 6' (加转置符) x = 1 4 2 5 3 6注意上面两式的区别。2) 下面三条命令可以产生一个行矢量1 a=linspace(x,y,n) 2 a=logspace(x,y,n) 3 a=x:n:y

33、第一条命令可以在线性空间产生一个值在10x至10y之间间隔点数为n的行矢量(一组数据)。第二条命令可以在对数空间产生一个值在x至y之间等间隔的行矢量(一组数据)。其行矢量的起始值是x, 终值为y, 点数为n。第三条命令产生X至y步长为n的行矢量。但是, 三个命令之间存在差别，下面的例子可以说明这一点。 例一 x=logspace(0,5,6) x = 1 10 100 1000 10000 100000例二 x=linspace(0,10,11) x = 0 1 2 3 4 5 6 7 8 9 10例三 x=0:1:10 x = 0 1 2 3 4 5 6 7 8 9 10 通过上面三个例子可

34、以看出例一, 例二中n代表选取的点数。而在例三 中n则表示步长. 我们应当注意它们的区别。3) 矩阵的加, 减, 乘, 除等, 和其它语言书写一样。但要注意的是在运算符 前面加有(.)则表示是元素对元素的操作. 4)以下是常用的运算命令运算命令名 功能Angle求复数的角Min求最小值Max求最大值Sum求和Roots求多项式的根Poly由多项式的根求多项式的系数Polyval求给定点多项式的值Polyder多项式求导在进行潮流分布计算时，实际上是由多个简单系统构成的复杂系统，在求节点导纳矩阵时的要用到反馈的指令，所以在MATLAB中有下面几种命令可以解决两个系统间的连接问题。1）系统的并联p

35、arallel 命令可以实现两个系统的并联。示意图如下：系统1系统2 u1 y1 u + y u2 y2 并联后的系统传递函数表示式为： 其中n1、d1和n2、d2分别为g1（s）、g2（s）的传递函数分子、分母系数行矢量。命令格式： n，d = paralltl（n1，d1，n2，d2）a，b，c，d = paralltl（a1，b1，c1，d1，a2，b2，c2，d2）2）系统的串联series命令实现两个系统的串联，示意图如下：g2（s）g1（s） u1 y u2 y2串联后系统的传递函数为 命令格式： n，d = series（n1，d1，n2，d2） a，b，c，d = series

36、（a1，b1，c1，d1，a2，b2，c2，d2）3）系统的反馈feedback命令实现两个系统的反馈连接，示意图如下： u1 + y1 ±g2（s）g1（s） y2 u2连接后系统的传递函数表示为： 命令格式： n，d = feedback（n1，d1，n2，d2）或： n，d = feedback（n1，d1，n2，d2，sign） a，b，c，d = feedback（a1，b1，c1，d1，a2，b2，c2，d2，sign）其中sign是指示y2到u1连接的符号，缺省时默认为负（即sign = -1）。4）系统的闭环cloop命令可以将系统的输出反馈到系统的输入构成闭环系统，

37、示意图如下：g1（s）u y ±正、负反馈后闭环系统为： 命令格式： n，d = cloop （n1，d1，sign） ac，bc，cc，dc = cloop（a，b，c，d，sign）通过以上对MATLAB基本指令的了解，我们就可以对所求的电力系统网络的节点导纳矩阵进行画编程框架图。6.1 潮流计算流程图在上一章节，我们对MATLAB的基本指令有了初步的了解，我们就可以画出程序编程框图，按照程序编程框图，我们可以通过MATLAB软件进行编程，实现我们所求的节点导纳矩阵。程序编程框图如下：启动输入原始数据形成节点导纳矩阵YB给定电压初值e(0)、f(0)置K=0利用式（4-46）（4

38、-47）计算P（K）、Q（K）、U2（K）计算线路功率及平衡节点功率|P（K）、Q（K）|<？输出利用式（4-49）（4-50）求雅可比矩阵各元素解修正式（4-48）求e(K)、f(K)修正各节点电压e(K+1)= e(K) +e(K)f(K+1) =f(K) +f(K)e(K+1)= e(K)、f(K+1) =f(K) K+K+16.2 潮流计算源程序图disp('电力系统极坐标下的牛顿-拉夫逊法潮流计算:');clearn=input('请输入结点数：n=');n1=input('请输入PV结点数：n1=');n2=input('

39、;请输入PQ结点数：n2=');isb=input('请输入平衡结点：isb=');pr=input('请输入精确度：pr=');K=input('请输入变比矩阵:K=');C=input('请输入支路阻抗矩阵：C=');y=input('请输入支路导纳矩阵：y=');U=input('请输入结点电压矩阵：U=');S=input('请输入各结点的功率：S=');Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros(1

40、,n1+n2);for m=1:n for R=1:n C(m,m)=C(m,m)+y(m,R); if K(m,R)=0 C(m,m)=C(m,m)+1/(C(m,R) /( K(m,R) * (K(m,R)-1) ; C(R,R)=C(R,R)+1/(C(m,R)/(1-K(m,R); C(m,R)=C(m,R)/K(m,R); C(R,m)=C(m,R); endendendfor m=1:n for R=1:n if m=R Z(m)=Z(m)+1/C(m,R); end endendfor m=1:n for R=1:n if m=R Y(m,m)=C(m,m)+Z(m); else

41、 Y(m,R)=-1/C(m,R); end endenddisp('结点导纳矩阵:');disp(Y);disp('迭代中的雅克比矩阵:');G=real(Y);B=imag(Y);O=angle(U);U1=abs(U);k=0;PR=1;P=real(S);Q=imag(S);while PR>pr for m=1:n2 UD(m)=U1(m); end for m=1:n1+n2 for R=1:n PT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R)+B(m,R)*sin(O(m)-O(R); end PT1(m)=su

42、m(PT); PP(m)=P(m)-PT1(m); PP1(k+1,m)=PP(m); end for m=1:n2 for R=1:n QT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R)-B(m,R)*cos(O(m)-O(R); end QT1(m)=sum(QT); QQ(m)=Q(m)-QT1(m); QQ1(k+1,m)=QQ(m); end PR1=max(abs(PP); PR2=max(abs(QQ); PR=max(PR1,PR2); for m=1:n1+n2 for R=1:n1+n2 if m=R H(m,m)=U1(m)2*B(m,m)+Q

43、T1(m); else H(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R)-B(m,R)*cos(O(m)-O(R); end end end for m=1:n1+n2 for R=1:n2 if m=R N(m,m)=-U1(m)2*G(m,m)-PT1(m); else N(m,R)=-U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R)+B(m,R)*sin(O(m)-O(R); end end end for m=1:n2 for R=1:n1+n2 if m=R J(m,m)=U1(m)2*G(m,m)-PT1(m); else J(m,R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R)+B(m,R)*sin(O(m)-O(R); end end end for m=1:n2 for R=1:n2 if m=R L(m,m)=U1(m)2*B(m,m)-QT1(m); else L(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R)-B(m,R)*cos(O(m)-O(R); end end end JJ=H N;J L; disp(JJ); PQ=PP'QQ' DA=-inv(JJ)*PQ; DA1=DA' for m=