潮流计算的快速分解法.doc

电力系统稳态分析——潮流计算的快速分解法
潮流计算的快速分解法

摘要：本文采用快速分解法进行潮流计算，分析其基本理论，并使用MATLAB软件进行编程设计。最后运用实例进行验证。结果表明快速分解法具有较好的迭代速度。
关键词：潮流计算 快速分解法 MATLAB编程，实例验证
1引言
潮流计算是电力系统分析最基本、最重要的计算，是电力系统运行、规划以及安全性、可靠性分析和优化的基础，也是各种电磁暂态和机电暂态分析的基础和出发点。潮流计算要求具有可靠的收敛性，占用内存少，计算速度快，调整和修改容易，使用灵活方便。各种算法的改进以及新算法的提出，很多都是为了使潮流计算能更好地满足计算要求。本文应用快速分解法进行潮流计算，并给出算例分析。
2潮流计算的快速分解法
研究表明，用牛顿-拉夫逊法计算潮流时，每次迭代都要重新形成雅可比矩阵，然后重新对它进行因子表分解并求解修正方程。为避免每次迭代重新形成雅可比矩阵及其因子表，人们研究用定雅可比矩阵取代随迭代过程不断变化的雅可比矩阵，这种方法叫定雅可比法。此外，人们还结合电力系统的物理特点，发展了各种版本的解耦潮流算法，20世纪70年代初提出的快速分解法是这一阶段的主要研究成果。
关于快速分解潮流算法，有三项里程碑意义的研究成果。其一是Stott在1974年发现的XB型算法；其二是Van Amerongen在1989年发现的BX型算法；其三是Monticelli等人在1990年所作的关于快速分解潮流算法收敛机理的理论阐述。这些研究工作不仅是电力系统计算方面的典范，也揭示了这样一个事实：工程上有效的方法一定有其深刻的理论来支持。
2.1 快速分解法的修正方程及迭代格式
将极坐标型定雅可比法的修正公式重写如下：
（2.1）
经验表明，电力系统中有功功率主要受电压相角的影响，而无功功率主要受电压幅值的影响，同时由于高压电网大部分线路的电阻比电抗小，因此在牛顿-拉夫逊迭代中可以忽略雅可比矩阵的非对角块，即将，设为零，从而实现有功和无功潮流修正方程的解耦。Stott通过大量的计算实践发现，为了获得最好的收敛性，还要对雅可比矩阵的对角块作特殊的常数化处理：对系数矩阵，忽略支路电阻和接地支路的影响，即用为支路电纳建立的节点电纳矩阵代替；对系数矩阵，用节点导纳矩阵中不包含节点的虚部代替；前的电压幅值用标幺值1代替。于是可得简化的修正方程式如下：
（2.2）
（2.3）
在潮流计算中，上述两个修正方程式依次交替迭代，Stott把在此基础发展起来的潮流算法称为快速分解法（fast decoupled load flow）。假定当前点为，则求解的连续迭代格式如下：
（2.4）
（2.5）
快速分解法公式的特点是：①和迭代分别交替进行；②功率偏差计算时使用最近修正过得电压值，且有功无功偏差都用电压幅值去除；③和的构成不同，应用建立，并忽略所有接地支路（对非标准变比变压器支路，变比可取为1），而就是导纳矩阵的虚部，不包括节点。在快速分解法的实施中，这些技术细节缺一不可，否则程序的收敛性将受到影响。
1989年，荷兰学者Van Amerongen通过大量仿真计算发现了另一版本的快速分解潮流算法，他把该算法称为BX型算法，而把Stott的算法称为XB型算法，用以区分二者。BX型算法与XB型算法的主要不同在于雅可比矩阵对角块的形成上。BX型算法的处理方式是：在对系数矩阵进行简化时，保留了支路电阻的影响，但忽略了接地支路项。BX型算法的迭代格式与XB型算法是相同的。计算经验表明，BX型和XB型两种快速分解潮流算法在大部分情况下性能接近，在某些情况下BX型算法收敛性略好。
快速分解法只对雅可比矩阵作了简化，但节点功率偏差量的计算及收敛条件仍是严格的，因此收敛后的潮流结果仍然是准确的。由于方程的维数减小了，且和是常数矩阵，只需在迭代计算之前形成一次，然后分解成因子表，并一直在迭代过程中使用，所以计算效率大幅提高。快速分解法是一种定雅可比法，虽然只具有线性收敛速度，但由于其鲁棒性好，适应性强，在电力工业界被广泛采用，特别适合在线计算。
2.2 快速分解法的理论基础
Stott的快速分解法提出时并没有任何理论解释，它是计算实践的产物。多年来，人们普遍认为在满足的系统中，快速分解法才能有较好的收敛性。但在许多实际应用中，当时，快速分解法也能很好收敛。因此，从理论上解释快速分解法的收敛机理，便成为一个有趣的研究课题。20世纪80年代末，Monticelli等人的研究工作对这一问题做了比较完整的解释，在一定程度上阐明了XB型和BX型快速分解潮流算法的收敛机理。
Monticelli等人的分析工作是以定雅可比牛顿-拉夫逊迭代方程为出发点的。具体过程如下：①通过高斯消去法，把牛顿-拉夫逊法的每一次迭代等价地细分为三步计算；②对每一步计算作详细分析，证明了在连续的两次牛顿-拉夫逊迭代中，上一次迭代的第三步和下一次迭代的第一步可以合并，从而导出等效的两步式分解算法；③论证了该两步式分解算法的系数矩阵与快速分解法的系数矩阵是一致的。推导过程并未引用任何解耦的假设。
为以后书写方便，将式（2.1）中的用代替，用代替，而用代替，则给出的定雅可比法的修正公式改写如下
（2.6）
式中
，
，
整个推导分为三步。
1) 将原问题分解成，子问题
首先，对式（2.6）用高斯消去法消去子块，有
（2.7）
记
，
并定义
，
则式（2.7）的解可以表示为

