两机五节点网络潮流计算方法牛拉法和pq法电力系统稳态分析课程设计报告书

电力系统稳态分析课程设计

题目名称两机五节点网络潮流计算方法牛拉法和pq法

目 录

摘

要…………………………………………………….............................................

第

一章 原理简

介………………………………………………………..................…3 1.1对潮流分析的简介…………………………………………………..……………3 1.1.1

潮

流

计

算

方

法

分

析

比

较……………………………………………...…….3

1 / 33

1.2 MATLAB简介…………………………………………………………….………4

1.2.1 矩阵的运算…………………………………………………………………5 1.3牛顿拉夫逊法计算潮流分布……………………………………………………..6 第二章 程序与结果…………………………………………………………..……..10 2.1 设计资料与参数………………………………………………………………...10 2.1.1 牛顿拉夫逊法的程序框图………………………………………………13

2.2 用Matlab设计程序……………………………………………………………..14

2.2.1程序的编写………………………………………………………………14 2.2.2程序运行结果…………………………………………………...………..19 2.2．3p_q法程序编写………………………………………………………………22 总结……………………………………………………………………..………32

参考文献………………………………………………………………….…….32

2 / 33

电力系统稳态分析课程设计

1.1对潮流分析的简介

潮流分析是研究电力系统的一种最基本和最重要的计算。最初，电力系统潮流计算是通过人工手算的，后来为了适应电力系统日益发展的需要，采用了交流计算台。随着电子数字计算机的出现，1956 年 Ward等人编制了实际可行的计算机潮流计算程序。这样，就为日趋复杂的大规模电力系统提供了极其有力的计算手段。经过几十年的时间，电力系统潮流计算已经发展得十分成熟。潮流计算是研究电力系统稳态运行情况的一种计算，是根据给定的运行条件与系统接线情况确定整个电力系统各个部分的运行状态，如各母线的电压、各元件中流过的功率、系统的功率损耗等等。电力系统潮流计算是计算系统动态稳定和静态稳定的基础。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用电力系统潮流计算来定量的比较供电方案或运行方式的合理性、可

3 / 33

靠性和经济性。

1.1.1潮流计算方法分析比较

高斯 - 赛德尔潮流计算法原理简单，编程实现容易，特别是对于配网潮流有其独特优势。但是高斯 - 赛德尔潮流计算法在牛顿法以与各种解耦法出现以后似乎成了一种边缘性的方法。牛顿 - 拉夫逊法的优点是收敛速度快，若初值选择较好，算法将具有平方收敛特性，一般迭代 4~5 次便可以收敛到一个非常精确的解，而且其迭代次数与所计算的网络规模基本无关。牛顿 - 拉夫逊法也具有良好的收敛可靠性 对于呈病态的系统，牛顿-拉夫逊法均能可靠地收敛。牛顿法的缺点是每次迭代的计算量和所需的存量较大。这是因为雅可比阵元素的数目约为 2(n- 1)×2(n- 1)个（直角坐标），且其数值在迭代过程中不断变化。不过，存占用量与每次迭代所需的时间与程序设计技巧密切相关。牛顿-拉夫逊法的可靠收敛取决于一个良好的启动初值，如果初值选择不当，算法有可能根本不收敛或收敛到一个无法运行的解点上。对于正常运行的系统，各节电电压一般均在额定值附近，偏移不会太大，并且各节电的相角差也不大，所以对各节电可以采用统一的电压初值。P- Q分解法是为了改进牛顿-拉夫逊法在存占用量与计算速度方面的不足，P- Q分解法根据电力系统实际运行状态的物理特点，对极坐标形式的牛顿 - 拉夫逊法修正方程式进行了合理的简化。它无论在存占用量还是计算速度方面都比牛顿-拉夫逊法有较大的改进，是目前计算速度最快的潮流算法。

1.2 Matlab 的简介

4 / 33

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

1.2.1 矩阵的运算

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

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

5 / 33

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

trace(a) 求矩阵a的迹（对角线元素之和） 我们在进行工程计算时常常遇到矩阵对应元素之间的运算。这种运算不同于前面讲的数算，为有所区别，我们称之为数组运算。数组的加、减与矩阵的加、减运算完全一样。而乘除法运算有相当大的区别，数组的乘除法是指两同维数组对应元素之间的乘除法，它们的运算符为“.*”和“./”或“.\\”。前面讲过常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定，数组与常数之间的除法运算没有任何。另外，矩阵的数组运算中还有幂运算（运算符为.^ ）、指数运算（exp）、对数运算（log）、和开方运算（sqrt）等。有了“对应元素”的规定，数组的运算实质上就是针对数组部的每个元素进行的。矩阵的幂运算与数组的幂运算有很大的区别。

