【文章內(nèi)容簡介】 折算到K所在側(cè)的值,K在i側(cè) 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。 %最大迭代次數(shù)eps1=。 %收斂精度eps2=。 fprintf(myf, 39。節(jié)點導納矩陣Y\n39。)。 %輸出節(jié)點導納矩陣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]。 %從bus矩陣中得到PQ PV節(jié)點的相角、PQ節(jié)點的電壓形成的列向量作為x Jac=form_jac(bus,Y)。 %得到雅可比矩陣 del=dPQ(Y,bus)。 %得到功率偏差列向量 dx=Jac\del。 %求得dx fprintf(myf,39。第%d次迭代結(jié)果\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))。 %輸出每次迭代的節(jié)點相角和電壓的偏差 fprintf(myf, 39。\n39。)。 end for I=nb:nb+nPQ1 dx(I,1)=dx(I,1)*x(I,1)。 %由于求得的是delU/U，故還需作相應處理 end x=xdx。 %修正 fprintf(myf,39。第%d次迭代的節(jié)點相角delta和電壓U\n39。,i)。 angl=x(1:nb1,1)。 %將合一起的x，拆成相角、電壓兩個列向量 u=x(nb:nb+nPQ1,1)。 for I=1:nb1 %輸出修正的bus矩陣中PQ、PV結(jié)點的相角 bus(I,3)=angl(I,1)。 fprintf(myf, 39。ang%d %10f \n39。, I,angl(I,1))。 end for I=1:nPQ %輸出修正的bus矩陣中PQ結(jié)點的電壓 bus(I,2)=u(I,1)。 fprintf(myf, 39。U%d %10f \n39。, I,u(I,1))。 end if (max(abs(dx))eps1)amp。(max(abs(dPQ(Y,bus)))eps2) %判斷是否達到所指定的精度要求 break endendif i==max1 %超出最大迭代次數(shù)。出錯指示 error(39。超過最大迭代次數(shù)，不收斂停機！39。)。end%*******************計算節(jié)點注入功率************function PQ global nb。global nPQ。global bus。global line。global Y。 for I=nPQ+1:nb %求PV結(jié)點的Q，平衡節(jié)點的P,Q if bus(I,6)==3 %對于平衡節(jié)點，求其P sum=0。 for J=1:nb ang=bus(I,3)bus(J,3)。 A=real(Y(I,J))*cos(ang)+imag(Y(I,J))*sin(ang)。 sum=sum+bus(I,2)*bus(J,2)*(A)。 end bus(I,4)=sum。 %求取的平衡節(jié)點的P存回bus矩陣 end sum=0。 for J=1:nb ang=bus(I,3)bus(J,3)。 B=real(Y(I,J))*sin(ang)imag(Y(I,J))*cos(ang)。 sum=sum+bus(I,2)*bus(J,2)*B。 end bus(I,5)=sum。 %求取的PV、平衡節(jié)點的Q存回bus矩陣end%***************計算線路功率及損耗******************function flow global nl。global lPQ。global bus。global line。 for k=1:nl I=line(k,1)。 %讀入線路參數(shù) J=line(k,2)。 lPQ(k,1)=I。 %用lPQ矩陣第一、二列保存線路兩端節(jié)點編號 lPQ(k,2)=J。 if J~=0 Zt=line(k,3)+j*line(k,4)。 Yt=1/Zt。 end Ym=line(k,5)+j*line(k,6)。 K=line(k,7)。 Ui=bus(I,2)*(cos(bus(I,3))+j*sin(bus(I,3)))。 %用復數(shù)表示Ui if J~=0 Uj=bus(J,2)*(cos(bus(J,3))+j*sin(bus(J,3)))。 %非接地節(jié)點，則表示Uj end if (K==0)amp。(J~=0) %普通線路非接地支路。 Iij=Ui*(Yt+Ym)Uj*Yt。 Iji=Uj*(Yt+Ym)Ui*Yt。 end if (K==0)amp。(J==0) %接地支路 Iij=Ui*Ym。 Iji=0。 end if K0 %變壓器線路: Zt和Ym為折算到i側(cè)的值,K在j側(cè) Iij=(UiUj)*Yt/K+Ui*(Ym+Yt*(K1)/K)。 Iji=(UjUi)*Yt/K+Uj*Yt*(1K)/K^2。 end if K0 %變壓器線路: Zt和Ym為折算到K側(cè)的值,K在i側(cè) K=K。 Iij=(UiUj)*Yt*K+Ui*(Ym+Yt*(1K))。 Iji=(UjUi)*Yt*K+Uj*Yt*K*(K1)。 end Sij=Ui*conj(Iij)。 %求取節(jié)點注入功率 Sji=Uj*conj(Iji)。delS=Sij+Sji。 %求線路損耗 lPQ(k,3:5)=[Sij Sji delS]。 %保存到lPQ矩陣end%******************恢