電力系統(tǒng)分析潮流計(jì)算的計(jì)算機(jī)算法

1、潮流計(jì)算的計(jì)算機(jī)算法實(shí)驗(yàn)報(bào)告姓名： 學(xué)號： 班級： 一、 實(shí)驗(yàn)?zāi)康恼莆粘绷饔?jì)算的計(jì)算機(jī)算法。熟悉MATLAB，并掌握MATLAB程序的基本調(diào)試方法。二、 實(shí)驗(yàn)準(zhǔn)備根據(jù)課程內(nèi)容，熟悉MATLAB軟件的使用方法，自行學(xué)習(xí)MATLAB程序的基礎(chǔ)語法，并根據(jù)所學(xué)知識編寫潮流計(jì)算牛頓拉夫遜法(或PQ分解法) 的計(jì)算程序，用相應(yīng)的算例在MATLAB上進(jìn)行計(jì)算、調(diào)試和驗(yàn)證。三、 實(shí)驗(yàn)要求每人一組，在實(shí)驗(yàn)課時(shí)內(nèi)，用MATLAB調(diào)試和修改運(yùn)行程序，用算例計(jì)算輸出潮流結(jié)果。四、 實(shí)驗(yàn)程序clear; %清空內(nèi)存n=input(請輸入節(jié)點(diǎn)數(shù):n=);n1=input(請輸入支路數(shù):n1=);isb=input(請輸

2、入平衡節(jié)點(diǎn)號:isb=);pr=input(請輸入誤差精度:pr=);B1=input(請輸入支路參數(shù):B1=);B2=input(請輸入節(jié)點(diǎn)參數(shù):B2=);X=input(節(jié)點(diǎn)號和對地參數(shù):X=);Y=zeros(n); Times=1; %一：創(chuàng)建節(jié)點(diǎn)導(dǎo)納矩陣for i=1:n1 if 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)

3、; 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); Y(q,q)=Y(q,q)+1/(B1(i,5)2*B1(i,3); endendY; %將OrgS、DetaS初始化 OrgS=zeros(2*n-2,1); DetaS=zeros(2*n-2,1); %二：創(chuàng)建OrgS，用于存儲初始功率參數(shù)h=0;j=0;for i=1:n %對PQ節(jié)點(diǎn)的處理 if i=isb&B2(i,6)=2 h=h+1; for j=1:n

4、OrgS(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)*I mag(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)*rea

5、l(B2(j,3); end endend%三：對PV節(jié)點(diǎn)的處理，注意這時(shí)不可再將h初始化為0for i=1:n if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(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,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*

6、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); end endendOrgS;%四：創(chuàng)建PVU 用于存儲PV節(jié)點(diǎn)的初始電壓PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)=3 t=t+1; PVU(t,1)=B2(i,3); endendPVU;%五：創(chuàng)建DetaS,用于存儲有功功率、無功功率和電壓幅值的不平衡量h=0;for i=1:n %對PQ節(jié)點(diǎn)的處理 if i=isb&B2(i,6)=2 h=

7、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節(jié)點(diǎn)的處理，注意這時(shí)不可再將h初始化為0 if i=isb&B2(i,6)=3 h=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)2+imag(PVU(t,1)2-real(B2(i,3)2-imag( B2(i,3)2; endendDetaS;%七：創(chuàng)建

8、I，用于存儲節(jié)點(diǎn)電流參數(shù)i=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endendI;%八：創(chuàng)建Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %對PQ節(jié)點(diǎn)的處理 if 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)+real(Y(

9、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)=-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); else %非對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacb

10、i(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); end if k=(n-1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行 k=0; end end end endendk=0;for i=1:n %對PV節(jié)點(diǎn)的處理 if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j %對角元素的處理 Jacbi(2*h-1,2*

11、k-1)=-imag(Y(i,j)*real(B2(i,3)+real(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); else %非對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); J

12、acbi(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; end if k=(n-1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行 k=0; end end end endendJacbi;%九：求解修正方程，獲取節(jié)點(diǎn)電壓的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU;%修正節(jié)點(diǎn)電壓j=0;for i=1:n %對PQ節(jié)點(diǎn)處理 if B2(i,6)=2 j=j+1; B2(

13、i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n %對PV節(jié)點(diǎn)的處理 if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendB2;%十：開始循環(huán)*while abs(max(DetaU)prOrgS=zeros(2*n-2,1); %!初始功率參數(shù)在迭代過程中是不累加的，所以在這里必須將其初始化為零矩陣h=0;j=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; for j=1:n

14、 OrgS(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,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)*rea

15、l(B2(j,3); end endendfor i=1:n if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(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,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)*i

16、mag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendOrgS;%創(chuàng)建DetaSh=0;for i=1:n 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 if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,

17、2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; endendDetaS;%創(chuàng)建Ii=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endendI;%創(chuàng)建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)=2 h=h+1; for j=1:n if j=isb

18、k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(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)=-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); else Jacbi(2*h-1,2*k-1)=

19、-imag(Y(i,j)*real(B2(i,3)+real(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); end if k=(n-1) k=0; end end end endendk=0;for i=1:n if B2(i,6)=3 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)+real(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); else Jacbi(2*h-1,2*k-1)=-imag(