%程序的功能是用牛頓——拉夫法進(jìn)行潮流計(jì)算%B1矩：1支路首端號；、末端號；、支路阻抗、支路對地電納%5、路的變比、支路首端處于K側(cè)為，1側(cè)%B2矩：1該節(jié)點(diǎn)發(fā)電機(jī)功率、該節(jié)點(diǎn)負(fù)荷功率；、節(jié)點(diǎn)電壓初始值%4節(jié)電壓V的定值、節(jié)點(diǎn)所接的無功補(bǔ)償設(shè)備的容量%6、點(diǎn)分類標(biāo)號：1為衡點(diǎn)（應(yīng)為節(jié)點(diǎn)為節(jié)；%3為PV點(diǎn)；請入節(jié)點(diǎn)數(shù)：n=');nl=10;%input('請入支路數(shù)：nl=');isb=1;%input('輸入平衡母線節(jié)點(diǎn)號：isb=');pr=0.00001;%input('請輸入誤差精度pr=');B1=[120.03512+0.08306i1301.0238140.05620+0.13289i150.00811+0.24549i1.0238160.05620+0.13289i60.04215+0.09967i701.0238180.02810+0.06645i100.00811+0.24549i01;90.03512+0.08306i0]B2=[001.11.11;102;102;1102;12;102;10.51.11.103;102];%input('請輸各節(jié)點(diǎn)參數(shù)形成的矩陣：B2=');%%%---------------------------------------------------fori=1:nl%路數(shù)if%節(jié)點(diǎn)處于側(cè)左點(diǎn)處于K側(cè)非對角元Y(q,p)=Y(p,q);%非角元Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;對角元K側(cè)對角元側(cè)%導(dǎo)納矩陣導(dǎo)納矩disp(Y)%----------------------------------------------------------分解出導(dǎo)納陣的實(shí)部虛部fori=1:n%定各節(jié)點(diǎn)初始電壓的實(shí)部和虛部(i)=imag(B2(i,3));V(i)=B2(i,4);%PV節(jié)電壓給定模值fori=1:n%定各節(jié)點(diǎn)注入功率%i點(diǎn)注入功率SG-SLB(i,i)=B(i,i)+B2(i,5);%i節(jié)無功補(bǔ)償量%===================================================================分解出各節(jié)點(diǎn)入的有功和無功功率迭代次數(shù)ICT1、；不滿足收斂要求的節(jié)點(diǎn)數(shù)IT2IT2~=0雅可比矩陣的階數(shù)擴(kuò)列IT2=0;a=a+1;fori=1:nifi~=isb%平衡節(jié)點(diǎn)C(i)=0;D(i)=0;forj1=1:n(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%(Gij*ej-Bij*fj)(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)P1=C(i)*e(i)+f(i)*D(i);%點(diǎn)功率P計(jì)(Gij*ej-Bij*fj)+fi(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%節(jié)功率計(jì)fiΣ(Gij*fj+Bij*ej)%i節(jié)有功和無功功率P',Q'的計(jì)算值V2=e(i)^2+f(i)^2;%電壓模平方%=========以下針對非節(jié)來求取功率差及矩陣元素ifB2(i,6)~=3非PV節(jié)點(diǎn)DP=P(i)-P1;節(jié)點(diǎn)有功功率差節(jié)無功功率差以為除平衡節(jié)點(diǎn)外其它節(jié)點(diǎn)的功率計(jì)算%=================求矩===================forj1=1:nif%平衡節(jié)&對角元X3=X2;X4=-X1;X1=dP/dep=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;%節(jié)點(diǎn)無功功率差J(m,q)=X1;J(m,N)=DP;q=q+1;%J(m,N)=DP節(jié)點(diǎn)有功功率差X2=dp/dfelseif非平衡節(jié)點(diǎn)&對角元dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/df%X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%dQ/df擴(kuò)列eq\o\ac(△,Q)m=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴(kuò)展列eq\o\ac(△,P)J(m,q)=X2;下是針對PV節(jié)來求取矩的元素===========DP=P(i)-P1;%PV節(jié)有功誤差%PV點(diǎn)電壓誤差forj1=1:nif%平衡節(jié)&對角元%dP/dedP/dfX5=0;X6=0;%PV點(diǎn)電壓誤差m=p+1;%點(diǎn)有功誤差J(m,q)=X2;elseif非平衡節(jié)點(diǎn)&對角元dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%dP/df%PV點(diǎn)電壓誤差m=p+1;%點(diǎn)有功誤差J(m,q)=X2;%=========以上為求雅可比矩陣的各個(gè)元素?cái)U(kuò)展列功率差或壓差=====================fork=3:N0%N0=2*n（從第三行開始，第一、二行是平衡節(jié)點(diǎn)）%N=N0+1即擴(kuò)列eq\o\ac(△,P)eq\o\ac(△,)eq\o\ac(△,、)Q或eq\o\ac(△,U)eq\o\ac(△,)fork2=k1:N1從列元素到擴(kuò)展列的eq\o\ac(△,P)eq\o\ac(△,)eq\o\ac(△,、)Q或eq\o\ac(△,U)J(k,k2)=J(k,k2)./J(k,k);%用K行K列角元素去除K行K列的非對角元素進(jìn)行規(guī)格化J(k,k)=1;%對角元規(guī)格化KK列角元素賦%====================

