潮流计算的快速分解法可编辑修改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%有变压器时,Yij为线路阻抗乘以KT后分之一再取负值Y(I(m),J(m))=-1/(KT(m)*ZT(m));
Y(J(m),I(m))=-1/(KT(m)*ZT(m));
end
x%求Y中的Yiiforn=1:
y
ifKT(n)==0%无变压器时Yii为Yij加上线路对地电导的一半乘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算到I(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