电力系统潮流分析计算的MATLAB仿真-周明亮01

第一篇：电力系统潮流分析计算的MATLAB仿真-周明亮01

成都理工大学工程技术学院毕业论文

电力系统潮流分析计算的MATLAB仿真

作者姓名：周明亮

专业名称：电气工程及其自动化

指导教师：赵

熹讲师电力系统潮流分析计算的MATLAB仿真

摘要

传统的潮流计算程序缺乏图形用户界面，结果显示不直观，难于与其他分析功能集成。网络原始数据输入工作量大且易于出错。随着计算机技术的飞速发展，MICROSOFT WINDOWS操作系统早已被大家所熟悉，其友好的图形用户界面已成为PC机的标准，而DOS操作系统下的应用程序因其界面不够友好，开发具有WINDOWS风格界面的电力系统分析软件已成为当前的主流趋势。另外，传统的程序设计方法是结构化程序设计方法，该方法基于功能分解，把整个软件工程看作是一个个对象的组合，由于对某个特定问题域来说，该对象组成基本不变，因此，这种基于对象分解方法设计的软件结构上比较稳定，易于维护和扩充。

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

关键词：电力系统

潮流仿真计算

牛顿—拉夫逊法潮流计算

MATLAB

-I电力系统潮流分析计算的MATLAB仿真

摘要............................................................................................................I Abstract.....................................................................................................II 目录.........................................................................................................III 前言...........................................................................................................1 1电力系统潮流计算概述........................................................................3

1.1电力系统概述..............................................................................3

1.2潮流计算介绍..............................................................................3 1.3国内用得较多的几种潮流计算软件简介..................................4 2潮流计算的数学模型............................................................................6

2.1导纳矩阵的原理及计算方法......................................................6 2.1.1自导纳和互导纳的确定方法................................................6 2.1.2 节点导纳矩阵的性质及意义...............................................7 2.1.3 非标准变比变压器等值电路...............................................9 2.2潮流计算的基本方程................................................................11 2.3 电力系统节点分类...................................................................14 2.4潮流计算的约束条件................................................................15 3牛顿－拉夫逊法概述..........................................................................17

3.1牛顿-拉夫逊法基本原理...........................................................17 3.2牛顿--拉夫逊法潮流求解过程.................................................18 3.3 牛顿—拉夫逊法的程序框图...................................................23 4潮流计算程序的实现..........................................................................25

4.1MATLAB软件简介....................................................................25 4.2矩阵的运算................................................................................26 4.3牛顿—拉夫逊法潮流计算实例................................................27 结束语.....................................................................................................34

-III电力系统潮流分析计算的MATLAB仿真

前言

潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况下确定系统稳态运行状态的一种基本方法，是电力系统规划和运营中不可缺少的一个重要组成部分。可以说，它是电力系统分析中最基本、最重要的计算，是系统安全、经济分析和实时控制与调度的基础。是电力系统研究人员长期研究的一个课题。MATLAB自1980年问世以来，它的强大的矩阵处理功能给电力系统的分析、计算带来许多方便。在处理潮流计算时，其计算机软件的速度已无法满足大电网模拟和实时控制的仿真要求，而高效的潮流问题相关软件的研究已成为大规模电力系统仿真计算的关键。随着计算机技术的不断发展和成熟，对MATLAB潮流计算的研究为快速、详细地解决大电网的计算问题开辟了新思路。

电力系统潮流计算是电力系统分析中的一种最基本的计算，是对复杂电力系统正常和故障条件下稳态运行状态的计算。潮流计算的目标是求取电力系统在给定运行状态的计算。即节点电压和功率分布，用以检查系统各元件是否过负荷。各点电压是否满足要求，功率的分布和分配是否合理以及功率损耗等。对现有电力系统的运行和扩建，对新的电力系统进行规划设计以及对电力系统进行静态和暂态稳定分析都是以潮流计算为基础。潮流计算结果可用如电力系统稳态研究，安全估计或最优潮流等对潮流计算的模型和方法有直接影响。在用数字计算机解电力系统潮流问题的开始阶段，普遍采取以节点导纳矩阵为基础的逐次代入法。这个方法的原理比较简单，要求的数字计算机内存量比较下，适应50年代电子计算机制造水平和当时电力系统理论水平。但它的收敛性较差，当系统规模变大时，迭代次数急剧上升，在计算中往往出现迭代不收敛的情况。这就迫使电力系统计算人员转向以阻抗矩阵为基础的逐次代入法。阻抗法改善了系统潮流计算问题的收敛性，解决了导纳法无法求解的一些系统的潮流计算，在60年代获得了广泛的应用。阻抗法的主要缺点是占用计算机内存大，每次迭代的计算量大。当系统不断扩大时，这些缺点就更加突出。为

-1-电力系统潮流分析计算的MATLAB仿真

了克服阻抗法在内存和速度方面的缺点，60年代中期发展了以阻抗矩阵为基础的分块阻抗法。这个方法把一个大系统分割为几个小的地区系统，在计算机内只需要存储各个地区系统的阻抗矩阵及它们之间联络线的阻抗，这样不仅大幅度地节省了内存容量，同时也提高了计算速度。克服阻抗法缺点的另一途径是采用牛顿-拉夫逊法。这是数学中解决非线性方程式的典型方法，有较好的收敛性。在解决电力系统潮流计算问题时，是以导纳矩阵为基础的，因此，只要我们能在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿法潮流程序的效率。自从60年代中期，在牛顿法中利用了最佳顺序消去法以后，牛顿法在收敛性。内存要求。速度方面都超过了阻抗法，成为60年代末期以后广泛采用的优秀方法。

-2-电力系统潮流分析计算的MATLAB仿真

1电力系统潮流计算概述

1.1电力系统概述

电力工业发展初期，电能是直接在用户附近的发电站(或称发电厂)中生产的，各发电站孤立运行。随着工农业生产和城市的发展，电能的需要量迅速增加，而热能资源(如煤田)和水能资源丰富的地区又往往远离用电比较集中的城市和工矿区，为了解决这个矛盾，就需要在动力资源丰富的地区建立大型发电站，然后将电能远距离输送给电力用户。同时，为了提高供电可靠性以及资源利用的综合经济性，又把许多分散的各种形式的发电站，通过送电线路和变电所联系起来。这种由发电机、升压和降压变电所，送电线路以及用电设备有机连接起来的整体，即称为电力系统。

现代电力系统提出了“灵活交流输电与新型直流输电”的概念。灵活交流输电技术是指运用固态电子器件与现代自动控制技术对交流电网的电压、相位角、阻抗、功率以及电路的通断进行实时闭环控制，从而提高高压输电线路的输送能力和电力系统的稳定水平。新型直流输电技术是指应用现电力电子技术的最新成果，改善和简化变流站的造价等。

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

1.2潮流计算介绍

电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行

-3-电力系统潮流分析计算的MATLAB仿真

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

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

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

(1)计算方法的可靠性或收敛性。(2)对计算机内存量的要求。(3)计算速度。

(4)计算的方便性和灵活性。

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

1.3国内用得较多的几种潮流计算软件简介

(1)BPA潮流计算程序

美国帮涅维尔电力局(BPA，Bonneville Power Administr-ation)开发，被中国电力科学院引进吸收，从1984年开始在中国得到推广应用。程序提供两种潮流计算方法：P_Q分解法和牛顿法。

(2)PSASP 潮流计算程序

中国电力科学院开发。程序提供五种潮流计算方法： P-Q分解

-4-电力系统潮流分析计算的MATLAB仿真

法、牛顿法(功率式)、最佳乘子法、牛顿法(电流式)、P-Q分解法转牛顿法(电流式)。

(3)PSS/E 潮流计算程序

美国PTI开发，70年代推向市场，目前已有40个国家200多家公司应用该程序。提供5种潮流计算方法：牛顿法、解耦牛顿法、快速牛顿法、高斯－塞德尔法、改进的高斯－塞德尔法。

-5-电力系统潮流分析计算的MATLAB仿真

2潮流计算的数学模型

2.1导纳矩阵的原理及计算方法 2.1.1自导纳和互导纳的确定方法

电力网络的节点电压方程：

IBYBUB

(2-1)式(2-1)IB为节点注入电流列向量，注入电流有正有负，注入网络的电流为正，流出网络的电流为负。根据这一规定，电源节点的注入电流为正，负荷节点为负。既无电源又无负荷的联络节点为零，带有地方负荷的电源节点为二者代数之和。

式(2-1)UB为节点电压列向量，由于节点电压是对称于参考节点而言的，因而需先选定参考节点。在电力系统中一般以地为参考节点。如整个网络无接地支路，则需要选定某一节点为参考。设网络中节点数为(不含参考节点)，则IB，UB均为n*n列向量。YB为n*n阶节点导纳矩阵。

节电导纳矩阵的节点电压方程：

IBYBUB 展开为：

I1Y11YI221IY313IYnn1YBY12Y22Y32Yn2Y13Y23Y33Yn3Y1nU1Y2nU2Y3nU3

(2-2)YnnUn是一个n*n阶节点导纳矩阵，其阶数就等于网络中除参考节点

-6-电力系统潮流分析计算的MATLAB仿真

外的节点数。节点导纳矩阵的对角元素Yii(i=1，2，n)成为自导纳。自导纳数Yii值上就等于在i节点施加单位电压，其他节点全部接地时，经节点i注入网络的电流，因此，它可以定义为：

/U(U0,ji)

(2-3)YiiIiij节点i的自导纳Yii数值上就等于与节点直接连接的所有支路导纳的总和。

节点导纳矩阵的非对角元素Yij(j =1，2，…，n;i =1，2，…，n;j = i)称互导纳，由此

可得互导纳Yij数值上就等于在节点i施加单位电压，其他节点全部接地时，经节点j注入网络的电流，因此可定义为：

