潮流计算的快速分解法程序

研究生课程设计 (论文)电力系统稳态分析Stability Analysis of Power Systems教学单位自动化学院姓名薛媛媛学号 *********年级2011 级专业电力系统及其自动化指导教师杨伟职称副教授2011年12月6日摘要：快速分解法是一种定雅克比法，形成系数矩阵'B 、''B 时忽略了支路电阻、对地导纳和理想变压器非标准变化，以及θP -迭代过程中节点电压的不同取值，从而避免每次迭代重新形成雅克比矩阵及其因子表，计算效率大幅提高。

本文采用快速分解法中的XB 型算法，并基于MATLAB 软件仿真分析具体实例，发现具有较好的收敛性，特别适合在线计算。

关键词：电力系统 潮流计算 快速分解法 雅可比矩阵0 引言用牛顿-拉夫逊法计算潮流时，每次迭代都需要重新形成雅克比矩阵并分解因子表。

为避免每次迭代重新形成雅克比矩阵及其因子表，20世纪70年代初Stott 提出了快速分解法。

快速分解法采用了一些假设：a 、电力系统有功功率主要受电压相角影响，无功功率主要受电压幅值影响。

b 、高压网线路的r<<x 。

c 、系数矩阵'B 、''B 忽略了支路电阻、对地导纳和理想变压器非标准变化。

d 、θP -迭代过程中的V Δθ前的电压幅值用标幺值1代替。

