但在许多实际应用中，当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
∂VT
整个推导分为三步。
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
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））
⎨~L
V（k+2）=V（k+1）+∆V（k+1）
⎩L
（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）
VM
（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）
=V（k）
+
∆V
（k）
H
⎧⎪∆（k）=-~-1
∆P（（k）,V
（k+1））
⎩
（2.16）
⎨⎪（k+1）=（k）+∆（k）
在式（2.6）处已说明，∆P实际是∆PV，∆Q实际是∆QV，∆实际是
H
V∆，式（2.15）和式（2.16）和快速分解法迭代格式相同。
显然，这种迭代算法是否与快速分解法等效，取决于系数矩阵L和~。
与XB型快速分解法的修正
H
方程相比，系数矩阵L是导纳矩阵的虚部，这与B''相同，所以关键要看~是否
与B'有相同或相似的关系。
3）
H
简化有功迭代矩阵~
由
~
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）
H
式中，b'为以-1x为支路电纳组成的对角线矩阵；B'为以-1x为支路电纳建立的节点电纳矩阵。
这说明对树形电网，~就是XB型快速分解法中的B'阵。
对于环形网络，如果电网是均一网，即对任一支路l有rlxl=，则得
gl=
rl
r2+x2
=xl
r2+x2
=-bl
llll
并有
-r2x1
b+gb1g=（1+2）b=（1+l）（-l）==-
llll
lx2r2+x2x
llll
所以
g=-b，（1+2）b=b'
故有
~
H=AbAT
+（AgAT
）（AbAT
）-1（AgAT）
=AbAT+2（AbAT）（AbAT）-1（AbAT）
=（1+2）AbAT
=Ab'AT
=B'（2.19）
网不是均一网，上述结论不再严格成立。
但~和B'相比，在B'的零
H
~的元素近似等于零；在B'的非零元素处，相应
~的元素近似和B'
如果电
元素处，相应HH
的非零元素相等。
这可以用下面的例子来说明。
图2.1四节点电力系统
H
以图2.1所示的四节点系统为例，图中给出了支路阻抗。
该例中H，~和
B'分别为
⎡
⎢
H=⎢
⎢
1.5
-1
0
-1
1.2
-0.2
0
-0.2
0.7
-0.5⎤
0
⎥
⎥
-0.5⎥
⎢-0.50-0.5⎥
⎡1.9
-0.8
-0.1-1⎤
⎡2-10-1⎤
~⎢-0.8
1.6
-0.80⎥
⎢-12-10⎥
H=⎢
⎢-0.1
-0.8
1.9
⎥
-1⎥
B'=-⎢
⎢0-1
⎥
2-1⎥
⎢-1
0-1⎥
⎢-1
0-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（:
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
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
end
form=1:
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;
end
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）节点
end
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））;
end
form=1:
x1
forn=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节点编号对应行置空
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;
form=1:
（x-1）%求deltaPisum1=0;
forn=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的不平衡量form=1:
（x-1）%求相角a的新迭代值矩阵
a（notph（m））=a（notph（m））+deltaa（m）;
end
max1=abs（deltaPi
（1）/U（notph
（1）））;%求deltaP/U绝对值的最大值form=1:
（x-2）
ifabs（deltaPi（m）/U（notph（m）））max1=abs（deltaPi（m+1）/U（notph（m+1）））;
end
end
ifmax1<=e%如果最大值满足要求,则kp置为"1",表示收敛kp=1;
end
form=1:
pqnum%求deltaQi
sum2=0;forn=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的不平衡量deltaUmax2=max（abs（deltaQi./Upq））;%求deltaQ/U绝对值的最大值ifmax2<=e%如果最大值满足要求,则kq置为"1",表示收敛
kq=1;
end
form=1:
pqnum%求U的迭代新值U
升级会员 