潮流计算的快速分解法可编辑修改word版文档格式.docx
《潮流计算的快速分解法可编辑修改word版文档格式.docx》由会员分享，可在线阅读，更多相关《潮流计算的快速分解法可编辑修改word版文档格式.docx（29页珍藏版）》请在冰豆网上搜索。
⎧∆V（k）=-B'
-1∆Q（（k）,V（k））V（k）
（2.4）
⎨（k+1）
⎩V
=V（k）
+
∆V
（k）
⎧∆（k）=-B'
-1∆P（（k）,V（k+1））V（k+1）
⎨（k+1）（k）
（2.5）
⎩=
+∆
快速分解法公式的特点是：
①P-和Q-V迭代分别交替进行；
②功率偏
差计算时使用最近修正过得电压值，且有功无功偏差都用电压幅值去除；
③B'
和B'
的构成不同，B'
应用-1x建立，并忽略所有接地支路（对非标准变比变压
器支路，变比可取为1），而B'
就是导纳矩阵的虚部，不包括PV节点。
在快速
分解法的实施中，这些技术细节缺一不可，否则程序的收敛性将受到影响。
1989年，荷兰学者VanAmerongen通过大量仿真计算发现了另一版本的快速分解潮流算法，他把该算法称为BX型算法，而把Stott的算法称为XB型算法，用以区分二者。
BX型算法与XB型算法的主要不同在于雅可比矩阵对角块的形成
上。
BX型算法的处理方式是：
在对系数矩阵BH进行简化时，保留了支路电阻的
影响，但忽略了接地支路项。
BX型算法的迭代格式与XB型算法是相同的。
计算经验表明，BX型和XB型两种快速分解潮流算法在大部分情况下性能接近，在某些情况下BX型算法收敛性略好。
快速分解法只对雅可比矩阵作了简化，但节点功率偏差量的计算及收敛条件仍是严格的，因此收敛后的潮流结果仍然是准确的。
由于方程的维数减小了，
且B'
和B'
是常数矩阵，只需在迭代计算之前形成一次，然后分解成因子表，并
一直在迭代过程中使用，所以计算效率大幅提高。
快速分解法是一种定雅可比法，虽然只具有线性收敛速度，但由于其鲁棒性好，适应性强，在电力工业界被广泛采用，特别适合在线计算。
2.2快速分解法的理论基础
Stott的快速分解法提出时并没有任何理论解释，它是计算实践的产物。
多年来，人们普遍认为在满足r<
<
x的系统中，快速分解法才能有较好的收敛性。
但在许多实际应用中，当r>
x时，快速分解法也能很好收敛。
因此，从理论上解释快速分解法的收敛机理，便成为一个有趣的研究课题。
20世纪80年代末，Monticelli等人的研究工作对这一问题做了比较完整的解释，在一定程度上阐明了XB型和BX型快速分解潮流算法的收敛机理。
Monticelli等人的分析工作是以定雅可比牛顿-拉夫逊迭代方程为出发点的。
具体过程如下：
①通过高斯消去法，把牛顿-拉夫逊法的每一次迭代等价地细分为三步计算；
②对每一步计算作详细分析，证明了在连续的两次牛顿-拉夫逊迭代中，上一次迭代的第三步和下一次迭代的第一步可以合并，从而导出等效的两步式分解算法；
③论证了该两步式分解算法的系数矩阵与快速分解法的系数
矩阵是一致的。
推导过程并未引用任何解耦的假设。
为以后书写方便，将式（2.1）中的∆PV用∆P代替，∆QV用∆Q代替，而V∆用∆代替，则给出的定雅可比法的修正公式改写如下
-⎡HN⎤⎡∆⎤=⎡∆P⎤
⎢ML⎥⎢∆V⎥⎢∆Q⎥
（2.6）
式中
⎣⎦⎣⎦⎣⎦
H=BH
≈∂∆P，N=-G
∂TN
≈∂∆P
∂VT
M=GM
≈∂∆Q，L=B
∂TL
≈∂∆Q
整个推导分为三步。
1）将原问题分解成P，Q子问题
首先，对式（2.6）用高斯消去法消去子块N，有
-⎡H-NL-1M0⎤⎡∆⎤=⎡∆P-NL-1∆Q⎤
⎢ML⎥⎢∆V⎥⎢∆Q⎥
（2.7）
记
~-1~-1
H=H-NLM，∆P=∆P-NL∆Q
并定义
LM
∆V=-L-1∆Q，∆V
=-L-1M∆
则式（2.7）的解可以表示为
H
⎧∆=-~-1∆P
⎨∆V=∆V+∆V
M
P
上式中对∆~的计算可以采用较简单的方法。
在给定的电压幅值和相角初值附近，保持电压相角不变，考虑只有电压幅值的变化∆VL时，有功功率的偏差量
为
∆P（,V+∆VL
）≈∆P（,V）+∂∆P∆V
∂VTL
=∆P（,V）-NL-1∆Q=∆~
（2.8）
综合上述结果，如果当前的迭代点为（（k）,V（k）），则第k次迭代对式（2.6）
的计算可以分解为以下三步。
①
⎧∆V（k）=-L-1∆Q（（k）,V（k））
⎨~L
V（k+1）=V（k）+∆V（k）
（2.9）
②
⎩L
⎧⎪∆（k）=-~-1（k）
~（k+1）
（2.10）
③
⎨⎪（k+1）=（k）+∆（k）
⎩
⎧∆V（k）=-L-1M∆（k）
⎨M~
V（k+1）=V（k+1）+∆V（k）
⎩M
（2.11）
2）简化无功迭代步骤
按①～③完成第k次迭代后，下面再考察第k+1次迭代的①，有
⎧∆V（k+1）=-L-1∆Q（（k+1）,V（k+1））
V（k+2）=V（k+1）+∆V（k+1）
（2.12）
利用式（2.11），上式中的无功功率偏差为
∆Q（（k+1）,V（k+1））=∆Q（（k+1）
~（k+1）+∆V（k））
VM
≈∆Q（（k+1）
~（k+1））+∂∆Q∆V（k）
V∂VTM
=∆Q（（k+1）
~（k+1））+L∆V（k）
（2.13）
代入式（2.12），经整理得
∆V（k+1）+∆V（k）=-L-1∆Q（（k+1）
~（k+1））
LM,V
（2.14）
V
式（2.14）说明，如果将第k次迭代的①计算出的~（k+1）和②计算出的（k+1），
用于计算第k+1次迭代的无功偏差量，即式（2.14）中的∆Q，则所求得的第k+1
次迭代的电压修正量将自动包含第k次迭代的③的式（2.11）与第k+1次迭代的①的式（2.12）合并，只需保留式（2.9）和式（2.10）。
因此，第k次迭代对式（2.6）的计算可以用以下两步计算完成：
⎧⎪∆V（k）=-L-1∆Q（（k）,V（k））
（2.15）
⎨
⎪⎩V
（k+1）
⎧⎪∆（k）=-~-1
∆P（（k）,V
（k+1））
（2.16）
在式（2.6）处已说明，∆P实际是∆PV，∆Q实际是∆QV，∆实际是
V∆，式（2.15）和式（2.16）和快速分解法迭代格式相同。
显然，这种迭代算法是否与快速分解法等效，取决于系数矩阵L和~。
与XB型快速分解法的修正
方程相比，系数矩阵L是导纳矩阵的虚部，这与B'
相同，所以关键要看~是否
与B'
有相同或相似的关系。
3）
简化有功迭代矩阵~
由
~
H的定义，有
~-1-1
H=H-NLM=BH+GNBLGM
（2.17）
对于一般的电网，~可能有较复杂的结构。
为了对~有直观的认识，假定
HH
网络中无PV接点，则式（2.17）中各矩阵的维数相等，并且接点导纳矩阵可用节点支路关联矩阵A和支路导纳对角矩阵（分别用的b和g表示电纳和电导）表
示。
下面将证明，对于树形电网或所有支路的r
x比值都相同的环形网络，~与
B'
相等。
如果网络是树状的，其关联矩阵A是方阵且非奇异，此时对式（2.17）有
H=AbAT
（AgAT
）（AbAT
）-1（AgAT）
=A（b+gATA-Tb-1A-1Ag）AT
=A（b+gb-1b）AT
=Ab'
AT
=B'
（2.18）
式中，b'
为以-1x为支路电纳组成的对角线矩阵；
B'
为以-1x为支路电纳建立的节点电纳矩阵。
这说明对树形电网，~就是XB型快速分解法中的B'
阵。
对于环形网络，如果电网是均一网，即对任一支路l有rlxl=，则得
gl=
rl
r2+x2
=xl
=-bl
llll
并有
-r2x1
b+gb1g=（1+2）b=（1+l）（-l）==-
llll
lx2r2+x2x
所以
g=-b，（1+2）b=b'
故有
+（AgAT
=AbAT+2（AbAT）（AbAT）-1（AbAT）
=（1+2）AbAT
（2.19）
网不是均一网，上述结论不再严格成立。
但~和B'
相比，在B'
的零
~的元素近似等于零；
在B'
的非零元素处，相应
~的元素近似和B'
如果电
元素处，相应HH
的非零元素相等。
这可以用下面的例子来说明。
图2.1四节点电力系统
以图2.1所示的四节点系统为例，图中给出了支路阻抗。
该例中H，~和
分别为
⎡
⎢
H=⎢
1.5
-1
1.2
-0.2
0.7
-0.5⎤
⎥
-0.5⎥
⎢-0.50-0.5⎥
⎡1.9
-0.8
-0.1-1⎤
⎡2-10-1⎤
~⎢-0.8
1.6
-0.80⎥
⎢-12-10⎥
⎢-0.1
1.9
-1⎥
=-⎢
⎢0-1
2-1⎥
⎢-1
0-1⎥
⎢-1
可见，B'
比H更接近于~，而用B'
代替~即得到XB型快速分解法。
3程序流程
计算时，只要有电路图及其各支路阻抗初始值、节点类型，即可进行迭代计算得出结果。
利用快速分解法编程流程图如下。
是
kp=1
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;
%支路数ye=0.00005;
%误差要求
point=[110-2-1;
21.0100.50.25;
31000];
%从exel中读取节点信息矩阵
branch=[120.010.20.02000;
130000.010.11.05;
230.020.2
0.04000];
%从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（:
J=branch（:
Rij=branch（:
Xij=branch（:
Zij=Rij+j*Xij;
B0=branch（:
RT=branch（:
6）;
XT=branch（:
7）;
KT=branch（:
8）;
%求节点导纳矩阵Y
form=1:
y%求Y中非对角元元素Yij
ifKT（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
x%求Y中的Ｙiiforn=1:
y
ifKT（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;
elseifI（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））;
elseifJ（n）==m%有变压器时，若支路终止节点为m,则变压器等值模型的对地支路的（KT-1）/KT算到J（m）节点
endY
Y（m,m）=Y（m,m）-Y（I（n）,J（n））+（KT（n）-1）/KT（n）*（1/ZT（n））;
elseY（m,m）=Y（m,m）;
end
G=real（Y）;
%求B'
矩阵及其逆矩阵B1
B=imag（Y）;
%求Y的虚部,节点电纳矩阵[y1]=[y-1];
[x1]=[x-1];
B11=zeros（x1,x1）;
form=1:
y1
B11（I（m）,J（m））=1/（Xij（m）+XT（m））;
B11（J（m）,I（m））=1/（Xij（m）+XT（m））;
x1
forn=1:
B11（m,m）=B11（m,m）-1/（Xij（n）+XT（n））;
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节点编号对应行置空
B22B2=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;
（x-1）%求deltaPisum1=0;
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）））;
deltaPi（m）=P（notph（m））-sum1;
Up=U;
%Up为中间变量Up（ph）=[];
%将平衡节点所在行置空
Unotph=Up;
%求得除平衡节点外的电压列向量
deltaa=（-B1*（deltaPi./Unotph））./Unotph;
%求相角a的不平衡量form=1:
（x-1）%求相角a的新迭代值矩阵
a（notph（m））=a（notph（m））+deltaa（m）;
max1=abs（deltaPi
（1）/U（notph
（1）））;
%求deltaP/U绝对值的最大值form=1:
（x-2）
ifabs（deltaPi（m）/U（notph（m）））<
abs（deltaPi（m+1）/U（notph（m+1）））
max1=abs（deltaPi（m+1）/U（notph（m+1）））;
ifmax1<
=e%如果最大值满足要求,则kp置为"
1"
表示收敛kp=1;
pqnum%求deltaQi
sum2=0;
forn=1:
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）））;
deltaQi（m）=Q（pq（m））-sum2;
Uq=U;
%Uq为中间变量
Uq（phpv）=[];
%将非PQ节点所在行置空Upq=Uq;
%求得包括PQ节点电压的电压列向量
deltaU=-B2*（deltaQi./Upq）;
%求U的不平衡量deltaUmax2=max（abs（deltaQi./Upq））;
%求deltaQ/U绝对值的最大值ifmax2<
=e%如果最大值满足要求,则kq置为"
表示收敛
kq=1;
pqnum%求U的迭代新值U
升级会员 