课 程 设 计 任 务 书

课程设计任务书

题目电力系统课程设计

题目C

学 院信息科学与电气工程学院 专业电气工程及其自动化 班级电气 学生姓名 学号

11月9日至11月20日共2周

指导教师(签字)

院 长(签字)

一、设计内容及要求 复杂网络牛顿—拉夫逊法潮流分析与计算的设计 电力系统潮流计算是电力系统中一项最基本的计算，设计内容为复杂网络潮流计算的计算机算法——牛顿-拉夫逊法。 首先，根据给定的电力系统简图，通过手算完成计算机算法的两次迭代过程，从而加深对牛顿-拉夫逊法的理解，有助于计算机编程的应用。 其次，利用计算机编程对电力系统稳态运行的各参数进行解析和计算；编程完成复杂网络的节点导纳矩阵的形成；电力系统支路改变、节点增减的程序变化；编程完成各元件的功率损耗、各段网络的电压损耗、各点电压、功率大小和方向的计算。 二、设计原始资料 试用牛顿-拉夫逊潮流计算方法求解教材P102例题3-5。求解精度为10e-6。

三、设计完成后提交的文件和图表 1．计算说明书部分 设计报告和手算潮流的步骤及结果 四、进程安排 第一天上午：选题，查资料，制定设计方案； 第一天下午——第三天下午：手算完成潮流计算的要求； 第四天上午——第五天上午：编程完成潮流计算，并对照手算结果，分析误差 第五天下午：答辩，交设计报告。 五、主要参考资料 《电力系统分析（第三版）》 于永源主编，中国电力出版社，2007年 《电力系统分析》，何仰赞 温增银编著，华中科技大学出版社，2002年版； 《电力系统分析》，韩桢祥主编，浙江大学出版社，2001年版； 《电力系统稳态分析》，陈珩 编，水利电力出版社；

课程设计成绩评定用表 平时成绩 答辩成绩 报告成绩 总成绩

目录

摘 要.............................................................. 2 1 任务提出与方案论证................................................ 3

1.1 潮流计算题目 ................................................ 3 1.2 对课题的分析及求解思路 ...................................... 4 1.3潮流计算要求................................................. 4 1.4直角牛顿-拉夫逊法潮流计算求解过程............................ 4 2 总体设计.......................................................... 9

2.1 潮流计算算法 ................................................ 9 2.2 潮流计算的手工计算 ......................................... 14 3设计过程 ......................................................... 16

3.1 Matlab ..................................................... 12 3.2 潮流计算流程图 ............................................. 13 3.3 潮流计算源程序图及结果 ..................................... 15 4 心得体会......................................................... 28

摘 要

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

关键字： 潮流计算 牛顿-拉夫逊 Matlab软件

1 任务提出与方案论证

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

1.1 潮流计算题目

1.2 对课题的分析及求解思路

此电力系统是一个5节点，5支路的电力网络。其中1节点为平衡节点，2、3节点为PV节点，4、5节点为PQ节点。综合比较牛顿拉夫逊法（直角坐标、极坐标）、PQ分解法等多种求解方法的特点，最后确定采用牛顿拉夫逊法（直角坐标）。因为此方法所需解的方程组通俗易懂。

1.3潮流计算要求

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

Uimin?Ui?Uimax(i?1,2,?n)

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

2.节点的有功功率和无功功率应满足

PGimin?PGi?PGimax?Q?Q?

Gimin?QGiGimax?PQ节点的有功功率和无功功率，以及PU节点的有功功率，在给定是就必须满足上述条件，因此，对平衡节点的P和Q以及PU节点的Q应按上述条件进行检验。

3.节点之间电压的相位差应满足

|?ij|?|?i??j|?|?i??j|max

因此，潮流计算可以归结为求解一组非线性方程组，并使其解答满足一定的约束条件。常用的方法是迭代法和牛顿法，在计算过程中，或得出结果之后用约束条件进行检验。如果不能满足要求，则应修改某些变量的给定值，甚至修改系统的运行方式，重新进行计算。 1.4直角牛顿-拉夫逊法潮流计算求解过程

以下是用直角坐标形式的牛顿—拉夫逊法潮流的求解过程。当采用直角坐标时，潮流问题的待求量为各节点电压的实部和虚部两个分量

