潮流计算的快速分解法可编辑修改word版.docx
- 文档编号:6933130
- 上传时间:2023-01-12
- 格式:DOCX
- 页数:29
- 大小:206.49KB
潮流计算的快速分解法可编辑修改word版.docx
《潮流计算的快速分解法可编辑修改word版.docx》由会员分享,可在线阅读,更多相关《潮流计算的快速分解法可编辑修改word版.docx(29页珍藏版)》请在冰豆网上搜索。
潮流计算的快速分解法可编辑修改word版
潮流计算的快速分解法
摘要:
本文采用快速分解法进行潮流计算,分析其基本理论,并使用MATLAB软件进行编程设计。
最后运用实例进行验证。
结果表明快速分解法具有较好的迭代速度。
关键词:
潮流计算快速分解法MATLAB编程,实例验证
1引言
潮流计算是电力系统分析最基本、最重要的计算,是电力系统运行、规划以及安全性、可靠性分析和优化的基础,也是各种电磁暂态和机电暂态分析的基础和出发点。
潮流计算要求具有可靠的收敛性,占用内存少,计算速度快,调整和修改容易,使用灵活方便。
各种算法的改进以及新算法的提出,很多都是为了使潮流计算能更好地满足计算要求。
本文应用快速分解法进行潮流计算,并给出算例分析。
2潮流计算的快速分解法
研究表明,用牛顿-拉夫逊法计算潮流时,每次迭代都要重新形成雅可比矩阵,然后重新对它进行因子表分解并求解修正方程。
为避免每次迭代重新形成雅可比矩阵及其因子表,人们研究用定雅可比矩阵取代随迭代过程不断变化的雅可比矩阵,这种方法叫定雅可比法。
此外,人们还结合电力系统的物理特点,发展了各种版本的解耦潮流算法,20世纪70年代初提出的快速分解法是这一阶段的主要研究成果。
关于快速分解潮流算法,有三项里程碑意义的研究成果。
其一是Stott在1974年发现的XB型算法;其二是VanAmerongen在1989年发现的BX型算法;其三是Monticelli等人在1990年所作的关于快速分解潮流算法收敛机理的理论阐述。
这些研究工作不仅是电力系统计算方面的典范,也揭示了这样一个事实:
工程上有效的方法一定有其深刻的理论来支持。
2.1快速分解法的修正方程及迭代格式
将极坐标型定雅可比法的修正公式重写如下:
-⎡BH-GN⎤⎡V∆⎤=⎡∆PV⎤
⎢GB
⎥⎢∆V⎥
⎢∆QV⎥
(2.1)
⎣ML⎦⎣⎦⎣⎦
经验表明,电力系统中有功功率主要受电压相角的影响,而无功功率主要受电压幅值的影响,同时由于高压电网大部分线路的电阻比电抗小,因此在牛顿-拉夫逊迭代中可以忽略雅可比矩阵的非对角块,即将GN,GM设为零,从而实
现有功和无功潮流修正方程的解耦。
Stott通过大量的计算实践发现,为了获得最好的收敛性,还要对雅可比矩阵的对角块作特殊的常数化处理:
对系数矩阵BH
,忽略支路电阻和接地支路的影响,即用-1x为支路电纳建立的节点电纳矩阵B'
HL
代替B;对系数矩阵B,用节点导纳矩阵中不包含PV节点的虚部B''代替;
V∆前的电压幅值用标幺值1代替。
于是可得简化的修正方程式如下:
-B'∆=∆PV
(2.2)
-B''∆V=∆QV
(2.3)
在潮流计算中,上述两个修正方程式依次交替迭代,Stott把在此基础发展起来的潮流算法称为快速分解法(fastdecoupledloadflow)。
假定当前点为((k),V(k)),则求解((k+1),V(k+1))的连续迭代格式如下:
⎧∆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)
(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< 但在许多实际应用中,当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%有变压器时,Yij为线路阻抗乘以KT后分之一再取负值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中的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; end 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)节点 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 潮流 计算 快速 解法 编辑 修改 word
![提示](https://static.bdocx.com/images/bang_tan.gif)