牛拉法潮流计算程序(附3机9节点结果对比).docx

摘要

电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各局部的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。潮流计算主要用于电网规划和静态平安分析，它可为扩建电力网络，以到达规划周期内所需要的输电能力提供依据；也可以对预想事故进行模拟和分析，校核预想事故下的电力系统平安性。

本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法，然后用具体的实例，利用MATLAB对牛顿-拉夫逊法的算法进行了验证。

关键词：电力系统潮流计算牛顿-拉夫逊法MATLAB

一、牛拉法的数学模型

对一个N节点的电力网路，列写节点电压方程，即

I=

式中，I为节点注入电流列相量，Y为节点导纳矩阵，V为节点电压列相量。

由于异地测量的两个电流缺少时间同步信息，以注入功率替换注入电流作为量。即

(1.2)

其中,Yij=

Pi

Qi

大局部情况下，PQ，求解Vθ。考虑到电网的功率平衡，至少选择一台发电机来平衡全网有功功率，即至少有一个平衡节点，常选择调频或出线较多的发电机作为平衡节点。具有无功补偿的母线能保持电压幅值恒定，这类节点可作为PV节点。

潮流计算中节点分类总结如下：

表SEQTable\*ARABIC1.1节点分类

PQ节点

负荷母线、变电站母线

绝大局部是PQ节点

PQ，求Vθ

PV节点

有无功储藏的发电机母线

很小局部是PV节点

PV〔V=1〕，求Qθ

平衡节点

选择一个发电机母线

Vθ〔V=1,θ=0〕,求PQ

电压相量作为参考轴θ=0

电力系统有m个PQ节点，r个PV节点和1个平衡节点，那么可以提取m+r个有功功率方程和m个无功功率方程，从而求解出m+r个θ和m个V，其余节点的有功和无功可通过式〔1.3〕、〔1.4〕求得，这样就完成了潮流计算。

二、潮流计算流程

用于潮流计算的有功和无功功率方程的阶数一般很高，常选择牛顿-拉夫逊迭代法进行求解高维的非线性方程组。

结束maxF(Xk)ε？XkXk+1=X

结束

maxF(

X

X

maxΔ

Δ

图2.1牛拉法迭代流程图

将功率方程改成Fx

ΔP

ΔQ

可以得到

Δθ

为了便于计算，将上式改为

Δθ

求取功率方程的雅克比矩阵

i≠j时，

?Δ

?Δ

?Δ

?Δ

i=j时，

?Δ

?Δ

?Δ

?Δ

牛拉法求解潮流的迭代过程如下

图2.2牛拉法潮流迭代流程图

三、MATLAB编程实例

实例采用王锡凡主编的《现代电力系统分析》p326-327的例题。

图3.1某小型电力网络图

例题中，1号发电机为平衡节点，2号和3号发电机为PV节点，其余为PQ节点。为了便于编程，可将PV节点和平衡节点放在节点描述图的最后，故要重新对节点进行编号。

。

图3.2修改节点编号后的电力网络图

支路表如下：

%首节点末节点电阻电抗容纳之半变比

12 0.010 0.0850.088 1

13 0.017 0.0920.079 1

24 0.032 0.1610.153 1

36 0.039 0.1700.179 1

45 0.0085 0.0720.0745 1

56 0.0119 0.10080.1045 1

91 00.057601

74 00.062501

86 00.058601

节点表如下：

%节点电压幅值电压相角节点有功节点无功节点类型/“1”为PQ节点，“2”为PV节点，“3”为平衡节点

11.00 00 1

21.00 -1.2500-0.5 1

31.00 -0.9-0.3 1

41.00 00 1

51.00 -1-0.351

61.00 00 1

71.0250 1.6300 0 2

81.0250 0.8500 0 2

91.0400 00 3

程序清单：取ε=10-5

a=load(zhilu.txt);

b=load(jiedian.txt);

N=size(a,1);%支路数

Nbus=size(b,1);%节点数

pq=0;

fork=1:Nbus%PQ节点个数

ifb(k,6)==1

pq=pq+1;

end

end

Y=zeros(Nbus);%节点导纳矩阵

fork=1:N

t1=a(k,1);t2=a(k,2);r=a(k,3);x=a(k,4);ban=a(k,5);K=a(k,6);

Y(t1,t1)=Y(t1,t1)+1/(r+j*x)+j*ban;

Y(t1,t2)=Y(t1,t2)-1/(K*(r+j*x));

Y(t2,t1)=Y(t2,t1)-1/(K*(r+j*x));

Y(t2,t2)=Y(t2,t2)+1/(K*K*(r+j*x))+j*ban;

end

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

precision=1