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

电力系统潮流计算的MATLAB辅助程序 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 潮流计算，通常指负荷潮流，是电力系统 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。现代电力系统潮流计算的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 主要:高斯法、牛顿法、快速解耦法和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程序清单%PowerflowsolutionbyGauss-SeidelmethodVm=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=[];fork=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);ifVm(n)<=0Vm(n)=1.0;V(n)=1+j*0;elsedelta(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);endDV(n)=0;endnum=0;AcurBus=0;converge=1;Vc=zeros(nbus,1)+j*zeros(nbus,1);Sc=zeros(nbus,1)+j*zeros(nbus,1);whileexist('accel')~=1accel=1.3;endwhileexist('accuracy')~=1accuracy=0.001;endwhileexist('basemva')~=1basemva=100;endwhileexist('maxiter')~=1maxiter=100;endmline=ones(nbr,1);fork=1:nbrform=k+1:nbrif((nl(k)==nl(m))&(nr(k)==nr(m)));mline(m)=2;elseif((nl(k)==nr(m))&(nr(k)==nl(m)));mline(m)=2;else,endendenditer=0;maxerror=10;whilemaxerror>=accuracy&iter<=maxiteriter=iter+1;forn=1:nbus;YV=0+j*0;forL=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);endendSc=conj(V(n))*(Ybus(n,n)*V(n)+YV);Sc=conj(Sc);DP(n)=P(n)-real(Sc);DQ(n)=Q(n)-imag(Sc);ifkb(n)==1S(n)=Sc;P(n)=real(Sc);Q(n)=imag(Sc);DP(n)=0;DQ(n)=0;Vc(n)=V(n);elseifkb(n)==2Q(n)=imag(Sc);S(n)=P(n)+j*Q(n);ifQmax(n)~=0Qgc=Q(n)*basemva+Qd(n)-Qsh(n);ifabs(DQ(n))<=.005&iter>=10ifDV(n)<=0.045ifQgcQmax(n),Vm(n)=Vm(n)-0.005;DV(n)=DV(n)+.005;endelse,endelse,endelse,endendifkb(n)~=1Vc(n)=(conj(S(n))/conj(V(n))-YV)/Ybus(n,n);else,endifkb(n)==0V(n)=V(n)+accel*(Vc(n)-V(n));elseifkb(n)==2VcI=imag(Vc(n));VcR=sqrt(Vm(n)^2-VcI^2);Vc(n)=VcR+j*VcI;V(n)=V(n)+accel*(Vc(n)-V(n));endendmaxerror=max(max(abs(real(DP))),max(abs(imag(DQ))));ifiter==maxiter&maxerror>accuracyfprintf('\nWARNING:Iterativesolutiondidnotconvergedafter')fprintf('%g',iter),fprintf('iterations.\n\n')fprintf('PressEntertoterminatetheiterationsandprinttheresults\n')converge=0;pause,else,endendifconverge~=1tech=('ITERATIVESOLUTIONDIDNOTCONVERGE');else,tech=('PowerFlowSolutionbyGauss-SeidelMethod');endk=0;forn=1:nbusVm(n)=abs(V(n));deltad(n)=angle(V(n))*180/pi;ifkb(n)==1S(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);elseifkb(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);Qsht=sum(Qsh);busdata(:,3)=Vm';busdata(:,4)=deltad';clearAcurBusDPDQDVLScVcVcIVcRYVconvergedelta2.lfybus.m程序清单%ThisprogramobtainstheBusAdmittanceMatrixforpowerflowsolutionj=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;%支路导纳forn=1:nbrifa(n)<=0a(n)=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%对角元素的数值forn=1:nbusfork=1:nbrifnl(k)==nYbus(n,n)=Ybus(n,n)+y(k)/(a(k)^2)+Bc(k);elseifnr(k)==nYbus(n,n)=Ybus(n,n)+y(k)+Bc(k);else,endendendclearPgg3.busout.m程序清单%Thisprogramprintsthepowerflowsolutioninatabulatedform%onthescreen.disp(tech)fprintf('MaximumPowerMismatch=%g\n',maxerror)fprintf('No.ofIterations=%g\n\n',iter)head=['BusVoltageAngle------Load---------Generation---Injected''No.Mag.DegreeMWMvarMWMvarMvar'''];disp(head)forn=1:nbusfprintf('%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))endfprintf('\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程序清单%ThisprogramisusedinconjunctionwithlfgaussorlfNewton%forthecomputationoflineflowandlinelosses.SLT=0;fprintf('\n')fprintf('LineFlowandLosses\n\n')fprintf('--Line--Poweratbus&lineflow--Lineloss--Transformer\n')fprintf('fromtoMWMvarMVAMWMvartap\n')forn=1:nbusbusprt=0;forL=1:nbr;ifbusprt==0fprintf('\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,endifnl(L)==nk=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;elseifnr(L)==nk=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,endifnl(L)==n|nr(L)==nfprintf('%12g',k),fprintf('%9.3f',real(Snk)),fprintf('%9.3f',imag(Snk))fprintf('%9.3f',abs(Snk)),fprintf('%9.3f',real(SL)),ifnl(L)==n&a(L)~=1fprintf('%9.3f',imag(SL)),fprintf('%9.3f\n',a(L))else,fprintf('%9.3f\n',imag(SL))endelse,endendendSLT=SLT/2;fprintf('\n'),fprintf('Totalloss')fprintf('%9.3f',real(SLT)),fprintf('%9.3f\n',imag(SLT))clearIkInSLSLTSknSnk5.lfnewton.m程序清单%PowerflowsolutionbyNewton-Raphsonmethodns=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=[];fork=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);ifVm(n)<=0Vm(n)=1.0;V(n)=1+j*0;elsedelta(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);endendfork=1:nbusifkb(k)==1,ns=ns+1;else,endifkb(k)==2ng=ng+1;else,endngs(k)=ng;nss(k)=ns;endYm=abs(Ybus);t=angle(Ybus);m=2*nbus-ng-2*ns;maxerror=1;converge=1;iter=0;mline=ones(nbr,1);fork=1:nbrform=k+1:nbrif((nl(k)==nl(m))&(nr(k)==nr(m)));mline(m)=2;elseif((nl(k)==nr(m))&(nr(k)==nl(m)));mline(m)=2;else,endendend%雅可比矩阵clearADCJDXwhilemaxerror>=accuracy&iter<=maxiterforii=1:mfork=1:mA(ii,k)=0;%初始化雅可比矩阵end,enditer=iter+1;forn=1:nbusnn=n-nss(n);lm=nbus+n-ngs(n)-nss(n)-ns;J11=0;J22=0;J33=0;J44=0;forii=1:nbrifmline(ii)==1ifnl(ii)==n|nr(ii)==nifnl(ii)==n,l=nr(ii);endifnr(ii)==n,l=nl(ii);endJ11=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));ifkb(n)~=1J22=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,endifkb(n)~=1&kb(l)~=1lk=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));ifkb(l)==0%J2的非对角元素A(nn,lk)=Vm(n)*Ym(n,l)*cos(t(n,l)-delta(n)+delta(l));endifkb(n)==0%J3的非对角元素A(lm,ll)=-Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)-delta(n)+delta(l));endifkb(n)==0&kb(l)==0%J4的非对角元素A(lm,lk)=-Vm(n)*Ym(n,l)*sin(t(n,l)-delta(n)+delta(l));endelseendelse,endelse,endendPk=Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;Qk=-Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;ifkb(n)==1P(n)=Pk;Q(n)=Qk;end%SwingbusPifkb(n)==2Q(n)=Qk;ifQmax(n)~=0Qgc=Q(n)*basemva+Qd(n)-Qsh(n);ifiter<=7ifiter>2ifQgcQmax(n),Vm(n)=Vm(n)-0.01;endelse,endelse,endelse,endendifkb(n)~=1A(nn,nn)=J11;%J1对角元素DC(nn)=P(n)-Pk;endifkb(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;endendDX=A\DC';forn=1:nbusnn=n-nss(n);lm=nbus+n-ngs(n)-nss(n)-ns;ifkb(n)~=1delta(n)=delta(n)+DX(nn);endifkb(n)==0Vm(n)=Vm(n)+DX(lm);endendmaxerror=max(abs(DC));ifiter==maxiter&maxerror>accuracyfprintf('\nWARNING:Iterativesolutiondidnotconvergedafter')fprintf('%g',iter),fprintf('iterations.\n\n')fprintf('PressEntertoterminatetheiterationsandprinttheresults\n')converge=0;pause,else,endendifconverge~=1tech=('ITERATIVESOLUTIONDIDNOTCONVERGE');else,tech=('PowerFlowSolutionbyNewton-RaphsonMethod');endV=Vm.*cos(delta)+j*Vm.*sin(delta);deltad=180/pi*delta;i=sqrt(-1);k=0;forn=1:nbusifkb(n)==1k=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);elseifkb(n)==2k=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);endyload(n)=(Pd(n)-j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);endbusdata(:,3)=Vm';busdata(:,4)=deltad';Pgt=sum(Pg);Qgt=sum(Qg);Pdt=sum(Pd);Qdt=sum(Qd);Qsht=sum(Qsh);6.decouple.m程序清单%FastDecoupledPowerFlowSolutionns=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=[];fork=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);ifVm(n)<=0Vm(n)=1.0;V(n)=1+j*0;elsedelta(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);endifkb(n)==1,ns=ns+1;else,endnss(n)=ns;endYm=abs(Ybus);t=angle(Ybus);ii=0;forib=1:nbusifkb(ib)==0|kb(ib)==2ii=ii+1;jj=0;forjb=1:nbusifkb(jb)==0|kb(jb)==2jj=jj+1;B1(ii,jj)=imag(Ybus(ib,jb));else,endendelse,endendii=0;forib=1:nbusifkb(ib)==0ii=ii+1;jj=0;forjb=1:nbusifkb(jb)==0jj=jj+1;B2(ii,jj)=imag(Ybus(ib,jb));else,endendelse,endendB1inv=inv(B1);B2inv=inv(B2);maxerror=1;converge=1;iter=0;mline=ones(nbr,1);fork=1:nbrform=k+1:nbrif((nl(k)==nl(m))&(nr(k)==nr(m)));mline(m)=2;elseif((nl(k)==nr(m))&(nr(k)==nl(m)));mline(m)=2;else,endendend%开始迭代whilemaxerror>=accuracy&iter<=maxiter%检验不平衡功率iter=iter+1;id=0;iv=0;forn=1:nbusnn=n-nss(n);J11=0;J33=0;forii=1:nbrifmline(ii)==1ifnl(ii)==n|nr(ii)==nifnl(ii)==n,l=nr(ii);endifnr(ii)==n,l=nl(ii);endJ11=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,endelse,endendPk=Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;Qk=-Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;ifkb(n)==1P(n)=Pk;Q(n)=Qk;end%SwingbusPifkb(n)==2Q(n)=Qk;Qgc=Q(n)*basemva+Qd(n)-Qsh(n);ifQmax(n)~=0ifiter<=20%Betweenthe1th&6thiterationsifiter>=10%theMvarofgeneratorbusesareifQgcQmax(n),%bringthegeneratorMvarwithinVm(n)=Vm(n)-0.005;end%thespecifiedlimits.else,endelse,endelse,endendifkb(n)~=1id=id+1;DP(id)=P(n)-Pk;DPV(id)=(P(n)-Pk)/Vm(n);endifkb(n)==0iv=iv+1;DQ(iv)=Q(n)-Qk;DQV(iv)=(Q(n)-Qk)/Vm(n);endendDd=-B1\DPV';DV=-B2\DQV';id=0;iv=0;forn=1:nbusifkb(n)~=1id=id+1;delta(n)=delta(n)+Dd(id);endifkb(n)==0iv=iv+1;Vm(n)=Vm(n)+DV(iv);endendmaxerror=max(max(abs(DP)),max(abs(DQ)));ifiter==maxiter&maxerror>accuracyfprintf('\nWARNING:Iterativesolutiondidnotconvergedafter')fprintf('%g',iter),fprintf('iterations.\n\n')fprintf('PressEntertoterminatetheiterationsandprinttheresults\n')converge=0;pause,else,endendifconverge~=1tech=('ITERATIVESOLUTIONDIDNOTCONVERGE');else,tech=('PowerFlowSolutionbyFastDecoupledMethod');endk=0;V=Vm.*cos(delta)+j*Vm.*sin(delta);deltad=180/pi*delta;clearA;clearDC;clearDXi=sqrt(-1);forn=1:nbusifkb(n)==1S(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);elseifkb(n)==2S(n)=P(n)+j*Q(n);Qg(n)=Q(n)*basemva+Qd(n)-Qsh(n);k=k+1;Pgg(k)=Pg(n);endyload(n)=(Pd(n)-j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);endbusdata(:,3)=Vm';busdata(:,4)=deltad';Pgt=sum(Pg);Qgt=sum(Qg);Pdt=sum(Pd);Qdt=sum(Qd);Qsht=sum(Qsh);clearPkQkDPDQJ11J33B1B1invB2B2invDPVDQVDddeltaibidiiivjbjj四、30节点电力系统计算实例潮流计算时，必须将前面的六个程序保存在MATLAB目录下格式为.m的文件，然后在MATLAB的命令窗口输入如下命令：clearbasemva=100;accuracy=0.001;accel=1.8;maxiter=100;%30节点电力系统%母线--母线电压相角负载发电机注入功率%编号节点幅值角度有功无功有功无功无功最小值无功最大值无功busdata=[111.060.00.00.00.00.0000221.0430.021.7012.740.00.0-40500301.00.02.41.20.00.0000401.060.07.61.60.00.0000521.010.094.219.00.00.0-40400601.00.00.00.00.00.0000701.00.022.810.90.00.0000821.010.030.030.00.00.0-30400901.00.00.00.00.00.00001001.00.05.82.00.00.0-624191121.0820.00.00.00.00.00001201.0011.27.5000001321.071000.000-6240140106.21.600000150108.22.500000160103.51.800000170109.05.800000180103.20.900000190109.53.400000200102.20.7000002101017.511.2000002201000.000000230103.21.600000240108.76.700004.32501000.000000260103.52.3000002701000.0000002801000.000000290102.40.9000003001010.61.900000];%线路数据%busbusRX1/2B1forlineslinedata=[120.01920.05750.026401130.04520.18520.020401240.05700.17370.018401340.01320.03790.004201250.04720.19830.020901260.05810.17630.018701460.01190.04140.004501570.04600.11600.010201670.02670.08200.008501680.01200.04200.004501690.00.20800.00.9786100.556000.9699110.2080019100.1100014120.256000.93212130.1400011214.1231.2559011215.0662.1304011216.0945.1987011415.2210.1997011617.0824.1923011518.1073.2185011819.0639.1292011920.0340.0680011020.0936.2090011017.0324.0845011021.0348.0749011022.0727.1499012122.0116.0236011523.1000.2020012224.1150.1790012324.1320.2700012425.1885.3292012526.2544.3800012527.1093.20870128270.396000.9682729.2198.4153012730.3202.6027012930.2399.453301828.0636.20000.02141628.0169.05990.0651];最后运行程序输入以下命令：lfybus%形成节点导纳矩阵lfgauss%高斯-赛德尔法潮流计算busout%屏幕显示潮流计算结果lineflow%计算并显示线路潮流和损耗将lfgauss变为lfnewton/decouple,即可使用牛顿-拉夫逊法/快速解耦法进行潮流计算，输入以上4个命令行后，即可得到潮流计算结果：