【正文】 第一篇：基于MATLAB的電力系統(tǒng)潮流計算基于MATLAB的電力系統(tǒng)潮流計算%簡單潮流計算的小程序，相關(guān)的原始數(shù)據(jù)數(shù)據(jù)數(shù)據(jù)輸入格式如下：%B1是支路參數(shù)矩陣，第一列和第二列是節(jié)點(diǎn)編號。節(jié)點(diǎn)編號由小到大編寫 %對于含有變壓器的支路，第一列為低壓側(cè)節(jié)點(diǎn)編號，第二列為高壓側(cè)節(jié)點(diǎn) %編號，將變壓器的串聯(lián)阻抗置于低壓側(cè)處理。%第三列為支路的串列阻抗參數(shù)。%第四列為支路的對地導(dǎo)納參數(shù)。%第五烈為含變壓器支路的變壓器的變比%第六列為變壓器是否是否含有變壓器的參數(shù)，其中“1”為含有變壓器，%“0”為不含有變壓器。%B2為節(jié)點(diǎn)參數(shù)矩陣，其中第一列為節(jié)點(diǎn)注入發(fā)電功率參數(shù)；第二列為節(jié)點(diǎn) %負(fù)荷功率參數(shù)；第三列為節(jié)點(diǎn)電壓參數(shù)；第六列為節(jié)點(diǎn)類型參數(shù)，其中 %“1”為平衡節(jié)點(diǎn)，“2”為PQ節(jié)點(diǎn)，“3”為PV節(jié)點(diǎn)參數(shù)。%X為節(jié)點(diǎn)號和對地參數(shù)矩陣。其中第一列為節(jié)點(diǎn)編號，第二列為節(jié)點(diǎn)對地 %參數(shù)。n=input(39。請輸入節(jié)點(diǎn)數(shù):n=39。)。n1=input(39。請輸入支路數(shù):n1=39。)。isb=input(39。請輸入平衡節(jié)點(diǎn)號:isb=39。)。pr=input(39。請輸入誤差精度:pr=39。)。B1=input(39。請輸入支路參數(shù):B1=39。)。B2=input(39。請輸入節(jié)點(diǎn)參數(shù):B2=39。)。X=input(39。節(jié)點(diǎn)號和對地參數(shù):X=39。)。Y=zeros(n)。Times=1。%置迭代次數(shù)為初始值 %創(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)+*B1(i,4)。Y(q,q)=Y(q,q)+1/B1(i,3)+*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)。Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3))。end end Y OrgS=zeros(2*n2,1)。DetaS=zeros(2*n2,1)。%將OrgS、DetaS初始化%創(chuàng)建OrgS，用于存儲初始功率參數(shù) h=0。j=0。for i=1:n %對PQ節(jié)點(diǎn)的處理 if i~=isbamp。B2(i,6)==2 h=h+1。for j=1:n OrgS(2*h1,1)=OrgS(2*h1,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))*real(B2(j,3)))。end end end for i=1:n %對PV節(jié)點(diǎn)的處理，注意這時不可再將h初始化為0 if i~=isbamp。B2(i,6)==3 h=h+1。for j=1:n OrgS(2*h1,1)=OrgS(2*h1,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))*real(B2(j,3)))。end end end OrgS %創(chuàng)建PVU 用于存儲PV節(jié)點(diǎn)的初始電壓 PVU=zeros(nh1,1)。t=0。for i=1:n if B2(i,6)==3 t=t+1。PVU(t,1)=B2(i,3)。end end PVU %創(chuàng)建DetaS，用于存儲有功功率、無功功率和電壓幅值的不平衡量 h=0。for i=1:n %對PQ節(jié)點(diǎn)的處理 if i~=isbamp。B2(i,6)==2 h=h+1。DetaS(2*h1,1)=real(B2(i,2))OrgS(2*h1,1)。DetaS(2*h,1)=imag(B2(i,2))OrgS(2*h,1)。end end t=0。for i=1:n %對PV節(jié)點(diǎn)的處理，注意這時不可再將h初始化為0 if i~=isbamp。B2(i,6)==3 h=h+1。t=t+1。DetaS(2*h1,1)=real(B2(i,2))OrgS(2*h1,1)。DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2real(B2(i,3))^2imag(B2(i,3))^2。end end DetaS %創(chuàng)建I，用于存儲節(jié)點(diǎn)電流參數(shù) i=zeros(n1,1)。h=0。for i=1:n if i~=isb h=h+1。I(h,1)=(OrgS(2*h1,1)OrgS(2*h,1)*sqrt(1))/conj(B2(i,3))。end end I %創(chuàng)建Jacbi(雅可比矩陣)Jacbi=zeros(2*n2)。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*h1,2*k1)=imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))。Jacbi(2*h1,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*k1)=Jacbi(2*h1,2*k)+2*real(I(h,1))。Jacbi(2*h,2*k)=Jacbi(2*h1,2*k1)2*imag(I(h,1))。else %非對角元素的處理Jacbi(2*h1,2*k1)=imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))。Jacbi(2*h1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))。Jacbi(2*h,2*k1)=Jacbi(2*h1,2*k)。Jacbi(2*h,2*k)=Jacbi(2*h1,2*k1)。end if k==(n1)%將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行k=0。end end end end end k=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*h1,2*k1)=imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))。Jacbi(2*h1,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*k1)=2*imag(B2(i,3))。Jacbi(2*h,2*k)=2*real(B2(i,3))。else %非對角元素的處理Jacbi(2*h1,2*k1)=imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))。Jacbi(2*h1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))。Jacbi(2*h,2*k1)=0。Jacbi(2*h,2*k)=0。end if k==(n1)%將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行k=0。end end end end end Jacbi %求解修正方程，獲取節(jié)點(diǎn)電壓的不平衡量 DetaU=zeros(2*n2,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(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j1,1)*sqrt(1)。end end for 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*j1,1)*sqrt(1)。end end B2 %開始循環(huán)********************************************************************** while abs(max(DetaU))pr OrgS=zeros(2*n2,1)。%！!初始功率參數(shù)在迭代過程中是不累加的，所以在這里必須將其初始化為零矩陣 h=0。j=0。for i=1:n if i~=isbamp。B2(i,6)==2 h=h+1。for j=1:n OrgS(2*h1,1)=OrgS(2*h1,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))*real(B2(j,3)))。end end end for i=1:n if i~=isbamp。B2(i,6)==3 h=h+1。for j=1:n OrgS(2*h1,1)=OrgS(2*h1,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))*real(B2(j,3)))。end end end OrgS %創(chuàng)建DetaS h=0。for i=1:n if i~=isbamp。B2(i,6)==2 h=h+1。DetaS(2*h1,1)=real(B2(i,2))OrgS(2*h1,1)。DetaS(2*h,1)=imag(B2(i,2))OrgS(2*h,1)。end end t=0。for i=1:n if i~=isbamp。B2(i,6)==3 h=h+1。t=t+1。DetaS(2*h1,1)=real(B2(i,2))OrgS(2*h1,1)。DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2real(B2(i,3))^2imag(B2(i,3))^2。end end DetaS %創(chuàng)建I i=zeros(n1,1)。h=0。for i=1:n if i~=isb h=h+1。I(h,1)=(OrgS(2*h1,1)OrgS(2*h,1)*sqrt(1))/conj(B2(i,3))。end end I %創(chuàng)建Jacbi Jacbi=zeros(2*n2)。h=0。k=0。for i=1:n if B2(i,6)==2 h=h+1。for j=1:n if j~=isb k=k+1。if i==j Jacbi(2*h1,2*k1)=imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))。Jacbi(2*h1,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*k1)=Jacbi(2*h1,2*k)+2*real(I(h,1))。Jacbi(2*h,2*k)=Jacbi(2*h1,2*k1)2*imag(I(h,1))。else Jacbi(2*h1,2*k1)=imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))。Jacbi(2*h1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))。Jacbi(2*h,2*k1)=Jacbi(2*h1,2*k)。Jacbi(2*h,2*k)=Jacbi(2*h1,2*k1)。end if k==(n1)k=0。end end end end end k=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*h1,2*k1)=imag(Y(i,j))*real(