1.3牛顿拉夫逊法计算潮流分布

一、牛顿拉夫逊法求解过程大致可以分为以下步骤： （1）形成节点导纳矩阵 （2）将各节点电压设初值U

（3）将节点初值代入相关求式，求出修正方程式的常数项向量 （4）将节点电压初值代入求式，求出雅可比矩阵元素

6 / 33

（5）求解修正方程，求修正向量 （6）求取节点电压的新值

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

（8）计算支路功率分布，PV节点无功功率和平衡节点柱入功率。 二、直角坐标系计算

(1)牛顿-拉夫逊法潮流计算的公式。把牛顿法用于潮流计算，采用直角坐标形式。其中电压和支路导纳可表示为：

VieijfiYijGijjBij（1-1） (2)PQ节点的输出有功功率和无功功率是给定的，则第i节点的给定功率设为P。假定系统中的第1、2、…、m节点为PQ节点，is和Qis（称为注入功率）对其中每一个节点的N-R法表达式:

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1nnQiQisQiQisfi(GijejBijfj)ei(GijfjBijej)0j1j1nn

i=（1、2、…、m） （1-2）

(3)PV节点的有功功率和节点电压幅值是给定的。假定系统中的第

m+1、m+2、…、n-1节点为PV节点，则对其中每一PV节点可以列写方程：

7 / 33

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j122U2UisUi2Uis(ei2fi2) （1-3）

i=（m+1、m+2、…、n-1）

nn(4)形成雅可比矩阵。 当j=i时，对角元素为

nPi(GijejBijfj)GiieiBiifiNiieij1nPi(GijfjBijej)BiieiGiifiHiifij1nQi(GijfjBijej)BiieiGiifiLiieij1nQi(GijejBijfj)GiieiBiifiiJii

fij12Ui2eieiUi22fifi（1-4）

当ji时,矩阵非对角元素为：

8 / 33

PiQi(GijeiBijfi)NijJijeifjPiQiBijeiGijfiHijLij （1-5） fjejUi2Ui20ejfj

三、极坐标计算

对于潮流计算中待于求出功率的节点功率方程组，在某个近似解附近用泰勒级数展开略去二阶与以上的高阶得到已矩阵表示的修正方程：

式中节点n为节点数，m为PV节点，雅阁比矩阵是（2n-2m-2）阶非奇异矩阵，雅可比矩阵各元素表示如下：

9 / 33

牛顿拉夫逊极坐标潮流计算的修正方程的迭代方程为：

10 / 33

2.1设计资料与参数：

课程名称 设计题目 指导教师 一、教学要求 电力系统稳态分析课程设计 两机五节点网络潮流计算—牛拉法 景霞 时间 1周 电力系统稳态分析课程设计以设计和优化电力系统的潮流分析为重点，提高学生综合能力为目标，尽可能结合实际工程进行。设计容的安排要充分考虑学校现有的设备，设计时间与工程实际需要，并使学生初步学会运用所学知识解决工程中的实际问题。 二、设计资料与参数 （一）设计原始资料 11 / 33

1、待设计电气设备系统图 2、电力系统网络各元件参数 3、电力系统电气元件的使用规 4、电力工程电气设计手册 （二）设计参考资料 1、《电力系统稳态分析》，珩，中国电力，2007，第三版 2、《电力系统分析》，祯祥，大学，2005，第三版 3、《电力系统分析课程实际设计与综合实验》，祝书萍，中国电力，2007，第一版 三、设计要求与成果 1.根据给定的参数或工程具体要求，收集和查阅资料；学习相关软件（软件自选）。 2.在给定的电力网络上画出等值电路图。 3.运用计算机进行潮流计算。 4.编写设计说明书。 基本要求： 1.编写潮流计算程序； 2.在计算机上调试通过（？）； 3.运行程序并计算出正确结果（？）； 4.写出课程设计报告（包括以下容）（1份） （1）程序框图；（2）源程序；（3）符号说明表；（4）算例与计算结果 5.编写计算说明书（1份）。 四、进度安排 根据给定的参数或工程具体要求，收集和查阅资料（半天） 学习软件（MATLAB或C语言等）（一天半） 编程计算复杂系统潮流计算（三天） 编写计算设计书（一天） 五、评分标准 课程设计成绩评定依据包括以下几点： 1）工作态度（占10%）； 2）基本技能的掌握程度（占20%）； 3）程序编写是否合理是否有运行结果(40%)； 4）课程设计说明书编写水平(占30%)。 5）分为优、良、中、合格、不合格五个等级。

