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

1、精选优质文档-倾情为你奉上电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法。高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和

2、快速解耦法的decouple程序联合使用。(读者可以到MATPOWER主页下载MATPOWER4.1,然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算)一、高斯-赛德尔法潮流计算使用的程序： 高斯-赛德法的具体使用方法读者可参考后面的实例，这里仅介绍各程序的编写格式:lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。对于PV节点，如发电机节点，要提供一个无功功率限定值。当给定电压过高

3、或过低时，无功功率可能超出功率限定值。在几次迭代之后(高斯-塞德尔迭代为10次)，需要检查一次发电机节点的无功出力，如果接近限定值，电压幅值进行上下5%的调整，使得无功保持在限定值内。lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。并将这些参数放在名为linedata的文件中。这个程序将阻抗转换为导纳，并得到节点导纳矩阵。busout：该程序以表格形式输出结果，节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。lineflow：该程序输出线路的相关数据，程序设计输出流入线路终端的有功和无功的功率、线损以及节点功率，还包含整

4、个系统的有功和无功损耗。lfnewton是牛顿-拉夫逊法对实际电力系统潮流计算开发的程序，数据准备和程序格式和高斯-赛德尔法一样，包括程序lfybus，busout和 lineflow。 decouple是快速解耦法对实际电力系统潮流计算开发的程序，同高斯法和牛顿法一样需要用到三个程序：lfybus、busout、lineflow。 二、数据准备为了在MATLAB环境下用高斯法进行潮流计算，必须定义下列变量：基准功率，功率允许误差，加速因子和最大迭代次数。上述变量命名(小写字母)为：basemva、accuracy、accel和maxiter，一般规定为：basemva=100； accura

5、cy=0.001；accel=1.6；maxiter=80；输入文件准备的第一步是给节点编号，节点号码必须是连续的，但节点数据输入不一定按顺序来编写。此外，还需要下列数据文件：1.节点数据文件busdata：节点信息输入格式为单行输入，输入的数据形成一个矩阵，叫做busdata矩阵。第一列为节点号；第二列为节点类型；第三列和第四列分别为节点电压幅值(标幺值)和相角(单位为度)；第五列和第六列分别为负荷的有功功率和无功功率；第七列到十列分别为发电机的有功功率、无功功率、最小无功出力和最大无功出力；最后一列为并联电容器注入无功功率。第二列的编码用0、1、2来区分PQ节点、平衡节点和PV节点：0表示

6、PQ节点，输入正的有功功率（MW）和无功功率（Mvar），并且要设定节点电压初始估计值，一般幅值和相角分别设为1和0，若已经给定初始值，则用其给定值来代替1和0。1表示平衡节点，且已知该节点的电压幅值和相角。2表示PV节点，要设定该节点的节点电压幅值和发电机的有功功率（MW），并设定发电机的无功最小出力和最大出力（Mvar）。2.线路数据文件linedata 线路数据用节点对的方法来确定，数据包含在称为linedata的矩阵中。第一列和第二列为节点号码，第三列到第五列为线路电阻、电抗及该线路电纳值的一半，以标幺值表示。最后一列为变压器分接头设定值，对线路来说，需要输入1。线路输入为无输入顺序，

7、对变压器来说，左侧的节点号设为分接头端。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=; fo

8、r 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 + j*0; else delta(n) =

9、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); endDV(n)=0;endnum = 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 ac

10、cel = 1.3;endwhile exist('accuracy')=1 accuracy = 0.001;endwhile exist('basemva')=1 basemva= 100;endwhile 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); ml

11、ine(m)=2; else, end end end iter=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); e

12、nd 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) -

13、 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);

14、 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('nWARNI

15、NG: 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&#

16、39;); 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

17、); 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.lf

18、ybus.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(

19、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 Y

20、bus(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

21、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('

22、 %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', Qsh

23、t)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

24、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 = n

25、r(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

26、)/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(S

27、L), 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

28、 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); de

29、lta(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)=

30、(Pg(n)-Pd(n)/basemva; Q(n)=(Qg(n)-Qd(n)+ Qsh(n)/basemva; S(n) = P(n) + j*Q(n); endendfor k=1:nbusif kb(k) = 1, ns = ns+1; else, endif kb(k) = 2 ng = ng+1; else, endngs(k) = ng;nss(k) = ns;endYm=abs(Ybus); t = angle(Ybus);m=2*nbus-ng-2*ns;maxerror = 1; converge=1;iter = 0; mline=ones(nbr,1);for k=1:n

31、br 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 DXwhile maxerror >= 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=n

32、bus+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(

33、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的非对角元素

34、 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

35、*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)

36、- 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

37、(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(&#

38、39; 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(d

39、elta)+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

40、); 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 = lengt

41、h(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(

42、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); endif kb(n) = 1, ns = ns+1; else, endnss(n) = ns;endYm = abs(Y

43、bus); t = angle(Ybus);ii=0;for ib=1:nbus if kb(ib) = 0 | kb(ib) = 2 ii = ii+1; jj=0; for jb=1:nbus if kb(jb) = 0 | kb(jb) = 2 jj = jj+1; B1(ii,jj)=imag(Ybus(ib,jb); else,end end else, endend ii=0;for ib=1:nbus if kb(ib) = 0 ii = ii+1; jj=0; for jb=1:nbus if kb(jb) = 0 jj = jj+1; B2(ii,jj)=imag(Ybus(

44、ib,jb); else,end end else, endendB1inv=inv(B1); B2inv = inv(B2); 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 endend % 开始迭代while maxerror >= accuracy &

45、amp; 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

46、,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 <= 20 % Between the 1th &

47、 6th iterations if iter >= 10 % the Mvar of generator buses are if Qgc < Qmin(n), % tested. If not within limits Vm(n) Vm(n) = Vm(n) + 0.005; % is changed in steps of 0.05 pu to elseif 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