(U0,jIj/iU

Yji

(2-4)ij节点j，i之间的互导纳Yij数值上就等于连接节点j，i支路到导纳的负值。显然，恒Yij等于Yji。互导纳的这些性质决定了节点导纳矩阵是一个对称稀疏矩阵。而且，由于每个节点所连接的支路数总有一个限度，随着网络中节点数的增加非

零元素相对愈来愈少，节点导纳矩阵的稀疏度，即零元素数与总元素的比值就愈来愈高。

2.1.2 节点导纳矩阵的性质及意义

节点导纳矩阵的性质：

(1)YB为对称矩阵，Yij=Yji。如网络中含有源元件，如移相变压器，则对称性不再成立。

(2)YB对无接地支路的节点，其所在行列的元素之和均为零，即

nnj1Yi,j0,i1Yj,i0。对于有接地支路的节点，其所在行列的元素之和等于该点接地支路的导纳。利用这一性质，可以检验所形成节点导纳矩阵的正确性。

(3)YB具有强对角性：对角元素的值不小于同一行或同一列中任

-7-电力系统潮流分析计算的MATLAB仿真

一元素。

(4)YB为稀疏矩阵，因节点i，j 之间无支路直接相连时Yij=0，这种情况在实际电力系统中非常普遍。矩阵的稀疏性用稀疏度表示，其定义为矩阵中的零元素与全部元素之比，即 SZ/n2，式中Z 为YB中的零元素。S 随节点数n 的增加而增加：n=50，S可达92%；n=100，S 可达90%；n=500，S可达99%，充分利用节点导纳矩阵的稀疏性可节省计算机内存，加快计算速度，这种技巧称为稀疏技术。

节点导纳矩阵的意义：

是n*n阶方阵，其对角元素 Yii(i=1，2，----n)称为自导纳，非对角元素Yij(i，j=1，2，n，ij)称为互导纳。将节点电压方程YBIBYBUB展开为：

I1Y11I2Y21InYn1可见，/U(UYiiIiijY12YY22n2UY1n1U2Y2n YnnUn0,i,j1,2,,n,ij)

(2-5)表明，自导纳Yii在数值上等于仅在节点i施加单位电压而其余节点电压均为零(即其余节点全部接地)时，经节点i注入网络的电流。其显然等于与节点i直接相连的所有支路的导纳之和。同时可见(U0,i,j1,2,n,ji)。表明，互导纳在数值上等于YijIi/Uji仅在节点j施加单位电压而其余节点电压均为零时，经节点i注入网络的电流，其显然等于(yij)即Yij=yij。yij为支路的导纳，负号表示该电流流出网络。如节点i j之间无支路直接相连，则该电流为0，从而Yij=0。

注意字母几种不写法的不同意义：粗体黑字表示导纳矩阵，大写字母Yij代矩阵YB中的第i行第j列元素，即节点i和节点j之间的互导

-8-电力系统潮流分析计算的MATLAB仿真

纳。小写字母i，j支路的导纳等于支路阻抗的倒数数，yij1/Zij。

根据定义直接求取节点导纳矩阵时，注意以下几点：

(1)节点导纳矩阵是方阵，其阶数就等于网络中除去参考节点外的节点数。参考节点一般取大地，编号为零。

(2)节点导纳矩阵是稀疏矩阵，其各行非零非对角元素就等于与该行相对应节点所连接的不接地支路数。

(3)节点导纳矩阵的对角元素就等于各该节点所连接导纳的总和。因此，与没有接地支路的节点对应的行或列中，对角元素为非对角元素之和的负值。

(4)节点导纳矩阵的非对角元素等于连接节点i，j支路导纳的负值。因此，一般情况下，节点导纳矩阵的对角元素往往大于非对角元素的负值。

(5)节点导纳矩阵一般是对称矩阵，这是网络的互易特性所决定的。从而，一般只要求求取这个矩阵的上三角或下三角部分。

2.1.3 非标准变比变压器等值电路

变压器型等值电路更便于计算机反复计算,更适宜于复杂网络的潮流计算.双绕组变压器可用阻抗与一个理想变压器串联的电路表示.理想变压器只是一个参数,那就是变比KU1/U2。现在变压器阻抗按实际变比归算到低压侧为例,推导出变压器型等值电路。

图2.1双绕组变压器原理图

图2.2变压器阻抗归算到低压侧等值模型

-9-电力系统潮流分析计算的MATLAB仿真

流入和流出理想变压器的功率相等

IUI/K U111

2I1I2/K

(2-6)式(2-6)中, KU1/U2是理想变压器的变比,U1和 U2分别为变压器高,低绕组的实际电压.从图2-2直接可得：

KUIZU122T

(2-7)

从而可得: I1U1KZ2TU2KZTYTU1K2YTU2K

UUYTU121

(2-8)I2=-=-YTU2KZTZTK式(2-8)中YT1/ZT,又因节点电流方程应具有如下形式:

=YU+YU I1111122+YU

(2-9)-I2Y21U1222将式(2-8)与(2-9)比较，得：Y11=YT/K2，Y12=-YT/K；Y21=-YT/K，Y22=YT。

因此可得各支路导纳为：

Y21=-Y21=YT/K1-KY10=Y11-Y12=YT

(2-10)2KK-1Y20=Y22-Y21=YTKY12=-Y12=YT/K由此可得用导纳表示的变压器型等值电路：

-10-电力系统潮流分析计算的MATLAB仿真

图2.3 变压器型等值电路

2.2潮流计算的基本方程

在潮流问题中，任何复杂的电力系统都可以归纳为以下元件(参数)组成。

(1)发电机(注入电流或功率)(2)负荷(注入负的电流或功率)(3)输电线支路(电阻，电抗)(4)变压器支路(电阻，电抗，变比)(5)母线上的对地支路(阻抗和导纳)(6)线路上的对地支路(一般为线路充电点容导纳)集中了以上各类型的元件的简单网络如图2-4。

图2.4计算用的电网结构图

-11-电力系统潮流分析计算的MATLAB仿真

图2.5 潮流计算等值网络

采用导纳矩阵时，节点注入电流和节点电压构成以下线性方程组：

I=YU

(2-11)

 I1U1I2U2 其中，I=

； U=。

MMUn In 可展开如下形式： Ii n Yj 1ij Uj(i=1,2,n)

(2-12)由于实际电网中测量的节点注入量一般不是电流而是功率，因此必须将式中的注入电流用节点注入功率来表示。

节点功率与节点电流之间的关系为：

I

(2-13)

Si=PijQiUii式中PiPGiPLDi，QiQGiQLDi

PijQiIS/U因此用导纳矩阵时，PQ节点可以表示为i iiUi把这个关系代入式中，得

-12-电力系统潮流分析计算的MATLAB仿真

PijQiUinYj1ijUj(i1,2,n)

(2-14)式(3-4)就是电力系统潮流计算的数学模型-----潮流方程。它具有如下特点：

(1)它是一组代数方程，因而表征的是电力系统的稳定运行特性。

(2)它是一组非线性方程，因而只能用迭代方法求其数值解。(3)由于方程中的电压和导纳既可以表为直角坐标，又可表为极坐标，因而潮流方程有多种表达形式---极坐标形式，直角坐标形式和混合坐标形式。

U，Y|y|，得到潮流方程的极坐标①：取 Uijijijiii形式：

nPijQiUiiYijUji

(2-15)

j1ejf，YGjB，得到潮流方程的直角②：取 Uijijijiii坐标形式：

Piei(GijejBijfj)fi(GijfjBijej)j1j1

(2-16)nnQifi(GijejBijfj)ei(GijfjBijej)j1j1nnU YGjB，得到潮流方程的混合坐标形③：取 Uijijijiii式：

nPiUiUj(GijcosijBijsinij)j1

(2-17)nQiUUj(GijsinijBijcosij)ij1不同坐标形式的潮流方程适用于不同的迭代解法。例如：利用牛

-13-电力系统潮流分析计算的MATLAB仿真

顿---拉夫逊迭代法求解，以直角坐标和混合坐标形式的潮流方程为方便；而P-Q解耦法是在混合坐标形式的基础上发展而成，故当然采用混合坐标形式。

(4)它是一组n个复数方程，因而实数方程数为2n个但方程中共含4n个变量：P，Q，U和，i=1，2，，n，故必须先指定2n个变量才能求解。

2.3 电力系统节点分类

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

PQ节点

对这一类点，事先给定的是节点功率(P，Q)，待求的未知量是节点电压向量(U，)，所以叫PQ节点。通常变电所母线都是PQ节点，当某些发电机的输出功率P。Q给定时，也作为PQ节点。PQ节点上的发电机称之为PQ机(或PQ给定型发电机)。在潮流计算中，系统大部分节点属于PQ节点。

PU节点

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

在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在计算中取该节点电压向量的方向作为参考轴，相当于给定该

-14-电力系统潮流分析计算的MATLAB仿真

点电压向量的角度为零。也就是说，对平衡节点给定的运行参数是U和，因此有城为U节点，而待求量是该节点的P。Q，整个系统的功率平衡由这一节点承担。

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

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

2.4潮流计算的约束条件

电力系统运行必须满足一定技术和经济上的要求。这些要求够成了潮流问题中某些变量的约束条件，常用的约束条件如下：

① 节点电压应满足

UiminUiUim(axi1,2,n)

(2-18)从保证电能质量和供电安全的要求来看，电力系统的所有电气设备都必须运行在额定电压附近。PU节点电压幅值必须按上述条件给定。因此，这一约束条件对PQ节点而言。

② 节点的有功功率和无功功率应满足

PGQGmiinmiinPGiGiPQ

(2-19)QGmiaxGmiaxPQ节点的有功功率和无功功率，以及PU节点的有功功率，在给定是就必须满足上述条件，因此，对平衡节点的P和Q以及PU节点的Q应按上述条件进行检验。

③ 节点之间电压的相位差应满足：

|ij||ij||ij|max

(2-20)为了保证系统运行的稳定性，要求某些输电线路两端的电压相位不超过一定的数值。这一约束的主要意义就在于此。

因此，潮流计算可以归结为求解一组非线性方程组，并使其解答满足一定的约束条件。常用的方法是迭代法和牛顿法，在计算过程

-15-电力系统潮流分析计算的MATLAB仿真

中，或得出结果之后用约束条件进行检验。如果不能满足要求，则应修改某些变量的给定值，甚至修改系统的运行方式，重新进行计算。

-16-电力系统潮流分析计算的MATLAB仿真

3牛顿－拉夫逊法概述

3.1牛顿-拉夫逊法基本原理

牛顿--拉夫逊法(简称牛顿法)在数学上是求解非线性代数方程式的有效方法。其要点是把非线性方程式的求解过程变成反复地对相应的线性方程式进行求解的过程。即通常所称的逐次线性化过程。

对于非线性代数方程组：

f(x)0 即fi(x1,x2,,xn)0

(i1,2,n,(3-1)在待求量x的某一个初始估计值x(0)附近，将上式展开成泰勒级数并略去二阶及以上的高阶项，得到如下的经线性化的方程组：

'(0)(0)x)x 0

f(x(0))f（(3-2)上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量

x(0)[f(x'(0))]1f(x(0))

(3-3)

(1)(1)将x(0)和x(0)相加，得到变量的第一次改进值x。接着就从x(0)出发，重复上述计算过程。因此从一定的初值x出发，应用牛顿法求解的迭代格式为：

-17-电力系统潮流分析计算的MATLAB仿真

f(x'(k))x(k)f(x(k))

(3-4)

x(k1)x(k)x(k)

(3-5)上两式中：f(x)是函数f(x)'对于变量x的一阶偏导数矩阵，即雅可比矩阵J；k为迭代次数。

有上式可见，牛顿法的核心便是反复形式并求解修正方程式。牛顿法当初始估计值x和方程的精确解足够接近时，收敛速度非常快，具有平方收敛特性。

牛顿潮流算法突出的优点是收敛速度快，若选择到一个较好的初值，算法将具有平方收敛特性，一般迭代4~5次便可以收敛到一个非常精确的解。而且其迭代次数与所计算网络的规模基本无关。牛顿法也具有良好的收敛可靠性，对于对以节点导纳矩阵为基础的高斯法呈病态的系统，牛顿法也能可靠收敛。牛顿法所需的内存量及每次迭代所需时间均较高斯法多。

牛顿法的可靠收敛取决于有一个良好的启动初值。如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的节点上。对于正常运行的系统，各节点电压一般均在额定值附近，偏移不会太大，并且各节点间的相位角差也不大，所以对各节点可以采用统一的电压初值(也称为平直电压)，如假定：

Ui(0)(0)1 i(0)0或ei(0)1 fi(0)0(iq1,2,,n;is)

(3-6)

这样一般能得到满意的结果。但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时，仍用上述初始电压就有可能出现问题。解决这个问题的办法可以用高斯法迭代1~2次，以此迭代结果作为牛顿法的初值。也可以先用直流法潮流求解一次以求得一个较好的角度初值，然后转入牛顿法迭代。

3.2牛顿--拉夫逊法潮流求解过程

以下讨论的是用直角坐标形式的牛顿—拉夫逊法潮流的求解过

-18-电力系统潮流分析计算的MATLAB仿真

程。当采用直角坐标时，潮流问题的待求量为各节点电压的实部和虚部两个分量e1,f,e12,f2...en,fn由于平衡节点的电压向量是给定的，因此待求两共2(n-1)需要2(n-1)个方程式。事实上，除了平衡节点的功率方程式在迭代过程中没有约束作用以外，其余每个节点都可以列出两个方程式。对PQ节点来说，P而可以写出

PiQiis和Q是给定的，因

ispQisei(GjiijeijjBfijj)f(GfjijjijjBeijjisf(GeijijBfijj)e(GfijjijBeij)0 )0j

(3-7)对PV节点来说，给定量是PiU2iPis和Q，因此可以列出

isjPisei(Gji2ijejBfij)f(GfiijjijBeijjU(eiis2f2i)0)0

(3-8)求解过程大致可以分为以下步骤：(1)形成节点导纳矩阵；(2)将各节点电压设初值U

(3)将节点初值代入相关求式，求出修正方程式的常数项向量；(4)将节点电压初值代入求式，求出雅可比矩阵元素；(5)求解修正方程，求修正向量；(6)求取节点电压的新值；

(7)检查是否收敛，如不收敛，则以各节点电压的新值作为初值自第3步重新开始进行狭义次迭代，否则转入下一步；

(8)计算支路功率分布，PV节点无功功率和平衡节点注入功率。

以直角坐标系形式表示：

1：迭代推算式

采用直角坐标时,节点电压相量及复数导纳可表示为:

ejfViiiYijGijjBij

(3-9)

-19-电力系统潮流分析计算的MATLAB仿真

将以上二关系式代入上式中,展开并分开实部和虚部;假定系统中的第1,2,,m号为P—Q节点,第m+1,m+2,,n-1为P—V节点,根据节点性质的不同,得到如下迭代推算式:(1)对于PQ节点

PiPiei(GijejBijfj)fi(GijfjBijej)j1j1

(3-10)nnQiQifi(GijejBijfj)ei(GijfjBijej)j1j1nni1,2,,m(2)对于PV节点

PiPiei(GijejBijfj)fi(GijfjBijej)j1j1

(3-11)2222UIUi(eifi)nnim1,m2,,n1

(3)对于平衡节点

平衡节点只设一个,电压为已知,不参见迭代,其电压为:

Unenjfn

(3-12)2：修正方程

两组迭代式中包括2(n-1)个方程.选定电压初值及变量修正量符号之后代入,并将其按泰勒级数展开,略去ei,fi二次方程及以后各项,得到修正方程如下:

WJU

(3-13)

-20-电力系统潮流分析计算的MATLAB仿真

P1Q1PmQmWPm1U2m1Pn12Un1e1f1emfmUem1fm1en1fn1其中，。

P1e1Q1e1Pme1Qme1JPm1e12Um1e1Pn1e12Un1e1P1f1Q1f1Pmf1Qmf1Pm1f1U2m1P1emQ1emPmemQmemPm1emU2m1P1fmQ1fmPmfmQmfmPm1fmU2m1P1em1Q1em1Pmem1Qmem1Pm1em1U2m1P1fm1Q1fm1Pmfm1Qmfm1Pm1fm1U2m1P1en1Q1en1Pmen1Qmen1Pm1en1U2m1f1Pn1f1U2n1emPn1emU2n1fmPn1fmU2n1em1Pn1em1U2n1fm1Pn1fm1U2n1en1Pn1en1U2n1f1emfmem1fm1en1fn1Q1fn1Pmfn1Qmfn1Pm1fn12Um1fn1Pn1fn12Un1fn1P1

(3-14)3：雅可比矩阵各元素的算式

式(3-14)中, 雅可比矩阵中的各元素可通过对式(3-10)和(3-11)进行偏导而求得.当ji时, 雅可比矩阵中非对角元素为

-21-电力系统潮流分析计算的MATLAB仿真

(GijeiBijfi)ejfjPiQiBijeiGijfi

(3-15)fjej22UU0ejfjPiQi当ji时,雅可比矩阵中对角元素为:

(GijejBijfj)GiieiBiifieij1nPi(GijfjBijej)GiifiBiieifjj1nQi(GijfjBijej)GiifiBiieieij1

nQi(GijejBijfj)GiieiBiififjj12Ui2eiej2Ui2fifiPin

(3-16)由式(3-15)和(3-16)看出,雅可比矩阵的特点：

(1)矩阵中各元素是节点电压的函数,在迭代过程中,这些元素随着节点电压的变化而变化；

(2)导纳矩阵中的某些非对角元素为零时,雅可比矩阵中对应的元素也是为零.若Yij0,则必有Jij0；

(3)雅可比矩阵不是对称矩阵;(iq1,2,,n;is)

雅可比矩阵各元素的表示如下：

HijPiej(GijeiBijfi)(ji)

(3-17)(GijejBijfj)GiieiBiifi(ji)ji-22-电力系统潮流分析计算的MATLAB仿真

BijeiGijfi)(ji)Pi Nij

(3-18)(GijfjBijej)BiieiGiifi(ji)fjjiBijeiGijfi)(ji)Qi

(3-19)Mijej(GijfjBijej)BiieiGiifi(ji)jiGijeiBijfi)LQi(ji)ijfj(GB

ijejijfj)GiieiBiifi(ji)ji2Ri0(ji)ijUej2e

i(ji)U2Siijf0(ji)j2f

i(ji)3.3 牛顿—拉夫逊法的程序框图

图3.1牛顿—拉夫逊法的程序框图

(3-20)

(3-21)(3-22)

电力系统潮流分析计算的MATLAB仿真

图3.1牛顿—拉夫逊法的程序框图

-24-电力系统潮流分析计算的MATLAB仿真

4潮流计算程序的实现

4.1MATLAB软件简介

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

MATLAB是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整，且具有优良的移植性，它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。

另外，MATLAB提供了一种特殊的工具：工具箱(TOOLBOXES).这些工具箱主要包括：信号处理(SIGNAL PROCESSING)、控制系统(CONTROL SYSTEMS)、神经网络(NEURAL NETWORKS)、模糊逻辑(FUZZY LOGIC)、小波(WAVELETS)和模拟(SIMULATION)等等。不同领域、不同层次的用户通过相应工具的学习和应用，可以方便地进行计算、分析及设计工作。

MATLAB设计中，原始数据的填写格式是很关键的一个环节，它与程序使用的方便性和灵活性有着直接的关系。原始数据输入格式的设计，主要应从使用的角度出发，原则是简单明了，便于修改。

-25-电力系统潮流分析计算的MATLAB仿真

4.2矩阵的运算

矩阵是MATLAB数据存储的基本单元，而矩阵的运算是MATLAB语言的核心，在MATLAB语言系统中几乎一切运算均是以对矩阵的操作为基础的。矩阵的基本数学运算包括矩阵的四则运算、与常数的运算、逆运算、行列式运算、秩运算、特征值运算等基本函数运算，这里进行简单介绍。

(1)四则运算

矩阵的加、减、乘运算符分别为“+，-，*”，用法与数字运算几乎相同，但计算时要满足其数学要求在MATLAB中矩阵的除法有两种形式：左除“”和右除“/”。在传统的MATLAB算法中，右除是先计算矩阵的逆再相乘，而左除则不需要计算逆矩阵直接进行除运算。通常右除要快一点，但左除可避免被除矩阵的奇异性所带来的麻烦。在MATLAB6中两者的区别不太大。

(2)与常数的运算

常数与矩阵的运算即是同该矩阵的每一元素进行运算。但需注意进行数除时，常数通常只能做除数。

(3)基本函数运算

矩阵的函数运算是矩阵运算中最实用的部分，常用的主要有以下几个：

det(a)

求矩阵a的行列式

eig(a)

求矩阵a的特征值 inv(a)或a ^(-1)

求矩阵a的逆矩阵 rank(a)

求矩阵a的秩

trace(a)

求矩阵a的迹(对角线元素之和)我们在进行工程计算时常常遇到矩阵对应元素之间的运算。这种运算不同于前面讲的数学运算，为有所区别，我们称之为数组运算。

(4)基本数学运算

数组的加、减与矩阵的加、减运算完全相同。而乘除法运算有相当大的区别，数组的乘除法是指两同维数组对应元素之间的乘除法，它们的运算符为“.*”和“./”或“.”。前面讲过常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定，数组与常数之

-26-电力系统潮流分析计算的MATLAB仿真

间的除法运算没有任何限制。

另外，矩阵的数组运算中还有幂运算(运算符为.^)、指数运算(exp)、对数运算(log)、和开方运算(sqrt)等。有了“对应元素”的规定，数组的运算实质上就是针对数组内部的每个元素进行的。矩阵的幂运算与数组的幂运算有很大的区别。

(5)逻辑关系运算

逻辑运算是MATLAB中数组运算所特有的一种运算形式，也是几乎所有的高级语言普遍适用的一种运算。

4.3牛顿—拉夫逊法潮流计算实例

潮流计算例题：根据给定的参数或工程具体要求(如图4-1)，图中节点1为平衡节点，节点2、3、4、5为PQ节点。

图4.1系统接线图

1.画出程序流程图。

2.运用MATLAB软件进行潮流计算。解：① 画程序流程图。参照图3.1

② 根据程序流程图，然后写出程序(详细见附录)，最后运行MATLAB软件进行计算，计算结果如下。

-27-电力系统潮流分析计算的MATLAB仿真

Y值：

图4.2计算结果截图

迭代过程：

图4.3计算结果截图

-28-电力系统潮流分析计算的MATLAB仿真

图4.4计算结果截图

图4.5计算结果截图

电力系统潮流分析计算的MATLAB仿真

图4.6计算结果截图

图4.7计算结果截图

-30-电力系统潮流分析计算的MATLAB仿真

图4.8计算结果截图

图4.9计算结果截图

-31-电力系统潮流分析计算的MATLAB仿真

图4.10计算结果截图

电压值：

图4.12计算结果截图

-32-电力系统潮流分析计算的MATLAB仿真

平衡节点注入功率及电流：

4.13计算结果截图

图电力系统潮流分析计算的MATLAB仿真

结论

在电力系统的潮流计算算法中牛顿一拉夫逊法是得到电力系统研究人员认可的算法之一，在本文中我们采用牛顿一拉夫逊法，主要是同时考虑到计算的准确和程序的运行速度。没有采用经常用的高斯迭代法，而是采用了传统的逆阵方法，是考虑到用MATLAB实现高斯迭代将会通过很多的循环迭代才能实现，而逆阵可以直接通过命令来求解，这必然可以大大节省时间。对于不能求逆的矩阵我们通过在电力系统中至少有一条支路上有接地支路来实现其求逆。对于P-Q分解法不能使用于有些R/X比较大的电力系统的缺点，可以通过在电力系统中并联补偿法或虚构节点来得以解决。通过实例计算分析，取得了比较满意的效果。基于MATLAB的电力系统潮流计算使计算机在计算、分析、研究复杂的电力系统潮流分布问题上又前进了一步。不管采用什么算法，所有的潮流计算都是基于矩阵的迭代运算。而MATLAB语言正是以处理矩阵见长，实践证明，MATLAB语言在电力系统潮流计算仿真研究中的应用是可行的，而且由于其强大的矩阵处理功能，完全可以应用于电力系统的其它分析计算中；用MATLAB语言编程效率高，程序调试十分方便，可大大缩减软件开发周期，如果像控制界一样开发出电力系统自己的专用工具箱，将系统分析用的一些基本计算以函数的形式直接调用，那么更高层次的系统软件也可以很容易地实现。

-34-电力系统潮流分析计算的MATLAB仿真

致谢

大学四年的学习生活即将结束了，本次毕业设计也已经接近尾声，作为一个本科生的毕业设计，由于经验的匮乏，难免有许多考虑不周全的地方，但在各位老师的督促指导下，以及同学们的关心和帮助下，我的各方面都有所提高，整体素质也有了很大的进步。这一切都是和我系的老师对我的指引和教育是分不开的。

首先我要衷心感谢我的指导老师赵熹老师，这篇论文是在赵老师的悉心指导下，结合相关的书籍、资料完成的，更是倾注了她大量的心血。赵老师在毕业设计过程中给予了我细心指导和无私帮助，从初稿到定稿，赵老师不厌其烦，一审再审，大到篇章布局的偏颇，小到语句格式的瑕疵，都一一予以指出，帮助我解决了论文中出现的许多难题。她的教诲，令我终身难忘。在论文完成之际，特向赵老师致以衷心的感谢！

同时感谢学校各位领导和老师，是他们在我学业和生活中给予了极大的关心和帮助。是他们的默默工作，为我们所有的学生创建了一个安定和有序的学习环境。

感谢我的室友们，从遥远的家来到这个陌生的城市里，是你们和我共同维系着彼此之间兄弟般的感情，维系着寝室那份家的融洽。四年的时光，仿佛就在昨天。

谨向我的父母和家人表示诚挚的谢意。他们是我生命中永远的依靠和支持，他们无微不至的关怀，是我前进的动力；他们的殷殷希望，激发我不断前行。没有他们就没有我，我的点滴成就都来自他们。

感谢老师们百忙中审阅我的论文！

在论文即将完成之际，我的心情非常激动，从开始进入课题到论文的顺利完成，有多少可敬的师长、同学、朋友给了我无言的帮助，在这里请接受我诚挚的谢意！

-35-电力系统潮流分析计算的MATLAB仿真

参考文献

[1] 于永源,杨绮雯.电力系统分析(第二版).中国电力出版社, 2004.[2] 刘从爱.电力工程基础.山东科学技术出版社, 1997.[3] 《电气工程师手册》第二版编辑委员会.电气工程师手册.机械工业出版社, 2000.[4] 何仰赞,温增银.电力系统分析.华中科技大学出版社, 2002.[5] 吴际舜,侯志捡.电力系统潮流计算的计算机方法.上海交通大学出版社, 2000.[6] 陈珩.电力系统稳态分析(第三版).中国电力出版社, 2007.[7] 丽华.电力工程基础.机械工业出版社, 2006.[8] 李维波.MATLAB在电气工程中的应用.中国电力出版社, 2007.[9] 陈恳.直角坐标牛顿一拉夫逊法潮流计算新解法[J].电力系统及其自动化学报, 1999.[10] 张铮.MATLAB程序设计与实例应用.中国铁道出版社, 2003.[11] 张志涌,徐彦琴.MATLAB教程-基于6.x版本.北京航空航天大学出版社, 2001.-36-电力系统潮流分析计算的MATLAB仿真

附件 MATLAB计算实例程序

clear;clc %重新编号，把原题中的节点1,2,3,4,5重新依次编号为5,1,2,3,4，其中1-4号为PQ节点，5号为平衡节点 y=0;%输入原始数据，求节点导纳矩阵

y(1,2)=1/(0.06+0.18i);y(1,3)=1/(0.06+0.18i);y(1,4)=1/(0.04+0.12i);y(1,5)=1/(0.02+0.06i);y(2,3)=1/(0.01+0.03i);y(2,5)=1/(0.08+0.24i);y(3,4)=1/(0.08+0.24i);y(4,5)=0;for i=1:5 for j=i:5 y(j,i)=y(i,j);end end Y=0;%求互导纳 for i=1:5 for j=1:5 if i~=j Y(i,j)=-y(i,j);end end end %求自导纳 for i=1:5 Y(i,i)=sum(y(i,:));end Y %Y 为导纳矩阵 G=real(Y);B=imag(Y);

-37-电力系统潮流分析计算的MATLAB仿真

%原始节点功率

S(1)=0.2+0.2i;S(2)=-0.45-0.15i;S(3)=-0.4-0.05i;S(4)=-0.6-0.1i;S(5)=0;P=real(S);Q=imag(S);%赋初值

U=ones(1,5);U(5)=1.06;e=zeros(1,5);ox=ones(8,1);fx=ones(8,1);count=0 %计算迭代次数 while max(fx)>1e-5 for i=1:4

for j=1:4

H(i,j)=0;N(i,j)=0;M(i,j)=0;L(i,j)=0;oP(i)=0;oQ(i)=0;

end end for i=1:4 for j=1:5 oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)));end oP(i)=oP(i)+P(i);oQ(i)=oQ(i)+Q(i);end fx=[oP,oQ]';%求雅克比矩阵