回代運(yùn)算=======================================if%不是第三行k3for用行從第三行開始到當(dāng)前行的前一行行消去fork2=k1:N1行后各行上三角元素消運(yùn)算（當(dāng)前k列素消為0）用前行列元素減去當(dāng)前行k列素乘以第k行K2列素%當(dāng)前行第k列素已消為if%已到最后一行break;==================================

前代運(yùn)算fork3=k1:N0%從行2*n最一行fork2=k1:N1從k+1列擴(kuò)展列消去后各行下三角元素消運(yùn)算用前行列元素減去當(dāng)前行k列素以第列素%當(dāng)前行第k列元素已消為0%第三行k=3第三行的前代運(yùn)算========================fork3=k1:N0從第四行到2n（最后一行）fork2=k1:N1從第四列到2n+1（即擴(kuò)展列）消運(yùn)算（當(dāng)前3列素消為）用前行列素減去當(dāng)前行3列元素乘以第三行K2列素%當(dāng)前行第3列素已消為0上面是用線性變換方式高消去法將Jacobi矩化成單位矩====fore(L)=e(L)-J(k,N);%修節(jié)電壓實(shí)部f(L)=f(L)-J(k1,N);修改節(jié)點(diǎn)電壓虛部%------修節(jié)點(diǎn)電壓-----------fork=3:N0if%壓偏差量是否滿足要求IT2=IT2+1;不滿足要求的節(jié)點(diǎn)數(shù)加1滿足要求的節(jié)點(diǎn)數(shù)ICT1=ICT1+1;迭代次%高斯消去法"w=-J*V"迭代次：沒有達(dá)精度要求的個(gè)數(shù)fork=1:nV(k)=sqrt(e(k)^2+f(k)^2);計(jì)算各節(jié)點(diǎn)電壓的模值計(jì)算各節(jié)點(diǎn)電壓的角度%將各節(jié)電壓用復(fù)數(shù)表示計(jì)各輸出量===========================各節(jié)點(diǎn)實(shí)際電壓標(biāo)幺值E為節(jié)號從小到大排)disp(E);顯示各節(jié)點(diǎn)的實(shí)際電壓標(biāo)幺值E用數(shù)表示各節(jié)點(diǎn)電壓大小V為(節(jié)點(diǎn)號從小到大排列)：');disp(V);顯示各節(jié)點(diǎn)的電壓大小V的模值各節(jié)點(diǎn)電壓相角sida為(節(jié)點(diǎn)號從小到大排列：disp(sida);顯示各節(jié)點(diǎn)的電壓相角forp=1:nforq=1:n%算各節(jié)點(diǎn)的注入電流的共軛值S(p)=E(p)*C(p);計(jì)各節(jié)點(diǎn)的功率S=電X注電流的軛值各節(jié)點(diǎn)功率S(節(jié)點(diǎn)號從小到大排列：%顯各節(jié)點(diǎn)的注入功率各條支的首端功率為順序同您輸入B1時(shí)致)：fori=1:nlifSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));disp(Si(p,q));SSi(p,q)=Si(p,q);disp(ZF);各條支的末端功率為順序同您輸入B1時(shí)致)：fori=1:nlifSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);disp(Sj(q,p));SSj(q,p)=Sj(q,p);disp(ZF);各條支的功率損耗DS為順序同您輸入B1時(shí)致：fori=1:nlDS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i);ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);figure(1);subplot(1,2,1);plot(V);xlabel('點(diǎn)號');ylabel('電標(biāo)幺值);gridon;subplot(1,2,2);plot(sida);xlabel('點(diǎn)號');ylabel('電角度'gridon;figure(2);subp