上式中对的计算可以采用较简单的方法。在给定的电压幅值和相角初值附近，保持电压相角不变，考虑只有电压幅值的变化时，有功功率的偏差量为
（2.8）
综合上述结果，如果当前的迭代点为，则第次迭代对式（2.6）的计算可以分解为以下三步。
①
（2.9）
②
（2.10）
③
（2.11）
2) 简化无功迭代步骤
按①～③完成第次迭代后，下面再考察第次迭代的①，有
（2.12）
利用式（2.11），上式中的无功功率偏差为

（2.13）
代入式（2.12），经整理得
（2.14）
式（2.14）说明，如果将第次迭代的①计算出的和②计算出的，用于计算第次迭代的无功偏差量，即式（2.14）中的，则所求得的第次迭代的电压修正量将自动包含第次迭代的③的式（2.11）与第次迭代的①的式（2.12）合并，只需保留式（2.9）和式（2.10）。因此，第次迭代对式（2.6）的计算可以用以下两步计算完成：
（2.15）
（2.16）
在式（2.6）处已说明，实际是，实际是，实际是，式（2.15）和式（2.16）和快速分解法迭代格式相同。显然，这种迭代算法是否与快速分解法等效，取决于系数矩阵和。与XB型快速分解法的修正方程相比，系数矩阵是导纳矩阵的虚部，这与相同，所以关键要看是否与有相同或相似的关系。
3) 简化有功迭代矩阵
由的定义，有
（2.17）
对于一般的电网，可能有较复杂的结构。为了对有直观的认识，假定网络中无接点，则式（2.17）中各矩阵的维数相等，并且接点导纳矩阵可用节点支路关联矩阵和支路导纳对角矩阵（分别用的和表示电纳和电导）表示。下面将证明，对于树形电网或所有支路的比值都相同的环形网络，与相等。
如果网络是树状的，其关联矩阵是方阵且非奇异，此时对式（2.17）有

（2.18）
式中，为以为支路电纳组成的对角线矩阵；为以为支路电纳建立的节点电纳矩阵。这说明对树形电网，就是XB型快速分解法中的阵。
对于环形网络，如果电网是均一网，即对任一支路有，则得

并有

所以
，
故有

（2.19）
如果电网不是均一网，上述结论不再严格成立。但和相比，在的零元素处，相应的元素近似等于零；在的非零元素处，相应的元素近似和的非零元素相等。这可以用下面的例子来说明。

图2.1 四节点电力系统
以图2.1所示的四节点系统为例，图中给出了支路阻抗。该例中，和分别为

可见，比更接近于，而用代替即得到XB型快速分解法。
3程序流程
计算时，只要有电路图及其各支路阻抗初始值、节点类型，即可进行迭代计算得出结果。利用快速分解法编程流程图如下。

程序开始
输入节点和支路数据数据
建立节点导纳矩阵
形成矩阵和
置迭代次数

计算无功不平衡量

计算相角修正量，求得
计算有功不平衡量