%当i~=j时候求H,N,M,L 如下： for i=1:4 for j=1:4 if i~=j H(i,j)=-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)));N(i,j)=-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));L(i,j)=H(i,j);M(i,j)=-N(i,j);end end

-38-电力系统潮流分析计算的MATLAB仿真

end H,N,M,L %当i=j 时H,N,M,L如下： for i=1:4 for j=1:5 if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i, j)*cos(e(i)-e(j)));

N(i,i)=N(i,i)-U(i)*U(j)*(G(i, j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)));end end N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);end J=[H,N;M,L] %J 为雅克比矩阵 ox=-((inv(J))*fx);for i=1:4 oe(i)=ox(i);oU(i)=ox(i+4)*U(i);end for i=1:4 e(i)=e(i)+oe(i);U(i)=U(i)+oU(i);end count=count+1;end ox,U,e,count %求节点注入的净功率 i=5;for j=1:5 P(i)=U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)))+P(i);Q(i)=U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)))+Q(i);end S(5)=P(5)+Q(5)*sqrt(-1);S %求节点注入电流 I=Y*U'

第二篇：用matlab电力系统潮流计算

题目：潮流计算与matlab

教学单位电气信息学院

姓名学号

年级

专业电气工程及其自动化

指导教师

职称副教授

摘要

