电力系统潮流计算的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 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) = 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) = 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, end endif 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); endyload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n)/(basemva*Vm(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 + j*X; 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)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.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=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 2 if 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; endendDX=ADC;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.nn) 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); endyload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n)/(basemva*Vm(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) = accuracy & iter = maxiter % 检验不平衡功率iter = iter+1;id=0; iv=0;for n=1:nbusnn=n-nss(n);J11=0; J33=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); 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; Qgc = Q(n)*basemva + Qd(n) - Qsh(n); if Qmax(n) = 0 if iter = 10 % the Mvar of generator buses are if Qgc Qmax(n), % bring the generator Mvar within Vm(n) = Vm(n) - 0.005;end % the specified limits. else, end else,end else,end end if kb(n) = 1 id = id+1; DP(id) = P(n)-Pk; DPV(id) = (P(n)-Pk)/Vm(n); end if kb(n) = 0 iv=iv+1; DQ(iv) = Q(n)-Qk; DQV(iv) = (Q(n)-Qk)/Vm(n); endendDd=-B1DPV;DV=-B2DQV;id=0;iv=0; for n=1:nbus if kb(n) = 1 id = id+1; delta(n) = delta(n)+Dd(id); end if kb(n) = 0 iv = iv+1; Vm(n)=Vm(n)+DV(iv); end end maxerror=max(max(abs(DP),max(abs(DQ); if iter = maxiter & maxerror accuracy fprintf(nWARNING: Iterative solution did not converged after ) fprintf(%g, iter), fprintf( iterations.nn) fprintf(Press Ent