潮流计算matlab-牛顿拉夫逊法法

用matlab潮流计算(牛顿拉夫逊法)%主程序：[dfile,pathname]=uigetfile('ieee14.m','SelectDataFile');ifpathname==0error('youmustselectavaliddatafile')elselfile=length(dfile);%stripoff.meval(dfile(1:lfile-2));%打开数据文件endglobaln;globalm;[nb,mb]=size(bus);%节点重新编号[nl,ml]=size(line);nSW=0;%平衡节点数目nPV=0;%PV节点数目nPQ=0;%PQ节点数目fori=1:nb,%nb为总节点数type=bus(i,6);iftype==3,nSW=nSW+1;%统计平衡节点数目SW(nSW,:)=bus(i,:);elseiftype==2,nPV=nPV+1;%统计PV节点数目PV(nPV,:)=bus(i,:);elsenPQ=nPQ+1;%统计PQ节点数目PQ(nPQ,:)=bus(i,:);endend;bus=[PQ;PV;SW];newbus=[1:nb]';f=bus(:,1);nodenum=[newbusbus(:,1)];bus(:,1)=newbus;fori=1:nlforj=1:2fork=1:nbifline(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);breakendendendendY=y(bus,line);%形成节点导纳矩阵K=0;%迭代次数初值Kmax=10;%最大迭代次数eps1=1.0e-10;eps2=1.0e-10;m=nPQ;n=nb;Um=eye(m,m);myf=fopen('output1.dat','w');forK=1:Kmaxfori=1:mforj=1:mifi==jUm(i,j)=bus(i,2);endendendb=dPQ(Y,bus);C=jac(bus,Y);dX=C\b';dx=dX';[nx,mx]=size(dx);fori=1:n-1%计算相角bus(i,3)=bus(i,3)-dX(i,1);endB=dx(nx,n:mx)*Um;%计算电压差bus(1:m,2)=bus(1:m,2)-B';%计算电压值dx(nx,n:mx)=B;fprintf(myf,'--第%d次迭代时雅可比矩阵--',K)fprintf(myf,'\n');fori=1:(n+m-1)forj=1:(n+m-1)fprintf(myf,'%8.6f',C(i,j));fprintf(myf,'');endfprintf(myf,'\n');endfprintf(myf,'--第%d次迭代时dPQ的误差--',K)fprintf(myf,'\n');fori=1:(n+m-1)fprintf(myf,'%8.6e',b(1,i));fprintf(myf,'\n');endfprintf(myf,'\n');fprintf(myf,'--第%d次迭代时dx(误差)--',K)fprintf(myf,'\n');fori=1:(n+m-1)fprintf(myf,'%8.6e',dX(i,1));fprintf(myf,'\n');endfprintf(myf,'\n');fprintf(myf,'第%d次迭代后节点电压(仅PQ节点)',K)fprintf(myf,'\n');fori=1:mfprintf(myf,'%8.6f',bus(i,2));fprintf(myf,'\n');endfprintf(myf,'\n');fprintf(myf,'第%d次迭代后相角(角度)',K)fprintf(myf,'\n');fori=1:nfprintf(myf,'%8.6f',bus(i,3)*180/pi);fprintf(myf,'\n');endfprintf(myf,'\n');if(max(abs(dx))eps1)&(max(abs(b))eps2)%判断是否达到计算精度break;endend%计算功率P1=0;T=0;%计算平衡节点的有功forj=1:nT=bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));P1=P1+T;endbus(n,4)=P1;fork=m+1:n%计算PV节点、平衡节点的无功Q1=0;T=0;forj=1:nT=bus(k,2)*bus(j,2)*(real(Y(k,j))*sin(bus(k,3)-bus(j,3))-imag(Y(k,j))*cos(bus(k,3)-bus(j,3)));Q1=Q1+T;endbus(k,5)=Q1;endbus(:,1)=f;%换回各节点、支路的初始编号r=zeros(1,mb);fort=1:nforl=t+1:nifbus(t,1)bus(l,1)r=bus(t,:);bus(t,:)=bus(l,:);bus(l,:)=r;endendendfori=1:nlforj=1:2fork=1:nbifline(i,j)==nodenum(k,1)line(i,j)=nodenum(k,2);breakendendendendfclose(myf);Pf=loss(bus,line);%计算支路潮流及损耗%将节点导纳矩阵、节点潮流计算结果写入文件output2myf=fopen('output2.dat','w');fprintf(myf,'--节点导纳矩阵--\n');fork=1:nforj=1:nfprintf(myf,'%8.6f',real(Y(k,j)));fprintf(myf,'+i*');fprintf(myf,'%8.6f',imag(Y(k,j)));fprintf(myf,'');endfprintf(myf,'\n');endfprintf(myf,'------------------牛顿－拉夫逊法潮流计算结果-------------\n');fprintf(myf,'------------节点计算结果------------\n');fprintf(myf,'--节点节点电压节点相角注入有功功率（P）注入无功功率(Q)类型--\n');forl=1:nbforj=1:mbifj==1|j==6fprintf(myf,'%8.1f',bus(l,j));elseifj==3fprintf(myf,'%8.6f',bus(l,j)*180/pi);elsefprintf(myf,'%8.6f',bus(l,j));endendfprintf(myf,'\n');endfprintf(myf,'--支路计算结果--\n');fprintf(myf,'--节点(I)节点(J)线路功率S(I,J)线路功率S(J,I)线路损耗dS(I,J)--\n');fork=1:nlforj=1:5ifj=2fprintf(myf,'%8.1f',Pf(k,j));fprintf(myf,'');elsefprintf(myf,'%8.6f',real(Pf(k,j)));fprintf(myf,'+i*');fprintf(myf,'%8.6f',imag(Pf(k,j)));fprintf(myf,'');endendfprintf(myf,'\n');endfclose(myf);%根据支路参数建立节点导纳矩阵程序：functionY=y(bus,line)%目的：根据支路参数建立节点导纳矩阵%输入：节点参数矩阵--bus；支路参数矩阵--line%输出：节点导纳矩阵--Y[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);fork=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4);ifZt==0disp('对地支路');Yt=inf;elseYt=1/Zt;endYm=line(k,5)+j*line(k,6);K=line(k,7);if(K==0)&(J~=0)%普通线路Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif(K==0)&(J==0)%对地支路K=J=0，R=X=0Y(I,I)=Y(I,I)+Ym;endifK0%K0时变压器支路Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/K^2;Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endifK0%K0时变压器支路Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt*K^2;Y(I,J)=Y(I,J)+Yt*K;Y(J,I)=Y(I,J);endend%形成雅克矩阵程序：functionJ=jac(bus,Y)%形成雅可比矩阵globaln;globalm;fori=1:(n-1)%计算PQ、PV节点的有功P1=0;T=0;forj=1:nT=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));P1=P1+T;endbus(i,4)=P1;endfori=1:n-1%计算PV、PQ节点的无功Q1=0;T=0;forj=1:nT=bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));Q1=Q1+T;endbus(i,5)=Q1;endfori=1:n-1%计算Hforj=1:n-1ifi~=jH(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));N(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));K(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));L(i,j)=H(i,j);endendendfori=1:n-1%计算Hforj=1:n-1ifj==iH(i,i)=bus(i,5)+imag(Y(i,i))*bus(i,2)^2;N(i,i)=-bus(i,4)-real(Y(i,i))*bus(i,2)^2;K(i,i)=N(i,i)+2*real(Y(i,i))*bus(i,2)^2;L(i,i)=-H(i,i)+2*imag(Y(i,i))*bus(i,2)^2;endendendN=N(1:n-1,1:m);K=K(1:m,1:n-1);L=L(1:m,1:m);J=[H,N;K,L];%计算dPQ的程序：functiondPQ=dPQ(Y,bus)%delP--有功偏差量%delQ--无功偏差量%Y--节点导纳矩阵%bus--节点参数（P,Q,U及相角）矩阵globaln;globalm;delP=zeros(1,n-1);delQ=zeros(1,m);fori=1:(n-1)%形成delP矩阵(PQ、PV节点)P1=0;T=0;forj=1:nT=bus(i,2)*bus(j,2)*(real(Y(i,j))*co