12 / 33

系统接线图

其中节点1为平衡节点，节点2、3、4、5为PQ节点。

2.1.1 牛顿拉夫逊法程序框图

13 / 33

电力系统稳态分析课程设计

14 / 33

2.2.1 程序的编写

%说明：为了使节点按照先PQ，再PV节点，最后平衡节点的次序编号，以便与公式对照，节点1与节点5对调。 %节点阻抗矩阵

Z=[0,0.04+0.12i,0,0.08+0.24i,0;0.04+0.12i,0,0.06+0.18i,0.06+0.18i,0.02+0.06i;0,0.06+0.18i,0,0.01+0.03i,0.08+0.24i;0.08+0.24i,0.06+0.18i,0.01+0.03i,0,0;0,0.02+0.06i,0.08+0.24i,0,0]; %求互导纳 for m=1:5 for n=1:5

if Z(m,n)==0 y(m,n)=0; else

y(m,n)=1/Z(m,n); end end end %求导纳 for m=1:5 for n=1:5 if m~=n

Y(m,n)=-y(m,n); else

Y(m,n)=sum(y(m,:)); end end end %导纳矩阵

G=real(Y);B=imag(Y);

%计算各节点功率的不平衡量设U=E+jF ；Y=G+Bj； Y

U=[1,1,1,1,1.06]; E=real(U); F=imag(U);

15 / 33

%设S=P+Qj；

S=[-0.60-0.10i;0.20+0.20i;-0.45-0.15i;-0.40-0.05i;0] P=real(S);Q=imag(S);

k=0;C=1;

while C > 0.00001

E(5)=1.06;F(5)=0; for m=1:4 for n=1:5

%计算Pi,Qi,设Pi=Pt;Qi=Qt,按照书上的公式：

Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n))+F(m)*(G(m,n)*F(n)+B(m,n)*E(n))); Qt(n)=(F(m)*(G(m,n)*E(n)-B(m,n)*F(n))-E(m)*(G(m,n)*F(n)+B(m,n)*E(n))); end

%设P,Q的改变量为dP,dQ dP(m)=P(m)-sum(Pt); dQ(m)=Q(m)-sum(Qt); end

%计算Hii,Nii,Jii,Lii,由公式4-41b 左侧公式实现,sum(Ai)，sum(Bi)用于实现公式中的sigerma从j到n的求和; for m=1:4 for n=1:5

Bi(n)=G(m,n)*F(n)+B(m,n)*E(n);

Ai(n)=G(m,n)*E(n)-B(m,n)*F(n); end

H(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))+2*G(m,m)*F(m); N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m))+2*G(m,m)*E(m); J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)); L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))); end

%设雅可比矩阵为JJ，以下语句用来实现雅可比矩阵中对角线上元素H N J L 的排列 for m=1:4

JJ(2*m-1,2*m-1)=H(m,m); JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m); JJ(2*m,2*m)=L(m,m);

16 / 33

end

%以下语句用于实现雅可比矩阵非对角线上元素的排列 for m=1:4 for n=1:4 if m==n else

H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m); N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m); L(m,n)=G(m,n)*F(m)-B(m,n)*E(m);

JJ(2*m-1,2*n-1)=H(m,n); JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n); JJ(2*m,2*n)=L(m,n); end end end JJ

%设由P,Q的改变量组成的8×1矩阵为PQ，由E,F的改变量组成的8×1矩阵为dU for m=1:4

PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m); end

dU=inv(JJ)*PQ';%求逆矩阵

C=max(abs(dU));%abs绝对值或复数绝对值 for n=1:4

F(n)=F(n)+dU(2*n-1); E(n)=E(n)+dU(2*n); end

for n=1:5

U(n)=E(n)+(F(n))*i; end

k=k+1;

k, dU=dU',PQ,U end

%计算S(5),也就是题目中的S1，即平衡节点功率

17 / 33

for m=1:5

I(m)=Y(5,m)*U(m); end

S(5)=U(5)*sum(conj(I)) Y =

3.7500 -11.2500i -2.5000 + 7.5000i 0 -1.2500 + 3.7500i 0

