1、电力系统潮流计算的MATLAB辅助程序设计潮流计算程序电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法。高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。lfybus、busout和lineflow程序也可与牛
2、顿法的lfnewton程序和快速解耦法的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=
5、100； accuracy=0.001；accel=1.6；maxiter=80；输入文件准备的第一步是给节点编号，节点号码必须是连续的，但节点数据输入不一定按顺序来编写。此外，还需要下列数据文件：1.节点数据文件busdata：节点信息输入格式为单行输入，输入的数据形成一个矩阵，叫做busdata矩阵。第一列为节点号；第二列为节点类型；第三列和第四列分别为节点电压幅值(标幺值)和相角(单位为度)；第五列和第六列分别为负荷的有功功率和无功功率；第七列到十列分别为发电机的有功功率、无功功率、最小无功出力和最大无功出力；最后一列为并联电容器注入无功功率。第二列的编码用0、1、2来区分PQ节点、平衡
6、节点和PV节点：0表示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
8、=; 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
9、= 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);endend 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;
10、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) = 10if DV(n) = 0.045 if Qgc Qmax(n), Vm(n) = Vm(n) - 0.005; DV(n)=DV(n)+.005; endelse, endelse,endelse,endendif kb(
11、n) = 1 Vc(n) = (conj(S(n)/conj(V(n) - YV )/ Ybus(n,n);else, endif 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);endend maxerror=max( max(abs(real(DP), max(abs(imag(DQ) );if iter = maxiter
12、& 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,
13、tech=( Power Flow Solution by Gauss-Seidel Method);endk=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
14、) = 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 AcurBusDPDQDVLScVcVcIVcRYVconvergedelta2.lfybus.m程序清单% This program obtains the Bus Admi
15、ttance 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)
16、 = 1; elseendYbus=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:nbusfor k=1:nbrif 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, endendendclear P
17、gg3. 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)
18、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), fpr
19、intf( %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
20、 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, endif 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)
21、; 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(
22、Ik)*basemva; SL = Snk + Skn; SLT = SLT + SL;else, endif 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)endelse, e
23、ndendendSLT = SLT/2;fprintf( n), fprintf( Total loss )fprintf(%9.3f, real(SLT), fprintf(%9.3fn, imag(SLT)clear IkInSLSLTSknSnk5.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
24、=; 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
25、 = 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:nbrif mline(ii)=1 if nl(ii) = n | nr(ii) = nif nl(ii) = n , l = nr(ii); endif 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, endif kb(n) = 1 & kb(l)
