电力系统潮流计算的MATLAB辅助程序设计潮流计算程序精编文档doc.docx
- 文档编号:3940203
- 上传时间:2022-11-26
- 格式:DOCX
- 页数:26
- 大小:347.10KB
电力系统潮流计算的MATLAB辅助程序设计潮流计算程序精编文档doc.docx
《电力系统潮流计算的MATLAB辅助程序设计潮流计算程序精编文档doc.docx》由会员分享,可在线阅读,更多相关《电力系统潮流计算的MATLAB辅助程序设计潮流计算程序精编文档doc.docx(26页珍藏版)》请在冰豆网上搜索。
电力系统潮流计算的MATLAB辅助程序设计潮流计算程序精编文档doc
【最新整理,下载后即可编辑】
电力系统潮流计算的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
ifQgc Vm(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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 电力系统 潮流 计算 MATLAB 辅助 程序设计 程序 精编 文档 doc