电力系统潮流计算的MATLAB辅助程序设计-潮流计算程序.docx

电力系统潮流计算的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程序清单
％ Power flow solution by Gauss—Seidel method
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=［］;
for k=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）;
if Vm（n) <= 0 Vm(n) = 1.0； V(n) = 1 + j＊0；
else delta（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）;
while exist（'accel’)～=1
accel = 1.3;
end
while exist(’accuracy’)~=1
accuracy = 0。001；
end
while exist('basemva’）~=1
basemva= 100;
end
while exist（’maxiter’)～=1
maxiter = 100;
end
mline=ones(nbr，1）;
for k=1：nbr
for m=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;
while maxerror 〉= accuracy ＆ iter 〈= maxiter
iter=iter+1；
for n = 1:nbus;
YV = 0+j*0;
for L = 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）；
if kb（n） == 1
S(n) =Sc; P（n） = real(Sc)； Q（n) = imag（Sc)； DP（n) =0； DQ（n）=0;
Vc（n） = V(n);
elseif kb（n) == 2
Q（n) = imag(Sc)； S(n) = P(n) + j*Q（n）;
if Qmax（n） ～= 0
Qgc = Q（n)*basemva + Qd（n) - Qsh（n)；
if abs(DQ（n)） 〈= 。005 & iter 〉= 10
if DV(n） 〈= 0.045
if Qgc < Qmin(n）,
Vm（n) = Vm（n) + 0。005;
DV(n） = DV(n)+.005；
elseif Qgc > Qmax（n)，
Vm(n) = Vm（n） — 0。005;
DV（n)=DV（n）+。005； end
else， end
else，end
else,end
end
if kb(n) ~= 1
Vc(n） = (conj（S(n）)/conj(V（n）) — YV )/ Ybus（n，n）;
else， end
if kb(n） == 0
V(n） = V（n) + accel*（Vc(n)—V(n));
elseif kb（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））） ）;
if iter == maxiter & maxerror 〉 accuracy
fprintf(’\nWARNING: Iterative solution did not converged after ’)
fprintf(’%g', iter）， fprintf（' iterations.\n\n'）
fprintf(’Press Enter to terminate the iterations and print the results \n'）
converge = 0; pause， else， end
end
if converge ～= 1
tech= （' ITERATIVE SOLUTION DID NOT CONVERGE'）； else，
tech=(’ Power Flow Solution by Gauss-Seidel Method’）;
end
k=0；
for n = 1：nbus
Vm(n） = abs(V（n))； deltad(n） = angle(V(n））*180/pi;
if kb(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）；
elseif kb（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'；
clear AcurBus DP DQ DV L Sc Vc VcI VcR YV converge delta
2。lfybus.m程序清单
% This program obtains the Bus Admittance Matrix for power flow solution
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; %支路导纳
for n = 1:nbr
if a（n) 〈= 0 a（n) = 1； else end
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
% 对角元素的数值
for n=1:nbus
for k=1：nbr
if nl(k）==n
Ybus(n,n) = Ybus（n，n）+y(k）/(a（k)^2) + Bc(k）；
elseif nr(k）==n
Ybus（n,n） = Ybus(n,n）+y(k) +Bc(k）；
else, end
end
end
clear Pgg
3。 busout.m程序清单
% This program prints the power flow solution in a tabulated form
% on the screen.
disp(tech）
fprintf（' Maximum Power Mismatch = %g \n’， maxerror)
fprintf（' No. of Iterations = ％g \n\n', iter）
head =［' Bus Voltage Angle -——-——Load—---—— —-—Generation—-— Injected'
’ No. Mag. Degree MW Mvar MW Mvar Mvar '
' ’]；
disp（head）
for n=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程序清单
％ This program is used in conjunction with lfgauss or lfNewton
％ for the computation of line flow and line losses。
SLT = 0;
fprintf(’\n'）
fprintf（' Line Flow and Losses \n\n'）
fprintf（’ —-Line-- Power at bus & line flow -—Line loss-- Transformer\n'）
fprintf(’from to MW Mvar MVA MW Mvar tap\n'）
for n = 1：nbus
busprt = 0;
for L = 1：nbr；
if busprt == 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
if nl（L）==n k = 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;
elseif nr(L）==n k = 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
if nl（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）)，
if nl（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（' Total loss ’)
fprintf(’%9。3f’， real（SLT)）, fprintf(’％9。3f\n'， imag(SLT）)
clear Ik In SL SLT Skn Snk
5.lfnewton。m程序清单
% Power flow solution by Newton-Raphson method
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=[］;
for k=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)；
if Vm(n） <= 0 Vm（n） = 1。0; V（n) = 1 + j＊0；
else delta(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
for k=1:nbus
if kb(k） == 1， ns = ns+1； else, end
if kb（k) == 2 ng = 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）；
for k=1：nbr
for m=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
%雅可比矩阵
clear A DC J DX
while maxerror >= accuracy & iter <= maxiter
for ii=1：m
for k=1：m
A(ii,k）=0; %初始化雅可比矩阵
end， end
iter = iter+1;
for n=1：nbus
nn=n—nss(n)；
lm=nbus+n-ngs(n)-nss（n）—ns;
J11=0; J22=0; J33=0； J44=0；
for ii=1：nbr
if mline（ii）==1
if nl（ii） == n | nr(ii） == n
if nl(ii） == n ， l = nr(ii）; end
if nr（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）)；
if kb（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）);
else, end
if kb(n) ～= 1 & kb(l） ～=1
lk = nbus+l—ngs(l）-nss(l)—ns；
ll = l —nss(l）;
％ J1的非对角元素
A（nn, ll) =—Vm(n）*Vm(l)＊Ym（n,l)＊sin（t（n,l）— delta（n） + delta（l））；
if kb(l) == 0 ％ J2的非对角元素
A（nn， lk) =Vm(n）＊Ym(n，l)*cos(t(n,l)— delta（n) + delta(l)）；end
if kb（n） == 0 % J3的非对角元素
A(lm， ll） =-Vm(n）＊Vm（l）*Ym(n,l）*cos（t(n,l）- delta(n)+delta（l)); end
if kb(n) == 0 & kb(l） == 0 ％ J4的非对角元素
A（lm, lk) =—Vm（n）*Ym(n,l)＊sin(t(n，l）- delta（n) + delta(l)）;end
else end
else , end
else， end
end
Pk = Vm（n)^2*Ym(n，n)＊cos（t（n,n))+J33;
Qk = -Vm(n）^2*Ym(n，n)*sin(t(n,n））—J11；
if kb(n） == 1 P(n）=Pk; Q（n） = Qk； end % Swing bus P
if kb（n） == 2 Q(n）=Qk;
if Qmax(n) ~= 0
Qgc = Q(n）＊basemva + Qd(n） — Qsh(n);
if iter <= 7
if iter > 2
if Qgc 〈 Qmin（n）,
Vm(n） = Vm（n) + 0。01;
elseif Qgc > Qmax(n),
Vm(n） = Vm（n) - 0。01;end
else， end
else，end
else，end
end
if kb（n) ～= 1
A(nn，nn） = J11; % J1对角元素
DC(nn） = P(n）—Pk;
end
if kb（n） == 0
A（nn，lm） = 2＊Vm（n）＊Ym（n，n)*cos(t(n，n)）+J22; % J2对角元素
A(lm，nn)= J33； % J3对角元素
A(lm,lm) =-2＊Vm（n)*Ym(n，n)＊sin（t（n,n)）—J44; % J4对角元素
DC（lm) = Q（n）—Qk；
end
end
DX=A\DC';
for n=1：nbus
nn=n-nss(n）;
lm=nbus+n—ngs(n)-nss（n)—ns;
if kb(n） ~= 1
delta（n） = delta（n)+DX(nn); end
if kb（n） == 0
Vm（n)=Vm(n)+DX（lm)； end
end
maxerror=max（abs（DC)）；
if iter == maxiter ＆ maxerror 〉 accuracy
fprintf（’\nWARNING: Iterative solution did not converged after ’）
fprintf（’％g’, iter）, fprintf(' iterations.\n\n'）
fprintf('Press Enter to terminate the iterations and print the results \n')
converge = 0； pause, else， end
end
if converge ～= 1
tech= （' ITERATIVE SOLUTION DID NOT CONVERGE'); else,
tech=（' Power Flow Solution by Newton-Raphson Method')；
end
V = Vm。＊cos（delta)+j＊Vm.*sin(delta）;
deltad=180/pi*delta；
i=sqrt(-1)；
k=0;
for n = 1：nbus
if kb（n) == 1
k=k+1；
S（n）= P(n)+j＊Q(n)；
Pg(n) = P(n）＊basemva + Pd（n);
Qg（n） = Q（n）＊basemva + Qd(n) — Qsh(n）；
Pgg（k）=Pg（n）;
Qgg（k）=Qg(n)；
elseif kb（n） ==2
k=k+1;
S（n)=P（n）+j*Q(n)；
Qg(n） = Q（n)*basemva + Qd(n） - Qsh（n）;
Pgg（k)=Pg(n）；
Qgg（k)=Qg（n)；
end
yload（n） = (Pd(n）— j＊Qd(n）+j*Qsh(n)）/（basemva＊Vm(n)^2)；
end
busdata(：，3)=Vm'； busdata(:,4）=deltad';
Pgt = sum（Pg); Qgt = sum（Qg）; Pdt = sum（Pd）; Qdt = sum(Qd); Qsht = sum(Qsh);
6。decouple.m程序清单
% Fast Decoupled Power Flow Solution
ns=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=[]；
for k=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）;
if Vm（n） 〈= 0 Vm(n） =