e,f,e,f...e,f1122nn由于平衡节点的电压向量是给定的，因此待求两共2(n?1)需要2(n-1)个方程式。事实上，除了平衡节点的功率方程式在迭代过程中没有

和Qis约束作用以外，其余每个节点都可以列出两个方程式。对PQ节点来说，Pis是给定的，因而可以写出

?Pi?ip?e?(Ge?Bfisij?iijjijisij?iijjijj)?jf?(Gfjj?iijj?ij?Q?Q?f?(Ge?Bf)?ej?(Gijf?Bijej)?0??? （2-2-1）

?Bijej)?0?j?对PV节点来说，给定量是Pis和Vis，因此可以列出

?Pi?Pis?ei?(Gijej?Bijj?i2fj)?f?(Gfij?iijj?Vi?Vis?(ei?22f2i)?0?Bijej)?0?? ? （2-2-2）

??

求解过程大致可以分为以下步骤： （1）形成节点导纳矩阵Y （2）将各节点电压设初值U，

（3）将节点初值代入式(2-2-1)或式(2-2-2)，求出修正方程式的常数项向量

（4）将节点电压初值代入求式，求出雅可比矩阵元素 （5）求解修正方程，求修正向量 （6）求取节点电压的新值

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

（8）计算支路功率分布，PV节点无功功率和平衡节点柱入功率。 以直角坐标系形式表示 迭代推算式

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

??e?jfViii (2-2-3)

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

⑴ 于PQ节点

nn?PG?i?Pi?ei?(ijej?Bijfj)?fi?(Gijfj?Bijej)j?1j?1??nn? (2-2-4) ?Qi?Qi?fi?(Gijej?Bijfj)?ei?(Gijfj?Bijej)?j?1j?1??i?1,2,?,m

⑵ 于PV节点

nn?P?P?(Ge?ii?eiijej?Bijfj)?fi?(Gijfj?Bijj)j?1j?1??2(2-2-5)

?V2V2?I?i?(ei?f2i)?i?m?1,m?2,?,n?1

⑶对于平衡节点

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

Vn?en?jfn (2-2-6)

①. 修正方程式迭代式共包括2(n-1)个方程.选定电压初值及变量

修正量符号,代入方程并按泰勒级数展开,略去?ei,?fi二次方程及以后各项,得到一组线性方程组或线性化了的方程组，常称修正方程组:

?W??J?U (2-2-7)

???P1???e1???Q??1???f?1????????P?????m???em??W????Qm???P??U????fm??m?1???e?m??U2m?1???1? ?????fm?1?????????P?n?1????e???U2?n?1n?1??????fn?1??

???P1??e1????Q1??e1???????Pm??e1????Qm??e1J????Pm?1??e1????U2m?1??e1???????Pn?1??e1?2???Un?1??e1???P1?f1??Q1?f1???Pm?f1??Qm?f1?Pm?1?f1????????P1?em??Q1?em???Pm?em??Qm?em?Pm?1?em??P1?fm??Q1?fm???Pm?fm??Qm?fm?Pm?1?fm??U2m?1?fm???Pn?1?fm??U2n?1?fm??P1?em?1??Q1?em?1???Pm?em?1??Qm?em?1?Pm?1?em?1??U2m?1?em?1???Pn?1?em?1??U2n?1?em?1??P1?fm?1??Q1?fm?1???Pm?fm?1??Qm?fm?1?Pm?1?fm?1????????P1?en?1??Q1?en?1???Pm?en?1??Qm?en?1?Pm?1?en?1??U2m?1??U2m?1??f1?em???Pn?1?f1?????Pn?1?em??U2m?1??U2m?1??fm?1?en?1???Pn?1?fm?1?????Pn?1?en?1??U2n?1??U2n?1??f1?em??U2n?1??U2n?1??fm?1?en?1??P?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????U2n?1??fn?1??(2-2-8)

③.雅可比矩阵各元素的算式

式(2-2-8)中, 雅可比矩阵中的各元素可通过对式(2-2-4)和(2-2-5)进行偏

导而求得.当j?i时, 雅可比矩阵中非对角元素为

???Pi??Qi????(Gijei?Bijfi)???ej??fj????Pi??Qi???Bijei?Gijfi? (2-2-9) ??fj??ej????U2??U2???0?ej?fj??当j?i时,雅可比矩阵中对角元素为:

