1、%簡單潮流計算的牛頓拉夫遜程序，相關的原始數據數據數據輸入格式如下： %B1是支路參數矩陣，第一列和第二列是節點編號。節點編號由小到大編寫 %對于含有變壓器的支路，第一列為低壓側節點編號，第二列為高壓側節點 %編號，將變壓器的串聯阻抗置于低壓側處理。%第三列為支路的串列阻抗參數。%第四列為支路的對地導納參數。%第五列為含變壓器支路的變壓器的變比 %第六列為變壓器是否含有變壓器的參數，其中“ 1”為含有變壓器， %“0”為不含有變壓器。%B2為節點參數矩陣，其中第一列為節點注入發電功率參數；第二列為節點負荷 功率參數；第三列為節點電壓參數；第六列為節點類型參數，其中“ 1”為平衡 節點，“2”為

2、PQ節點，“3”為PV節點參數。%X為節點號和對地參數矩陣。其中第一列為節點編號，第二列為節點對地 %參數。n=input(' 請輸入節點數 :n='); n1=input(' 請輸入支路數 :n1='); isb=input(' 請輸入平衡節點號 :isb='); pr=input(' 請輸入誤差精度 :pr=');B1=input(' 請輸入支路參數 :B1=');B2=input(' 請輸入節點參數 :B2=');X=input(' 節點號和對地參數 :X=');Y=zeros

3、(n);Times=1;%置迭代次數為初始值%創建節點導納矩陣for i=1:n1if B1(i,6)=0%不含變壓器的支路p=B1(i,1);q=B1(i,2); Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else%含有變壓器的支路p=B1(i,1);q=B1(i,2); Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5);Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3);

4、丫(q,q)=Y(q,q)+1/(B1(i,5F2*B1(i,3);endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %將 OrgS、DetaS 初始化%創建OrgS,用于存儲初始功率參數h=0;j=0;for i=1:n%對PQ節點的處理if i=isb&B2(i,6)=2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*r eal(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*

5、 imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real( B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag (B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendh 初始化為 0endfor i=1:n %對PV節點的處理，注意這時不可再將if i=isb&B2(i,6)=3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,

6、1)+real(B2(i,3)*(real(Y(i,j)*r eal(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)* imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag (B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendOrgS%

7、創建PVU用于存儲PV節點的初始電壓PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)=3t=t+1;PVU(t,1)=B2(i,3);endend PVU %創建DetaS,用于存儲有功功率、無功功率和電壓幅值的不平衡量h=0;for i=1:n %對PQ節點的處理if i=isb&B2(i,6)=2 h=h+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);endendt=0;for i=1:n%對PV節點的處理，注意這時不可再將 h

8、初始化為0if i=isb&B2(i,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1)F2+imag(PVU(t,1)F2-real(B2(i,3)F 2-imag(B2(i,3)F2;endendDetaS%創建 I ，用于存儲節點電流參數i=zeros(n-1,1);h=0;for i=1:nif i=isbh=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3);endendI%創建 Jacbi

9、( 雅可比矩陣 )Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n % 對PQ節點的處理if B2(i,6)=2h=h+1;for j=1:n if j=isb k=k+1;if i=j % 對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=-Jac

10、bi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1);else % 非對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1

11、);endif k=(n-1) % 將用于內循環的指針置于初始值，以確保雅可 比矩陣換行k=0;endendendend end k=0;for i=1:n % 對PV節點的處理if B2(i,6)=3h=h+1;for j=1:n if j=isbk=k+1;if i=j %對角元素的處理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2(i,3)+rea

12、l(I(h,1);Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3);else %非對角元素的處理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0;endif k=(n-1) % 將用于內循環的指針置于初始值，以確保 雅可比

13、矩陣換行k=0;endendendendendJacbi%求解修正方程，獲取節點電壓的不平衡量 DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正節點電壓j=0;for i=1:n% 對PQ節點處理if B2(i,6)=2j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n% 對PV節點的處理if B2(i,6)=3j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);enden

14、dB2%開始循環*while abs(max(DetaU)>prOrgS=zeros(2*n-2,1); %!初始功率參數在迭代過程中是不累加的，以在這里必須將其初始化為零矩陣h=0;j=0;for i=1:nif i=isb&B2(i,6)=2 h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)* imag(B2(j,3)+imag(Y(i,j)*real(B2(j,

15、3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag (B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendfor i=1:nif i=isb&B2(i,6)=3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*r eal(B2(j,3)-imag(Y(i,j)*imag(B2(j

16、,3)+imag(B2(i,3)*(real(Y(i,j)* imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag (B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendOrgS%創建 DetaSh=0;for i=1:nif i=isb&B2(i,6)=2h=h+1;DetaS(2*h-1,1)=

17、real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);endendt=0;for i=1:nif i=isb&B2(i,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1)F2+imag(PVU(t,1)F2-real(B2(i,3)F2-imag(B2(i,3)F2;endendDetaS%創建 Ii=zeros(n-1,1);h=0;for i=1:nif i=isbh=h+1; I(h,1)

18、=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3);endendI%創建 JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2

19、(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(

20、2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k=(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)=3 h=h+1;for j=1:n if j=isb k=k+1; if i=jJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3)+imag(I(h,1);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag( Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=2*imag(B2(i,3); Jacbi(2*h,2*k)=2*real(B2(i,3); elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3