《潮流计算的快速分解法程序》由会员分享，可在线阅读，更多相关《潮流计算的快速分解法程序（22页珍藏版）》请在金锄头文库上搜索。
1、研 究 生 课 程 设 计 (论 文)电力系统稳态分析Stability Analysis of Power Systems 教学单位 自 动 化 学 院 姓 名 薛 媛 媛 学 号 111101130 年 级 2011 级 专 业 电力系统及其自动化 指导教师 杨 伟 职 称 副 教 授 2011 年 12 月 6 日摘要：快速分解法是一种定雅克比法，形成系数矩阵 、 时忽略了支路电阻、B对地导纳和理想变压器非标准变化，以及 迭代过程中节点电压的不同取值，从而避免每次迭代重新形成雅克比矩阵及其因子表，计算效率大幅提高。本文采用快速分解法中的 XB 型算法，并基于 MATLAB 软件仿真分析具
2、体实例，发现具有较好的收敛性，特别适合在线计算。 关键词：电力系统 潮流计算 快速分解法 雅可比矩阵0 引言用牛顿-拉夫逊法计算潮流时，每次迭代都需要重新形成雅克比矩阵并分解因子表。为避免每次迭代重新形成雅克比矩阵及其因子表，20 世纪 70 年代初Stott 提出了快速分解法。快速分解法采用了一些假设：a、电力系统有功功率主要受电压相角影响，无功功率主要受电压幅值影响。b、高压网线路的 r1);%找出非 PQ 节点的编号 B22=B; %BB 矩阵为中间变量B22(:,npq)=;%非 PQ 节点编号对应行置空B22(npq,:)=;%非 PQ 节点编号对应行置空B22B2=B22; B2=
3、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(
4、pqn,1);%dQi 为 pqn*1 阶矩阵while(eq|ep)&(k100) 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(np(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) %求相角
5、 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),
6、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;end for m=1:pqn %求 U 的迭代新值U(pq(m)=U(pq(m)+dU(m);endendsum3=0+j*0;%求平衡节点功率 Sphfor m
7、=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);1 1.5 2 2.5 3 3.5 4 4.5 5-0.200.20.40.60.811.2 代代代代代代代代图 3 电压幅值经过 23 次迭代，得到正常运行下的参数，见表 3。表 3 节点电压、负荷值图 4 节点导纳矩阵、B和 B矩阵2.程序二-主程序-function Uo,Oo=powerflow(x,y);x= 1,1,2,2,
8、0.040,0.250,0.250,0; 1,1,3,3,0.100,0.350,0.000,0;2,2,3,3,0.080,0.300,0.250,0;2,2,4,4,0.000,0.015,0,1.050;3,3,5,5,0.000,0.030,0,1.050;y= 0,5,0,1.05,1.05,0; 1,4,5.0,1.05,1.05,0; 2,3,-3.7,-1.3,1,0;2,2,-2.0,-1.0,1,0;2,1,-1.6,-0.8,1,0;Pointstyle=y(:,1);PointNumber=y(:,2); Ps=y(:,3);Qs=y(:,4);Uk=y(:,5);Ok
9、=y(:,6);Y=Yn(x,y);Ba=B1(x,y);Bb=B11(x,y);G=real(Y);B=imag(Y);Uo=Uk;Oo=Ok;L=length(PointNumber);P=zeros(L,1);Q=zeros(L,1);dP=zeros(L,1);dQ=zeros(L,1);number=1;i=1;k=1;while Pointstyle(i)=1; i=i+1;k=k+1;end;while number23;for i=2:L;sum1=0;for j=1:L;sum1=sum1+Uo(j)*(G(i,j)*cos(Oo(i)-Oo(j)+B(i,j)*sin(Oo
10、(i)-Oo(j);end;dP(i)=Ps(i)-Uo(i)*sum1;end;for i=k:L;sum2=0;for j=1:L;sum2=sum2+Uo(j)*(G(i,j)*sin(Oo(i)-Oo(j)-B(i,j)*cos(Oo(i)-Oo(j);end; dQ(i)=Qs(i)-Uo(i)*sum2;end;dP1=dP(2:L)./Uo(2:L);dQ1=dQ(k:L)./Uo(k:L);a=max(norm(dP1,inf);b=max(norm(dQ1,inf);if max(a,b)0.00001break;disp(迭代)disp(k);disp(次后收敛);else
11、dO=-inv(Ba)*dP1;dU=-inv(Bb)*dQ1; zeros1=zeros(k-1,1);zeros2=0;DU=zeros1;dU;DO=zeros2;dO;Uo=Uo+DU;Oo=Oo+DO;number=number+1;end;if number=23;disp(迭代 23 次结束);end;end;disp(Yn);-形成节点导纳 Yn 矩阵-function Y=Yn(x,y)x=1,1,2,2,0.040,0.250,0.250,0;1,1,3,3,0.100,0.350,0.000,0;2,2,3,3,0.080,0.300,0.250,0;2,2,4,4,0.
12、000,0.015,0,1.050;3,3,5,5,0.000,0.030,0,1.050;y= 0,5,0,1.05,1.05,0; 1,4,5.0,1.05,1.05,0; 2,3,-3.7,-1.3,1,0;2,2,-2.0,-1.0,1,0;2,1,-1.6,-0.8,1,0;HeadPoint=x(:,1);HeadNumber=x(:,2);EndPoint=x(:,3);EndNumber=x(:,4); R=x(:,5);X=x(:,6);B=x(:,7);k=x(:,8);L=length(HeadNumber);Y=zeros(L,L);for i=1:L;m=HeadNu
13、mber(i);n=EndNumber(i);if n=0;Y(m,m)=Y(m,m)+1j*B(i)+1/(R(i)+1j*X(i);Y(n,n)=Y(n,n)+1j*B(i)+1/(R(i)+1j*X(i);Y(m,n)=Y(m,n)-1/(R(i)+1j*X(i);Y(n,m)=Y(n,m)-1/(R(i)+1j*X(i);elseY(m,m)=Y(m,m)+R(i)+1j*X(i);end;end;maxm=HeadNumber(1);for i=1:L;if maxm=HeadNumber(i);maxm=HeadNumber(i);end;end;maxn=EndNumber(1)
14、;for i=1:L;if maxn=EndNumber(i);maxn=EndNumber(i);end;end;Y=Y(1:max(maxm,maxn),1:max(maxm,maxn);-形成 矩阵-Bfunction Ba=B1(x,y)x=1,1,2,2,0.040,0.250,0.250,0;1,1,3,3,0.100,0.350,0.000,0;2,2,3,3,0.080,0.300,0.250,0;2,2,4,4,0.000,0.015,0,1.050;3,3,5,5,0.000,0.030,0,1.050;y= 0,5,0,1.05,1.05,0; 1,4,5.0,1.05,1.05,0; 2,3,-3.7,-1.3,1,0;2,2,-2.0,-1.0,1,0;2,1,-1.6,-0.8,1,0;X=x(:,6);for i=1:4Ba(i)=-1/X(i);i=i+1;if i=4Ba(i)=-1/X(i+1);end;end;Bc=zeros(4,4);for m=1:4;for n=1:4;if m=n;Bc(m,m)=Ba(m);elseBc(m,n)=0;end;end;end;Ba=Bc;-