??P?ei?ni??(G?ijej?Bijfj)?Giiei?Biifij?1????Pni?f???(G?ijfj?Bijej)?Giifi?Biiei?jj?1???Qni?e??(G?ijfj?Bijej)?Giifi?Biiei?ij?1????Qn? (2-2-10) i?f???(Gij?ej?Bijfj)?Giiei?Biifi?jj?1???U2?i?e??2e?ij????U2i?f??2f?i?i??由式(2-2-9)和(2-2-10)看出,雅可比矩阵的特点:

⒈矩阵中各元素是节点电压的函数,在迭代过程中,这些元素随着节点电压

的变化而变化;

⒉导纳矩阵中的某些非对角元素为零时,雅可比矩阵中对应的元素也是为零.若Yij?0,则必有Jij?0;

⒊雅可比矩阵不是对称矩阵;(i?q1,2,?,n;i?s) 雅可比矩阵各元素的表示如下式(2-2-11):

??(Gijei?Bijfi)H??Pi??e?(j?i)ij??j??(Ge?Bf)?Ge?B

??ijjijjiiiiifi(j?i)j?i

?Bijei?Gijfi)N??Pi?(j?i)ij??f??j??(Gf?Be)?B

??ijjijjiiei?Giifi(j?i)j?i

BijeiM???Q??Gijfi)i?e??(j?i)ij?j?(Gf?Be)?

??ijjijjBiiei?Giifi(j?i)j?iGijeL??Q?i?Bijfi)i?(j?i)ij??f??j??(Ge?Bf)?Ge?Bf

??ijjijjiiiiii(j?i)j?i

??U2i?0(j?i)??U2i?0(j?i)Rij???Sij???

?ej?fj??2ei(j?i)??2fi(j?i)

2 总体设计

在电力系统潮流计算中，往往要计算节点导纳矩阵，而我们计算节点导纳矩阵采用节点电压法来实现，在利用电路知识列节点电压方程，从而导出所需的导

纳矩阵。 2.1节点导纳矩阵

节点导纳矩阵既可根据自导纳和互导纳的定义直接求取，也可根据电路知识中找出改网络的关联矩阵，在节点电压方程的矩阵形式进行求解。本章节我们主要讨论的是直接求解导纳矩阵。根据节点电压方程章节我们知道，在利用电子数字计算机计算电力系统运行情况是，多采用I?YV形式的节点方程式。其中阶数等于电力网络的节点数。从而可以得到n个节点时的节点导纳矩阵方程组（2-11）如下：

?Y11V1?Y12V2?......Y1nVn?I1?....?Y21V1?Y22V2?......Y23Vn?I2??.......? (2-11） ....?Yn1V1?Yn2V2?......Yn3Vn?In??

由此可以得到n个节点导纳矩阵：

....?Y11Y12...Y1n??Y?Y22...Y2n?21?Y??? ......???Yn1Yn2...Ynn?它反映了网络的参数及接线情况，因此导纳矩阵可以看成是对电力网络电气特性的一 种数学抽象。由导纳短阵所联系的节点方程式是电力网络广泛应用的一种数学模

型。

通过上面的讨论，可以看出节点导纳矩阵的有以下特点： （1）导纳矩阵的元素很容易根据网络接线图和支路参数直观地求得，形成节点导纳矩阵的程序比较简单。

（2）导纳矩阵为对称矩阵。由网络的互易特性易知Yij?Yji。

（3）导纳矩阵是稀疏矩阵。它的对角线元素一般不为零，但在非对角线元素中则存在不少零元素。在电力系统的接线图中，一般每个节点与平均不超过3~4个其他节点有直接的支路连接。因此，在导纳矩阵的非对角线元素中每行仅有3~4个非零元素，其余的都是零元素，而且网络的规模越大，这种现象越显著。

导纳矩阵的对称性和稀疏性对于应用计算机求解电力系统问题有很大的影响。如果能充分地利用这两个特点，如在程序设计中储存导纳矩阵的对角元素和上三角元素（或下三角元素），排除零元素的储存和运算，就可以大大地节省储存单元和提高计算速度。

节点导纳矩阵的形式可归纳如下：

（1）导纳矩阵的阶数等于电力网络的节点数。

（2）导纳矩阵各行非对角元素中非零元素的个数等于对应节点所连得不接地支路数。

（3）导纳矩阵各对角元素，即节点的自导纳等于相应节点之间的支路导纳之和。

（4）导纳矩阵非对角元素，即节点之间的互导纳等于相应节点之间的支路导纳的负值。

而在电力系统中进行潮流计算时，往往要计算不同接线下的运行状况，例如，改变变压器主抽头时，潮流分布也随之变化，以及改变其他设备参数进行计算潮流分布，此时就需要导出变化时的导纳矩阵就需要对所设计的程序进行参数设定，而不需要重复上述步骤去导出所求的导纳矩阵。 2.2 潮流计算的手工计算

1 所求电力系统网络的各个节点的相关参数如下： y11=1/(0.0194+0.0592j)+1/(0.054+0.223j); y12=y21=-1/(0.0194+0.0592j);

y22=1/（0.04699+0.198j)+0.067j+1/(0.0581+0.1763j)); y23=y32=1/(0.04699+0.198j); y24=y42=1/(0.0581+0.1763j); y33=1/(0.04699+0.198j)+0.022j; y44=1/（0.0581+j0.1763）+j0.0187；