-2.5000 + 7.5000i 10.8333 -32.5000i -1.6667 + 5.0000i -1.6667 + 5.0000i -5.0000 +15.0000i

0 -1.6667 + 5.0000i 12.9167 -38.7500i -10.0000 +30.0000i -1.2500 + 3.7500i

-1.2500 + 3.7500i -1.6667 + 5.0000i -10.0000 +30.0000i 12.9167 -38.7500i 0

0 -5.0000 +15.0000i -1.2500 + 3.7500i 0 6.2500 -18.7500i S =

-0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.1500i -0.4000 - 0.0500i 0 JJ =

11.2500 3.7500 -7.5000 -2.5000 0 0 -3.7500 -1.2500

-3.7500 11.2500 2.5000 -7.5000 0 0 1.2500 -3.7500

-7.5000 -2.5000 33.4000 10.5333 -5.0000 -1.6667 -5.0000 -1.6667

2.5000 -7.5000 -11.1333 31.6000 1.6667 -5.0000 1.6667 -5.0000

0 0 -5.0000 -1.6667 38.9750 12.8417 -30.0000

18 / 33

-10.0000

0 0 1.6667 -5.0000 -12.9917 38.5250 10.0000 -30.0000

-3.7500 -1.2500 -5.0000 -1.6667 -30.0000 -10.0000 38.7500 12.9167

1.2500 -3.7500 1.6667 -5.0000 10.0000 -30.0000 -12.9167 38.7500 k =

1 dU =

-0.1076 0.0093 -0.0473 0.0430 -0.0863 0.01 -0.0922 0.0141 PQ =

-0.6000 -0.1000 0.5000 1.1000 -0.3750 0.0750 -0.4000 -0.0500 U =

1.0093 - 0.1076i 1.0430 - 0.0473i 1.01 - 0.0863i 1.0141 - 0.0922i 1.0600 JJ =

11.0516 4.3956 -7.3011 -3.3304 0 0 -3.6505 -1.6652

-5.5956 10.8516 3.3304 -7.3011 0 0 1.6652 -3.6505

-7.7040 -2.9621 33.1594 13.0913 -5.1360 -1.9747 -5.1360 -1.9747

19 / 33

2.9621 -7.7040 -12.5804 33.6083 1.9747 -5.1360 1.9747 -5.1360

0 0 -4.9331 -2.1238 38.3848 16.0299 -29.5988 -12.7427

0 0 2.1238 -4.9331 -16.8888 38.0788 12.7427 -29.5988

-3.6876 -1.6135 -4.9168 -2.1513 -29.5009 -12.9079 38.1553 16.2727

1.6135 -3.6876 2.1513 -4.9168 12.9079 -29.5009 -17.0727 38.0553 k =

2 dU =

0.0032 -0.0131 -0.0004 -0.0075 0.0017 -0.0101 0.0021 -0.0108 PQ =

0.01 -0.0636 -0.0770 -0.0220 -0.0008 -0.0317 0.0103 -0.0362 U =

0.9963 - 0.1044i 1.0355 - 0.0477i 1.0053 - 0.0846i 1.0033 - 0.0902i 1.0600 JJ =

10.9774 4.3253 -7.2110 -3.2738 0 0 -3.6055 -1.6369

-5.4960 10.6556 3.2738 -7.2110 0 0 1.6369

20 / 33

-3.6055

-7.66 -2.9466 32.9334 12.9530 -5.0977 -1.94 -5.0977 -1.94

2.9466 -7.66 -12.5843 33.3371 1.94 -5.0977 1.94 -5.0977

0 0 -4.8857 -2.0983 38.0494 15.8299 -29.3142 -12.58

0 0 2.0983 -4.8857 -16.6938 37.6790 12.58 -29.3142

-3.99 -1.5922 -4.8665 -2.1230 -29.1988 -12.7380 37.7997 16.0621

1.5922 -3.99 2.1230 -4.8665 12.7380 -29.1988 -16.8443 37.6306 k =

3 dU =

1.0e-003 *

0.0241 -0.1731 -0.0025 -0.02 0.0098 -0.1231 0.0128 -0.1327 PQ =

1.0e-003 *

-0.0034 -0.8423 -0.5236 -0.2008 -0.0815 -0.3198 0.0241 -0.3906 U =

0.9961 - 0.1044i 1.03 - 0.0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600

21 / 33

JJ =

10.97 4.3243 -7.2098 -3.2731 0 0 -3.6049 -1.6366