电力系统稳态分析包括潮流计算和静态安全分析。本文主要运用的事潮流计算，潮流计算是电力网络设计与运行中最基本的运算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中的各元件的电力损耗，进而求得电能损耗。本位就是运用潮流计算具体分析，并有MATLAB仿真。

关键词：电力系统潮流计算 MATLAB

Abstract Electric power system steady flow calculation and analysis of the static safety analysis.This paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power.The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation.Key words: Power system;Flow calculation;MATLAB simulation

目录任务提出与方案论证....................................................................................................................................2 2 总体设计........................................................................................................................................................3 2.1潮流计算等值电路.............................................................................................................................3 2.2建立电力系统模型.............................................................................................................................3 2.3模型的调试与运行.............................................................................................................................3 3 详细设计........................................................................................................................................................4 3.1 计算前提............................................................................................................................................4 3.2手工计算.............................................................................................................................................7 4设计图及源程序...........................................................................................................................................11 4.1MATLAB仿真.......................................................................................................................................11 4.2潮流计算源程序...............................................................................................................................11 5 总结.............................................................................................................................................................31 参考文献..........................................................................................................................................................32 任务提出与方案论证

潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况下确定系统稳态运行状态的一种基本方法,是电力系统规划和运营中不可缺少的一个重要组成部分。可以说,它是电力系统分析中最基本、最重要的计算,是系统安全、经济分析和实时控制与调度的基础。常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。是电力系统研究人员长期研究的一个课题。它既是对电力系统规划设计和运行方式的合理性、可靠性及经济性进行定量分析的依据，又是电力系统静态和暂态稳定计算的基础。

潮流计算经历了一个由手工到应用数字电子计算机的发展过程，现在的潮流算法都以计算机的应用为前提用计算机进行潮流计算主要步骤在于编制计算机程序，这是一项非常复杂的工作。对系统进行潮流分析，本文利用 MATLAB中的SimpowerSystems工具箱设计电力系统，在simulink 环境下，不仅可以仿真系统的动态过程，还可以对系统进行稳态潮流分析。

总体设计

SimpowerSystems使用Simulink环境，可以将该系统中的发电机、变压器，线路等模型联结起来，形成电力系统仿真模拟图。在加人测量模块，并对各元件的参数进行设置后，用measurement和sink中的仪器可以观察各元件的电压、电流、功率的大小。

2.1潮流计算等值电路