y55=1/(0.054+0.233j)+j0.0246;

2 设定各节点的电压

取U1=1.06+j0 U2=1.045+j0 U3=1.01+j0 U4=1.00+j0 U5=1.00+j0

??Pi?Pi?ei?(Gijej?Bijfj)?fi?(Gijfj?Bijej)?j?1j?1??nn?Qi?Qi?fi?(Gijej?Bijfj)?ei?(Gijfj?Bijej)??j?1j?1?

nn3.计算雅可比矩阵中各元素

当j?i时, 雅可比矩阵中非对角元素为

???Pi??Qi????(Gijei?Bijfi)???ej??fj????Pi??Qi???Bijei?Gijfi???fj??ej????U2??U2???0?ej?fj??

当j?i时,雅可比矩阵中对角元素为:

n??P?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??

3设计过程

3.1 Matlab

MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB不仅仅可以实现计算潮流分布，还可以模拟仿真各式各样的数值系统，工程。这里将借助MATLAB来完成用直角牛顿-拉夫逊法计算电力系统潮流分布。 与常数的运算

常数与矩阵的运算即是同该矩阵的每一元素进行计算。但需注意进行数除时，常数通常只能做除量。 基本函数运算中，矩阵的函数运算是矩阵预算中最实用的部分，常用的主要有以下几个：

det(a) 求矩阵a的行列式 eig(a) 求矩阵a的特征值 inv(a)或a^(-1) 求矩阵a的逆矩阵 rank(a) 求矩阵a的秩

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

数组的加、减与矩阵的加、减运算完全相同。而乘除法运算有相当大的区别，数组的乘除法是指两同维数组对应元素之间的乘除法，它们的运算符为“.*”和“./”或“.\\。”前面讲过常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定，数组与常数之间的除法运算没有任何限制。 另外，矩阵的数组运算中还有幂运算（运算符.^）、指数运算（exp）、对数运算(log)、和开方运算（sqrt）、等，有了“对应元素”的规定，数组的运算实质上就是针对数组内部的每个元素进行的。矩阵的幂运算与数组的幂运算有很大的区别。

逻辑关系运算

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

1）一般MATLAB命令格式为 [输出参数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矩阵 a=[1 2 3;4 5 6;7 8 9] a = 1 2 3 4 5 6 7 8 9 3.2 潮流计算流程图

在上一节，我们对MATLAB的基本指令有了初步的了解，我们就可以画出程序编程框图，按照程序编程框图，我们可以通过MATLAB软件进行编程，实现我们所求的节点导纳矩阵。

程序编程框图如下：

启动 输入原始数形成节点导纳矩分解各节点初始电压的实部和虚部 迭代次数K=0 求PQ节点的?Pi(k) ,?Qi(k)，求PV节点的?Pi(k)，?Ui(k) 置节点号i=0 是 雅克比矩阵是否形成，i>n 否 求得雅克比矩阵个元素?Hij(k)?Nij(k)?Jij(k)?Lij(k)?Rij(k)?Sij(k) 把雅克比矩阵单位增大节点号i=i+1 求解修正方程，得?ei(k)，?fi(k) 回带各电压新值，K=K+1 求解最大修正量|?ei(k)|max，|?fi(k)|max 否 是 是否收敛 计算输出电压大小及相角，节电功率及支路损耗 停止

3.3 潮流计算源程序图及结果

%本程序 的功能是用牛顿——拉夫逊法进行潮流计算

% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳 % 5、支路的变比；6、支路首端处于K侧为1，1侧为0不含变压器；

% 7、支路末段对地电纳

% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值

% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量

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

clear;

n=5;%input('请输入节点数：n='); nl=4;%input('请输入支路数：nl=');

isb=1;%input('请输入平衡母线节点号：isb='); pr=0.000001;%input('请输入误差精度：pr=');

B1=[1 2 0.0194+0.0592i 0 1 0 0.067i;

2 3 0.04699+0.198i 0 1 0 0.022i;

1 5 0.054+0.223i 0 1 0 0.0246i;

2 4 0.0581+0.1763i 0 1 0 0.0187i]

B2=[0 0 1.06 1.06 0 1; 0 0.217+0.121i 1.045 1.045 0 3; 0 0.749+0.19i 1.01 1.01 0 3; 0 0.658-0.039i 1 0 0 2; 0 0.076+0.016i 1 0 0 2] ;%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,7); %对角元K侧 Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4); %对角元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); end

