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 和 ma
5、xiter，一般规定为：basemva=100； accuracy=0.001；accel=1.6；maxiter=80；输入文件准备的第一步是给节点编号，节点号码必须是连续的，但节点数据输入不一定按顺序来编写。此外，还需要下列数据文件：1.节点数据文件 busdata：节点信息输入格式为单行输入，输入的数据形成一个矩阵，叫做 busdata 矩阵。第一列为节点号；第二列为节点类型；第三列和第四列分别为节点电压幅值(标幺值)和相角(单位为度)；第五列和第六列分别为负荷的有功功率和无功功率；第七列到十列分别为发电机的有功功率、无功功率、最小无功出力和最大无功出力；最后一列为并联电容器注入无功功率
6、。第二列的编码用 0、1、2 来区分 PQ 节点、平衡节点和 PV 节点：0 表示 PQ 节点，输入正的有功功率（MW）和无功功率（Mvar） ，并且要设定节点电压初始估计值，一般幅值和相角分别设为 1 和 0，若已经给定初始值，则用其给定值来代替 1 和 0。1 表示平衡节点，且已知该节点的电压幅值和相角。2 表示 PV 节点，要设定该节点的节点电压幅值和发电机的有功功率（MW） ，并设定发电机的无功最小出力和最大出力（Mvar） 。2.线路数据文件 linedata 线路数据用节点对的方法来确定，数据包含在称为linedata 的矩阵中。第一列和第二列为节点号码，第三列到第五列为线路电阻、
7、电抗及该线路电纳值的一半，以标幺值表示。最后一列为变压器分接头设定值，对线路来说，需要输入 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=;
8、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);i
9、f Vm(n) = accuracy DV(n)=DV(n)+.005; endelse, endelse,endelse,endendif kb(n) = 1Vc(n) = (conj(S(n)/conj(V(n) - YV )/ Ybus(n,n);else, endif kb(n) = 0V(n) = V(n) + accel*(Vc(n)-V(n);elseif kb(n) = 2VcI = imag(Vc(n);VcR = sqrt(Vm(n)2 - VcI2);Vc(n) = VcR + j*VcI;V(n) = V(n) + accel*(Vc(n) -V(n);endendma
10、xerror=max( max(abs(real(DP), max(abs(imag(DQ) );if iter = maxiter pause, else, endendif converge = 1tech= ( ITERATIVE SOLUTION DID NOT CONVERGE); else, tech=( Power Flow Solution by Gauss-Seidel Method);end k=0;for n = 1:nbusVm(n) = abs(V(n); deltad(n) = angle(V(n)*180/pi;if kb(n) = 1S(n)=P(n)+j*Q(
11、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) =2k=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); Qsh
12、t = 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 =
13、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) = accuracy end else, endelse,endelse,endendif kb(n) = 1A(nn,nn) = J11; % J1对角元素DC(nn) = P(n)-Pk;endif kb(n) = 0A(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;endend