10MWYN,d11463MWVA15MWGGGG120MW10kVp015.7kWps73kWI%0.50Vs%10.5YN,d1116MWVA463MW“xd0.134x20.161x0.060cosN0.8510kVp011kWps50kWI%0.550Vs%10.5YN,d11210MWVA35kV32km25MW110kV80km25MW110kV70km110kVYN,d11220MWVA20MWGGG415MW”xd0.136x20.16x0.0730cosN0.8p018.6kWps89kWI%0.530MW0Vs%10.510kVp015.7kWps73kWI%0.5GG0V%10.5sYN,d11216MWVA63MWVAp044kWps121kW10kVI%0.350Vs%10.5GG35MWYN,Y,d11210MVA10kVGGG312MW150MW“xd0.128x20.154x0.0540cosN0.85225MW”xd0.128x20.157x0.05910cosN0.8x0.136x20.161x0.0750cosN0.8“d

2.2建立电力系统模型

在Simulink中按照电力系统原型选择元件进行建模。所建立的模型和建立的方法在详细设计中详述。

在电力系统模型的建立工程中主要涉及到的是：元器件的选择及其参数的设置；发电机选型；变压器选择；线路的选择；负荷模型的选择；母线选择。

2.3模型的调试与运行

建立系统模型，并设置好参数以后，就可以在Simulink环境下进行仿真运行。运行的具体结果和分析也在详细设计中详述。

30km35kV0km31p044kWps121kWI%0.350Vs%10.5p013.2kWps63kWI0%0.55V%10.5s(12)80MWVs(23)%6.55MWVs(13)%17.53 详细设计

3.1 计算前提

首先是发电机的参数计算，先对5个发电厂简化为5台发电机来计算。发电机G1：

P141560MWQ160tan(arccos0.8)45MVar发电机G2：

P2463252MWQ2252tan(arccos0.85)156MVar发电机G3：

P331236MWQ336tan(arccos0.8)27MVar发电机G4：

P415050MWQ450tan(arccos0.85)31MVar发电机G5：

P522550MWQ550tan(arccos0.8)37.5MVar

其次是变电站的参数计算，我们还是对7个变电站简化为7台变压器来计算。变压器T1：

2psVN7311023RT1101033.450232SN(1610)2Vs%VN10.51102XT1101079.406SN16103S01p0j变压器T2：（双并联）

I0%SN(0.0157j0.0800)MVA 100RT2XT221psVN18911023101031.3462322SN2(2010)21Vs%VN110.51102101031.7625 32SN22010S022(p0j变压器T3：（四并联）

I0%SN)(0.0372j0.2000)MVA100 1psVN2112111023RT3101030.0922324SN4(6310)XT31Vs%VN2110.5110210105.042 4SN463103I0%SN)(0.1760j0.8820)MVA100S034(p0j变压器T4：（双并联）

1RT11.725021 XT4XT139.70302S042S01(0.0314j0.1600)MVART4变压器T5：

RT54RT30.3680XT54XT320.168S051S03(0.0440j0.2205)MVA41633523100.386 322(1010)

变压器T6：（两个三绕组变压器并联）

RT61RT62RT631[Vs(12)%Vs(13)%Vs(23)%]10.7521Vs2%[Vs(12)%Vs(23)%Vs(13)%]0.25

21Vs3%[Vs(13)%Vs(23)%Vs(12)%]6.752Vs1%21Vs1%VNXT61106.5842SNXT62XT6321Vs2%VN100.1532SN21Vs3%VN104.134 2SNS062(P06j变压器T7：（双并联）

I0%10)(0.0264j0.1100)MVA 100 RT7XT721psVN1503523101030.3062322SN2(1010)21Vs%VN110.535210106.4312SN210103

S072(p0jI0%SN)(0.0220j0.1100)MVA100再次是传输线参数计算，5条传输线的具体计算如下。

根据教材查得r00.21/km x00.4/km b02.810S/km 6线路L1：

线路L2：

线路L3：（双回路）

线路L4：

线路L5：（双回路）RL1r0l10.21408.4XL1x0l10.44016B6L1b0l12.810401.12104S Q1L12BL1V2N0.6776MVarRL2r0l20.2113027.3XL2x0l20.413052B6L2b0l22.8101303.64104S Q1L22BL2V2N2.2022MVarR12rl1L30320.21707.35X11L32x0l320.47014 BL32b40l322.8106703.9210SQ1L32B2L3VN2.3716MVarRL4r0l40.216012.6XL4x0l40.46024BL4b0l42.8106601.68104S Q12L42BL4VN1.0164MVar

11RL5r0l50.21202.12211XL5x0l50.4204 22BL52b0l522.8106201.12104S1QL5BL3VN20.0686MVar23.2手工计算

FLR1：

P2102ST12(RT1jXT1)(3.450j74.406)(0.0285j0.6562)MVA2VN110Sa10MWST1S01jQL1(10.0442j0.1142)MVAP2Q210.044220.11422SL1(RL1jXL1)(8.4j16)(0.070j0.1334)MVAVN21102ST2P2Q2402452(RT2jXT2)(1.346j31.7625)(0.4032j9.5156)MVAVN21102FLR2SbSG120ST260j45200.4032j9.5156(39.5968j35.4844)MVAScSbSa25jQL1SL1(4.4826j35.9144)MVA：

ST3P2Q225221562(RT3jXT3)(0.092j5.042)(0.6679j36.6024)MVA22VN110PQ4.493134.1048(RjX)(27.3j52)(2.67j5.0854)MVAL2L222VN1102222Sc'(4.4931j34.1048)MVASL2FLRSdSG2Sc'120ST3S03jQL2SL2(132.9792j149.229)MVA3：

ST4P2Q262272(RT4jXT4)(1.725j39.703)(0.1091j2.5101)MVAVN21102P2Q2133.59552149.99562(RL3jXL3)(7.35j14)(24.51j46.682)MVA22VN110'Sd(133.5955j149.9956)MVASL3'SeSG3Sd3025ST4S04jQL3SL3(89.945j130.0151)MVAFLR4： ST5P2Q2502312(RT5jXT5)(0.368j20.168)(0.1052j5.7687)MVA22VN110P2Q292.74872133.99372(RL4jXL4)(12.6j24)(27.654j52.674)MVA22VN110Se'(92.7481j133.9937)MVASL4SfSG4Se'80ST5S05jQL4SL4(34.9449j107.3469)MVAFLR5: 152ST72(0.306j6.431)(0.0562j1.1812)MVA35Sh15ST7S07jQL5(15.0782j0.3422)MVA15.078220.34222SL5(2.1j4)(0.3899j0.743)MVA352SiShSL5S06jQL55(20.4945j1.1266)MVA 15237.52ST63(0.386j4.34)(0.514j5.7793)MVA35220.650520.54512ST62(0.386j0.153)(0.1345j0.0533)MVA23526.336298.73692ST61(0.386j6.584)(3.2905j56.1256)MVA352SgSfST61SG5ST62ST63Si35(25.5114j194.12)MVA计算每一个FLR的功率分布和电压分布计算如下： FLR1：

VT2PRQX401.3464531.762512.8970kVVN115 Vb115VT2102.1030kVPRQX10.04428.40.144216VL10.8489kVVb102.1030VaVbVL1101.2541kVFLR2：

功率分布：

SL2ZZZ*T3*L2*SdT3(VbVN)ZZL2****VN(0.092j5.042)(132.9792j149.229)1418.6727.392j57.042T3(4.8812j13.8097)MVAST3ZZZ*L2*L2*SdT3(VbVN)ZZL2VN(27.3j52)(132.9792j149.229)1418.6727.392j57.042T3(108.687j122.62)MVA 电压分布：

Sc1SL2SL2(4.8812j13.8097)(2.67j5.0854)(7.5512j8.7243)MVA7.551227.38.7243522.424kV102.1030VdVbVL2102.103(2.424)104.527kVVL2功率分布：

FLR3：

SL3ZZZ*T4*L3*SeT4(VG3Vd)ZZL3****VN(1.725j39.703)(89.945j130.0151)1037.9279.075j53.73T4(59.444j16.846)MVAST4ZZZ*L3*L3*SeT4(VbVN)ZZL3VN(7.35j14)(89.945j130.0151)1037.9279.075j53.73T4(31.811j60.1256)MVA 电压分布：

Se1SL3SL3(59.444j19.846)(24.51j46.682)(83.954j26.836)MVA83.9547.3526.836149.404kV105.5643VeVdVL396.16kVVL3功率分布：

FLR4：

SL4ZZZ*T5*L4*SfT5(VG3Vd)ZZL4****VN=(0.368j20.168)(34.9449j107.3469)1037.92712.968j44.168T5(20.843j19.689)MVAST4ZZZ*L4*L4*SfT5(VG3Vd)ZZL4VN=(12.6j24)(34.9449j107.3469)1037.92712.968j44.168T5(1.398j44.389)MVA 电压分布：

Se1SL3SL3(59.444j16.846)(24.51j46.682)(83.954j63.528)MVA83.95412.663.5282424.464kV105.5643VeVdVL381.10kVVL4FLR5: 这里我们先将f点和发电机G5当做电源，经过ZT61和ZT63构成两端供电网络以g点作为运算负荷进行计算。ST6ST4(0.386j4.134)(20.2656j70.9293)(22.093837)35(3.900j25.1175)MVA0.772j10.718(0.386j6.584)(20.2656j70.9293)(22.093837)35(16.5061j91.7905)MVA0.772j10.718电压分布：

ST631ST63ST63(16.6421j97.5698)MVA16.64210.38697.56984.13410.9186kV37Vg37VT6326.0814VVT6320.26560.38670.9293(0.153)0.1162kV

26.0814ViVgVT6226.1976VT6220.49452.11.126641.815kV26.1976VhViVL524.3826VL5

4设计图及源程序

4.1MATLAB仿真

相关的原始数据输入格式如下：

1、B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写。

2、对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗置于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为含变压器支路的变压器的变比，第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器。

3、B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

4、X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地参数。

4.2潮流计算源程序

%本程序的功能是用牛顿——拉夫逊法进行11节点潮流计算 clear;n=11;%input('请输入节点数：n=');nl=11;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1

0.03512+0.08306i

0.13455i

0.0068+0.18375i

1.02381

0.05620+0.13289i

0.05382i

0.00811+0.24549i

1.02381

0.05620+0.13289i

0.05382i

0.04215+0.09967i

0.04037i

0.0068+0.18375i

1.02381

0.02810+0.06645i

0.10764i

0.05620+0.13289i

0.05382i

0;0.00811+0.24549i

0.03512+0.08306i

0.13455i

0] B2=[0

1.1

0.343+0.21256i

0.204+0.12638i

0.306+0.18962i

0.5

1.1

0.343+0.21256i

2];% B1矩阵：

1、支路首端号；

2、末端号；

3、支路阻抗；

4、支路对地电纳 %

5、支路的变比；

6、支路首端处于K侧为1，1侧为0 % B2矩阵：

1、该节点发电机功率；

2、该节点负荷功率；

3、节点电压初始值 %

4、PV节点电压V的给定值；