-5.4951 10.6529 3.2731 -7.2098 0 0 1.6366 -3.6049

-7.59 -2.94 32.9307 12.9517 -5.0973 -1.93 -1.93

2.94 -7.59 -12.5839 33.3340 1.93 -5.0973 -5.0973

0 0 -4.8851 -2.0980 38.0451 15.8278 -12.5883

0 0 2.0980 -4.8851 -16.6919 37.6740 -29.3106

-3.94 -1.5920 -4.8658 -2.1227 -29.1950 -12.7363 16.0599

1.5920 -3.94 2.1227 -4.8658 12.7363 -29.1950 37.6252 k =

4 dU =

1.0e-007 *

0.0232 -0.28 -0.0013 -0.1352 0.00 -0.1876 -0.2038 PQ =

1.0e-006 *

-0.0240 -0.1486 -0.0424 -0.01 -0.0165 -0.0415 -0.0042 -0.0525

22 / 33

-5.0973 1.93 -29.3106 12.5883 37.7952 -16.8421 0.0087

U =

0.9961 - 0.1044i 1.03 - 0.0477i 1.0052 - 0.0845i 1.0032 - 0.0901i 1.0600 S =

-0.6000 - 0.1000i 0.2000 + 0.2000i -0.4500 - 0.1500i -0.4000 - 0.0500i 1.2982 + 0.2445i

2.2.3p—q法程序编写

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

%本程序的功能是用PQ分解法进行潮流计算 n=input(‘请输入节点数：n=’); nl==input(‘请输入支路数：nl==’);

isb==input(‘请输入平衡母线节点号：isb=’); pr==input(‘请输入误差精度:pr=’);

B1==input(‘请输入由支路参数形成的矩阵:B1=’); B2==input(‘请输入由节点参数形成的矩阵:B2=’);

X==input(‘请输入由节点号与其对地阻抗形成的矩阵:X=’); na==input(‘请输入PQ节点号:na=’);

Y==zeros(n);YI= zeros(n);e= zeros(1,n);f= zeros(1,n);V== zeros(1,n);O= zeros(1,n); for i=1:n if X(i,2)~=0; p=x(i,1);

Y(p,p)=1./X(i,2); end end

for i=1:nl

23 / 33

if B1(i,6)= =0

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

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

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

Y(p,p)= Y(p,p) +1./B1(i,3)+ B1(i,4)./2; YI(p,p)= YI(p,p) +1./B1(i,3); end %求导纳矩阵

G=real(Y);B=imag(YI);BI=imag(Y); for i=1:n

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

BI(i,i)= BI(i,i)+ B2(i,5); end

P=real(S);Q=imag(S); for i=1:n

e(i)=real(B2(i,3)); f(i)=imag(B2(i,3)); V(i)=B2(i,4); end

for i=1:n

if B2(i,6)= =2

V(i)=sqrt(e(i)~2+f(i)~2); O(i)=atan(f(i)./e(i)); end end

for i=2:n if i= =n

B(i,i)=1./B(i,i); else IC1=i+1; for j1=IC1:n

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

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

24 / 33

B(k,j1)=B(k,j1)-B(k,i)*B(i,j1); end end end end

p=0;q=0; for i=1:n

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

if B2(i1,6)= =2 k=k+1;

A(p,k)=BI(i,j1); end end end end

for i=1:na if i= =na

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

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

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

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

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

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

25 / 33

end

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

Np(K)=ICT2; if ICT2~=0 for i=2:n

DP(i)=B(i,i)*DP(i); if i~ =n IC1=i+1;

for k=ic1:n

DP(k)=DP(k)-B(k,i)*DP(i); end else

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

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

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

for i=2:n

O(i)=O(i)-DP(i); end

kq=1;L=0; for i=1:n

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

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

DQ1(i)=Q(i)-V(i)*C(i);

26 / 33

DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i)); if DET>=pr ICT3=ICT3+1; end end end end

else kp=0; if kq~=0; L=0;

for i=1:n

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

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

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

Nq(K)=ICT3; if ICT3~=0 L=0;

for i=1:na

DQ(i)=A(i,i)*DQ(i); if i= =na for LZ=2:i L=i+2-LZ; IC4=L-1; for MZ=1:IC4 I=IC4+1-MZ;

DQ(I)=DQ(I)-A(I,L)*DQ(L); end end else

IC1=i+1;

27 / 33

for k=IC1:na

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

for i=1:n

