電力系統(tǒng)分析課程設(shè)計(jì)——牛頓拉夫遜潮流計(jì)算

1、電力系統(tǒng)編程潮流計(jì)算1 設(shè)計(jì)任務(wù)及初步分析1.1 設(shè)計(jì)任務(wù)條件：節(jié)點(diǎn)數(shù)：3 支路數(shù)：3 計(jì)算精度：0.00010支路1： 0.0300+j0.0900 12支路2： 0.0200+j0.0900 23支路3： 0.0300+j0.0900 31節(jié)點(diǎn)1：PQ節(jié)點(diǎn)，S（1）節(jié)點(diǎn)2：PQ節(jié)點(diǎn)，S（2）節(jié)點(diǎn)3：平衡節(jié)點(diǎn)，U（3）=1.00000.0000要求：編寫程序計(jì)算潮流1.2 初步分析潮流計(jì)算在數(shù)學(xué)上可歸結(jié)為求解非線性方程組，其數(shù)學(xué)模型簡(jiǎn)寫如下： 2 牛頓-拉夫遜法簡(jiǎn)介2.1概述牛頓拉夫遜法是目前求解非線性方程最好的一種方法。這種方法的特點(diǎn)就是把對(duì)非線性方程的求解過(guò)程變成反復(fù)對(duì)相應(yīng)的線性方程求解

2、的過(guò)程，通常稱為逐次線性化過(guò)程，就是牛頓拉夫遜法的核心。牛頓-拉夫遜法的基本原理是在解的某一鄰域內(nèi)的某一初始點(diǎn)出發(fā)，沿著該點(diǎn)的一階偏導(dǎo)數(shù)雅可比矩陣，朝減小方程的殘差的方向前進(jìn)一步，在新的點(diǎn)上再計(jì)算殘差和雅可矩陣?yán)^續(xù)前進(jìn)，重復(fù)這一過(guò)程直到殘差達(dá)到收斂標(biāo)準(zhǔn)，即得到了非線性方程組的解。因?yàn)樵娇拷?，偏?dǎo)數(shù)的方向越準(zhǔn)，收斂速度也越快，所以牛頓法具有二階收斂特性。而所謂“某一鄰域”是指雅可比方向均指向解的范圍，否則可能走向非線性函數(shù)的其它極值點(diǎn)，一般來(lái)說(shuō)潮流由平電壓即各母線電壓(相角為0，幅值為1)啟動(dòng)即在此鄰域內(nèi)。22 一般概念對(duì)于非線性代數(shù)方程組即 (21)在待求量的某一個(gè)初始計(jì)算值附件，將上式展開(kāi)

3、泰勒級(jí)數(shù)并略去二階及以上的高階項(xiàng)，得到如下的線性化的方程組 (22)上式稱之為牛頓法的修正方程式。由此可以求得第一次迭代的修正量 (23)將和相加，得到變量的第一次改進(jìn)值。接著再?gòu)某霭l(fā)，重復(fù)上述計(jì)算過(guò)程。因此從一定的初值出發(fā)，應(yīng)用牛頓法求解的迭代格式為 (24) (25)上兩式中：是函數(shù)對(duì)于變量的一階偏導(dǎo)數(shù)矩陣，即雅可比矩陣；為迭代次數(shù)。由式（24）和式子（25）可見(jiàn)，牛頓法的核心便是反復(fù)形成求解修正方程式。牛頓法當(dāng)初始估計(jì)值和方程的精確解足夠接近時(shí)，收斂速度非常快，具有平方收斂特性。2.3 潮流計(jì)算的修正方程運(yùn)用牛頓拉夫遜法計(jì)算潮流分布時(shí)，首先要找出描述電力系統(tǒng)的非線性方程。這里仍從節(jié)點(diǎn)電壓