5、节点所接的无功补偿设备的容量 %

6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %

3为PV节点；

%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl

%支路数

if B1(i,6)==0

%左节点处于1侧

p=B1(i,1);q=B1(i,2);

else

%左节点处于K侧

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元

Y(q,p)=Y(p,q);

%非对角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

%对角元1侧

end %求导纳矩阵

disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);

%分解出导纳阵的实部和虚部

for i=1:n

%给定各节点初始电压的实部和虚部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i,4);

%PV节点电压给定模值

end for i=1:n

%给定各节点注入功率

S(i)=B2(i,1)-B2(i,2);

%i节点注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量

end %===================== P=real(S);Q=imag(S);

%分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0

% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列

IT2=0;a=a+1;

for i=1:n

if i~=isb

%非平衡节点

C(i)=0;D(i)=0;

for j1=1:n

C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)

D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)

end

P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)

Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值

V2=e(i)^2+f(i)^2;

%电压模平方

%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========

if B2(i,6)~=3

%非PV节点

DP=P(i)-P1;

%节点有功功率差

DQ=Q(i)-Q1;

%节点无功功率差

%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================

for j1=1:n

if j1~=isb&j1~=i

%非平衡节点&非对角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df

X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de

X3=X2;

% X2=dp/df X3=dQ/de

X4=-X1;

% X1=dP/de X4=dQ/df

p=2*i-1;q=2*j1-1;

J(p,q)=X3;J(p,N)=DQ;m=p+1;

% X3=dQ/de J(p,N)=DQ节点无功功率差

J(m,q)=X1;J(m,N)=DP;q=q+1;

% X1=dP/de J(m,N)=DP节点有功功率差

J(p,q)=X4;J(m,q)=X2;

% X4=dQ/df X2=dp/df

elseif j1==i&j1~=isb %非平衡节点&对角元

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df

X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de

X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q

m=p+1;

J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P

J(m,q)=X2;

end

else

%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========

DP=P(i)-P1;

% PV节点有功误差

DV=V(i)^2-V2;

% PV节点电压误差

for j1=1:n

if j1~=isb&j1~=i

%非平衡节点&非对角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

% dP/de

X2=B(i,j1)*e(i)-G(i,j1)*f(i);

% dP/df

X5=0;X6=0;

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV节点有功误差

J(m,q)=X2;

elseif j1==i&j1~=isb %非平衡节点&对角元

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df

X5=-2*e(i);

X6=-2*f(i);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV节点有功误差

J(m,q)=X2;

end

end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================

for k=3:N0

% N0=2*n（从第三行开始，第一、二行是平衡节点）

k1=k+1;N1=N;

% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U

for k2=k1:N1

% 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U

J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化

end

J(k,k)=1;

% 对角元规格化K行K列对角元素赋1

%==================== 回代运算

=======================================

if k~=3

% 不是第三行

k > 3

k4=k-1;

for k3=3:k4

% 用k3行从第三行开始到当前行的前一行k4行消去

for k2=k1:N1 %

k3行后各行上三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）

end

%用当前行K2列元素减去当前行k列元素乘以第k行K14 列元素

J(k3,k)=0;%当前行第k列元素已消为0

end

if k==N0

%若已到最后一行

break;

end

%================== 前代运算

==================================

for k3=k1:N0

% 从k+1行到2*n最后一行

for k2=k1:N1

% 从k+1列到扩展列消去k+1行后各行下三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算

end

%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素

J(k3,k)=0;%当前行第k列元素已消为0

end

else

%是第三行k=3

%====================== 第三行k=3的前代运算 ========================

for k3=k1:N0

%从第四行到2n行（最后一行）

for k2=k1:N1

%从第四列到2n+1列（即扩展列）

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）

end

%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素

J(k3,k)=0;%当前行第3列元素已消为0

end

end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====

for k=3:2:N0-1

L=(k+1)./2;

e(L)=e(L)-J(k,N);

%修改节点电压实部

k1=k+1;

f(L)=f(L)-J(k1,N);

%修改节点电压虚部

end

%------修改节点电压-----------

for k=3:N0

DET=abs(J(k,N));

if DET>=pr

%电压偏差量是否满足要求

IT2=IT2+1;%不满足要求的节点数加1

end

ICT2(a)=IT2;

%不满足要求的节点数

ICT1=ICT1+1;

%迭代次数 end %用高斯消去法解”w=-J*V“

disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n

V(k)=sqrt(e(k)^2+f(k)^2);

%计算各节点电压的模值

sida(k)=atan(f(k)./e(k))*180./pi;

%计算各节点电压的角度

E(k)=e(k)+f(k)*j;

%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);

%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);

%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);

%显示各节点的电压相角 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);

%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);

%显示各节点的注入功率

disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

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

Siz(i)=Si(p,q);

else

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

Siz(i)=Si(p,q);

end

disp(Si(p,q));

SSi(p,q)=Si(p,q);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];

disp(ZF);

disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');

for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

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

Sjy(i)=Sj(q,p);

else

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

Sjy(i)=Sj(q,p);

end

disp(Sj(q,p));

SSj(q,p)=Sj(q,p);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];

disp(ZF);

disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p);

disp(DS(i));

DDS(i)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];

disp(ZF);

disp('----------------------');end

%本程序的功能是用牛顿——拉夫逊法进行10节点潮流计算 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 clear;n=10;%input('请输入节点数：n=');nl=10;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1

0.03512+0.08306i

0.13455i

0.0068+0.18375i

1.02381

0.05620+0.13289i

0.05382i

0.00811+0.24549i

1.02381

0.05620+0.13289i

0.05382i

0.04215+0.09967i

0.04037i

0.0068+0.18375i

1.02381

0.02810+0.06645i

0.10764i

0;0.00811+0.24549i

0.03512+0.08306i

0.13455i

0] B2=[0

1.1

0.343+0.21256i

0.204+0.12638i

0.306+0.18962i

0.5

1.1

0.343+0.21256i

2];% B1矩阵：

1、支路首端号；

2、末端号；

3、支路阻抗；

4、支路对地电纳 %

5、支路的变比；

6、支路首端处于K侧为1，1侧为0 % B2矩阵：

1、该节点发电机功率；

2、该节点负荷功率；

3、节点电压初始值 %

4、PV节点电压V的给定值；

5、节点所接的无功补偿设备的容量 %

6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %

3为PV节点；

%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl

%支路数

if B1(i,6)==0

%左节点处于1侧

p=B1(i,1);q=B1(i,2);

else

%左节点处于K侧

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元

Y(q,p)=Y(p,q);

%非对角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

%对角元1侧

end %求导纳矩阵

disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);

%分解出导纳阵的实部和虚部

for i=1:n

%给定各节点初始电压的实部和虚部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i,4);

%PV节点电压给定模值

end for i=1:n

%给定各节点注入功率

S(i)=B2(i,1)-B2(i,2);

%i节点注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量

end %===================== P=real(S);Q=imag(S);

%分解出各节点注入的有功和无功功率

ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0

% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列

IT2=0;a=a+1;

for i=1:n

if i~=isb

%非平衡节点

C(i)=0;D(i)=0;

for j1=1:n

C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)

D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)

end

P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)

Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值

V2=e(i)^2+f(i)^2;

%电压模平方

%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========

if B2(i,6)~=3

%非PV节点

DP=P(i)-P1;

%节点有功功率差

DQ=Q(i)-Q1;

%节点无功功率差

%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================

for j1=1:n

if j1~=isb&j1~=i

%非平衡节点&非对角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df

X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de

X3=X2;

% X2=dp/df X3=dQ/de

X4=-X1;

% X1=dP/de X4=dQ/df

p=2*i-1;q=2*j1-1;

J(p,q)=X3;J(p,N)=DQ;m=p+1;

% X3=dQ/de J(p,N)=DQ节点无功功率差

J(m,q)=X1;J(m,N)=DP;q=q+1;

% X1=dP/de J(m,N)=DP节点有功功率差

J(p,q)=X4;J(m,q)=X2;

% X4=dQ/df X2=dp/df

elseif j1==i&j1~=isb %非平衡节点&对角元

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df

X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de

X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q

m=p+1;

J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P

J(m,q)=X2;

end

else

%=============== 下面是针对PV节点来求取Jacobi矩阵的元素

===========

DP=P(i)-P1;

% PV节点有功误差

DV=V(i)^2-V2;

% PV节点电压误差

for j1=1:n

if j1~=isb&j1~=i

%非平衡节点&非对角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

% dP/de

X2=B(i,j1)*e(i)-G(i,j1)*f(i);

% dP/df

X5=0;X6=0;

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV节点有功误差

J(m,q)=X2;

elseif j1==i&j1~=isb %非平衡节点&对角元

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df

X5=-2*e(i);

X6=-2*f(i);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV节点有功误差

J(m,q)=X2;

end

end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================

for k=3:N0

% N0=2*n（从第三行开始，第一、二行是平衡节点）

k1=k+1;N1=N;

% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U

for k2=k1:N1

% 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U

J(k,k2)=J(k,k2)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化

end

J(k,k)=1;

% 对角元规格化K行K列对角元素赋1

%==================== 回代运算

=======================================

if k~=3

% 不是第三行

k > 3

k4=k-1;

for k3=3:k4

% 用k3行从第三行开始到当前行的前一行k4行消

去

for k2=k1:N1 %

k3行后各行上三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）

end

%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素

J(k3,k)=0;%当前行第k列元素已消为0

end

if k==N0

%若已到最后一行

break;

end

%================== 前代运算

==================================

for k3=k1:N0

% 从k+1行到2*n最后一行

for k2=k1:N1

% 从k+1列到扩展列消去k+1行后各行下三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算

end

%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素

J(k3,k)=0;%当前行第k列元素已消为0

end

else

%是第三行k=3

%====================== 第三行k=3的前代运算 ========================

for k3=k1:N0

%从第四行到2n行（最后一行）

for k2=k1:N1

%从第四列到2n+1列（即扩展列）

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）

end

%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素

J(k3,k)=0;%当前行第3列元素已消为0

end

end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====

for k=3:2:N0-1

L=(k+1)./2;

e(L)=e(L)-J(k,N);

%修改节点电压实部

k1=k+1;

f(L)=f(L)-J(k1,N);

%修改节点电压虚部

end

%------修改节点电压-----------

for k=3:N0

DET=abs(J(k,N));

if DET>=pr

%电压偏差量是否满足要求

IT2=IT2+1;%不满足要求的节点数加1

end

ICT2(a)=IT2;

%不满足要求的节点数

ICT1=ICT1+1;

%迭代次数 end %用高斯消去法解”w=-J*V“ disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n

V(k)=sqrt(e(k)^2+f(k)^2);

%计算各节点电压的模值

sida(k)=atan(f(k)./e(k))*180./pi;

%计算各节点电压的角度

E(k)=e(k)+f(k)*j;