计算电压修正量，求得
计算平衡节点功率

？
？
输出数据
是
否
是
是
否
否

4程序
function [Y,U,P,Q,deltaSij,a,S,Sij,Sji,sumdeltaS]=PQ()
%电力系统潮流计算程序；
%输出：U——节点电压，P--节点有功，Q--节点无功，deltaSij--支路功率损耗，
%Sij--从节点i流向节点j的功率,S--节点复功率,sumdeltaS--网络总损耗
%输入参数：point为节点信息矩阵，branch为支路信息矩阵；
[x]=3;%节点数x
[y]=3;%支路数y
e=0.00005;%误差要求
point= [1 1 0 -2 -1;2 1.01 0 0.5 0.25;3 1 0 0 0];%从exel中读取节点信息矩阵
branch=[1 2 0.01 0.2 0.02 0 0 0;1 3 0 0 0 0.01 0.1 1.05;2 3 0.02 0.2 0.04 0 0 0];%从exel中读取支路信息矩阵
TYPE=zeros(x,1);%TYPE为节点类型矩阵
U=zeros(x,1);%U为节点电压矩阵
a=zeros(x,1);%a为节点电压相角矩阵
P=zeros(x,1);%P为节点有功功率
Q=zeros(x,1);%Q为节点无功功率
I=zeros(y,1);%I为起始节点编号矩阵
J=zeros(y,1);%J为终止节点编号矩阵
Rij=zeros(y,1);%R为线路电阻
Xij=zeros(y,1);%X为线路电抗
Zij=Rij+j*Xij;%Yij为线路阻抗
Y=zeros(x);%Y为n阶节点导纳方阵
G=zeros(x);%G为n阶节点电导方阵
B=zeros(x);%B为n阶节点电纳方阵
B0=zeros(y,1);%B0为n*1阶线路对地电纳值
RT=zeros(y,1);%RT为ij支路y( 矩阵branch的行数)*1阶变压器电阻
XT=zeros(y,1);%XT为ij支路y*1阶变压器电抗
ZT=RT+j*XT;%求变压器阻抗
KT=zeros(y,1); %K为ij支路y*1阶变压器变比,若k=0表示无变压器,K=1则为标准变比,k不等于1为非标准变比
%------------------------------矩阵赋初值：
TYPE=point(:,1);%将point矩阵的第一列赋给TYPE,以下类似
U=point(:,2);
a=point(:,3);
P=point(:,4);
Q=point(:,5);
I=branch(:,1);
J=branch(:,2);
Rij=branch(:,3);
Xij=branch(:,4);
Zij=Rij+j*Xij;
B0=branch(:,5);
RT=branch(:,6);
XT=branch(:,7);
ZT=RT+j*XT;
KT=branch(:,8);
%------------------------------求节点导纳矩阵Y
for m=1:y %求Y中非对角元元素Yij
if KT(m)==0%若无变压器,则Yij直接为线路阻抗分之一取负值．
Y(I(m),J(m))=-1/Zij(m);
Y(J(m),I(m))=-1/Zij(m);
else %有变压器时，Ｙij为线路阻抗乘以ＫＴ后分之一再取负值
Y(I(m),J(m))=-1/(KT(m)*ZT(m));
Y(J(m),I(m))=-1/(KT(m)*ZT(m));
end
end
for m=1:x %求Y中的Ｙii
for n=1:y
if KT(n)==0 %无变压器时Ｙii为Ｙij加上线路对地电导的一半乘j
if(I(n)==m|J(n)==m)
Y(m,m)=Y(m,m)-Y(I(n),J(n))+j*B0(n)/2;
end
elseif I(n)==m %有变压器时，若支路起始节点为m,则变压器等值模型的对地支路的(1-KT)/KT^2算到Ｉ(m)节点
Y(m,m)=Y(m,m)-Y(I(n),J(n))+(1-KT(n))/(KT(n)^2)*(1/ZT(n));
elseif J(n)==m %有变压器时，若支路终止节点为m,则变压器等值模型的对地支路的(KT-1)/KT算到J(m)节点
Y(m,m)=Y(m,m)-Y(I(n),J(n))+(KT(n)-1)/KT(n)*(1/ZT(n));
else Y(m,m)=Y(m,m);
end
end
end
Y
G=real(Y);
%-----------------------求B'矩阵及其逆矩阵B1
B=imag(Y);%求Y的虚部,节点电纳矩阵
[y1]=[y-1];
[x1]=[x-1];
B11=zeros(x1,x1);
for m=1:y1
B11(I(m),J(m))=1/(Xij(m)+XT(m));
B11(J(m),I(m))=1/(Xij(m)+XT(m));
end
for m=1:x1
for n=1:y
if(I(n)==m|J(n)==m)
B11(m,m)=B11(m,m)-1/(Xij(n)+XT(n));
end
end
end
ph=find(TYPE(:,1)==3);%找出平衡节点编号
B11(:,ph)=[];%平衡节点编号对应行置空
B11(ph,:)=[];%平衡节点编号对应列置空
B11
B1=B11;
B1=inv(B1);%B1求逆后得B1矩阵
%-----------------------%求B''及其逆矩阵B2
phpv=find(TYPE(:,1)>1);%找出非PQ节点的编号
B22=B; %BB矩阵为中间变量
B22(:,phpv)=[];%非PQ节点编号对应行置空
B22(phpv,:)=[];%非PQ节点编号对应行置空
B22
B2=B22;
B2=inv(B2);%求得B2矩阵
%-------------计算各节点有功功率不平衡量deltaPi
k=0; %k为迭代次数
kp=0; %计算P不平衡量deltaPi的收敛标志(0表示不收敛,1表示收敛)
kq=0; %计算U不平衡量deltaQi的收敛标志(0表示不收敛,1表示收敛)
notph=find(TYPE(:,1)<3);%找出非平衡节点编号
deltaPi=zeros(x-1,1);%deltaPi为x*1阶矩阵,x即为节点数
pq=find(TYPE(:,1)==1);%找出PQ节点编号
pqnum=size(B2);
pqnum=pqnum(1);%求PQ节点的个数(因B1矩阵的维数等于PQ节点数)
deltaQi=zeros(pqnum,1);%deltaQi为pqnum*1阶矩阵
while((~kq|~kp)&(k<100))
k=k+1;
for m=1:(x-1)%求deltaPi
sum1=0;
for n=1:x
sum1=sum1+U(notph(m))*U(n)*(G(notph(m),n)*cos(a(notph(m))-a(n))+B(notph(m),n)*sin(a(notph(m))-a(n)));
end
deltaPi(m)=P(notph(m))-sum1;
end
Up=U; %Up为中间变量
Up(ph)=[];%将平衡节点所在行置空
Unotph=Up;%求得除平衡节点外的电压列向量
deltaa=(-B1*(deltaPi./Unotph))./Unotph;%求相角a的不平衡量
for m=1:(x-1) %求相角a的新迭代值矩阵
a(notph(m))=a(notph(m))+deltaa(m);
end
max1=abs(deltaPi(1)/U(notph(1)));%求deltaP/U绝对值的最大值
for m=1:(x-2)
if abs(deltaPi(m)/U(notph(m)))<abs(deltaPi(m+1)/U(notph(m+1)))
max1=abs(deltaPi(m+1)/U(notph(m+1)));
end
end
if max1<=e %如果最大值满足要求,则kp置为"1",表示收敛
kp=1;
end
for m=1:pqnum %求deltaQi
sum2=0;
for n=1:x
sum2=sum2+U(pq(m))*U(n)*(G(pq(m),n)*sin(a(pq(m))-a(n))-B(pq(m),n)*cos(a(pq(m))-a(n)));
end
deltaQi(m)=Q(pq(m))-sum2;
end
Uq=U;%Uq为中间变量
Uq(phpv)=[];%将非PQ节点所在行置空
Upq=Uq;%求得包括PQ节点电压的电压列向量
deltaU=-B2*(deltaQi./Upq);%求U的不平衡量deltaU
max2=max(abs(deltaQi./Upq)); %求deltaQ/U绝对值的最大值
if max2<=e %如果最大值满足要求,则kq置为"1",表示收敛
kq=1;
end
for m=1:pqnum %求U的迭代新值
U(pq(m))=U(pq(m))+deltaU(m);
end

