第一篇：电力系统潮流计算的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)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)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))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.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 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);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 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: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.3fn', abs(S(n)*basemva))
busprt = 1;else, end
if nl(L)==n k = nr(L);In =(V(n)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)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 end end
SLT = SLT/2;fprintf(' n'), fprintf(' Total loss ')fprintf('%9.3f', real(SLT)), fprintf('%9.3fn', 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)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=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);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)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=-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 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);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个命令行后，即可得到潮流计算结果：
第二篇：电力系统潮流计算程序设计
电力系统潮流计算程序设计
姓名：韦应顺
学号：2011021052 电力工程学院
牛顿—拉夫逊潮流计算方法具有能够将非线性方程线性化的特点，而使用MATLAB语言是由于MATLAB语言的数学逻辑强，易编译。
【】【】1.MATLAB程序12
Function tisco %这是一个电力系统潮流计算的程序 n=input（‘n请输入节点数：n=’）； m=input(‘请输入支路数：m=’);ph=input(‘n请输入平衡母线的节点号：ph=’)； B1=input（‘n请输入支路信号：B1=’）;%它以矩阵形式存贮支路的情况，每行存贮一条支路 %第一列存贮支路的一个端点 %第二列存贮支路的另一个端点 %第三列存贮支路阻抗
%第四列存贮支路的对地导纳
%第五列存贮变压器的变比，注意支路为1 %第六列存贮支路的序号
B2=input(‘n请输入节点信息：B2=’)； %第一列为电源侧的功率 %第二列为负荷侧的功率 %第三列为该点的电压值
%第四列为该点的类型：1为PQ，2为PV节点，3为平衡节点 A=input（‘n请输入节点号及对地阻抗：A=’）； ip=input（‘n请输入修正值：ip=’）； %ip为修正值);Y=zeros（n）；
Y(p,q)=Y(p,q)-1./(B1(i3)*B1(i5);e=zeros（1，n）；
Y(p,q)=Y(p,q);f=zeros（1，n）；
no=2*ph=1； Y(q,q)=Y(q,q)+1./B1(i3)+B1(i4)/2；
End for i=1：n
G=real(Y);if A（i2）=0
B=imag(Y);p=A(i1)；
Y(p p)=1./A(i2);for i=1：n End e(i)=real(B2(i3));End f(i)=imag(B2(i3));For i=1:m S(i)=B2(i1)-B2(i2);p=B1(i1);V(i)=B2(i3);p=B1(i2);end Y(p,p)=Y(p,p)+1./(B1(i3)*B1(i5)^2+B1(i4)./2P=real(S);Q=imag(S);[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);[De,Di]=hxf(J,D,F,ph,n,no);t=0;while
max(abs(De))>ip&max(abs(Dfi)>ip
t=t+1;
e=e+De；
f=f+Df;
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)；
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)；
[De,Df]=hxf(J,Df,ph,n,no);end v=e+f*j;for i=1:n hh(i)=conj(Y(ph,i)*v(i));end S(ph)=sum(hh)*v(ph);B2(ph,1)=S(ph);V=abs(v);
jd=angle(v)*180/p;resulte1=[A(:,1),real(v),imag(v),V,jd,real(S’)，imag(S’),real(B2(:1)),imag(B2(:1)),real(B2(:2)),imag(B2(:,2))];for i=1:m
a(i)=conj((v(B1(i1))/B1(i5)-v(B1(i2))/B1(i3));
b(i)=v(B1(i1))*a(i)-j*B1(i4)*v(B1(i))^2/2;
c(i)=-v(B1(i2))*a(i)-j*B1(i4)*v(B1(i2))^2/2;end result2=[B1(:,6),B1(:,1),B1(:,2),real(b’),imag(b’),real(c’),imag(c’), real(b’+c’),imag(b’+c’)];printcut(result1,S,b,c,result2);type resultm function [C,D,Df]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)%该子程序是用来求取Df for i=1:n
If
i=ph
C(i)=0;
D(i)=0;
For j=i:n
C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);end
P1=C(i)*e(i)+D(i)*f(i);Q1=C(i)*f(i)-D(i)*e(i);V1=e(i)^2+f(i)^2;If
B2(i4)=2 p=2*i-1;
Df(p)=P(i)-P1;p=p+1;else p=2*i-1;
Df(p)=P(i)-P1;p=p+1;
Df(p)=Q(i)-Q1;end end end Df=Df’;If ph=n Df(no=[];end
function [De,Df]=hxf(J,Df,ph,n,no)%该子函数是为求取De Df DX=JDf;DX1=DX;
x1=length(DX1);if ph=n DX(no)=0;DX(no+1)=0;
For i=(no+2):(x1+2)DX(i)=DX1(i-2);End Else
DX=[DX1,0,0];End k=0;
[x,y]=size(DX);For i=1:2:x K=k+1;
Df(k)=DX(i);De(k)=DX(i+1);End End case 2 Function for j=1:n J=jacci(Y,G,B,PQ,e,f,V,C,D,B2,n,ph,no)X1=G(i,j)*f(i)-B(i,j)*e(i);
X2=G(i,j)*e(i)+B(i,j)*f(i);%该子程序是用来求取jacci矩阵
for i=1:n X3=0;switch B2(i4)X4=0;case 3 P=2*i-1;continue q=2*j-1;case 1 J(p,q)=X1;for j=1:n m=p+1;if
J=&J=ph J(m,q)=X3;X1=G(i)*f(i)-B(i,j)*e(i);q=q+1;X2=G(i,j)*e(i)+B(i,j)*f(i);J(p,q)=X2;X3=-X2;J(m,q)=X4;X4=X1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);p=2*i-1;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=2*j-1;X3=0;J(p,q)=X1;X4=0;m=p+1;P=2*i-1;J(p,q)=X2;q=2*j-1;J(m,q)=X4;J(p,q)=X1;Else if j=&j=jph m=p+1;X1=D(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X3;X2=C(i)+G(i,j)*e(i)+B(i,j)*f(i);q=q+1;X3= C(i)+G(i,j)*e(i)-B(i,j)*f(i);J(p,q)=X2;X4= C(i)+G(i,j)*f(i)-B(i,j)*e(i);J(m,q)=X4;P=2*i-1;end q=2*j-1;end J(p,q)=X1;end m=p+1;end J(m,q)=X3;if ph=n q=q+1;J(no:)=[];J(p,q)=X2;J(no:)=[];J(m,q)=X4;J(:,no)=[];End J(:,no)=[];End
2实例验证 【例题】设有一系统网络结线见图1，各支路阻抗和各节点功率均已以标幺值标示于图1中，其中节点2连接的是发电厂，设节点1电压保持U1＝1.06定值，试计算其中的潮流分布，请输入节点数：n=5 请输入支路数：m=7 请输入平衡母线的节点号：ph=l 请输入支路信息：
BI=[ l 2 0．02+0．06i O l 1；1 3 0．08+0．24i 0 1 2；2 3 0．06+0．18i 0 l 3： 2 4 0．06+0．18i O l 4： 2 5 0．04+0．12i 0 l 5： 3 4 0．01+0．03i 0 l 6： 4 5 0．08+0．24i O 1 7] 请输入节点信息：
B2=[ 0 0 1．06 3；0．2+0．20i 0 1 1；一O．45一O．15i 0 l l；一0．4-0．05i 0 l 1；一0．6—0．1i 0 1 l] 请输入节点号及对地阻抗： A=[l 0；2 0；3 0；4 0；5 O ] 请输入修正值：ip=0．000 0l
参考文献
[1]陈珩．电力系统稳定分析[M]．北京：中国电力出版社，2002：139—187．
[2]郑阿奇．MATLAB实用教程[M]．北京：电子工业出版社，2005：1-243．
[3] 束洪春，孙士云，等．云电送粤交商流混联系统全过 程动态电压研究[J】．中国电力，2008，4l(10)：l-4． SHU Hong—ch吼，SUN Shi-yun，et a1．Research on fun prc'cess dyn锄ic Voltage stabil时of hybrid AC／DC poWer tmnsmission System舶m Yu衄an proVince to G啪gdong province【J】．Electric Power，2008，4l(10)： l-4．
[4] 朱新立，汤涌，等．大电网安全分析的全过程动态仿 真技术[J】．电网技术，2008，32(22)：23—28． SONG Xin—Ii，TANG Yof唱，et a1． Full dyn锄ic simulation for the stabilhy a眦lysis of large power system【J】．Power System融IlrIolo影，2008，32(22)： 23．28．
[5]Roytelm锄I，Shallidehpour S M．A comprehcnsivc long teml dynaIIlic simulation for powcr system recoVery【J】． IEEE Transactions 0n Power Systems，1994，9(3)． [6] 石雩梅，汪志宏，等．发电机励磁系统数学模型及参 数对电网动态稳定性分析结果影响的研究[J】．继电 器，2007，35(21)：22-27．
SHI Xue．mei，WANG Zlli-hon舀et a1．Iksearch on the innuence of g锄e翰to璐baScd ∞de诅iled excitation system models柚d parameterS t0 power铲id dyn锄ic stabil时【J】．Relay，2007，35(2 1)：22-27．
[7] 方思立，朱方．快速励磁系统对系统稳定的影响[J】．中 国电机工程学报，1986，6(1)：20．28．
FANG Si．1i，ZHU Fang．The effbct of f弧t．respon∞
excitation system on the stability of power netwofk【J】． Proceedings ofthe CSEE，1986，6(1)：20-28．
[8] 刘取．电力系统稳定性及发电机励磁控制[M】．北京： 中国电力出版社，2007．
LIU Qu．Power system S诅bility锄d generator excitation control【M】．BeUing：ChiIla Electric Powef Press，2007． [9] Dallachy J L，Anderson T．EXperience with rcplacing ro诅ting exciters wim static exciters【J】．1k InStitution of Electrical Engineers，1 996．
[10] 陈利芳，陈天禄．浅谈自并励励磁系统在大容量机组 中的应用【J】．继电器，2007，35(1)：8l培4． CHEN Li-f抽岛CHEN Tian—lIL Application of 辩l仁exci组tion mode in large capacity髫memtor unit【J】． ReIay'2007，35(1)：81-84．
[11] 方思立，刘增煌，孟庆和．大型汽轮发电机自并励励 磁系统的应用条件【J]．中国电力，1994，27(12)：61．63． FANG Si．Ii，LIU Zeng-hu锄g，MENG Qin争hc．m application conditions of large turbine generator self-excitation system【J】．Electric Powef，1994，27(12)： 61．63．
[12]梁小冰，黄方能．利用EMTDC进行长持续时间过程 的仿真研究【J】．电网技术，2002，26(9)：55．57． LIANG Xiao-bing，HUANG Fan争眦ng．How to cany out simulalion of long dul‘ation processes by use of EMTDC【J】．Power System 11echnology，2002，26(9)： 55-57．
[13]王卉，陈楷，彭哲，等．数字仿真技术在电力系统中 的应用及常用的几种数字仿真工具【J】．继电器，2004，32(21)：7l一75．
wANG Hui，CHEN Kai，PENG zhe，et a1．Application of digital simulation眦hniques棚d severaJ simulation tools in power system[J】．Relay，2004，32(21)：71·75．
[14]IEEE Power Engmeering Socie哆．IEEE std 421．5．2005 IEEE玎ccOmmended practice for excitation system models for power system stabiI时studies【s】．
第三篇：用matlab电力系统潮流计算
题目：潮流计算与matlab
教学单位 电气信息学院
姓 名 学 号
年 级
专 业 电气工程及其自动化
指导教师
职 称 副教授
摘 要
电力系统稳态分析包括潮流计算和静态安全分析。本文主要运用的事潮流计算，潮流计算是电力网络设计与运行中最基本的运算，对电力网络的各种设计方案及各种运行方式进行潮流计算，可以得到各种电网各节点的电压，并求得网络的潮流及网络中的各元件的电力损耗，进而求得电能损耗。本位就是运用潮流计算具体分析，并有MATLAB仿真。
关键词： 电力系统 潮流计算 MATLAB
Abstract Electric power system steady flow calculation and analysis of the static safety analysis.This paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power.The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation.Key words: Power system;Flow calculation;MATLAB simulation
目 录 任务提出与方案论证....................................................................................................................................2 2 总体设计........................................................................................................................................................3 2.1潮流计算等值电路.............................................................................................................................3 2.2建立电力系统模型.............................................................................................................................3 2.3模型的调试与运行.............................................................................................................................3 3 详细设计........................................................................................................................................................4 3.1 计算前提............................................................................................................................................4 3.2手工计算.............................................................................................................................................7 4设计图及源程序...........................................................................................................................................11 4.1MATLAB仿真.......................................................................................................................................11 4.2潮流计算源程序...............................................................................................................................11 5 总结.............................................................................................................................................................31 参考文献..........................................................................................................................................................32 任务提出与方案论证
潮流计算是在给定电力系统网络结构、参数和决定系统运行状态的边界条件的情况下确定系统稳态运行状态的一种基本方法,是电力系统规划和运营中不可缺少的一个重要组成部分。可以说,它是电力系统分析中最基本、最重要的计算,是系统安全、经济分析和实时控制与调度的基础。常规潮流计算的任务是根据给定的运行条件和网路结构确定整个系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。潮流计算的结果是电力系统稳定计算和故障分析的基础。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。是电力系统研究人员长期研究的一个课题。它既是对电力系统规划设计和运行方式的合理性、可靠性及经济性进行定量分析的依据，又是电力系统静态和暂态稳定计算的基础。
潮流计算经历了一个由手工到应用数字电子计算机的发展过程，现在的潮流算法都以计算机的应用为前提用计算机进行潮流计算主要步骤在于编制计算机程序，这是一项非常复杂的工作。对系统进行潮流分析，本文利用 MATLAB中的SimpowerSystems工具箱设计电力系统，在simulink 环境下，不仅可以仿真系统的动态过程，还可以对系统进行稳态潮流分析。
总体设计
SimpowerSystems使用Simulink环境，可以将该系统中的发电机、变压器，线路等模型联结起来，形成电力系统仿真模拟图。在加人测量模块，并对各元件的参数进行设置后，用measurement和sink中的仪器可以观察各元件的电压、电流、功率的大小。
2.1潮流计算等值电路
10MWYN,d11463MWVA15MWGGGG120MW10kVp015.7kWps73kWI%0.50Vs%10.5YN,d1116MWVA463MW“xd0.134x20.161x0.060cosN0.8510kVp011kWps50kWI%0.550Vs%10.5YN,d11210MWVA35kV32km25MW110kV80km25MW110kV70km110kVYN,d11220MWVA20MWGGG415MW”xd0.136x20.16x0.0730cosN0.8p018.6kWps89kWI%0.530MW0Vs%10.510kVp015.7kWps73kWI%0.5GG0V%10.5sYN,d11216MWVA63MWVAp044kWps121kW10kVI%0.350Vs%10.5GG35MWYN,Y,d11210MVA10kVGGG312MW150MW“xd0.128x20.154x0.0540cosN0.85225MW”xd0.128x20.157x0.05910cosN0.8x0.136x20.161x0.0750cosN0.8“d
2.2建立电力系统模型
在Simulink中按照电力系统原型选择元件进行建模。所建立的模型和建立的方法在详细设计中详述。
在电力系统模型的建立工程中主要涉及到的是：元器件的选择及其参数的设置；发电机选型；变压器选择；线路的选择；负荷模型的选择；母线选择。
2.3模型的调试与运行
建立系统模型，并设置好参数以后，就可以在Simulink环境下进行仿真运行。运行的具体结果和分析也在详细设计中详述。
30km35kV0km31p044kWps121kWI%0.350Vs%10.5p013.2kWps63kWI0%0.55V%10.5s(12)80MWVs(23)%6.55MWVs(13)%17.53 详细设计
3.1 计算前提
首先是发电机的参数计算，先对5个发电厂简化为5台发电机来计算。发电机G1：
P141560MWQ160tan(arccos0.8)45MVar发电机G2：
P2463252MWQ2252tan(arccos0.85)156MVar发电机G3：
P331236MWQ336tan(arccos0.8)27MVar发电机G4：
P415050MWQ450tan(arccos0.85)31MVar发电机G5：
P522550MWQ550tan(arccos0.8)37.5MVar
其次是变电站的参数计算，我们还是对7个变电站简化为7台变压器来计算。变压器T1：
2psVN7311023RT1101033.450232SN(1610)2Vs%VN10.51102XT1101079.406SN16103S01p0j变压器T2：（双并联）
I0%SN(0.0157j0.0800)MVA 100RT2XT221psVN18911023101031.3462322SN2(2010)21Vs%VN110.51102101031.7625 32SN22010S022(p0j变压器T3：（四并联）
I0%SN)(0.0372j0.2000)MVA100 1psVN2112111023RT3101030.0922324SN4(6310)XT31Vs%VN2110.5110210105.042 4SN463103I0%SN)(0.1760j0.8820)MVA100S034(p0j变压器T4：（双并联）
1RT11.725021 XT4XT139.70302S042S01(0.0314j0.1600)MVART4变压器T5：
RT54RT30.3680XT54XT320.168S051S03(0.0440j0.2205)MVA41633523100.386 322(1010)
变压器T6：（两个三绕组变压器并联）
RT61RT62RT631[Vs(12)%Vs(13)%Vs(23)%]10.7521Vs2%[Vs(12)%Vs(23)%Vs(13)%]0.25
21Vs3%[Vs(13)%Vs(23)%Vs(12)%]6.752Vs1%21Vs1%VNXT61106.5842SNXT62XT6321Vs2%VN100.1532SN21Vs3%VN104.134 2SNS062(P06j变压器T7：（双并联）
I0%10)(0.0264j0.1100)MVA 100 RT7XT721psVN1503523101030.3062322SN2(1010)21Vs%VN110.535210106.4312SN210103
S072(p0jI0%SN)(0.0220j0.1100)MVA100再次是传输线参数计算，5条传输线的具体计算如下。
根据教材查得r00.21/km x00.4/km b02.810S/km 6线路L1：
线路L2：
线路L3：（双回路）
线路L4：
线路L5：（双回路）RL1r0l10.21408.4XL1x0l10.44016B6L1b0l12.810401.12104S Q1L12BL1V2N0.6776MVarRL2r0l20.2113027.3XL2x0l20.413052B6L2b0l22.8101303.64104S Q1L22BL2V2N2.2022MVarR12rl1L30320.21707.35X11L32x0l320.47014 BL32b40l322.8106703.9210SQ1L32B2L3VN2.3716MVarRL4r0l40.216012.6XL4x0l40.46024BL4b0l42.8106601.68104S Q12L42BL4VN1.0164MVar
11RL5r0l50.21202.12211XL5x0l50.4204 22BL52b0l522.8106201.12104S1QL5BL3VN20.0686MVar23.2手工计算
FLR1：
P2102ST12(RT1jXT1)(3.450j74.406)(0.0285j0.6562)MVA2VN110Sa10MWST1S01jQL1(10.0442j0.1142)MVAP2Q210.044220.11422SL1(RL1jXL1)(8.4j16)(0.070j0.1334)MVAVN21102ST2P2Q2402452(RT2jXT2)(1.346j31.7625)(0.4032j9.5156)MVAVN21102FLR2SbSG120ST260j45200.4032j9.5156(39.5968j35.4844)MVAScSbSa25jQL1SL1(4.4826j35.9144)MVA：
ST3P2Q225221562(RT3jXT3)(0.092j5.042)(0.6679j36.6024)MVA22VN110PQ4.493134.1048(RjX)(27.3j52)(2.67j5.0854)MVAL2L222VN1102222Sc'(4.4931j34.1048)MVASL2FLRSdSG2Sc'120ST3S03jQL2SL2(132.9792j149.229)MVA3：
ST4P2Q262272(RT4jXT4)(1.725j39.703)(0.1091j2.5101)MVAVN21102P2Q2133.59552149.99562(RL3jXL3)(7.35j14)(24.51j46.682)MVA22VN110'Sd(133.5955j149.9956)MVASL3'SeSG3Sd3025ST4S04jQL3SL3(89.945j130.0151)MVAFLR4： ST5P2Q2502312(RT5jXT5)(0.368j20.168)(0.1052j5.7687)MVA22VN110P2Q292.74872133.99372(RL4jXL4)(12.6j24)(27.654j52.674)MVA22VN110Se'(92.7481j133.9937)MVASL4SfSG4Se'80ST5S05jQL4SL4(34.9449j107.3469)MVAFLR5: 152ST72(0.306j6.431)(0.0562j1.1812)MVA35Sh15ST7S07jQL5(15.0782j0.3422)MVA15.078220.34222SL5(2.1j4)(0.3899j0.743)MVA352SiShSL5S06jQL55(20.4945j1.1266)MVA 15237.52ST63(0.386j4.34)(0.514j5.7793)MVA35220.650520.54512ST62(0.386j0.153)(0.1345j0.0533)MVA23526.336298.73692ST61(0.386j6.584)(3.2905j56.1256)MVA352SgSfST61SG5ST62ST63Si35(25.5114j194.12)MVA计算每一个FLR的功率分布和电压分布计算如下： FLR1：
VT2PRQX401.3464531.762512.8970kVVN115 Vb115VT2102.1030kVPRQX10.04428.40.144216VL10.8489kVVb102.1030VaVbVL1101.2541kVFLR2：
功率分布：
SL2ZZZ*T3*L2*SdT3(VbVN)ZZL2****VN(0.092j5.042)(132.9792j149.229)1418.6727.392j57.042T3(4.8812j13.8097)MVAST3ZZZ*L2*L2*SdT3(VbVN)ZZL2VN(27.3j52)(132.9792j149.229)1418.6727.392j57.042T3(108.687j122.62)MVA 电压分布：
Sc1SL2SL2(4.8812j13.8097)(2.67j5.0854)(7.5512j8.7243)MVA7.551227.38.7243522.424kV102.1030VdVbVL2102.103(2.424)104.527kVVL2功率分布：
FLR3：
SL3ZZZ*T4*L3*SeT4(VG3Vd)ZZL3****VN(1.725j39.703)(89.945j130.0151)1037.9279.075j53.73T4(59.444j16.846)MVAST4ZZZ*L3*L3*SeT4(VbVN)ZZL3VN(7.35j14)(89.945j130.0151)1037.9279.075j53.73T4(31.811j60.1256)MVA 电压分布：
Se1SL3SL3(59.444j19.846)(24.51j46.682)(83.954j26.836)MVA83.9547.3526.836149.404kV105.5643VeVdVL396.16kVVL3功率分布：
FLR4：
SL4ZZZ*T5*L4*SfT5(VG3Vd)ZZL4****VN=(0.368j20.168)(34.9449j107.3469)1037.92712.968j44.168T5(20.843j19.689)MVAST4ZZZ*L4*L4*SfT5(VG3Vd)ZZL4VN=(12.6j24)(34.9449j107.3469)1037.92712.968j44.168T5(1.398j44.389)MVA 电压分布：
Se1SL3SL3(59.444j16.846)(24.51j46.682)(83.954j63.528)MVA83.95412.663.5282424.464kV105.5643VeVdVL381.10kVVL4FLR5: 这里我们先将f点和发电机G5当做电源，经过ZT61和ZT63构成两端供电网络以g点作为运算负荷进行计算。ST6ST4(0.386j4.134)(20.2656j70.9293)(22.093837)35(3.900j25.1175)MVA0.772j10.718(0.386j6.584)(20.2656j70.9293)(22.093837)35(16.5061j91.7905)MVA0.772j10.718电压分布：
ST631ST63ST63(16.6421j97.5698)MVA16.64210.38697.56984.13410.9186kV37Vg37VT6326.0814VVT6320.26560.38670.9293(0.153)0.1162kV
26.0814ViVgVT6226.1976VT6220.49452.11.126641.815kV26.1976VhViVL524.3826VL5
4设计图及源程序
4.1MATLAB仿真
相关的原始数据输入格式如下：
1、B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写。
2、对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点编号，将变压器的串联阻抗置于低压侧处理，第三列为支路的串列阻抗参数，第四列为支路的对地导纳参数，第五烈为含变压器支路的变压器的变比，第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，“0”为不含有变压器。
3、B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
4、X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地参数。
4.2潮流计算源程序
%本程序的功能是用牛顿——拉夫逊法进行11节点潮流计算 clear;n=11;%input('请输入节点数：n=');nl=11;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;
0.05620+0.13289i
0.05382i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2];% B1矩阵：
1、支路首端号；
2、末端号；
3、支路阻抗；
4、支路对地电纳 %
5、支路的变比；
6、支路首端处于K侧为1，1侧为0 % B2矩阵：
1、该节点发电机功率；
2、该节点负荷功率；
3、节点电压初始值 %
4、PV节点电压V的给定值；
5、节点所接的无功补偿设备的容量 %
6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %
3为PV节点；
%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路数
if B1(i,6)==0
%左节点处于1侧
p=B1(i,1);q=B1(i,2);
else
%左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元
Y(q,p)=Y(p,q);
%非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
end %求导纳矩阵
disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出导纳阵的实部和虚部
for i=1:n
%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
end for i=1:n
%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量
end %===================== P=real(S);Q=imag(S);
%分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0
% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2;
%电压模平方
%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========
if B2(i,6)~=3
%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========
DP=P(i)-P1;
% PV节点有功误差
DV=V(i)^2-V2;
% PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
end
end
end
end
end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================
for k=3:N0
% N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U
for k2=k1:N1
% 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U
J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1;
% 对角元规格化K行K列对角元素赋1
%==================== 回代运算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行从第三行开始到当前行的前一行k4行消去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end
%用当前行K2列元素减去当前行k列元素乘以第k行K14 列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代运算
==================================
for k3=k1:N0
% 从k+1行到2*n最后一行
for k2=k1:N1
% 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
else
%是第三行k=3
%====================== 第三行k=3的前代运算 ========================
for k3=k1:N0
%从第四行到2n行（最后一行）
for k2=k1:N1
%从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end
%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0;%当前行第3列元素已消为0
end
end
end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
end
%------修改节点电压-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2(a)=IT2;
%不满足要求的节点数
ICT1=ICT1+1;
%迭代次数 end %用高斯消去法解”w=-J*V“
disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%计算各节点电压的角度
E(k)=e(k)+f(k)*j;
%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);
%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);
%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);
%显示各节点的电压相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);
%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);
%显示各节点的注入功率
disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');
for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
%本程序的功能是用牛顿——拉夫逊法进行10节点潮流计算 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 clear;n=10;%input('请输入节点数：n=');nl=10;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2];% B1矩阵：
1、支路首端号；
2、末端号；
3、支路阻抗；
4、支路对地电纳 %
5、支路的变比；
6、支路首端处于K侧为1，1侧为0 % B2矩阵：
1、该节点发电机功率；
2、该节点负荷功率；
3、节点电压初始值 %
4、PV节点电压V的给定值；
5、节点所接的无功补偿设备的容量 %
6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %
3为PV节点；
%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路数
if B1(i,6)==0
%左节点处于1侧
p=B1(i,1);q=B1(i,2);
else
%左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元
Y(q,p)=Y(p,q);
%非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
end %求导纳矩阵
disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出导纳阵的实部和虚部
for i=1:n
%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
end for i=1:n
%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量
end %===================== P=real(S);Q=imag(S);
%分解出各节点注入的有功和无功功率
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0
% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2;
%电压模平方
%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========
if B2(i,6)~=3
%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素
===========
DP=P(i)-P1;
% PV节点有功误差
DV=V(i)^2-V2;
% PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
end
end
end
end
end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================
for k=3:N0
% N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U
for k2=k1:N1
% 从k+1列的Jacobi元素到扩展列的△P、△Q 或 △U
J(k,k2)=J(k,k2)./J(k,k);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1;
% 对角元规格化K行K列对角元素赋1
%==================== 回代运算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行从第三行开始到当前行的前一行k4行消
去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代运算
==================================
for k3=k1:N0
% 从k+1行到2*n最后一行
for k2=k1:N1
% 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
else
%是第三行k=3
%====================== 第三行k=3的前代运算 ========================
for k3=k1:N0
%从第四行到2n行（最后一行）
for k2=k1:N1
%从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end
%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0;%当前行第3列元素已消为0
end
end
end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
end
%------修改节点电压-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2(a)=IT2;
%不满足要求的节点数
ICT1=ICT1+1;
%迭代次数 end %用高斯消去法解”w=-J*V“ disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%计算各节点电压的角度
E(k)=e(k)+f(k)*j;
%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);
%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);
%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);
%显示各节点的电压相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);
%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);
%显示各节点的注入功率
disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
%本程序的功能是用牛顿——拉夫逊法进行12节点潮流计算 %本程序的功能是用牛顿——拉夫逊法进行潮流计算 clear;n=12;%input('请输入节点数：n=');nl=12;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[1
0.03512+0.08306i
0.13455i
0;
0.0068+0.18375i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.00811+0.24549i
0
1.02381
1;
0.05620+0.13289i
0.05382i
0;
0.04215+0.09967i
0.04037i
0;
0.0068+0.18375i
0
1.02381
1;
0.02810+0.06645i
0.10764i
0;
0.05620+0.13289i
0.05382i
0;0.00811+0.24549i
0
1;
0.03512+0.08306i
0.13455i
0;
0.03512+0.08306i
0.13455i
0] B2=[0
0
1.1
1.1
0
1;
0
0
0
0
2;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0.204+0.12638i
0
0
2;
0
0
0
0
2;
0
0.306+0.18962i
0
0
2;
0
0
0
0
2;
0.5
0
1.1
1.1
0
3;
0
0.343+0.21256i
0
0
2;
0
0
0
0
2;
0
0
0
0
2];% B1矩阵：
1、支路首端号；
2、末端号；
3、支路阻抗；
4、支路对地电纳 %
5、支路的变比；
6、支路首端处于K侧为1，1侧为0 % B2矩阵：
1、该节点发电机功率；
2、该节点负荷功率；
3、节点电压初始值 %
4、PV节点电压V的给定值；
5、节点所接的无功补偿设备的容量 %
6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点； %
3为PV节点；
%input('请输入各节点参数形成的矩阵： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl
%支路数
if B1(i,6)==0
%左节点处于1侧
p=B1(i,1);q=B1(i,2);
else
%左节点处于K侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非对角元
Y(q,p)=Y(p,q);
%非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%对角元K侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
%对角元1侧
end %求导纳矩阵
disp('导纳矩阵 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);
%分解出导纳阵的实部和虚部
for i=1:n
%给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
%PV节点电压给定模值
end for i=1:n
%给定各节点注入功率
S(i)=B2(i,1)-B2(i,2);
%i节点注入功率SG-SL
B(i,i)=B(i,i)+B2(i,5);%i节点无功补偿量
end %===================== P=real(S);Q=imag(S);
%分解出各节点注入的有功和无功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次数ICT1、a；不满足收敛要求的节点数IT2 while IT2~=0
% N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列
IT2=0;a=a+1;
for i=1:n
if i~=isb
%非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)
end
P1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)
Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i节点有功和无功功率P',Q'的计算值
V2=e(i)^2+f(i)^2;
%电压模平方
%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========
if B2(i,6)~=3
%非PV节点
DP=P(i)-P1;
%节点有功功率差
DQ=Q(i)-Q1;
%节点无功功率差
%=============== 以上为除平衡节点外其它节点的功率计算 ================= %================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df
X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de
X3=X2;
% X2=dp/df X3=dQ/de
X4=-X1;
% X1=dP/de X4=dQ/df
p=2*i-1;q=2*j1-1;
J(p,q)=X3;J(p,N)=DQ;m=p+1;
% X3=dQ/de J(p,N)=DQ节点无功功率差
J(m,q)=X1;J(m,N)=DP;q=q+1;
% X1=dP/de J(m,N)=DP节点有功功率差
J(p,q)=X4;J(m,q)=X2;
% X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dQ/de
X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df
p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
m=p+1;
J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========
DP=P(i)-P1;
% PV节点有功误差
DV=V(i)^2-V2;
% PV节点电压误差
for j1=1:n
if j1~=isb&j1~=i
%非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i);
% dP/de
X2=B(i,j1)*e(i)-G(i,j1)*f(i);
% dP/df
X5=0;X6=0;
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df
X5=-2*e(i);
X6=-2*f(i);
p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV节点电压误差
m=p+1;
J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;
% PV节点有功误差
J(m,q)=X2;
end
end
end
end
end %========= 以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 =====================
for k=3:N0
% N0=2*n（从第三行开始，第一、二行是平衡节点）
k1=k+1;N1=N;
% N=N0+1 即 N=2*n+1扩展列△P、△Q 或 △U
for k2=k1:N1
% 从k+1列的Jacobi元素到扩展列的△P、△Q
或 △U
J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列对角元素去除K行K列后的非对角元素进行规格化
end
J(k,k)=1;
% 对角元规格化K行K列对角元素赋1
%==================== 回代运算
=======================================
if k~=3
% 不是第三行
k > 3
k4=k-1;
for k3=3:k4
% 用k3行从第三行开始到当前行的前一行k4行消去
for k2=k1:N1 %
k3行后各行上三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行k列元素消为0）
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
if k==N0
%若已到最后一行
break;
end
%================== 前代运算
==================================
for k3=k1:N0
% 从k+1行到2*n最后一行
for k2=k1:N1
% 从k+1列到扩展列消去k+1行后各行下三角元素
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
end
%用当前行K2列元素减去当前行k列元素乘以第k行K2列元素
J(k3,k)=0;%当前行第k列元素已消为0
end
else
%是第三行k=3
%====================== 第三行k=3的前代运算 ========================
for k3=k1:N0
%从第四行到2n行（最后一行）
for k2=k1:N1
%从第四列到2n+1列（即扩展列）
J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算（当前行3列元素消为0）
end
%用当前行K2列元素减去当前行3列元素乘以第三行K2列元素
J(k3,k)=0;%当前行第3列元素已消为0
end
end
end %====上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=====
for k=3:2:N0-1
L=(k+1)./2;
e(L)=e(L)-J(k,N);
%修改节点电压实部
k1=k+1;
f(L)=f(L)-J(k1,N);
%修改节点电压虚部
end
%------修改节点电压-----------
for k=3:N0
DET=abs(J(k,N));
if DET>=pr
%电压偏差量是否满足要求
IT2=IT2+1;%不满足要求的节点数加1
end
end
ICT2(a)=IT2;
%不满足要求的节点数
ICT1=ICT1+1;
%迭代次数 end %用高斯消去法解”w=-J*V" disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:n
V(k)=sqrt(e(k)^2+f(k)^2);
%计算各节点电压的模值
sida(k)=atan(f(k)./e(k))*180./pi;
%计算各节点电压的角度
E(k)=e(k)+f(k)*j;
%将各节点电压用复数表示 end %=============== 计算各输出量 =========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);
%显示各节点的实际电压标幺值E用复数表示 disp('----------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);
%显示各节点的电压大小V的模值 disp('----------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);
%显示各节点的电压相角 for p=1:n
C(p)=0;
for q=1:n
C(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算各节点的注入电流的共轭值
end
S(p)=E(p)*C(p);
%计算各节点的功率 S = 电压 X 注入电流的共轭值 end disp('各节点的功率S为(节点号从小到大排列)：');disp(S);
%显示各节点的注入功率
disp('----------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
else
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
Siz(i)=Si(p,q);
end
disp(Si(p,q));
SSi(p,q)=Si(p,q);
ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
disp(ZF);
disp('----------------------');end disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
if B1(i,6)==0
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
else
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
Sjy(i)=Sj(q,p);
end
disp(Sj(q,p));
SSj(q,p)=Sj(q,p);
ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
disp(ZF);
disp('----------------------');end disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nl
p=B1(i,1);q=B1(i,2);
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
DDS(i)=DS(i);
ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
disp(ZF);
disp('----------------------');end
如果源程序的运行结果需要作图可用下面的程序
figure(1);subplot(1,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(1,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('节点号');ylabel('节点注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('节点号');ylabel('节点注入无功');grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路号');ylabel('支路首端注入有功');grid on;subplot(2,2,4);bar(Q1);xlabel('支路号');ylabel('支路首端注入无功');grid on;
总结
通过本次课程设计让我有复习了一次潮流计算的相关知识，跟家清晰了什么事潮流计算以及潮流计算的在电力系统的重要性。电力系统的稳定运行状况即是正常运行状况，是指电力系统在稳定运行条件下电压、功率的分布，也称为潮流分布。电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。同时潮流分布也是电力系统规划设计的一项基础工作。
整个计算过程的模型建立并不是十分复杂，但计算过程十分繁琐、计算量相当的大，而且由于枝节太多很容易算错。不过在计算潮流计算的过程中却对以往学过的电力系统分析的相关知识进行了一次较为深入的复习。而且整个计算对计算量的要求很大，锻炼了我们的计算能力。而且对细节的把握也得到了锻炼，做题的精细程度得到了提高。
参考文献
[1]何仰赞, 温增银《电力系统分析》（第三版）［M].华中科技大学，2002 [2]http://baike.baidu.com/view/627420.[3]王守相，刘玉田 电力系统潮流计算研究现状--《山东电力技术》1996年05期
[4]何仰赞，温增银.电力系统分析（上册）第三版［M].湖北：华中科技大学出版社，2002 [5] 刘同娟.MATLAB在电路分析中的应用.电气电子教学学报.2002 [6] 西安交通大学等.电力系统计算[M].北京：水利电力出版社，1993.12 [7] 李光琦.电力系统暂态分析[M].北京： 水利电力出版社，2002.5 [8]何仰赞，温增银.电力系统分析（下册）第三版［M].湖北：华中科技大学出版社，2002 [9]韦化，李滨，杭乃善，等．大规模水一火电力系统最优潮流的现代内点算法实现[J]．中国电机工程学报，2003．23(6)：13一l8．
[10]Chen Luo—nan，Suzuki Hideki，Katou Kazuo．Mean fieldtheory for optimal power flow[J]．IEEE Transactions OilPower Systems，1997，12(4)：1481·1486
第四篇：电力系统潮流计算
南 京 理 工 大 学
《电力系统稳态分析》
课程报告
姓名
XX
学 号： 5*** 自动化学院 电气工程
基于牛顿-拉夫逊法的潮流计算例题编程报学院(系): 专
业: 题
目: 任课教师 硕士导师 告
杨伟 XX
2015年6月10号
基于牛顿-拉夫逊法的潮流计算例题编程报告
摘要：电力系统潮流计算的目的在于：确定电力系统的运行方式、检查系统中各元件是否过压或者过载、为电力系统继电保护的整定提供依据、为电力系统的稳定计算提供初值、为电力系统规划和经济运行提供分析的基础。潮流计算的计算机算法包含高斯—赛德尔迭代法、牛顿-拉夫逊法和P—Q分解法等，其中牛拉法计算原理较简单、计算过程也不复杂，而且由于人们引入泰勒级数和非线性代数方程等在算法里从而进一步提高了算法的收敛性和计算速度。同时基于MATLAB的计算机算法以双精度类型进行数据的存储和运算, 数据精确度高，能进行潮流计算中的各种矩阵运算，使得传统潮流计算方法更加优化。
一 研究内容
通过一道例题来认真分析牛顿-拉夫逊法的原理和方法（采用极坐标形式的牛拉法），同时掌握潮流计算计算机算法的相关知识，能看懂并初步使用MATLAB软件进行编程，培养自己电力系统潮流计算机算法编程能力。
例题如下：用牛顿-拉夫逊法计算下图所示系统的潮流分布，其中系统中5为平衡节点，节点5电压保持U=1.05为定值，其他四个节点分别为PQ节点，给定的注入功率如图所示。计算精度要求各节点电压修正量不大于10-6。
二 牛顿-拉夫逊法潮流计算 1 基本原理
牛顿法是取近似解x(k)之后，在这个基础上，找到比x(k)更接近的方程的根，一步步地迭代，找到尽可能接近方程根的近似根。牛顿迭代法其最大优点是在方程f(x)=0的单根附近时误差将呈平方减少，而且该法还可以用来求方程的重根、复根。电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点的电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成节点电压新的初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。2 基本步骤和设计流程图
形成了雅克比矩阵并建立了修正方程式，运用牛顿-拉夫逊法计算潮流的核心问题已经解决，已有可能列出基本计算步骤并编制流程图。由课本总结基本步骤如下：
1)形成节点导纳矩阵Y；
2)设各节点电压的初值，如果是直角坐标的话设电压的实部e和虚部f；如果是极坐标的话则设电压的幅值U和相角a；
3)将各个节点电压的初值代入公式求修正方程中的不平衡量以及修正方程的系数矩阵的雅克比矩阵；
4)解修正方程式，求各节点电压的变化量，即修正量； 5)计算各个节点电压的新值，即修正后的值；
6)利用新值从第（3）步开始进入下一次迭代，直至达到精度退出循环； 7)计算平衡节点的功率和线路功率，输出最后计算结果； ① 公式推导
② 流程图
三
matlab编程代码
clear;
% 如图所示1，2，3，4为PQ节点,5为平衡节点
y=0;
% 输入原始数据，求节点导纳矩阵
y(1,2)=1/(0.07+0.21j);
y(4,5)=0;y(1,3)=1/(0.06+0.18j);
y(1,4)=1/(0.05+0.10j);
y(1,5)=1/(0.04+0.12j);
y(2,3)=1/(0.05+0.10j);
y(2,5)=1/(0.08+0.24j);
y(3,4)=1/(0.06+0.18j);
for i=1:5
for j=i:5
y(j,i)=y(i,j);
end
end
Y=0;
% 求节点导纳矩阵中互导纳
for i=1:5
for j=1:5
if i~=j
Y(i,j)=-y(i,j);
end
end
end
% 求节点导纳矩阵中自导纳
for i=1:5
Y(i,i)=sum(y(i,:));
end
Y
% Y为导纳矩阵
G=real(Y);
B=imag(Y);% 输入原始节点的给定注入功率
S(1)=0.3+0.3j;
S(2)=-0.5-0.15j;
S(3)=-0.6-0.25j;
S(4)=-0.7-0.2j;
S(5)=0;
P=real(S);
Q=imag(S);
% 赋初值，U为节点电压的幅值，a为节点电压的相位角
U=ones(1,5);
U(5)=1.05;
a=zeros(1,5);
x1=ones(8,1);
x2=ones(8,1);
k=0;
while max(x2)>1e-6
for i=1:4
for j=1:4
H(i,j)=0;
N(i,j)=0;
M(i,j)=0;
L(i,j)=0;
oP(i)=0;
oQ(i)=0;
end
end
% 求有功、无功功率不平衡量
for i=1:4
for j=1:5
oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));
end
oP(i)=oP(i)+P(i);
oQ(i)=oQ(i)+Q(i);
end
x2=[oP,oQ]';
% x2为不平衡量列向量
% 求雅克比矩阵
% 当i~=j时，求H,N,M,L
for i=1:4
for j=1:4
if i~=j
H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));
N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
L(i,j)=H(i,j);
M(i,j)=-N(i,j);
end
end
end
% 当i=j时，求H,N,M,L
for i=1:4
for j=1:5
if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));
L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))
end
end
N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);
L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);
end
J=[H,N;M,L]
% J为雅克比矩阵
x1=-((inv(J))*x2);
% x1为所求△x的列向量
% 求节点电压新值，准备下一次迭代
for i=1:4
oa(i)=x1(i);
oU(i)=x1(i+4)*U(i);
end
for i=1:4
a(i)=a(i)+oa(i);
U(i)=U(i)+oU(i);
end
k=k+1;
end
k,U,a
% 求节点注入功率
i=5;
for j=1:5
P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i);
Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i);
end
S(5)=P(5)+Q(5)*sqrt(-1);
S
% 求节点注入电流
I=Y*U'
四
运行结果
节点导纳矩阵
经过五次迭代后的雅克比矩阵
迭代次数以及节点电压的幅值和相角（弧度数）
节点注入功率和电流
五 结果分析
在这次学习和实际操作过程里：首先，对电力系统分析中潮流计算的部分特别是潮流计算的计算机算法中的牛顿-拉夫逊法进行深入的研读，弄明白了其原理、计算过程、公式推导以及设计流程。牛顿-拉夫逊法是求解非线性方程的迭代过程，其计算公式为FJX，式中J为所求函数的雅可比矩阵；X为需要求的修正值；F为不平衡的列向量。利用x(*)=x(k+1)+X(k+1)进行多次迭代，通过迭代判据得到所需要的精度值即准确值x(*)。六 结论
通过这个任务，自己在matlab编程，潮流计算，word文档的编辑功能等方面均有提高，但也暴漏出一些问题：理论知识储备不足，对matlab的性能和特点还不能有一个全面的把握，对word软件也不是很熟练，相信通过以后的学习能弥补这些不足，达到一个新的层次。
第五篇：基于MATLAB的电力系统潮流计算
基于MATLAB的电力系统潮流计算
%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：
%B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写 %对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点 %编号，将变压器的串联阻抗置于低压侧处理。%第三列为支路的串列阻抗参数。%第四列为支路的对地导纳参数。
%第五烈为含变压器支路的变压器的变比
%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。
%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点 %负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中 %“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地 %参数。
n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');Y=zeros(n);Times=1;%置迭代次数为初始值 %创建节点导纳矩阵 for i=1:n1 if B1(i,6)==0 %不含变压器的支路 p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else %含有变压器的支路 p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));end end Y OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%将OrgS、DetaS初始化
%创建OrgS，用于存储初始功率参数 h=0;j=0;for i=1:n %对PQ节点的处理 if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end OrgS %创建PVU 用于存储PV节点的初始电压 PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)==3 t=t+1;PVU(t,1)=B2(i,3);end end PVU %创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1:n %对PQ节点的处理 if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end end t=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3 h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;end end DetaS %创建I，用于存储节点电流参数 i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理 if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==j %对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);end if k==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=0;end end end end end k=0;for i=1:n %对PV节点的处理 if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==j %对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理
Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k==(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行
k=0;end end end end end Jacbi %求解修正方程，获取节点电压的不平衡量 DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正节点电压 j=0;for i=1:n %对PQ节点处理 if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n %对PV节点的处理 if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 %开始循环********************************************************************** while abs(max(DetaU))>pr OrgS=zeros(2*n-2,1);%！!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵 h=0;j=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end for i=1:n if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));end end end OrgS %创建DetaS h=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);end end t=0;for i=1:n if i~=isb&B2(i,6)==3 h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;end end DetaS %创建I i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %创建Jacbi Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);end if k==(n-1)k=0;end end end end end k=0;for i=1:n if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k==(n-1)k=0;end end end end end Jacbi DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正节点电压 j=0;for i=1:n if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 Times=Times+1;%迭代次数加1 end Times 一个原始数据的例子 节点数 5 支路数 5平衡节点编号 5 精度pr 0.000001 B1(支路参数矩阵)[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1] B2(节点参数矩阵)[0-1.6-0.8i 1 0 0 2;0-2-1i 1 0 0 2;0-3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1] X(节点号和对地参数)[1 0;2 0;3 0;4 0;5 0]
电力系统潮流计算
——9结点算例-PQ法
原始数据录入data.txt文档:
标号,起始结点,终止结点,支路电阻参数,支路电抗参数,支路对地导纳参数 1,2,5,0.0,0.063,0.0, 2,5,9,0.019,0.072,0.075, 3,6,9,0.012,0.101,0.105, 4,3,6,0.0,0.059,0.0, 5,6,8,0.039,0.17,0.179, 6,4,8,0.017,0.092,0.079, 7,5,7,0.032,0.161,0.153, 8,4,7,0.01,0.085,0.088, 9,1,4,0.0,0.058,0.0, 潮流程序chaoliu.txt文档: #include ; float x1[N-1],x2 ;for(i=1;i guass(1,N-1,y1,B,x1);for(i=1;i guass(N-M,M,y2,B,x2);for(i=N-M;i else { kp=0;if(kq==0)val(u,g,b,r,ku,kr,h);else goto top;} } } void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]){ float ps=0,pv1=0,pv2=0;float qs=0,qv1=0,qv2=0;float p[N][N]={0};float q[N][N]={0};float s[N][N];float dp[N][N]={0};float dq[N][N]={0};float ds[N][N];float dSp=0,dSq=0;int i,j;FILE *fp1;printf(”n=====ping heng jie dian gong lv =====n“);getch();for(i=0;i printf(”n=======shu ju bao cun=====n“);fp1=fopen(”jieguo.txt“,”w+“);{ fprintf(fp1,”xian lu cao liu:n“);for(i=0;i
点击咨询 