电力系统潮流计算的MATLAB辅助程序设计潮流计算程序Word下载.docx
《电力系统潮流计算的MATLAB辅助程序设计潮流计算程序Word下载.docx》由会员分享，可在线阅读，更多相关《电力系统潮流计算的MATLAB辅助程序设计潮流计算程序Word下载.docx（26页珍藏版）》请在冰豆网上搜索。
输入文件准备的第一步是给节点编号，节点号码必须是连续的，但节点数据输入不一定按顺序来编写。
此外，还需要下列数据文件：
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;
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;
accuracy'
accuracy=0.001;
basemva'
basemva=100;
maxiter'
maxiter=100;
mline=ones（nbr,1）;
nbr
form=k+1:
if（（nl（k）==nl（m））&
（nr（k）==nr（m）））;
mline（m）=2;
elseif（（nl（k）==nr（m））&
（nr（k）==nl（m）））;
else,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）;
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
ifQmax（n）~=0
Qgc=Q（n）*basemva+Qd（n）-Qsh（n）;
ifabs（DQ（n））<
=.005&
iter>
=10
ifDV（n）<
=0.045
ifQgc<
Qmin（n）,
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
ifkb（n）~=1
Vc（n）=（conj（S（n））/conj（V（n））-YV）/Ybus（n,n）;
ifkb（n）==0
V（n）=V（n）+accel*（Vc（n）-V（n））;
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））;
maxerror=max（max（abs（real（DP）））,max（abs（imag（DQ））））;
ifiter==maxiter&
maxerror>
accuracy
fprintf（'
\nWARNING:
Iterativesolutiondidnotconvergedafter'
）
%g'
iter）,fprintf（'
iterations.\n\n'
PressEntertoterminatetheiterationsandprinttheresults\n'
converge=0;
pause,else,end
ifconverge~=1
tech=（'
ITERATIVESOLUTIONDIDNOTCONVERGE'
）;
else,
tech=（'
PowerFlowSolutionbyGauss-SeidelMethod'
k=0;
Vm（n）=abs（V（n））;
deltad（n）=angle（V（n））*180/pi;
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
yload（n）=（Pd（n）-j*Qd（n）+j*Qsh（n））/（basemva*Vm（n）^2）;
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（:
nbus=max（max（nl）,max（nr））;
Z=R+j*X;
y=ones（nbr,1）./Z;
%支路导纳
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））;
%对角元素的数值
forn=1:
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）;
clearPgg
3.busout.m程序清单
%Thisprogramprintsthepowerflowsolutioninatabulatedform
%onthescreen.
disp（tech）
fprintf（'
MaximumPowerMismatch=%g\n'
maxerror）
No.ofIterations=%g\n\n'
iter）
head=['
BusVoltageAngle------Load---------Generation---Injected'
'
No.Mag.DegreeMWMvarMWMvarMvar'
'
];
disp（head）
%5g'
n）,fprintf（'
%7.3f'
Vm（n））,
%8.3f'
deltad（n））,fprintf（'
%9.3f'
Pd（n））,
Qd（n））,fprintf（'
Pg（n））,
%9.3f'
Qg（n））,fprintf（'
%8.3f\n'
Qsh（n））
\n'
）,fprintf（'
Total'
Pdt）,fprintf（'
Qdt）,
Pgt）,fprintf（'
Qgt）,fprintf（'
%9.3f\n\n'
Qsht）
4.lineflow.m程序清单
%ThisprogramisusedinconjunctionwithlfgaussorlfNewton
%forthecomputationoflineflowandlinelosses.
SLT=0;
\n'
LineFlowandLosses\n\n'
--Line--Poweratbus&
lineflow--Lineloss--Transformer\n'
fromtoMWMvarMVAMWMvartap\n'
busprt=0;
ifbusprt==0
%6g'
P（n）*basemva）
%9.3f'
Q（n）*basemva）,fprintf（'
%9.3f\n'
abs（S（n）*basemva））
busprt=1;
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）;
ifnl（L）==n|nr（L）==n
%12g'
k）,
real（Snk））,fprintf（'
imag（Snk））
abs（Snk））,
real（SL））,
ifnl（L）==n&
a（L）~=1
imag（SL））,fprintf（'
a（L））
else,fprintf（'
imag（SL））
SLT=SLT/2;
Totalloss'
real（SLT））,fprintf（'
imag（SLT））
clearIkInSLSLTSknSnk
5.lfnewton.m程序清单
%PowerflowsolutionbyNewton-Raphsonmethod
ns=0;
ng=0;
Vm=0;
deltad=0;
ifkb（k）==1,ns=ns+1;
else,end
ifkb（k）==2ng=ng+1;
ngs（k）=ng;
nss（k）=ns;
Ym=abs（Ybus）;
t=angle（Ybus）;
m=2*nbus-ng-2*ns;
maxerror=1;
converge=1;
iter=0;
%雅可比矩阵
clearADCJDX
=maxiter
forii=1:
m
A（ii,k）=0;
%初始化雅可比矩阵
end,end
iter=iter+1;
nn=n-nss（n）;
lm=nbus+n-ngs（n）-nss（n）-ns;
J11=0;
J22=0;
J33=0;
J44=0;
ifmline（ii）==1
ifnl（ii）==n|nr（ii）==n
ifnl（ii）==n,l=nr（ii）;
ifnr（ii）==n,l=nl（ii）;
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）*cos（t（n,l）-delta（n）+delta（l））;
J44=J44+Vm（l）*Ym（n,l）*sin（t（n,l）-delta（n）+delta（l））;
ifkb（n）~=1&
kb（l）
升级会员 