end
sum3=0+j*0;%求平衡节点功率Sph
for m=1:x
sum3=sum3+conj(Y(ph,m))*(U(m)*cos(a(m))-i*U(m)*sin(a(m)));
end
Sph=(U(ph)*cos(a(ph))+j*U(ph)*sin(a(ph)))*sum3;%求平衡节点功率
fid=fopen('PQj.txt','wt'); %打开输出文件,PQ.txt,向其写入结果数据
fprintf(fid,'潮流计算输出结果\n');
fprintf(fid,'迭代次数k为: %d \n',k);
fprintf(fid,'节点电压: \n');
for m=1:x
fprintf(fid,'第%d个节点: %f\n',m,U(m));
end
fprintf(fid,'节点相角: \n');
for m=1:x
fprintf(fid,'第%d个节点: %f\n',m,a(m));
end
fprintf(fid,'节点有功功率: \n');
for m=1:x
fprintf(fid,'第%d个节点: %f\n',m,P(m));
end
fprintf(fid,'节点无功功率: \n');
for m=1:x
fprintf(fid,'第%d个节点: %f\n',m,Q(m));
end
fclose(fid); %关闭文件,程序结束
5算例验证
5.1 求解系统的节点导纳矩阵
如图5.1（a）所示的一个三母线电力系统，在母线①和母线③之间的输电线的母线③端连接有一个纵向串联加压器，可在同一电压等级改变电压幅值。该系统的网络原件用图（b）所示的等值电路表示，串联支路用电阻和电抗表示，并联支路用电纳表示。支路（1,3）用一个变比可调的等值变压器支路表示，非标准变比，在节点①侧。求解该网络的节点导纳矩阵。