%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);

%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);

%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);

%显示各节点的电压相角 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);

%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);

%显示各节点的注入功率

disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

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

Siz(i)=Si(p,q);

else

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

Siz(i)=Si(p,q);

end

disp(Si(p,q));

SSi(p,q)=Si(p,q);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];

disp(ZF);

disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

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

Sjy(i)=Sj(q,p);

else

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

Sjy(i)=Sj(q,p);

end

disp(Sj(q,p));

SSj(q,p)=Sj(q,p);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];

disp(ZF);

disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p);

disp(DS(i));

DDS(i)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];

disp(ZF);

disp('----------------------');end

%本程序的功能是用牛顿——拉夫逊法进行12节点潮流计算 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 clear;n=12;%input('请输入节点数：n=');nl=12;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1

0.03512+0.08306i

0.13455i

0.0068+0.18375i

1.02381

0.05620+0.13289i

0.05382i

0.00811+0.24549i

1.02381

0.05620+0.13289i

0.05382i

0.04215+0.09967i

0.04037i

0.0068+0.18375i

1.02381

0.02810+0.06645i

0.10764i

0.05620+0.13289i

0.05382i

0;0.00811+0.24549i

0.03512+0.08306i

0.13455i

0.03512+0.08306i

0.13455i

0] B2=[0

1.1

0.343+0.21256i

0.204+0.12638i

0.306+0.18962i

0.5

1.1

0.343+0.21256i

2];% B1矩阵：

1、支路首端号；

2、末端号；

3、支路阻抗；

4、支路对地电纳 %

5、支路的变比；

6、支路首端处于K侧为1，1侧为0 % B2矩阵：

1、该节点发电机功率；

2、该节点负荷功率；

3、节点电压初始值 %

4、PV节点电压V的给定值；

5、节点所接的无功补偿设备的容量 %

6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %

3为PV节点；

%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl

%支路数

if B1(i,6)==0

%左节点处于1侧

p=B1(i,1);q=B1(i,2);

else

%左节点处于K侧

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元

Y(q,p)=Y(p,q);

%非对角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

%对角元1侧

end %求导纳矩阵

disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);

%分解出导纳阵的实部和虚部

for i=1:n

%给定各节点初始电压的实部和虚部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i,4);

%PV节点电压给定模值

end for i=1:n

%给定各节点注入功率

S(i)=B2(i,1)-B2(i,2);

%i节点注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量

end %===================== P=real(S);Q=imag(S);

%分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0

% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列

IT2=0;a=a+1;

for i=1:n

if i~=isb

%非平衡节点

C(i)=0;D(i)=0;

for j1=1:n

C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)

D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)

end

P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)

Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值

V2=e(i)^2+f(i)^2;

%电压模平方

%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========

if B2(i,6)~=3

%非PV节点

DP=P(i)-P1;

%节点有功功率差

DQ=Q(i)-Q1;

%节点无功功率差

%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================

for j1=1:n

if j1~=isb&j1~=i

%非平衡节点&非对角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df

X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de

X3=X2;

% X2=dp/df X3=dQ/de

X4=-X1;

% X1=dP/de X4=dQ/df

p=2*i-1;q=2*j1-1;

J(p,q)=X3;J(p,N)=DQ;m=p+1;

% X3=dQ/de J(p,N)=DQ节点无功功率差

J(m,q)=X1;J(m,N)=DP;q=q+1;

% X1=dP/de J(m,N)=DP节点有功功率差

J(p,q)=X4;J(m,q)=X2;

% X4=dQ/df X2=dp/df

elseif j1==i&j1~=isb %非平衡节点&对角元

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df

X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de

X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df

p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q

m=p+1;

J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P

J(m,q)=X2;

end

else

%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========

DP=P(i)-P1;

% PV节点有功误差

DV=V(i)^2-V2;

% PV节点电压误差

for j1=1:n

if j1~=isb&j1~=i

%非平衡节点&非对角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

% dP/de

X2=B(i,j1)*e(i)-G(i,j1)*f(i);

% dP/df

X5=0;X6=0;

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV节点有功误差

J(m,q)=X2;

elseif j1==i&j1~=isb %非平衡节点&对角元

X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df

X5=-2*e(i);

X6=-2*f(i);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV节点有功误差

J(m,q)=X2;

end

end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================

for k=3:N0

% N0=2*n（从第三行开始，第一、二行是平衡节点）

k1=k+1;N1=N;

% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U

for k2=k1:N1

% 从k+1列的Jacobi元素到扩展列的△P、△Q

或 △U

J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化

end

J(k,k)=1;

% 对角元规格化K行K列对角元素赋1

%==================== 回代运算

=======================================

if k~=3

% 不是第三行

k > 3

k4=k-1;

for k3=3:k4

% 用k3行从第三行开始到当前行的前一行k4行消去

for k2=k1:N1 %

k3行后各行上三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）

end

%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素

J(k3,k)=0;%当前行第k列元素已消为0

end

if k==N0

%若已到最后一行

break;

end

%================== 前代运算

==================================

for k3=k1:N0

% 从k+1行到2*n最后一行

for k2=k1:N1

% 从k+1列到扩展列消去k+1行后各行下三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算

end

%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素

J(k3,k)=0;%当前行第k列元素已消为0

end

else

%是第三行k=3

%====================== 第三行k=3的前代运算 ========================

for k3=k1:N0

%从第四行到2n行（最后一行）

for k2=k1:N1

%从第四列到2n+1列（即扩展列）

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）

end

%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素

J(k3,k)=0;%当前行第3列元素已消为0

end

end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====

for k=3:2:N0-1

L=(k+1)./2;

e(L)=e(L)-J(k,N);

%修改节点电压实部

k1=k+1;

f(L)=f(L)-J(k1,N);

%修改节点电压虚部

end

%------修改节点电压-----------

for k=3:N0

DET=abs(J(k,N));

if DET>=pr

%电压偏差量是否满足要求

IT2=IT2+1;%不满足要求的节点数加1

end

ICT2(a)=IT2;

%不满足要求的节点数

ICT1=ICT1+1;

%迭代次数 end %用高斯消去法解”w=-J*V" disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n

V(k)=sqrt(e(k)^2+f(k)^2);

%计算各节点电压的模值

sida(k)=atan(f(k)./e(k))*180./pi;

%计算各节点电压的角度

E(k)=e(k)+f(k)*j;

%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);

%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);

%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);

%显示各节点的电压相角 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);

%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);

%显示各节点的注入功率

disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

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

Siz(i)=Si(p,q);

else

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

Siz(i)=Si(p,q);

end

disp(Si(p,q));

SSi(p,q)=Si(p,q);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];

disp(ZF);

disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

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

Sjy(i)=Sj(q,p);

else

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

Sjy(i)=Sj(q,p);

end

disp(Sj(q,p));

SSj(q,p)=Sj(q,p);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];

disp(ZF);

disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p);

disp(DS(i));

DDS(i)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];

disp(ZF);

disp('----------------------');end

如果源程序的运行结果需要作图可用下面的程序

figure(1);subplot(1,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(1,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('节点号');ylabel('节点注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('节点号');ylabel('节点注入无功');grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路号');ylabel('支路首端注入有功');grid on;subplot(2,2,4);bar(Q1);xlabel('支路号');ylabel('支路首端注入无功');grid on;

总结

通过本次课程设计让我有复习了一次潮流计算的相关知识，跟家清晰了什么事潮流计算以及潮流计算的在电力系统的重要性。电力系统的稳定运行状况即是正常运行状况，是指电力系统在稳定运行条件下电压、功率的分布，也称为潮流分布。电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。同时潮流分布也是电力系统规划设计的一项基础工作。

整个计算过程的模型建立并不是十分复杂，但计算过程十分繁琐、计算量相当的大，而且由于枝节太多很容易算错。不过在计算潮流计算的过程中却对以往学过的电力系统分析的相关知识进行了一次较为深入的复习。而且整个计算对计算量的要求很大，锻炼了我们的计算能力。而且对细节的把握也得到了锻炼，做题的精细程度得到了提高。

参考文献

[1]何仰赞, 温增银《电力系统分析》（第三版）［M].华中科技大学，2002 [2]http://baike.baidu.com/view/627420.[3]王守相，刘玉田电力系统潮流计算研究现状--《山东电力技术》1996年05期

[4]何仰赞，温增银.电力系统分析（上册）第三版［M].湖北：华中科技大学出版社，2002 [5] 刘同娟.MATLAB在电路分析中的应用.电气电子教学学报.2002 [6] 西安交通大学等.电力系统计算[M].北京：水利电力出版社，1993.12 [7] 李光琦.电力系统暂态分析[M].北京：水利电力出版社，2002.5 [8]何仰赞，温增银.电力系统分析（下册）第三版［M].湖北：华中科技大学出版社，2002 [9]韦化，李滨，杭乃善，等．大规模水一火电力系统最优潮流的现代内点算法实现[J]．中国电机工程学报，2003．23(6)：13一l8．

[10]Chen Luo—nan，Suzuki Hideki，Katou Kazuo．Mean fieldtheory for optimal power flow[J]．IEEE Transactions OilPower Systems，1997，12(4)：1481·1486

第三篇：基于MATLAB的电力系统潮流计算

基于MATLAB的电力系统潮流计算

%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：

%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写 %对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点 %编号，将变压器的串联阻抗置于低压侧处理。%第三列为支路的串列阻抗参数。%第四列为支路的对地导纳参数。

%第五烈为含变压器支路的变压器的变比

%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。

%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点 %负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中 %“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地 %参数。

n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');Y=zeros(n);Times=1;%置迭代次数为初始值 %创建节点导纳矩阵 for i=1:n1 if B1(i,6)==0 %不含变压器的支路 p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else %含有变压器的支路 p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));end end Y OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%将OrgS、DetaS初始化

%创建OrgS，用于存储初始功率参数 h=0;j=0;for i=1:n %对PQ节点的处理 if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end OrgS %创建PVU 用于存储PV节点的初始电压 PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)==3 t=t+1;PVU(t,1)=B2(i,3);end end PVU %创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1:n %对PQ节点的处理 if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end end t=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3 h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;end end DetaS %创建I，用于存储节点电流参数 i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理 if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==j %对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);end if k==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行

k=0;end end end end end k=0;for i=1:n %对PV节点的处理 if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==j %对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理

Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行

