【文章內容簡介】 global lPQ。 global myf。 openfile。 change。 ybus。 NR。 PQ。 flow。 ret。 output。 打開bus line的數(shù)據(jù)文件的子程序：function openfile global nb nl。global line。global bus。global myf。 [dfile,pathname]=uigetfile(39。*.m39。,39。Select Data File39。)。 if pathname == 0 error(39。 you must select a valid data file39。)else lfile =length(dfile)。 % strip off .m eval(dfile(1:lfile2))。 end [nl,ml]=size(line)。 [nb,mb]=size(bus)。 outfile=39。39。for I=1:length(dfile) if dfile(I)==39。.39。 break。 end outfile=strcat(outfile,dfile(I))。 endoutfile=strcat(outfile,39。.txt39。)。 myf=fopen(outfile,39。w39。)。 節(jié)點重新編號：function change global nb。global nl。global nPQ。global bus。global line。global nodenum。 nSW = 0。 nPV = 0。 nPQ = 0。 for I = 1:nb, type= bus(I,6)。if type == 3, nSW = nSW + 1。 SW(nSW,:)=bus(I,:)。elseif type == 2, nPV = nPV +1。 PV(nPV,:)=bus(I,:)。else nPQ = nPQ + 1。 PQ(nPQ,:)=bus(I,:)。 endend bus=[PQ。PV。SW]。 newbus=[1:nb]39。nodenum=[newbus bus(:,1)]。 bus(:,1)=newbus。 for I=1:nl for J=1:2 for k=1:nb if line(I,J)==nodenum(k,2) line(I,J)=nodenum(k,1)。 break end end endend 建立節(jié)點導納矩陣：function Y = ybus(bus,line) global nl。global nb。global bus。global line。global Y。global myf。 Y=zeros(nb,nb)。 for k=1:nlI=line(k,1)。 J=line(k,2)。 Zt=line(k,3)+j*line(k,4)。if J~=0 Yt=1/Zt。 end%Ym=line(k,5)+j*line(k,6)。 K=line(k,7)。 if (K==0)amp。(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)amp。(J==0) Y(I,I)=Y(I,I)+Ym。 endif K0 Y(I,I)=Y(I,I)+Yt+Ym。 Y(J,J)=Y(J,J)+Yt/K/K。 Y(I,J)=Y(I,J)Yt/K。 Y(J,I)=Y(I,J)。 endif K0 Y(I,I)=Y(I,I)+Yt+Ym。 Y(J,J)=Y(J,J)+K*K*Yt。 Y(I,J)=Y(I,J)+K*Yt。 Y(J,I)=Y(I,J)。 endend 牛拉法解方程組：function NR global nb。global nPQ。global bus。global line。global Y。global myf。 max1=100。 eps1=。 eps2=。 fprintf(myf, 39。節(jié)點導納矩陣Y\n39。)。 for I=1:nb for J=1:nb fprintf(myf, 39。%10f+j*(%10f) 39。, real(Y(I,J)),imag(Y(I,J)))。 end fprintf(myf, 39。\n39。)。end for i=1:max1 angl(:,1)=bus(1:nb1,3)。 u(:,1)=bus(1:nPQ,2)。 x=[angl。u]。 Jac=form_jac(bus,Y)。 del=dPQ(Y,bus)。 dx=Jac\del。 fprintf(myf,39。第%d次迭代結果\n39。,i)。 fprintf(myf,39。第%d次迭代的雅比矩陣J\n39。,i)。 for I=1:nb+nPQ1 for J=1:nb+nPQ1 fprintf(myf, 39。%10f 39。, Jac(I,J))。 end fprintf(myf, 39。\n39。)。 end fprintf(myf,39。第%d次迭代的功率偏差dP和dQ\n39。,i)。 for I=1:nb+nPQ1 fprintf(myf, 39。%10e 39。, del(I,1))。 fprintf(myf, 39。\n39。)。 end fprintf(myf,39。第%d次迭代的節(jié)點相角和電壓的偏差dx\n39。,i)。 for I=1:nb+nPQ1 fprintf(myf, 39。%10e 39。, dx(I,1))。 fprintf(myf, 39。\n39。)。 end for I=nb:nb+nPQ1 dx(I,1)=dx(I,1)*x(I,1)。 end x=xdx。 fprintf(myf,39。第%d次迭代的節(jié)點相角delta和電壓U\n39。,i)。 angl=x(1:nb1,1)。