if B2(i,6)= =2 L=L+1;

V(i)=V(i)-DQ(L); end end kp=1; K=K+1; else kq=0; if kp~=0 K=K+1; end end

for i=1:n

Dy(K-1,i)=V(i); end end

disp(‘迭代次数‘)； disp（K）； disp（‘每次没有达到精度要求的有功功率个数为‘）； disp（NP）； disp（‘每次没有达到精度要求的无功功率个数为‘）； disp（Nq）; for k=1:n

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

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

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

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

28 / 33

for p=1:n C(p)=0; for q=1:n

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

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

disp(‘各节点的功率S（节点号从小到大排）：‘)为； disp（S）；

disp(‘各条支路的首端功率Si为（顺序同您输入B1时的一样）：‘）为； for i=1:nl

if B1(i,6)= =0 p=B1(i,1);q=B1(i,2);

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

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

disp(‘各条支路的末端功率Sj（顺序同您输入B1时的一样）：‘）为； for i=1:nl

if B1(i,6)= =0 p=B1(i,1);q=B1(i,2);

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

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

disp(Sj(q,p)); end

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

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

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

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

for i=1:K

29 / 33

Cs(i)=i; for j=1:n

Dy(K,j)=Dy(K-1,j); end end

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

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

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

（1）n为节点数、nl为支路数、isb为平衡母线节点号（固定为 1）、pr为误差精度。 （2）输入由支路参数形成的矩阵 B1 矩阵 B1的每行是由下列参数构成的： ①某支路的首端号 P； ②末端号 Q，且 P ⑥折算到哪一侧的标志（如果支路的首端 P 处于高压侧则请输入“1”，否则请输入“0”）。

Figure 4.1 schematic

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

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

⑤节点所接的无功补偿设备的容量； ⑥节点分类标号 igl。 输入数据：

请输入节点数：n=5 请输入支路数：nl=5

请输入平衡母线节点号：isb=1 请输入误差精度：pr=0.00001

请输入由支路参数形成的矩阵：B1=[1 2 0.03i 0 1.050;2 3 0.08+0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1] 请输入由节点参数形成的矩阵：B2=[0 0 1.05 1.05 0 1;0 3.7+1.3i 1.05 0 0 2;0 2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0 0 2;5 0 1.05 1.05 0 3]

30 / 33

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

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

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

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

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

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

各节点的电压U大小（节点号从小到大排）为： 1.0500 1.03 1.0779 0.8662 1.0500 各节点的电压相角O（节点号从小到大排）为： 0 -3.519 17.8535 -4.7785 21.8433 各节点的功率S（节点号从小到大排）为：

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

各条支路的首端功率Si（顺序同您输入B1时的一样）为： 2.5794+2.2994i -1.2774+0.2032i 0.1568+0.4713i 1.5845+0.6725i 5.0000+1.8131i

各条支路的末端功率Sj（顺序同您输入B1时的一样）为： -2.8794-1.9745i 1.41-0.2443i -0.1338-0.3909i -1.4662-0.4091i -5.0000-1.4282i

各条支路的功率损耗DS（顺序同您输入B1时的一样）为： 0.0000+0.3249i 0.1381-0.0412i 0.0230+0.0804i 0.1184+0.2635i -0.0000+0.3849i

31 / 33

Figure 4.3 Simulation diagram

经过手动计算，得1至5号节点功率分别为： 2.5794+2.2994i -3.7000-1.3000i -2.0000-1.0000i -1.6000-0.8000i 5.0000+1.8131i

举例验算结论：通过手动运算，发现运算结果与手动结果运算结果一样，故软件计算正确。

32 / 33

总结

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

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

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

参考文献与资料： 1、《电力系统稳态分析》，珩，中国电力，2007，第三版 2、《电力系统分析》，祯祥，大学，2005，第三版 3、《电力系统分析》，小林 阎晓霞 ，中国电力 4、《电力系统分析课程实际设计与综合实验》，祝书萍，中国电力，2007，第一版 5、清华大学等校合编.电子数字计算机的应用—电力系统计算[M].:水利电力，1978 6、朱仁峰.精通Matlab7[M].：清华大学，2006 7、《火车发电厂、变电所二次接线技术规定（强电部分）》SDGJ─78 8、《MALAB 6.0 科算完整解决法案》，精锐创组编，人民邮电大学 9、《MALAB 7.0 应用集锦》，林雪松 周婧 林德新 10、朱仁峰.精通Matlab7[M].：清华大学，2006