k=0;end end end end end Jacbi %求解修正方程，获取节点电压的不平衡量 DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正节点电压 j=0;for i=1:n %对PQ节点处理 if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n %对PV节点的处理 if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 %开始循环********************************************************************** while abs(max(DetaU))>pr OrgS=zeros(2*n-2,1);%！!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵 h=0;j=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end for i=1:n if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end OrgS %创建DetaS h=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end end t=0;for i=1:n if i~=isb&B2(i,6)==3 h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;end end DetaS %创建I i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %创建Jacbi Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);end if k==(n-1)k=0;end end end end end k=0;for i=1:n if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k==(n-1)k=0;end end end end end Jacbi DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正节点电压 j=0;for i=1:n if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 Times=Times+1;%迭代次数加1 end Times 一个原始数据的例子节点数 5 支路数 5平衡节点编号 5 精度pr 0.000001 B1(支路参数矩阵)[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1] B2(节点参数矩阵)[0-1.6-0.8i 1 0 0 2;0-2-1i 1 0 0 2;0-3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1] X(节点号和对地参数)[1 0;2 0;3 0;4 0;5 0]

电力系统潮流计算

——9结点算例-PQ法

原始数据录入data.txt文档:

标号,起始结点,终止结点,支路电阻参数,支路电抗参数,支路对地导纳参数 1,2,5,0.0,0.063,0.0, 2,5,9,0.019,0.072,0.075, 3,6,9,0.012,0.101,0.105, 4,3,6,0.0,0.059,0.0, 5,6,8,0.039,0.17,0.179, 6,4,8,0.017,0.092,0.079, 7,5,7,0.032,0.161,0.153, 8,4,7,0.01,0.085,0.088, 9,1,4,0.0,0.058,0.0, 潮流程序chaoliu.txt文档: #include #include #define N 9 /*总结点数*/ #define M 6 /*PQ结点数*/ #define K 9 /*线路数*/ #define eps 1e-4 void guass(int n,int m,float c[],float b[][N],float x[])/*高斯函数*/ { float a[N][N],y[N];int i,j,k;for(i=0;i=0;i--){ for(j=i+1;jHeadnode-1;j=t->Endnode-1;lr=t->R;lx=t->X;lb1=t->b;lg=lr/(lr*lr+lx*lx);lb=-lx/(lr*lr+lx*lx);g[i][i]+=lg;g[j][j]+=lg;b[i][i]+=lb+lb1;b[j][j]+=lb+lb1;h[i][j]=h[j][i]=-lb1;g[i][j]-=lg;g[j][i]-=lg;b[i][j]-=lb;b[j][i]-=lb;} getch();printf(“n=====jie dian dao na ju zhen=====n”);for(i=0;i

;

float x1[N-1],x2

;for(i=1;imax)max=fabs(dpu[i]);} if(max>=eps){ for(i=0;i

guass(1,N-1,y1,B,x1);for(i=1;imax)max=fabs(dqu[i]);} if(max>=eps){ for(i=0;i

guass(N-M,M,y2,B,x2);for(i=N-M;i

else { kp=0;if(kq==0)val(u,g,b,r,ku,kr,h);else goto top;} } } void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]){ float ps=0,pv1=0,pv2=0;float qs=0,qv1=0,qv2=0;float p[N][N]={0};float q[N][N]={0};float s[N][N];float dp[N][N]={0};float dq[N][N]={0};float ds[N][N];float dSp=0,dSq=0;int i,j;FILE *fp1;printf(”n=====ping heng jie dian gong lv =====n“);getch();for(i=0;i

printf(”n=======shu ju bao cun=====n“);fp1=fopen(”jieguo.txt“,”w+“);{ fprintf(fp1,”xian lu cao liu:n“);for(i=0;i

第四篇：电力系统仿真MATPOWER潮流计算

IEEE30节点潮流计算

宁夏大学新华学院马智

潮流计算，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。它是基于配电网络特有的层次结构特性，论文提出了一种新颖的分层前推回代算法。该算法将网络支路按层次进行分类，并分层并行计算各层次的支路功率损耗和电压损耗，因而可大幅度提高配电网潮流的计算速度。论文在MATLAB环境下，利用其快速的复数矩阵运算功能，实现了文中所提的分层前推回代算法，并取得了非常明显的速度效益。另外，论文还讨论发现，当变压器支路阻抗过小时，利用Π型模型会产生数值巨大的对地导纳，由此会导致潮流不收敛。为此，论文根据理想变压器对功率和电压的变换原理，提出了一种有效的电压变换模型来处理变压器支路，从而改善了潮流算法的收敛特性。

关键词：电力系统；潮流分析；MATLAB

潮流计算的目的

电力系统的潮流计算最主要的目的是为了让电力系统能够安全稳定运行的同时做到经济运行。所以考留到经及调度、电网规划、电力系统可靠性分析。

具体表现在以下方面：

①在电网规划阶段,通过潮流计算,合理规划电源容量及接入点,合理规划网架,选择无功补偿方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。

②在编制年运行方式时,在预计负荷增长及新设备投运基础上,选择典型方式进行潮流计算,发现电网中薄弱环节,供调度员日常调度控制参考,并对规划、基建部门提出改进网架结构,加快基建进度的建议。

③正常检修及特殊运行方式下的潮流计算,用于日运行方式的编制,指导发电厂开机方式,有功、无功调整方案及负荷调整方案,满足线路、变压器热稳定要求及电压质量要求。

④预想事故、设备退出运行对静态安全的影响分析及作出预想的运行方式调整方案。

总结为在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。

MATLAB软件的应用

MATLAB Compiler是一种编译工具，它能够将M编写的函数文件生成函数库或者可执行文件COM组件等，以提供给其他高级语言如C++、C#等进行调用由此扩展MATLAB的应用范围，将MATLAB的开发效率与其他高级语言的运行结合起来，取长补短，丰富程序开发的手段。

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

MATLAB是一种交互式、面向对象的程序设计语言，广泛应用于工业界与学术界，主要用于矩阵运算，同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。

MATLAB程序设计语言结构完整，且具有优良的移植性，它的基本数据元素

是不需要定义的数组。它可以高效率地解决工业计算问题，特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言，可以用类似数学公式的方式来编写算法，大大降低了程序所需的难度并节省了时间，从而可把主要的精力集中在算法的构思而不是编程上。

另外，MATLAB提供了一种特殊的工具：工具箱（TOOLBOXES）.这些工具箱主要包括：信号处理（SIGNAL PROCESSING）、控制系统（CONTROL SYSTEMS）、神经网络（NEURAL NETWORKS）、模糊逻辑(FUZZY LOGIC)、小波(WAVELETS)和模拟（SIMULATION）等等。不同领域、不同层次的用户通过相应工具的学习和应用，可以方便地进行计算、分析及设计工作。

14611121416***25783***9202422302526

图1 IEEE-30节点系统接线图

总结及感想

通过这次的课程设计，我知道了潮流计算的基本步骤和方法，明白了潮流计算对于电力系统的重要性，准确的潮流计算对于工农业的生产有着十分重要的意义。这次实习忙碌但是充实，在其中我发现了自己的不足，自己知识的很多漏洞，和基础知识不扎实，课外知识知之甚少。看到了自己理论联系实际的能力还需提高，也知道了自己以后学习的方向和目的。这次课程设计对自己意义很大，自己从中获得很多东西。

第五篇：基于MATLAB的电力系统潮流计算设计

关键词：电力系统分析；潮流计算；matlab仿真

中图分类号：tm744 文献标识码：a 文章编号：1006-4311（2016）21-0185-03

0 引言

潮流计算是电力系统稳态运行中的基本计算方法中的一种计算方法，也是电力系统稳态运行中最重要的运算。潮流计算是保证电力系统安全、经济运行的根本。在新电网建设的初期规划中，有了潮流计算，可规划出电源的容量及其接入点，可计算出无功补偿的容量，选择合适的补偿方式，以满足在电网潮流的控制、调压、调相、调峰的交换要求。潮流计算可以选择电力系统的运行方式，便于定期对电力系统中的元件进行检修。

潮流计算的过程

1.1 原始资料

①系统图：两个发电厂分别通过变压器和输电线路与四个变电所相连。（图1）

②发电厂资料：

③变电所资料：

1）变电所1、2、3、4低压母线的电压等级分别为：10kv，35kv，10kv，35kv。

3）每个变电所的功率因数均为cosφ=0.9。

④输电线路资料：

发电厂和变电所之间的输电线路的电压等级及长度标于图中，单位长度的电阻为0.17ω，单位长度的电抗为0.402ω，单位长度的电纳为2.78*10-6s。

1.2 基本要求

①对给定的网络查找潮流计算所需的各元件等值参数，画出等值电路图。

②输入各支路数据，各节点数据，利用simulink搭建仿真模型等方法，进行在变电所的某一负荷情况下的潮流计算及仿真，并对计算结果进行分析。

③如果各母线电压不满足要求，进行电压的调整。（变电所低压母线电压10kv要求调整范围在9.5-10.5之间；电压35kv要求调整范围在35-36之间）。

④利用matlab软件，进行上述各种情况潮流的计算及仿真。

1.3 节点设置及分析

由上述系统图可知，该系统图为双端供电网络。将母线1，2设为节点1，10，将变电所1、2、3、4的高低压侧分别设为节点2、3、4、5、6、7、8、9。并且，将节点1设为平衡节点，将节点10设为pv节点，其余节点设为pq节点。

1.4 参数求取

将参数整理如表

1、表2所示。

1.5 进行潮流计算

图2为仿真模型图。

从潮流计算的结果可得到，系统的各个节点电压的标幺值可归纳为表3。

由matlab编程调节后，可得到表4的发电厂电压和变压器分接头电压得标幺值。

在得到了上述调节后的电压标幺值，对电机模型和变压器模型进行更改。表5为调节前后各节点的电压标幺值。

由题意可知，变电所低压母线电压10kv要求调整范围在9.5-10.5之间；电压35kv要求调整范围在35-36之间。因此我们可以看出，经过调节后，节点3、5、7、9点电压已经满足了系统的要求。表6是电压调节前后对线路损耗进行分析的记录。

由表6的电压调节前后功率损耗对比，可以看出有功功率随着变压器分接头变比的增大而逐渐增大，使得变压器的低压侧的电压处于允许范围内，符合其要求。

表7为调节后的各支路电压首末端的功率整理表

表8为各节点功率s的标幺值。

1.6 对比

由上面的三种方法简单地比较，我们可以看出，在同一个电力系统中，用不同的方法进行潮流计算，所得到的结果是大致相同的。

结束语

电力系统潮流分析计算的MATLAB仿真-周明亮01

第一篇：电力系统潮流分析计算的MATLAB仿真-周明亮01

第二篇：用matlab电力系统潮流计算

第三篇：基于MATLAB的电力系统潮流计算

第四篇：电力系统仿真MATPOWER潮流计算

第五篇：基于MATLAB的电力系统潮流计算设计

相关范文推荐

基于MATLAB的电力系统潮流计算_张宁

基于PSAT的电力系统潮流计算仿真

电力系统潮流计算

电力系统潮流计算的MATLAB辅助程序设计,潮流计算程序

电力系统潮流计算程序设计

电力系统潮流计算程序[本站推荐]

电力系统潮流计算发展史

电力系统潮流计算[大全五篇]