电力系统潮流计算的MATLAB辅助程序设计潮流计算程序文件.docx
《电力系统潮流计算的MATLAB辅助程序设计潮流计算程序文件.docx》由会员分享,可在线阅读,更多相关《电力系统潮流计算的MATLAB辅助程序设计潮流计算程序文件.docx(31页珍藏版)》请在冰豆网上搜索。
电力系统潮流计算的MATLAB辅助程序设计潮流计算程序文件
电力系统潮流计算的MATLAB辅助程序设计
潮流计算,通常指负荷潮流,是电力系统分析和设计的主要组成部分,对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。
此外,潮流计算还是其它电力系统分析的基础,比如暂态稳定,突发事件处理等。
现代电力系统潮流计算的方法主要:
高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1,这里主要介绍高斯法、牛顿法和快速解耦法。
高斯法的程序是lfgauss,其与lfybus、busout和lineflow程序联合使用求解潮流功率。
lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。
(读者可以到MATPOWER主页下载MATPOWER4.1,然后将其解压到MATLAB目录下,即可使用该软件进行潮流计算)
一、高斯-赛德尔法潮流计算使用的程序:
高斯-赛德法的具体使用方法读者可参考后面的实例,这里仅介绍各程序的编写格式:
lfgauss:
该程序是用高斯法对实际电力系统进行潮流计算,需要用到busdata和linedata两个文件。
程序设计为输入负荷和发电机的有功MW和无功Mvar,以及节点电压标幺值和相角的角度值。
根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。
对于PV节点,如发电机节点,要提供一个无功功率限定值。
当给定电压过高或过低时,无功功率可能超出功率限定值。
在几次迭代之后(高斯-塞德尔迭代为10次),需要检查一次发电机节点的无功出力,如果接近限定值,电压幅值进行上下5%的调整,使得无功保持在限定值内。
lfybus:
这个程序需要输入线路参数、变压器参数以及变压器分接头参数。
并将这些参数放在名为linedata的文件中。
这个程序将阻抗转换为导纳,并得到节点导纳矩阵。
busout:
该程序以表格形式输出结果,节点输出包括电压幅值和相角,发电机和负荷的有功和无功功率,以及并联电容器或电抗器的有功和无功功率。
lineflow:
该程序输出线路的相关数据,程序设计输出流入线路终端的有功和无功的功率、线损以及节点功率,还包含整个系统的有功和无功损耗。
lfnewton是牛顿-拉夫逊法对实际电力系统潮流计算开发的程序,数据准备和程序格式和高斯-赛德尔法一样,包括程序lfybus,busout和lineflow。
decouple是快速解耦法对实际电力系统潮流计算开发的程序,同高斯法和牛顿法一样需要用到三个程序:
lfybus、busout、lineflow。
二、数据准备
为了在MATLAB环境下用高斯法进行潮流计算,必须定义下列变量:
基准功率,功率允许误差,加速因子和最大迭代次数。
上述变量命名(小写字母)为:
basemva、accuracy、accel和maxiter,一般规定为:
basemva=100;accuracy=0.001;accel=1.6;maxiter=80;输入文件准备的第一步是给节点编号,节点号码必须是连续的,但节点数据输入不一定按顺序来编写。
此外,还需要下列数据文件:
1.节点数据文件busdata:
节点信息输入格式为单行输入,输入的数据形成一个矩阵,叫做busdata矩阵。
第一列为节点号;第二列为节点类型;第三列和第四列分别为节点电压幅值(标幺值)和相角(单位为度);第五列和第六列分别为负荷的有功功率和无功功率;第七列到十列分别为发电机的有功功率、无功功率、最小无功出力和最大无功出力;最后一列为并联电容器注入无功功率。
第二列的编码用0、1、2来区分PQ节点、平衡节点和PV节点:
0表示PQ节点,输入正的有功功率(MW)和无功功率(Mvar),并且要设定节点电压初始估计值,一般幅值和相角分别设为1和0,若已经给定初始值,则用其给定值来代替1和0。
1表示平衡节点,且已知该节点的电压幅值和相角。
2表示PV节点,要设定该节点的节点电压幅值和发电机的有功功率(MW),并设定发电机的无功最小出力和最大出力(Mvar)。
2.线路数据文件linedata线路数据用节点对的方法来确定,数据包含在称为linedata的矩阵中。
第一列和第二列为节点号码,第三列到第五列为线路电阻、电抗及该线路电纳值的一半,以标幺值表示。
最后一列为变压器分接头设定值,对线路来说,需要输入1。
线路输入为无输入顺序,对变压器来说,左侧的节点号设为分接头端。
3.zdata是线路数据输入变量,包括四项,前两项是节点编号,后两项是线路电阻和电抗,均以标幺值表示,函数返回节点导纳矩阵。
三、潮流计算的MATLAB程序清单
1.lfgauss.m程序清单
%PowerflowsolutionbyGauss-Seidelmethod
Vm=0;delta=0;yload=0;deltad=0;
nbus=length(busdata(:
1));
kb=[];Vm=[];delta=[];Pd=[];Qd=[];Pg=[];Qg=[];Qmin=[];Qmax=[];
Pk=[];P=[];Qk=[];Q=[];S=[];V=[];
fork=1:
nbus
n=busdata(k,1);
kb(n)=busdata(k,2);Vm(n)=busdata(k,3);delta(n)=busdata(k,4);
Pd(n)=busdata(k,5);Qd(n)=busdata(k,6);Pg(n)=busdata(k,7);Qg(n)=busdata(k,8);
Qmin(n)=busdata(k,9);Qmax(n)=busdata(k,10);
Qsh(n)=busdata(k,11);
ifVm(n)<=0Vm(n)=1.0;V(n)=1+j*0;
elsedelta(n)=pi/180*delta(n);
V(n)=Vm(n)*(cos(delta(n))+j*sin(delta(n)));
P(n)=(Pg(n)-Pd(n))/basemva;
Q(n)=(Qg(n)-Qd(n)+Qsh(n))/basemva;
S(n)=P(n)+j*Q(n);
end
DV(n)=0;
end
num=0;AcurBus=0;converge=1;
Vc=zeros(nbus,1)+j*zeros(nbus,1);Sc=zeros(nbus,1)+j*zeros(nbus,1);
whileexist('accel')~=1
accel=1.3;
end
whileexist('accuracy')~=1
accuracy=0.001;
end
whileexist('basemva')~=1
basemva=100;
end
whileexist('maxiter')~=1
maxiter=100;
end
mline=ones(nbr,1);
fork=1:
nbr
form=k+1:
nbr
if((nl(k)==nl(m))&(nr(k)==nr(m)));
mline(m)=2;
elseif((nl(k)==nr(m))&(nr(k)==nl(m)));
mline(m)=2;
else,end
end
end
iter=0;
maxerror=10;
whilemaxerror>=accuracy&iter<=maxiter
iter=iter+1;
forn=1:
nbus;
YV=0+j*0;
forL=1:
nbr;
if(nl(L)==n&mline(L)==1),k=nr(L);
YV=YV+Ybus(n,k)*V(k);
elseif(nr(L)==n&mline(L)==1),k=nl(L);
YV=YV+Ybus(n,k)*V(k);
end
end
Sc=conj(V(n))*(Ybus(n,n)*V(n)+YV);
Sc=conj(Sc);
DP(n)=P(n)-real(Sc);
DQ(n)=Q(n)-imag(Sc);
ifkb(n)==1
S(n)=Sc;P(n)=real(Sc);Q(n)=imag(Sc);DP(n)=0;DQ(n)=0;
Vc(n)=V(n);
elseifkb(n)==2
Q(n)=imag(Sc);S(n)=P(n)+j*Q(n);
ifQmax(n)~=0
Qgc=Q(n)*basemva+Qd(n)-Qsh(n);
ifabs(DQ(n))<=.005&iter>=10
ifDV(n)<=0.045
ifQgcVm(n)=Vm(n)+0.005;
DV(n)=DV(n)+.005;
elseifQgc>Qmax(n),
Vm(n)=Vm(n)-0.005;
DV(n)=DV(n)+.005;end
else,end
else,end
else,end
end
ifkb(n)~=1
Vc(n)=(conj(S(n))/conj(V(n))-YV)/Ybus(n,n);
else,end
ifkb(n)==0
V(n)=V(n)+accel*(Vc(n)-V(n));
elseifkb(n)==2
VcI=imag(Vc(n));
VcR=sqrt(Vm(n)^2-VcI^2);
Vc(n)=VcR+j*VcI;
V(n)=V(n)+accel*(Vc(n)-V(n));
end
end
maxerror=max(max(abs(real(DP))),max(abs(imag(DQ))));
ifiter==maxiter&maxerror>accuracy
fprintf('\nWARNING:
Iterativesolutiondidnotconvergedafter')
fprintf('%g',iter),fprintf('iterations.\n\n')
fprintf('PressEntertoterminatetheiterationsandprinttheresults\n')
converge=0;pause,else,end
end
ifconverge~=1
tech=('ITERATIVESOLUTIONDIDNOTCONVERGE');else,
tech=('PowerFlowSolutionbyGauss-SeidelMethod');
end
k=0;
forn=1:
nbus
Vm(n)=abs(V(n));deltad(n)=angle(V(n))*180/pi;
ifkb(n)==1
S(n)=P(n)+j*Q(n);
Pg(n)=P(n)*basemva+Pd(n);
Qg(n)=Q(n)*basemva+Qd(n)-Qsh(n);
k=k+1;
Pgg(k)=Pg(n);
elseifkb(n)==2
k=k+1;
Pgg(k)=Pg(n);
S(n)=P(n)+j*Q(n);
Qg(n)=Q(n)*basemva+Qd(n)-Qsh(n);
end
yload(n)=(Pd(n)-j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);
end
Pgt=sum(Pg);Qgt=sum(Qg);Pdt=sum(Pd);Qdt=sum(Qd);Qsht=sum(Qsh);
busdata(:
3)=Vm';busdata(:
4)=deltad';
clearAcurBusDPDQDVLScVcVcIVcRYVconvergedelta
2.lfybus.m程序清单
%ThisprogramobtainstheBusAdmittanceMatrixforpowerflowsolution
j=sqrt(-1);i=sqrt(-1);
nl=linedata(:
1);nr=linedata(:
2);R=linedata(:
3);
X=linedata(:
4);Bc=j*linedata(:
5);a=linedata(:
6);
nbr=length(linedata(:
1));nbus=max(max(nl),max(nr));
Z=R+j*X;y=ones(nbr,1)./Z;%支路导纳
forn=1:
nbr
ifa(n)<=0a(n)=1;elseend
Ybus=zeros(nbus,nbus);%将Ybus初始化为0
%非对角元素的数值
Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);
Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));
end
end
%对角元素的数值
forn=1:
nbus
fork=1:
nbr
ifnl(k)==n
Ybus(n,n)=Ybus(n,n)+y(k)/(a(k)^2)+Bc(k);
elseifnr(k)==n
Ybus(n,n)=Ybus(n,n)+y(k)+Bc(k);
else,end
end
end
clearPgg
3.busout.m程序清单
%Thisprogramprintsthepowerflowsolutioninatabulatedform
%onthescreen.
disp(tech)
fprintf('MaximumPowerMismatch=%g\n',maxerror)
fprintf('No.ofIterations=%g\n\n',iter)
head=['BusVoltageAngle------Load---------Generation---Injected'
'No.Mag.DegreeMWMvarMWMvarMvar'
''];
disp(head)
forn=1:
nbus
fprintf('%5g',n),fprintf('%7.3f',Vm(n)),
fprintf('%8.3f',deltad(n)),fprintf('%9.3f',Pd(n)),
fprintf('%9.3f',Qd(n)),fprintf('%9.3f',Pg(n)),
fprintf('%9.3f',Qg(n)),fprintf('%8.3f\n',Qsh(n))
end
fprintf('\n'),fprintf('Total')
fprintf('%9.3f',Pdt),fprintf('%9.3f',Qdt),
fprintf('%9.3f',Pgt),fprintf('%9.3f',Qgt),fprintf('%9.3f\n\n',Qsht)
4.lineflow.m程序清单
%ThisprogramisusedinconjunctionwithlfgaussorlfNewton
%forthecomputationoflineflowandlinelosses.
SLT=0;
fprintf('\n')
fprintf('LineFlowandLosses\n\n')
fprintf('--Line--Poweratbus&lineflow--Lineloss--Transformer\n')
fprintf('fromtoMWMvarMVAMWMvartap\n')
forn=1:
nbus
busprt=0;
forL=1:
nbr;
ifbusprt==0
fprintf('\n'),fprintf('%6g',n),fprintf('%9.3f',P(n)*basemva)
fprintf('%9.3f',Q(n)*basemva),fprintf('%9.3f\n',abs(S(n)*basemva))
busprt=1;
else,end
ifnl(L)==nk=nr(L);
In=(V(n)-a(L)*V(k))*y(L)/a(L)^2+Bc(L)/a(L)^2*V(n);
Ik=(V(k)-V(n)/a(L))*y(L)+Bc(L)*V(k);
Snk=V(n)*conj(In)*basemva;
Skn=V(k)*conj(Ik)*basemva;
SL=Snk+Skn;
SLT=SLT+SL;
elseifnr(L)==nk=nl(L);
In=(V(n)-V(k)/a(L))*y(L)+Bc(L)*V(n);
Ik=(V(k)-a(L)*V(n))*y(L)/a(L)^2+Bc(L)/a(L)^2*V(k);
Snk=V(n)*conj(In)*basemva;
Skn=V(k)*conj(Ik)*basemva;
SL=Snk+Skn;
SLT=SLT+SL;
else,end
ifnl(L)==n|nr(L)==n
fprintf('%12g',k),
fprintf('%9.3f',real(Snk)),fprintf('%9.3f',imag(Snk))
fprintf('%9.3f',abs(Snk)),
fprintf('%9.3f',real(SL)),
ifnl(L)==n&a(L)~=1
fprintf('%9.3f',imag(SL)),fprintf('%9.3f\n',a(L))
else,fprintf('%9.3f\n',imag(SL))
end
else,end
end
end
SLT=SLT/2;
fprintf('\n'),fprintf('Totalloss')
fprintf('%9.3f',real(SLT)),fprintf('%9.3f\n',imag(SLT))
clearIkInSLSLTSknSnk
5.lfnewton.m程序清单
%PowerflowsolutionbyNewton-Raphsonmethod
ns=0;ng=0;Vm=0;delta=0;yload=0;deltad=0;
nbus=length(busdata(:
1));
kb=[];Vm=[];delta=[];Pd=[];Qd=[];Pg=[];Qg=[];Qmin=[];Qmax=[];
Pk=[];P=[];Qk=[];Q=[];S=[];V=[];
fork=1:
nbus
n=busdata(k,1);
kb(n)=busdata(k,2);Vm(n)=busdata(k,3);delta(n)=busdata(k,4);
Pd(n)=busdata(k,5);Qd(n)=busdata(k,6);Pg(n)=busdata(k,7);Qg(n)=busdata(k,8);
Qmin(n)=busdata(k,9);Qmax(n)=busdata(k,10);
Qsh(n)=busdata(k,11);
ifVm(n)<=0Vm(n)=1.0;V(n)=1+j*0;
elsedelta(n)=pi/180*delta(n);
V(n)=Vm(n)*(cos(delta(n))+j*sin(delta(n)));
P(n)=(Pg(n)-Pd(n))/basemva;
Q(n)=(Qg(n)-Qd(n)+Qsh(n))/basemva;
S(n)=P(n)+j*Q(n);
end
end
fork=1:
nbus
ifkb(k)==1,ns=ns+1;else,end
ifkb(k)==2ng=ng+1;else,end
ngs(k)=ng;
nss(k)=ns;
end
Ym=abs(Ybus);t=angle(Ybus);
m=2*nbus-ng-2*ns;
maxerror=1;converge=1;
iter=0;
mline=ones(nbr,1);
fork=1:
nbr
form=k+1:
nbr
if((nl(k)==nl(m))&(nr(k)==nr(m)));
mline(m)=2;
elseif((nl(k)==nr(m))&(nr(k)==nl(m)));
mline(m)=2;
else,end
end
end
%雅可比矩阵
clearADCJDX
whilemaxerror>=accuracy&iter<=maxiter
forii=1:
m
fork=1:
m
A(ii,k)=0;%初始化雅可比矩阵
end,end
iter=iter+1;
forn=1:
nbus
nn=n-nss(n);
lm=nbus+n-ngs(n)-nss(n)-ns;
J11=0;J22=0;J33=0;J44=0;
forii=1:
nbr
ifmline(ii)==1
ifnl(ii)==n|nr(ii)==n
ifnl(ii)==n,l=nr(ii);end
ifnr(ii)==n,l=nl(ii);end
J11=J11+Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)-delta(n)+delta(l));
J33=J33+Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)-delta(n)+delta(l));
ifkb(n)~=1
J22=J22+Vm(l)*Ym(n,l