（a）

（b）

（c）
图5.1 三母线电力系统
解 首先对支路编号并规定串联支路的正方向如图5.1（c）所示，则广义节点支路关联矩阵是

A中行与节点对应，列与支路对应。
支路导纳是

建立节点导纳矩阵如下：
于是，Y的各元素是

导纳矩阵是对称矩阵，只需求出上三角部分。最后求得导纳矩阵为

4.2 用快速分解法计算潮流
如图5.2所示的电力系统，其导纳矩阵在上面已经求出。如果给定各节点的发电和负荷功率以及节点电压，写出极坐标形式的潮流方程，并用快速分解法计算潮流。已知，，，，。

图5.2 三母线电力系统
解 如图5.2所示节点①是节点，节点②是节点，节点③是节点。待求的状态变量是，共有两个有功潮流方程和一个无功潮流方程：

其具体表达式如下：

这个系统的节点导纳矩阵在5.1中已求出，为
快速分解法中的和分别为
,
其中是由建立所得，即通过
求出。
而则直接从导纳矩阵虚部中取。将矩阵的具体数值代入功率偏差方程有

从平启动开始，所有节点电压为，这时电压初值为

代入有功潮流方程中计算第一次迭代时的有功不平衡量

所以

计算相角修正量

式中，前得电压幅值项取值为1。
然后进行无功迭代。
先计算节点无功不平衡量。这时节点电压相角用最新修正后的值
，
只有一个节点，故只有一个无功平衡方程

再计算电压修正值：

再利用新计算出的电压转入第二次有功迭代。
整个迭代过程中功率不平衡量、电压幅值和角度的修正量以及电压幅值和角度的变化量列在表5.1中，其中、、的列用于判断收敛。
表5.1 快速分解法潮流计算迭代过程（）
迭代次数

1
-1.95260
0.49248
0.13670
0.01959
-0.13670
-0.01959
-0.77151
0.05527
0.94473
2
-0.13518
0.01872
0,.01071
0.00350
-0.14741
-0.02309
-0.08405
0.00637
0.93835
3
-0.01486
0.00431
0.00110
0.00012
-0.14851
-0.02321
-0.00904
0.00076
0.93759
4
-0.00167
0.00060
0.00012
0.00000
-0.014862
-0.02321
-0.00176
0.00009
0.93750
5
-0.00019
0.00007
0.00001
0.00000
-0.14854
-0.02321
-0.00014
0.00001
0.93749
6
-0.00002
0.00001

-0.00001

由表5.1可知第六次迭代时用功功率和无功功率不平衡量已小于预先制定的收敛阙值，因此不用继续做电压相角和幅值的修正。
潮流计算程序结果截图如图5.3所示：

（a）

（b）
图5.3 P-Q潮流计算程序输出结果截图

6参考文献
[1] 张伯明，陈寿孙，严正.高等电力网络分析[M].北京：清华大学出版社，2007.
[2] 韩平.潮流计算方法的实用化改进[D].北京：华北电力大学,2009.
[3] 刘国栋.电力网络潮流计算的模块开发[D].吉林：吉林大学,2009.
[4] 王忠礼，段慧达，高玉峰.MATLAB应用技术——在电气工程与自动化专业中的应用[M].北京：清华大学出版社，2007.
[5] 龚纯，王正林.MATLAB语言常用算法程序集[M].北京：电子工业出版社，2008.
- 22 -