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

电力系统潮流计算的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 method 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:nbus n=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) = 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); end

DV(n)=0; end

num = 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 accel = 1.3; end

while exist('accuracy')~=1 accuracy = 0.001; end

while exist('basemva')~=1 basemva= 100; end

while 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))); mline(m)=2; else, end end end

iter=0;

maxerror=10;

while maxerror >= accuracy & iter <= maxiter iter=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 < 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); 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 - VcI^2);

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.\\n\\n')

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 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); end

yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2); end

Pgt = sum(Pg); Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh); busdata(:,3)=Vm'; busdata(:,4)=deltad';

clear AcurBusDPDQDVLScVcVcIVcRYVconvergedelta

2.lfybus.m程序清单

% This program obtains the Bus Admittance Matrix for power flow solution j=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:nbr

if a(n) <= 0 a(n) = 1; elseend

Ybus=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)); end end

% 对角元素的数值 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 end end

clear Pgg

3. 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 \\n\\n', 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.3f\\n', 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.3f\\n\\n', 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 \\n\\n')

fprintf(' --Line-- Power at bus & line flow --Line loss-- Transformer\\n') fprintf('from to MW Mvar MVA MW Mvar tap\\n') for n = 1:nbus busprt = 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.3f\\n', 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('g', 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.3f\\n', a(L)) else, fprintf('%9.3f\\n', imag(SL)) end

else, end end end

SLT = SLT/2;

fprintf(' \\n'), fprintf(' Total loss ') fprintf('%9.3f', real(SLT)), fprintf('%9.3f\\n', imag(SLT)) clear IkInSLSLTSknSnk

5.lfnewton.m程序清单

%Power flow solution by Newton-Raphson method ns=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:nbus n=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) = 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); end end

for k=1:nbus

if kb(k) == 1, ns = ns+1; else, end if kb(k) == 2 ng = ng+1; else, end ngs(k) = ng; nss(k) = ns; end

Ym=abs(Ybus); t = angle(Ybus); m=2*nbus-ng-2*ns;

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 end end

%雅可比矩阵 clear ADCJDX

while maxerror >= accuracy & iter <= maxiter for ii=1:m for k=1:m

A(ii,k)=0; %初始化雅可比矩阵 end, end

iter = iter+1; for n=1:nbus nn=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 elseend 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 <= 7 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) = 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; end end

DX=A\\DC'; 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.\\n\\n')

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); end

yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2); end

busdata(:,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 Solution ns=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:nbus n=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) = 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); end

if kb(n) == 1, ns = ns+1; else, end nss(n) = ns; end

Ym = abs(Ybus); 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, end end

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(ib,jb)); else,end end

else, end end

B1inv=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 end end

% 开始迭代

while maxerror >= accuracy & iter <= maxiter % 检验不平衡功率 iter = iter+1; id=0; iv=0; for n=1:nbus nn=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 <= 20 % Between the 1th & 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); end end

Dd=-B1\\DPV'; DV=-B2\\DQV'; 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.\\n\\n')

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 Fast Decoupled Method');

end k=0;

V = Vm.*cos(delta)+j*Vm.*sin(delta); deltad=180/pi*delta;

clear A; clear DC; clear DX i=sqrt(-1); for n = 1:nbus 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

S(n)=P(n)+j*Q(n);

Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n); k=k+1;

Pgg(k)=Pg(n); end

yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2); end

busdata(:,3)=Vm'; busdata(:,4)=deltad';

Pgt = sum(Pg); Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh); clear PkQkDPDQJ11J33B1B1invB2B2invDPVDQVDddeltaibidiiivjbjj

四、30节点电力系统计算实例

潮流计算时，必须将前面的六个程序保存在MATLAB目录下格式为.m的文件，然后在MATLAB的命令窗口输入如下命令： clear

basemva = 100; accuracy = 0.001; accel = 1.8; maxiter = 100; % 30节点电力系统

% 母线--母线 电压 相角 负载 发电机 注入功率

