【正文】 s(k,1)。 S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))bus(I,2)*bus(J,2)*(conj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtYm(k,3))。 else S_res(k,4)=bus(J,2)^2*conj(YtYm(k,5))。 toc。 end a=0:pi/12:2*pi。 eps1=。 for i=1:max1, A=Fd(x)。 sprintf(39。 2*x(1), 162*(x(2)+),cos(x(3))。 y(3)=exp(x(2)*x(3))+20*x(3)+(10*pi3)/3。 end format long max1=100。 for i=1:max1, A=Fd(x)。 sprintf(39。 y(2)=3*x(1)*x(2)+x(1)*x(3)11*x(1)。y(3)]。 end 。 3*x(2)+x(3)11, 3*x(1),x(1)。 y=[y(1)。,i) [x dx F(x)] %在屏幕上輸出每次的 x(i),dx(i),F(x(i)) if (max(abs(dx))eps)amp。 dx=A\b。 x0=[ ]39。y(2)。 end 基于牛頓拉夫遜法潮流計(jì)算的 matlab 實(shí)驗(yàn)報(bào)告 function y=F(x) y(1)=3*x(1)cos(x(2)*x(3))。,i) [x dx F(x)] %在屏幕上輸出每次的 x(i),dx(i),F(x(i)) if (max(abs(dx))eps1)amp。 dx=A\b。 x0=[ ]39。 y=3*sin(a)。 for i=1:n y=y+(1)^(i1)*1/(2*i1)。 % 利用公式計(jì)算接地支路的潮流 end end end end 作業(yè) y=pifun(10000) 基于牛頓拉夫遜法潮流計(jì)算的 matlab 實(shí)驗(yàn)報(bào)告 y = clear all tic。 % 利用公式計(jì)算非接地支路的潮流 else if(J==0 ) S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4))。 if (J~=0)amp。 S_res = zeros(nl,5)。 % bus_res矩陣儲(chǔ)存著節(jié)點(diǎn)計(jì)算結(jié)果 bus_res(:,1:2) = bus(:,1:2)。 YtYm(k,4) = Ym+Yt*(1+K)。 end if K0 % 變壓器線路 : Zt和 Ym為折算到 i側(cè)的值 ,K在 j側(cè) YtYm(k,3) = Yt/K。 YtYm(k,4) = Ym。 end Ym=line(k,5)+j*line(k,6)。 for k=1:nl I=line(k,1)。 break end end end end % 恢復(fù) line的編號(hào) 作用為計(jì)算線路的等效 Yt和 Ym，以計(jì)算線路潮流 function YtYm = YtYm(line) [nl,ml]=size(line)。 for i = 1 :nb for j = 1 : nb if bus(j,1) == k bus_temp(k,:) = bus(j,:)。 end % 計(jì)算平衡節(jié)點(diǎn)的無功及有功注入 作用為對(duì)節(jié)點(diǎn)和線路數(shù)據(jù)恢復(fù)編號(hào) function [bus,line] = ReNum(bus,line,nodenum) [nb,mb]=size(bus)。 % 生成雅克比矩陣 作用為計(jì)算每個(gè)節(jié)點(diǎn)的功率注入 function bus = PQ(bus,Y,nPQ,nPV) n = nPQ+nPV+1。 end % 分別計(jì)算 K矩陣的對(duì)角及非對(duì)角元素 if j nPQ+1 if i~=j L(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)))。 end % 分別計(jì)算 H矩陣的對(duì)角及非對(duì)角元素 if j nPQ+1 if i~=j 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)))。 for i = 1:nb1 基于牛頓拉夫遜法潮流計(jì)算的 matlab 實(shí)驗(yàn)報(bào)告 for j = 1:nb Pi(i,1)=Pi(i,1)+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 = zeros(nPQ,nPQ)。 end end end % 利用循環(huán)計(jì)算求取 dP和 dQ 作用為計(jì)算雅克比矩陣 function J = Jac(bus,Y,nPQ) [nb,mb]=size(bus)。 % 總節(jié)點(diǎn)個(gè)數(shù) dP = bus(1:n1,4)。 Y(J,J)=Y(J,J)+K*K*Yt。 Y(J,J)=Y(J,J)+Yt/K/K。 end if (K==0)amp。 Y(I,I)=Y(I,I)+Yt+Ym。 Ym=line(k,5)+j*line(k,6)。 for k=1:nl I=line(k,1)。 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)。SW]。 % increment PQ bus counter PQ(nPQ,:)=bus(i,:)。 % increment swing bus counter 基于牛頓拉夫遜法潮流計(jì)算的 matlab 實(shí)驗(yàn)報(bào)告 SW(nSW,:)=bus(i,:)。 % number of swing bus counter nPV = 0。 end fclose(myf)。,real(S_res(i,4)),imag(S_res(i,4)))。% + j % 39。 fprintf(myf,39。)。% + j %\n39。 fprintf(myf,39。,bus_res(i,1))。牛頓－拉夫遜法潮流計(jì)算結(jié)果 節(jié)點(diǎn)計(jì)算結(jié)果： n節(jié)點(diǎn) 節(jié)點(diǎn)電壓 節(jié)點(diǎn)相角（角度） 節(jié)點(diǎn)注入功率 \n39。,39。 % 計(jì)算線路的等效 Yt和 Ym的子程序，以計(jì)算線路潮流 bus_res = bus_res(bus)。 % 修正相角 if (max(abs(dU))EPS)amp。 dAng = dAngU(1:nb1,1)。 % 計(jì)算雅克比矩陣的子程序 UD = zeros(nPQ,nPQ)。 fclose(myf)。)。 % 對(duì)節(jié)點(diǎn)重新排序的子程序 Y = y(bus,line) % 計(jì)算節(jié)點(diǎn)導(dǎo)納矩陣的子程序 myf = fopen(39。 % strip off .m eval(dfile(1:lfile2))。)。 心得：本程序分模塊進(jìn)行，先是排序，再是求導(dǎo)納陣，然后求雅閣比，再進(jìn)行迭代運(yùn)算，程序本身很簡(jiǎn)潔明了，運(yùn)行的時(shí)候只需要在 matlab 里輸入 main 就行了，然后打開 BUS 和 line 所在的 .m文件，結(jié)果就會(huì)自動(dòng)存在 result 文件中了，通過編寫牛頓拉夫遜法 matlab 潮流計(jì)算程序復(fù)習(xí)了潮流計(jì)算的知識(shí)，也實(shí)現(xiàn)了計(jì)算機(jī)算法 附錄： 實(shí)驗(yàn)源程序： Main函數(shù)： clear [dfile,pathname]=uigetfile(39。 38 29 0 0 。 基于牛頓拉夫遜法潮流計(jì)算的 matlab 實(shí)驗(yàn)報(bào)告 33 19 0 0 。 30 2 0 0 。 4 0 0 0 0 0。 27 26 0 0。 22 21 0 0。 21 16 0 0。 15 14 0 0。 39 9 0 0。 7 6 0 0。 5 4 0 0。 3 2 0 0。 二、程序流程圖 所用公式 112 2 2 2[ ( ) ( ) ][ ( ) ( ) ]()jni i i ij j ij j i ij j ij jjjni i i ij j ij j i ij j ij jji i i iP P e G e B f f G f B eQ Q f G e B f e G f B eU U e f?????? ? ? ? ? ?????? ? ? ? ? ???? ? ? ? ?