电力系统潮流计算的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 GaussSeidel methodVm=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：nbusn=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 + j0； else delta（n） = pi/180delta(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)； endDV（n）=0；endnum = 0； AcurBus = 0； converge = 1；Vc = zeros(nbus，1)+jzeros（nbus，1)； Sc = zeros（nbus，1)+j*zeros（nbus，1）;while exist（accel)=1 accel = 1.3;endwhile exist(accuracy)=1 accuracy = 0。001；endwhile exist(basemva）=1 basemva= 100;endwhile exist（maxiter)=1 maxiter = 100;endmline=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 enditer=0;maxerror=10;while maxerror = accuracy iter = maxiteriter=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 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 - VcI2）; 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.nn） fprintf(Press Enter to terminate the iterations and print the results n） converge = 0; pause， else， endendif 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)+jQ（n）; Qg（n） = Q（n)*basemva + Qd（n） - Qsh(n）； endyload(n） = (Pd(n)- jQd(n）+j*Qsh(n）/(basemvaVm（n）2）;endPgt = 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 delta2。lfybus.m程序清单% This program obtains the Bus Admittance Matrix for power flow solutionj=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 + jX; y= ones（nbr，1）。/Z; %支路导纳for n = 1:nbrif a（n) = 0 a（n) = 1； else endYbus=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)）; endend% 对角元素的数值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 endendclear Pgg3。 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 nn, 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.3fn， 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.3fnn, 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 nn）fprintf（ -Line- Power at bus & line flow -Line loss- Transformern）fprintf(from to MW Mvar MVA MW Mvar tapn）for n = 1：nbusbusprt = 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.3fn， 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)2V(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）2V(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。3fn, a（L） else, fprintf(9.3fn， imag(SL) end else， end endendSLT = SLT/2;fprintf( n）, fprintf（ Total loss )fprintf(%9。3f， real（SLT)）, fprintf(9。3fn， imag(SLT）)clear Ik In SL SLT Skn Snk5.lfnewton。m程序清单% Power flow solution by Newton-Raphson methodns=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：nbusn=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） = accuracy & iter = maxiter for ii=1：mfor k=1：m A(ii,k）=0; %初始化雅可比矩阵end， enditer = iter+1;for n=1：nbusnn=nnss(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+lngs(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 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） = 2Vm（n）Ym（n，n)*cos(t(n，n)）+J22; % J2对角元素 A(lm，nn)= J33； % J3对角元素 A(lm,lm) =-2Vm（n)*Ym(n，n)sin（t（n,n)）J44; % J4对角元素 DC（lm) = Q（n）Qk； endendDX=ADC;for n=1：nbus nn=n-nss(n）; lm=nbus+nngs(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.nn） fprintf(Press Enter to terminate the iterations and print the results n) converge = 0； pause, else， endendif converge = 1 tech= （ ITERATIVE SOLUTION DID NOT CONVERGE); else, tech=（ Power Flow Solution by Newton-Raphson Method)；end V = Vm。cos（delta)+jVm.*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)+jQ(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)； endyload（n） = (Pd(n） jQd(n）+j*Qsh(n)）/（basemvaVm(n)2)；endbusdata(：，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 Solutionns=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:nbusn=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） =