for i=1:n

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

%PV节点电压给定模值

%给定各节点注入功率

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;

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

while IT2~=0 %继续运行 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 -Nij

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

X3=X2; % X2=dp/df X3=dQ/de -Lij

X4=-X1; % X1=dP/de X4=dQ/df -Jij

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

J(p,q)=X1;J(p,N)=DP;m=p+1; % X1=dP/de J(m,N)=DP节点有功功率差

J(m,q)=X3;J(m,N)=DQ;q=q+1; % X3=dQ/de J(p,N)=DQ节点无功功率差

J(p,q)=X2;J(m,q)=X4; % 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)=X1;J(p,N)=DP; %扩展列△P=-Nij m=p+1;

J(m,q)=X3;q=q+1;J(p,q)=X2;J(m,N)=DQ; %扩展列△Q=-Hij

J(m,q)=X4;

end 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)=X1;J(p,N)=DP; % PV节点有功误差 m=p+1;

J(m,q)=X5;J(m,N)=DV;q=q+1;J(p,q)=X2; % PV节点电压误差 J(m,q)=X6; 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)=X1;J(p,N)=DP; %PV节点有功误差 m=p+1;

J(m,q)=X5;J(m,N)=DV;q=q+1;J(p,q)=X2; % PV节点电压误差 J(m,q)=X6; end end end end

end

%求JACOB矩阵

disp('JACOB矩阵 J='); disp(J)

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

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 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 end

ICT2(a)=IT2; %不满足要求的节点数 ICT1=ICT1+1; %迭代次数 end

%用高斯消去法解\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))+(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))+(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,7))+(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,7))+(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;

4 心得体会

潮流计算时电力系统的最基本、最常用的计算。根据系统给定的运行条件、网络接

线及元件参数，通过潮流计算可以确定各母线的电压（相角及幅值）、各元件中流过的功率、整个系统的功率损耗等。潮流计算是实现电力系统安全经济发展的必要手段和重要环节。 实际的情况远比我们计算的情况复杂，这让我深刻了解了潮流计算的重要性。精准的潮流计算不仅可以节省能源并充分发挥能源的作用，还可以使电网处于稳定且平衡的状态，这对电网的安全运行起到关键性的作用。所以我认为学好电力系统的关键在于学好潮流计算，电力系统的安全运行不是儿戏，一个很小的失误就可能造成整个电网的崩溃，这给国家和社会带来的危害和损失之大是可想而知的。学好电力系统这门课、学好潮流计算是我们电专业学生的当务之急。没有坚实的基础就不可能建造出高楼大厦，所以，要想在以后的工作岗位上有所成就，就必须打好基础，努力提高自己的专业素养。

参考文献

[1].陈珩.《电力系统稳态分析》[M],中国电力出版社，2007，第三版 [2].韩祯祥.《电力系统分析》[M]，浙江大学出版社，2005，第三版

[3].祝书萍.《电力系统分析课程实际设计与综合实验》[M]，中国电力出版社，2007，第一版

[4].《电力系统稳态分析》，百度文库

[5]. 陶海英. 基于MATLAB的电力系统稳态仿真分析[J],江西电力职业技术学院学报, 2010年 03期

本文来源：https://www.bwwdw.com/article/m2dg.html