4、方程入手，設(shè)電力系統(tǒng)導(dǎo)納矩陣已知，則系統(tǒng)中某節(jié)點(diǎn)（節(jié)點(diǎn)）電壓方程為從而得 進(jìn)而有 （26）式（26）中，左邊第一項(xiàng)為給定的節(jié)點(diǎn)注入功率，第二項(xiàng)為由節(jié)點(diǎn)電壓求得的節(jié)點(diǎn)注入功率。他們二者之差就是節(jié)點(diǎn)功率的不平衡量。現(xiàn)在有待解決的問(wèn)題就是各節(jié)點(diǎn)功率的不平衡量都趨近于零時(shí)，各節(jié)點(diǎn)電壓應(yīng)具有的價(jià)值。由此可見(jiàn)，如將式（26）作為牛頓拉夫遜中的非線性函數(shù)，其中節(jié)點(diǎn)電壓就相當(dāng)于變量。建立了這種對(duì)應(yīng)關(guān)系，就可列出修正方程式，并迭代求解。但由于節(jié)點(diǎn)電壓可有兩種表示方式以直角做表或者極坐標(biāo)表示，因而列出的迭代方程相應(yīng)地也有兩種，下面分別討論。2.4 直角坐標(biāo)表示的修正方程節(jié)點(diǎn)電壓以直角坐標(biāo)表示時(shí)，令、，且將導(dǎo)納矩陣

5、中元素表示為，則式（27）改變?yōu)?（27）再將實(shí)部和虛部分開(kāi)，可得 （28）這就是直角坐標(biāo)下的功率方程?？梢?jiàn)，一個(gè)節(jié)點(diǎn)列出了有功和無(wú)功兩個(gè)方程。對(duì)于節(jié)點(diǎn)（），給定量為節(jié)點(diǎn)注入功率，記為、，則由式（28）可得功率的不平衡量，作為非線性方程（29）式中、分別表示第節(jié)點(diǎn)的有功功率的不平衡量和無(wú)功功率的不平衡量。對(duì)于節(jié)點(diǎn)（），給定量為節(jié)點(diǎn)注入有功功率及電壓數(shù)值，記為、，因此，可以利用有功功率的不平衡量和電壓的不平衡量表示出非線性方程，即有（210）式中為電壓的不平衡量。對(duì)于平衡節(jié)點(diǎn)（），因?yàn)殡妷簲?shù)值及相位角給定，所以也確定，不需要參加迭代求節(jié)點(diǎn)電壓。因此，對(duì)于個(gè)節(jié)點(diǎn)的系統(tǒng)只能列出個(gè)方程，其中有功功率方

6、程個(gè)，無(wú)功功率方程個(gè)，電壓方程個(gè)。將式（29）、式（210） 非線性方程聯(lián)立，稱為個(gè)節(jié)點(diǎn)系統(tǒng)的非線性方程組，且按泰勒級(jí)數(shù)在、（）展開(kāi)，并略去高次項(xiàng)，得到以矩陣形式表示的修正方程如下。 （211）上式中雅可比矩陣的各個(gè)元素則分別為 將（211）寫成縮寫形式 （212）對(duì)雅可比矩陣各元素可做如下討論：當(dāng)時(shí)，對(duì)于特定的，只有該特定點(diǎn)的和是變量，于是雅可比矩陣中各非對(duì)角元素表示為 當(dāng)時(shí)，雅可比矩陣中各對(duì)角元素的表示式為由上述表達(dá)式可知，直角坐標(biāo)的雅可比矩陣有以下特點(diǎn)：1) 雅可比矩陣是階方陣，由于、等等，所以它是一個(gè)不對(duì)稱的方陣。2) 雅可比矩陣中諸元素是節(jié)點(diǎn)電壓的函數(shù)，在迭代過(guò)程中隨電壓的變化而不斷

7、地改變。3) 雅可比矩陣的非對(duì)角元素與節(jié)點(diǎn)導(dǎo)納矩陣中對(duì)應(yīng)的非對(duì)角元素有關(guān)，當(dāng)中的為零時(shí)，雅可比矩陣中相應(yīng)的、也都為零，因此，雅可比矩陣也是一個(gè)稀疏矩陣。3 程序設(shè)計(jì)3.1 程序流程圖圖3-1 程序流程圖3.2 源程序n=3n1=3isb=3pr=0.0001B1=1 2 0.03+0.09i 0 1 0;1 3 0.03+0.09i 0 1 0;2 3 0.02+0.09i 0 1 0B2=0 -0.5-0.2i 1 0 0 2;0 -0.6-0.25i 1 0 0 2;0 0 1 1 0 1X=1 0;2 0;3 0%X=input('節(jié)點(diǎn)號(hào)和對(duì)地參數(shù):X='); Y=zer

8、os(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)+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)=

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

10、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 %對(duì)PV節(jié)點(diǎn)的處理，注意這時(shí)不可再將h初始化為0 if i=isb&B2(i,6)=3 h=h+1; fo

11、r 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)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i

