电力系统潮流计算的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）*
升级会员 