对极坐标型定雅克比法的修正公式HN M L ⎡⎤-⎡⎤⎡⎤⎢⎥-=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎣⎦ΔP B G V ΔθV ΔQ G B ΔV V （1） 处理得到简化的修方程式'''==--ΔPB ΔθVΔQ B ΔV V（2）在潮流计算中，（2）式中两个修正方程式交替迭代，Stott 把在此基础发展起来的潮流算法称为快速分解法。

快速分解法是目前电力系统进行潮流计算的主要方法。

由于该方法将节点电压的相位和幅值的迭代过程分开进行，使修正方程的系数矩阵维数降低很多，系数矩阵为常数且对称，因而快速分解法具有计算速度快,占用内存少等特点。

1 理论基础以定雅克比牛顿-拉夫逊迭代方程为出发点，具体过程如下：通过高斯消去法，把牛顿-拉夫逊法的每一次迭代等价地细分为三步计算；对每一步运算作详细分析，证明了在连续的两次牛顿-拉夫逊迭代中，上一次迭代的第三步和下一次迭代的第一步可以合并，从而导出等效的两步式分解算法；论证了该两步式分解算法的系数矩阵与快速分解算法的系数矩阵是一致的。

推导过程并未因用任何解耦的假设。

将极坐标型定雅克比法的修正公式中ΔP V和ΔQ V 用ΔP 和ΔQ 代替，V Δθ用Δθ代替。

⎡⎤⎡⎤⎡⎤-=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦H N ΔθΔP M L ΔV ΔQ （3） 其中,,,H N M L T T T T∂∆∂∆∂∆∂∆=≈=-≈=≈=≈∂∂∂∂P P Q QH B N G M G L B θV θV1.1将原问题分解为P ,Q 子问题用高斯消去法消去子块N11--⎡⎤⎡⎤--⎡⎤-=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ΔθH NL M 0ΔP NL ΔQ ΔV M L ΔQ（4） 记11,--=-=-H H NL M ΔP ΔP NL ΔQ因-M ΔθL ΔV =ΔQ -，得到11--=--ΔV L ΔQ L M Δθ 于是定义11,L M --=-=-ΔV L ΔQ ΔV L M Δθ解为 ⎧-⎨⎩-1L MΔθ=H ΔPΔV =ΔV +ΔV （5）在给定的电压幅值和相角初值附近，保持电压相角不变，考虑只有电压幅值的变化L ΔV 时，有功功率的偏差量为()()()1+T-∂∆∆≈∆=∆-=∆∂L L P P θV +ΔV P θV ΔV P θV NL ΔQ P V，，， （6） 综合上述结果，如果当前的迭代点为()()()k kθV ，，则第k 次迭代可以分为三步()()()()()()()11=k k k L k k k L -+⎧∆-∆⎪⎨⎪=+∆⎩V L Q θV V V V ， （7）()()()()()()()111,k k k k k k+-+⎧∆=-∆⎪⎨⎪=+∆⎩θH P θV θθθ （8）()()()()()111k k M k k k M -++⎧∆=-∆⎪⎨=+∆⎪⎩V L M θVV V （9） 1.2简化无功迭代步骤按式（7）、式（8）和式（9）完成第k 次迭代后，下面考虑第k +1次迭代，有()()()()()()()+1+1+112+1+1=k k k L k k k L -+⎧∆-∆⎪⎨⎪=+∆⎩V L Q θV V V V ， （10）利用式（9），上式中的无功功率偏差为()()()()()()()()()()()()()()()+1+1+1+1+1+1+1+1=++=+k k k k k Mk k k M Tk k k M∆∆∆∂∆≈∆∆∂∆∆Q θV Q θV V Q Q θV V VQ θV L V ，，，， （11） 代入式（10）中，经整理得()()()()()+1+1+11+=k k k k LM -∆∆-∆V V L Q θV ， （12）式（12）说明，如果将第k 次迭代计算出的()1k +V 和()1k +θ用于计算第1k +迭代的无功功率偏差量，则所求得的第1k +次迭代的电压修正量将自动包含第k 次迭代计算出的()k M ∆V 。

所以，()k M ∆V 的计算可以省略，相当于将式（9）与式（10）合并，因此，第k 次迭代可以两步完成()()()()()()()11=k k k k k k -+⎧∆-∆⎪⎨⎪=+∆⎩V L Q θV V V V ， （13）()()()()()()()111,k k k k k k +-+⎧∆=-∆⎪⎨⎪=+∆⎩θH P θV θθθ （14）整理式（13）和式（14）得到简化的修正方程式=-∆=∆-∆∆H θP L V Q（15）因ΔP 和ΔQ 即为极坐标型定雅克比法的修正公式中ΔPV和ΔQ V ，所以式（15）与式（2）格式相同。

与XB 型快速分解法的修正公式相比，系数矩阵L 是导纳矩阵的虚部，与''B 相同。

由H 的定义，有11H N L M ---H =H NL M =B +G B G （16）对于一般的电网，H 可能有较复杂的结构。

为了对H 有直观的认识，假定网络中无PV 接点，则式（16）中各矩阵的维数相等，并且接点导纳矩阵可用节点支路关联矩阵A 和支路导纳对角矩阵（分别用的b 和g 表示电纳和电导）表示。

下面将证明，对于树形电网或所有支路的x r 比值都相同的环形网络，H 与'B 相等。

如果网络是树状的，其关联矩阵A 是方阵且非奇异，此时对式（16）有1-T T T T H =AbA +(AgA )(AbA )(AgA ) --T -T 11T =A(b +gA A b A Ag)A -1T =A(b +gb b)A'T '=Ab A =B （17）式中，'b 为以x 1-为支路电纳组成的对角线矩阵；'B 为以x 1-为支路电纳建立的节点电纳矩阵。

这说明对树形电网，H 就是XB 型快速分解法中的'B 阵。

对于环形网络，如果电网是均一网，即对任一支路l 有α=l l x r ，则得l ll ll l l l b x r x x r r g αα-=+=+=2222 并有l ll l l l l l ll l x x r x x r b g b g b 1))(1()1(222221-==+-+=+=+-α所以g =-αb ，2(1)α+='b b故有-T T T 1T H =AbA +(AgA )(AbA )(AgA ) 21α-=+T T T T AbA (AbA )(AbA )(AbA )2(1)α=+T AbA'T '=Ab A =B（18） 如果电网不是均一网，上述结论不再严格成立。

但H 和'B 相比，在'B 的零元素处，相应H 的元素近似等于零；在'B 的非零元素处，相应H 的元素近似和'B 的非零元素相等。

这可以用下面的例子来说明。

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

该例中H ，H 和'B 分别为1.5100.51 1.20.2000.20.70.50.500.51--⎡⎤⎢⎥--⎢⎥=⎢⎥--⎢⎥--⎣⎦H 1.90.80.110.8 1.60.800.10.8 1.911012---⎡⎤⎢⎥--⎢⎥=⎢⎥---⎢⎥--⎣⎦H 2101121001211012--⎡⎤⎢⎥--⎢⎥=-⎢⎥--⎢⎥--⎣⎦'B 可见，'B 比H 更接近于H ，而用'B 代替H 即得到XB 型快速分解法。

以上推导过程中，只在对有功功率偏差量和无功功率偏差量的计算处做了线性化近似，既没有r<<x 的假设，也没有引用PQ 解耦的假设。

2 快速分解法在地震系统安全分析中的应用某地区电力供应系统由2个发电厂，4个变电站组成，系统模型见图2，其中节点1、2、3为PQ节点，4为PV节点，5为平衡节点。

以100MVA、220KV为基准的节点初始标幺值和网络参数见表3、4。

利用快速分解法，当迭代功率偏差标准为 =0.0001时，迭代结束。

图2 电力系统模型节点 1 2 3 4 5P -1.60 -2.00 -3.70 5.00 --Q -0.80 -1.00 -1.30 1.05 1.05表1 节点负荷值支路r x b 变比1-2 0.040 0.250 0.250 -1-3 0.10 0.350 0．000 -2-3 0.080 0.300 0.300 -2-4 0.000 0.015 - 1.0503-5 0.000 0.030 - 1.050表2 线路参数3．程序流程图4．程序本文采用两种程序来仿真地震电力系统的潮流，两个程序得到的仿真结果是相同的。

1.程序一function [Y,U,P,Q,dSij, O,S,Sij,Sji,sumdS]=FDLFpowerflow()%输出：U——节点电压，P--节点有功，Q--节点无功，dSij--支路功率损耗，O 为相角；Sij--从节点i流向节点j的功率,S--节点复功率,sumdS--网络总损耗；point 为节点矩阵，branch为支路矩阵；[x]=5;%节点数x[y]=5;%支路数ye=0.00001;%误差要求point=[1 1 0 -1.6 -0.8;1 1 0 -2.0 -1.0;1 1 0 -3.7 -1.3;3 1.05 0 5.0 1.05;2 1.05 0 0 1.05];%节点矩阵branch=[1 2 0.04 0.250 0.250 0;1 3 0.10 0.350 0 0;2 3 0.080 0.300 0.300 0;2 4 0.0 0.015 0 1.05;3 5 0.0 0.030 0 1.05;];%支路矩阵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阶线路对地电纳值K=zeros(y,1); %K为ij支路y*1阶变压器变比,若K=0表示无变压器,K=1则为标准变比,K不等于1为非标准变比%------------------------------矩阵赋初值：TYPE=point(:,1);%将point矩阵的第一列赋给TYPE,以下类似U=point(:,2);O=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);KT=branch(:,6);%------------------------------求节点导纳矩阵Yfor m=1:y %求Y中非对角元元素YijY(I(m),J(m))=-1/Zij(m);Y(J(m),I(m))=-1/Zij(m);end;for m=1:x %求Y中的Yiifor n=1:yif(I(n)==m|J(n)==m)Y(m,m)=Y(m,m)-Y(I(n),J(n))+j*B0(n)/2; % Yii为Yij加上线路对地电导的一半乘j end;end;end;YG=real(Y);%-----------------------求B'矩阵及其逆矩阵B1B=imag(Y);%求Y的虚部,节点电纳矩阵B11=zeros(x,x);for m=1:yB11(I(m),J(m))=1/(Xij(m));B11(J(m),I(m))=1/(Xij(m));end;for m=1:xfor n=1:yif(I(n)==m|J(n)==m)B11(m,m)=B11(m,m)-1/(Xij(n));end;end;end;p=find(TYPE(:,1)==3);%找出平衡节点编号B11(:,p)=[];%平衡节点编号对应行置空B11(p,:)=[];%平衡节点编号对应列置空B11B1=B11;B1=inv(B1);%B1求逆后得B1矩阵%-----------------------%求B''及其逆矩阵B2npq=find(TYPE(:,1)>1);%找出非PQ节点的编号B22=B; %BB矩阵为中间变量B22(:,npq)=[];%非PQ节点编号对应行置空B22(npq,:)=[];%非PQ节点编号对应行置空B22B2=B22;B2=inv(B2);%求得B2矩阵%-------------计算各节点有功功率不平衡量dPik=0; %k为迭代次数ep=0; %计算P不平衡量dPi的收敛标志(0表示不收敛,1表示收敛)eq=0; %计算U不平衡量dQi的收敛标志(0表示不收敛,1表示收敛)np=find(TYPE(:,1)<3);%找出非平衡节点编号dPi=zeros(x-1,1);%dPi为x*1阶矩阵,x即为节点数pq=find(TYPE(:,1)==1);%找出PQ节点编号pqn=size(B2);pqn=pqn(1);%求PQ节点的个数(因B1矩阵的维数等于PQ节点数)dQi=zeros(pqn,1);%dQi为pqn*1阶矩阵while((~eq|~ep)&(k<100))k=k+1;for m=1:(x-1)%求dPisum1=0;for n=1:xsum1=sum1+U(np(m))*U(n)*(G(np(m),n)*cos(O(np(m))-O(n))+B(np(m),n)*sin(O(n p(m))-O(n)));enddPi(m)=P(np(m))-sum1;endUp=U; %Up为中间变量Up(p)=[];%将平衡节点所在行置空Unp=Up;%求得除平衡节点外的电压列向量dO=(-B1*(dPi./Unp))./Unp;%求相角O的不平衡量for m=1:(x-1) %求相角O的新迭代值矩阵O(np(m))=O(np(m))+dO(m);endmax1=abs(dPi(1)/U(np(1)));%求dP/U绝对值的最大值for m=1:(x-2)if abs(dPi(m)/U(np(m)))<abs(dPi(m+1)/U(np(m+1)))max1=abs(dPi(m+1)/U(np(m+1)));endendif max1<=e %如果最大值满足要求,则ep置为"1",表示收敛ep=1;endfor m=1:pqn %求dQisum2=0;for n=1:xsum2=sum2+U(pq(m))*U(n)*(G(pq(m),n)*sin(O(pq(m))-O(n))-B(pq(m),n)*cos(O(pq (m))-O(n)));enddQi(m)=Q(pq(m))-sum2;endUq=U;%Uq为中间变量Uq(npq)=[];%将非PQ节点所在行置空Upq=Uq;%求得包括PQ节点电压的电压列向量dU=-B2*(dQi./Upq);%求U的不平衡量dUmax2=max(abs(dQi./Upq)); %求dQ/U绝对值的最大值if max2<=e %如果最大值满足要求,则eq置为"1",表示收敛eq=1;endfor m=1:pqn %求U的迭代新值U(pq(m))=U(pq(m))+dU(m);endendsum3=0+j*0;%求平衡节点功率Sphfor m=1:xsum3=sum3+conj(Y(p,m))*(U(m)*cos(O(m))-i*U(m)*sin(O(m)));endSph=(U(p)*cos(O(p))+j*U(p)*sin(O(p)))*sum3;%求平衡节点功率disp(U);figure;plot(U);图3 电压幅值经过23次迭代，得到正常运行下的参数，见表3。