12、,j)*real(B2(j,3); end end end OrgS %創(chuàng)建PVU 用于存儲(chǔ)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); end end PVU %創(chuàng)建DetaS，用于存儲(chǔ)有功功率、無(wú)功功率和電壓幅值的不平衡量 h=0; for i=1:n %對(duì)PQ節(jié)點(diǎ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)-Or

13、gS(2*h,1); end end t=0; for i=1:n %對(duì)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; end end DetaS %創(chuàng)建I，用于存儲(chǔ)節(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

14、(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); end end I %創(chuàng)建Jacbi(雅可比矩陣) Jacbi=zeros(2*n-2); h=0; k=0; for i=1:n %對(duì)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 %對(duì)角元素的處理 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(

15、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 %非對(duì)角元素的處理 Jacbi(2*h-1,2*k-1)=-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); J

16、acbi(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 end end k=0; for i=1:n %對(duì)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 %對(duì)角元素的處理 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

17、(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 %非對(duì)角元素的處理 Jacbi(2*h-1,2*k-1)=-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

18、,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 end end Jacbi %求解修正方程，獲取節(jié)點(diǎn)電壓的不平衡量 DetaU=zeros(2*n-2,1); DetaU=inv(Jacbi)*DetaS; DetaU %修正節(jié)點(diǎn)電壓 j=0; for i=1:n %對(duì)PQ節(jié)點(diǎn)處理 if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end en

19、d for i=1:n %對(duì)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); end end B2 while abs(max(DetaS)>pr OrgS=zeros(2*n-2,1); %!初始功率參數(shù)在迭代過(guò)程中是不累加的，所以在這里必須將其初始化為零矩陣 h=0; j=0; for i=1:n if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y

20、(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=isb&B2(

21、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)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(

22、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=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); end end t=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,2)-OrgS(2*h-1,1); DetaS(2*h,

23、1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; end end DetaS %創(chuàng)建I 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); end end I %創(chuàng)建Jacbi Jacbi=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 k=k+1; if i=j Jac

24、bi(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)=-imag(Y(i,j)*real

25、(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 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*h-1,2*k-1)

26、=-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); Jacbi(2*h-1,2*k

27、)=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) k=0; end end end end end Jacbi DetaU=zeros(2*n-2,1); DetaU=inv(Jacbi)*DetaS; DetaU %修正節(jié)點(diǎn)電壓 j=0; for i=1:n if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); end end for i=

28、1: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); end end B2 Times=Times+1; %迭代次數(shù)加1 end disp('迭代次數(shù)為');disp(Times)for no=1:n Vn(no)=B2(no,3);str1 = sprintf('節(jié)點(diǎn)%d的電壓是%d',no);disp(str);disp(Vn(no);end3.3 結(jié)果及分析節(jié)點(diǎn)導(dǎo)納矩陣為Y= 6.6667 -j20.0000 -3.3333 +j10.0000 -3.33

29、33 +j10.0000-3.3333 +j10.0000 5.6863 -j20.5882 -2.3529 +j10.5882-3.3333 +j10.0000 -2.3529 +j10.5882 5.6863 -j20.5882迭代運(yùn)行了3次，每次功率的不平衡量表3-1所示表3-1 迭代過(guò)程中各節(jié)點(diǎn)的功率不平衡量kP1(k)+jQ1(k)P2(k)+jQ2(k)01-0.0083-j0.0279-0.0091-j0.03602-0.0000297-j0.000095530.0000+j0.00000.0000+j0.0000表3-2 迭代過(guò)程中各節(jié)點(diǎn)電壓ke1(k)+jf1(k)e2(k)

30、+jf2(k)10.9664 - j0.04230.9658 - j0.045720.9632 - j0.04230.9623 - j0.045730.9632 - j0.04230.9623 - j0.0457平衡節(jié)點(diǎn)功率為P3+jQ3= 1.1188 +j0.5188下面計(jì)算網(wǎng)絡(luò)中的功率分布，根據(jù)所有線路功率計(jì)算如下:S12=(0.9632-j0.0423)(0.9632+j0.0423)×0+(0.9632+j0.0423)-(0.9623+j0.0457)×(3.33+j10)=0.0355 - j0.0042 S13=0.9632-j0.04230.9632+j0.0423*0+0.9632+j0.0423-1×-5.6863-j20.5882=-0.5355 - j0.1958同樣可以得出 S21= -0.0354 + j0.0043 S23=-0.5