% 编号节点 幅值 角度有功 无功有功 无功 无功最小值 无功最大值 无功 busdata=[1 1 1.06 0.0 0.0 0.0 0.0 0.0 0 0 0 2 2 1.043 0.0 21.70 12.7 40.0 0.0 -40 50 0 3 0 1.0 0.0 2.4 1.2 0.0 0.0 0 0 0 4 0 1.06 0.0 7.6 1.6 0.0 0.0 0 0 0 5 2 1.01 0.0 94.2 19.0 0.0 0.0 -40 40 0 6 0 1.0 0.0 0.0 0.0 0.0 0.0 0 0 0

7 0 1.0 0.0 22.8 10.9 0.0 0.0 0 0 0 8 2 1.01 0.0 30.0 30.0 0.0 0.0 -30 40 0 9 0 1.0 0.0 0.0 0.0 0.0 0.0 0 0 0 10 0 1.0 0.0 5.8 2.0 0.0 0.0 -6 24 19 11 2 1.082 0.0 0.0 0.0 0.0 0.0 0 0 0 12 0 1.0 0 11.2 7.5 0 0 0 0 0 13 2 1.071 0 0 0.0 0 0 -6 24 0 14 0 1 0 6.2 1.6 0 0 0 0 0 15 0 1 0 8.2 2.5 0 0 0 0 0 16 0 1 0 3.5 1.8 0 0 0 0 0 17 0 1 0 9.0 5.8 0 0 0 0 0 18 0 1 0 3.2 0.9 0 0 0 0 0 19 0 1 0 9.5 3.4 0 0 0 0 0 20 0 1 0 2.2 0.7 0 0 0 0 0 21 0 1 0 17.5 11.2 0 0 0 0 0 22 0 1 0 0 0.0 0 0 0 0 0 23 0 1 0 3.2 1.6 0 0 0 0 0 24 0 1 0 8.7 6.7 0 0 0 0 4.3 25 0 1 0 0 0.0 0 0 0 0 0 26 0 1 0 3.5 2.3 0 0 0 0 0 27 0 1 0 0 0.0 0 0 0 0 0 28 0 1 0 0 0.0 0 0 0 0 0 29 0 1 0 2.4 0.9 0 0 0 0 0 30 0 1 0 10.6 1.9 0 0 0 0 0]; % 线路数据

% bus bus R X 1/2 B 1 for lines linedata=[1 2 0.0192 0.0575 0.02640 1 1 3 0.0452 0.1852 0.02040 1 2 4 0.0570 0.1737 0.01840 1 3 4 0.0132 0.0379 0.00420 1 2 5 0.0472 0.1983 0.02090 1 2 6 0.0581 0.1763 0.01870 1 4 6 0.0119 0.0414 0.00450 1 5 7 0.0460 0.1160 0.01020 1

6 7 0.0267 0.0820 0.00850 1 6 8 0.0120 0.0420 0.00450 1 6 9 0.0 0.2080 0.0 0.978 6 10 0 .5560 0 0.969 9 11 0 .2080 0 1 9 10 0 .1100 0 1 4 12 0 .2560 0 0.932 12 13 0 .1400 0 1 12 14 .1231 .2559 0 1 12 15 .0662 .1304 0 1 12 16 .0945 .1987 0 1 14 15 .2210 .1997 0 1 16 17 .0824 .1923 0 1 15 18 .1073 .2185 0 1 18 19 .0639 .1292 0 1 19 20 .0340 .0680 0 1 10 20 .0936 .2090 0 1 10 17 .0324 .0845 0 1 10 21 .0348 .0749 0 1 10 22 .0727 .1499 0 1 21 22 .0116 .0236 0 1 15 23 .1000 .2020 0 1 22 24 .1150 .1790 0 1 23 24 .1320 .2700 0 1 24 25 .1885 .3292 0 1 25 26 .2544 .3800 0 1 25 27 .1093 .2087 0 1 28 27 0 .3960 0 0.968 27 29 .2198 .4153 0 1 27 30 .3202 .6027 0 1 29 30 .2399 .4533 0 1 8 28 .0636 .2000 0.0214 1 6 28 .0169 .0599 0.065 1]; 最后运行程序输入以下命令：

lfybus % 形成节点导纳矩阵 lfgauss % 高斯-赛德尔法潮流计算 busout % 屏幕显示潮流计算结果 lineflow % 计算并显示线路潮流和损耗

将lfgauss变为lfnewton/decouple,即可使用牛顿-拉夫逊法/快速解耦法进行潮流计算，输入以上4个命令行后，即可得到潮流计算结果：

本文来源：https://www.bwwdw.com/article/1wz7.html