電力系統(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.090012支路2: 0.0200+j0.090023支路3: 0.0300+j0.090031節(jié)點(diǎn)1:PQ 節(jié)點(diǎn),S (1=-0.5000-j0.2000 節(jié)點(diǎn)2:PQ 節(jié)點(diǎn),S (2=-0.6000-j0.2500 節(jié)點(diǎn)3:平衡節(jié)點(diǎn),U (3=1.00000.0000要求:編寫程序計(jì)算潮流1.2 初步分析潮流計(jì)算在數(shù)學(xué)上可歸結(jié)為求解非線性方程組,其數(shù)學(xué)模型簡(jiǎn)寫如下:=0(0(0(21212211n n n n x x x f x x x f x x

2、x f ,2 牛頓-拉夫遜法簡(jiǎn)介2.1概述牛頓-拉夫遜法是目前求解非線性方程最好的一種方法。這種方法的特點(diǎn)就是把對(duì)非線性方程的求解過(guò)程變成反復(fù)對(duì)相應(yīng)的線性方程求解的過(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ù)的其它極值

3、點(diǎn),一般來(lái)說(shuō)潮流由平電壓即各母線電壓(相角為0,幅值為1啟動(dòng)即在此鄰域內(nèi)。2.2 一般概念對(duì)于非線性代數(shù)方程組(0=x f即 (0,21=n i x x x f (n i ,2,1= (2-1在待求量x 的某一個(gè)初始計(jì)算值(0x 附件,將上式展開(kāi)泰勒級(jí)數(shù)并略去二階及以上的高階項(xiàng),得到如下的線性化的方程組(000=+x x f xf (2-2上式稱之為牛頓法的修正方程式。由此可以求得第一次迭代的修正量(0100x f x f x-= (2-3將(0x 和(0x 相加,得到變量的第一次改進(jìn)值(1x 。接著再?gòu)?1x 出發(fā),重復(fù)上述計(jì)算過(guò)程。因此從一定的初值(0x 出發(fā),應(yīng)用牛頓法求解的迭代格式為(

4、k k k xf x x f -= (2-4(k k k xxx+=+1 (2-5上兩式中:(x f 是函數(shù)(x f 對(duì)于變量x 的一階偏導(dǎo)數(shù)矩陣,即雅可比矩陣J ;k為迭代次數(shù)。由式(2-4和式子(2-5可見(jiàn),牛頓法的核心便是反復(fù)形成求解修正方程式。牛頓法當(dāng)初始估計(jì)值(0x 和方程的精確解足夠接近時(shí),收斂速度非?？?具有平方收斂特性。2.3 潮流計(jì)算的修正方程運(yùn)用牛頓-拉夫遜法計(jì)算潮流分布時(shí),首先要找出描述電力系統(tǒng)的非線性方程。這里仍從節(jié)點(diǎn)電壓方程入手,設(shè)電力系統(tǒng)導(dǎo)納矩陣已知,則系統(tǒng)中某節(jié)點(diǎn)(i 節(jié)點(diǎn)電壓方程為=* =nj iijijU S UY1從而得=*=nj j ij i i U Y

5、U S 1進(jìn)而有(01=-+*=*j nj ij i i i U Y U jQ P(2-6式(2-6中,左邊第一項(xiàng)為給定的節(jié)點(diǎn)注入功率,第二項(xiàng)為由節(jié)點(diǎn)電壓求得的節(jié)點(diǎn)注入功率。他們二者之差就是節(jié)點(diǎn)功率的不平衡量?，F(xiàn)在有待解決的問(wèn)題就是各節(jié)點(diǎn)功率的不平衡量都趨近于零時(shí),各節(jié)點(diǎn)電壓應(yīng)具有的價(jià)值。由此可見(jiàn),如將式(2-6作為牛頓-拉夫遜中的非線性函數(shù)(0=X F ,其中節(jié)點(diǎn)電壓就相當(dāng)于變量X 。建立了這種對(duì)應(yīng)關(guān)系,就可列出修正方程式,并迭代求解。但由于節(jié)點(diǎn)電壓可有兩種表示方式以直角做表或者極坐標(biāo)表示,因而列出的迭代方程相應(yīng)地也有兩種,下面分別討論。2.4 直角坐標(biāo)表示的修正方程節(jié)點(diǎn)電壓以直角坐標(biāo)表示時(shí)

6、,令i i i jf e U +=、j j j jf e U +=,且將導(dǎo)納矩陣中元素表示為ij ij ij jB G Y +=,則式(2-7改變?yōu)?(01=-+-+=nj j j ij ij i i i i jf e jB G jf e jQ P(2-7再將實(shí)部和虛部分開(kāi),可得(=+-=+-=0011nj j ij j ij i j ij j ij i i nj j ij j ij i j ij j ij i i e B f G e f B e G f Q e B f G f f B e Ge P(2-8這就是直角坐標(biāo)下的功率方程。可見(jiàn),一個(gè)節(jié)點(diǎn)列出了有功和無(wú)功兩個(gè)方程。對(duì)于PQ 節(jié)點(diǎn)(1,

7、21-=m i ,給定量為節(jié)點(diǎn)注入功率,記為i P 、i Q ,則由式(2-8可得功率的不平衡量,作為非線性方程(+-=+-=nj j ij j ij i j ij j ij i i i nj j ij j ij i j ij j ij i i i e B f G e f B e G f Q Q e B f G f f B e Ge P P 11(2-9式中i P 、i Q 分別表示第i 節(jié)點(diǎn)的有功功率的不平衡量和無(wú)功功率的不平衡量。對(duì)于PV 節(jié)點(diǎn)(n m m i ,2,1 +=,給定量為節(jié)點(diǎn)注入有功功率及電壓數(shù)值,記為i P 、i U ,因此,可以利用有功功率的不平衡量和電壓的不平衡量表示出

8、非線性方程,即有(+-=+-=22221ii i inj j ij j ij i j ij j ij i i i f e U U e B f G f f B e Ge P P (2-10式中i U 為電壓的不平衡量。對(duì)于平衡節(jié)點(diǎn)(m i =,因?yàn)殡妷簲?shù)值及相位角給定,所以S s S jf e U +=也確定,不需要參加迭代求節(jié)點(diǎn)電壓。因此,對(duì)于n 個(gè)節(jié)點(diǎn)的系統(tǒng)只能列出(12-n 個(gè)方程,其中有功功率方程(1-n 個(gè),無(wú)功功率方程(1-m 個(gè),電壓方程(m n -個(gè)。將式(2-9、式(2-10 非線性方程聯(lián)立,稱為n 個(gè)節(jié)點(diǎn)系統(tǒng)的非線性方程組,且按泰勒級(jí)數(shù)在(0i f 、(0i e (m i n

9、 i =,2,1 展開(kāi),并略去高次項(xiàng),得到以矩陣形式表示的修正方程如下。=n n p p nn nnnpnpn n n n nn nn np np n n n n pn pn pp pp p p p p pn pnpp pp p p p p n n p p n n p p n n p p n n p p n n pp e f e f e f e f S R S R S R S R N H N H N H N H S R S R S R S RN HN HN HN H L J L J L J L J N H N H N H N H L J L J L J L J N H N H N H N H

10、U P U P Q P Q P 22112211221122211 (2-11 上式中雅可比矩陣的各個(gè)元素則分別為j i ij f P H =j i ij e P N =ji ij f Q J =ji ij e Q L =ji ij f U R =2ji ij e U S =2將(2-11寫成縮寫形式=e f e f S RL J N H U Q PJ 2 (2-12對(duì)雅可比矩陣各元素可做如下討論:當(dāng)i j 時(shí),對(duì)于特定的j ,只有該特定點(diǎn)的i f 和i e 是變量,于是雅可比矩陣中各非對(duì)角元素表示為i ij i ij j i ij f G e B f P H -=i ij i ij j i

11、ij f B e G e P N -=i ij i ij ji ij e G f B f Q J +=i ij i ij ji ij e B f G e Q L +-=02=ji ij f U R02=ji ij e U S當(dāng)i j =時(shí),雅可比矩陣中各對(duì)角元素的表示式為(iii i ii nj j ij j ij j i ij e B f G e B f G f P H +-+-=1(i ii i ii nj j ij j ij j i ij f B e G f B e G e P N -=1(fB e G f B e G f Q J ii i ii nj j ij j ij j i ij

12、+-=1(i ii i ii nj j ij j ij ji ij e B f G e B f G e Q L +-+=1i j i ij f f U R 22-=i ji ij e e U S 22-=由上述表達(dá)式可知,直角坐標(biāo)的雅可比矩陣有以下特點(diǎn):1 雅可比矩陣是(12-n 階方陣,由于ji ij H H 、ji ij N N 等等,所以它是一個(gè)不對(duì)稱的方陣。2 雅可比矩陣中諸元素是節(jié)點(diǎn)電壓的函數(shù),在迭代過(guò)程中隨電壓的變化而不斷地改變。3 雅可比矩陣的非對(duì)角元素與節(jié)點(diǎn)導(dǎo)納矩陣B Y 中對(duì)應(yīng)的非對(duì)角元素有關(guān),當(dāng)B Y 中的ij Y 為零時(shí),雅可比矩陣中相應(yīng)的ij H 、ij N 、ij J

13、 、ij L 也都為零,因此,雅可比矩陣也是一個(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=zeros(n;Times=1; %置迭代次數(shù)為初始值%創(chuàng)建節(jié)點(diǎn)導(dǎo)納矩陣for i=1:n1if B1(i,6=0 %不含

14、變壓器的支路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;Y(q,q=Y(q,q+1/(B1(i,52*B1(i,3;endendYOrgS=zeros(2*n-2,1;DetaS=zeros(2*n-2,1; %將OrgS、D

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

16、m ag(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:n %對(duì)PV節(jié)點(diǎn)的處理,注意這時(shí)不可再將h初始化為0if 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*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*

17、h,1=OrgS(2*h,1+imag(B2(i,3*(real(Y(i,j*real(B2(j,3-imag(Y(i,j*im ag(B2(j,3-real(B2(i,3*(real(Y(i,j*imag(B2(j,3+imag(Y(i,j*real(B2(j,3; endendendOrgS%創(chuàng)建PVU 用于存儲(chǔ)PV節(jié)點(diǎn)的初始電壓PVU=zeros(n-h-1,1;t=0;for i=1:nif B2(i,6=3t=t+1;PVU(t,1=B2(i,3;endendPVU%創(chuàng)建DetaS,用于存儲(chǔ)有功功率、無(wú)功功率和電壓幅值的不平衡量h=0;for i=1:n %對(duì)PQ節(jié)點(diǎn)的處理if i=

18、isb&B2(i,6=2h=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 %對(duì)PV節(jié)點(diǎn)的處理,注意這時(shí)不可再將h初始化為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,12+imag(PVU(t,12-real(B2(i,32-imag(B2(i,32; endendDetaS%創(chuàng)建I,用于存儲(chǔ)節(jié)點(diǎn)電流參數(shù)i=z

19、eros(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%創(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=2h=h+1;for j=1:nif j=isbk=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

20、*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 %非對(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;Jacbi(2

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

22、al(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,3;Jacbi(2*h,2*k-1=0;Jacbi(2*h,2*k=0;endif k=(n-1 %將用于內(nèi)循環(huán)的指

23、針置于初始值,以確保雅可比矩陣換行k=0;endendendendendJacbi%求解修正方程,獲取節(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=2j=j+1;B2(i,3=B2(i,3+DetaU(2*j,1+DetaU(2*j-1,1*sqrt(-1;endendfor i=1:n %對(duì)PV節(jié)點(diǎn)的處理if B2(i,6=3j=j+1;B2(i,3=B2(i,3+DetaU(2*j,1+DetaU(2*j-1,1*sqrt(-1;endendB2w

24、hile abs(max(DetaSprOrgS=zeros(2*n-2,1; %!初始功率參數(shù)在迭代過(guò)程中是不累加的,所以在這里必須將其初始化為零矩陣h=0;j=0;for i=1:nif 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*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

25、(i,3*(real(Y(i,j*real(B2(j,3-imag(Y(i,j*im ag(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*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

26、(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*im ag(B2(j,3-real(B2(i,3*(real(Y(i,j*imag(B2(j,3+imag(Y(i,j*real(B2(j,3; endendendOrgS%創(chuàng)建DetaSh=0;for i=1:nif i=isb&B2(i,6=2h=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:

27、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,12+imag(PVU(t,12-real(B2(i,32-imag(B2(i,32; endendDetaSi=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%創(chuàng)建JacbiJacbi=zeros(2*n-2;h=0;k=0;for i=1:nif B2(i,6=2

28、h=h+1;for j=1:nif j=isbk=k+1;if i=jJacbi(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;elseJacbi(2*h-1,2*k-1=-

29、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;endif k=(n-1k=0;endendendendendk=0;for i=1:nif B2(i,6=3h=h+1;for j=1:nk=k+1;if i=jJacbi(2*h-1,2*k-1=-imag(Y(i,j*real(B2(i,3+rea

30、l(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+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

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

32、%迭代次數(shù)加1enddisp(迭代次數(shù)為;disp(Timesfor no=1:nVn(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.3333 +j10.00003.3333 +j10.0000 5.6863 j20.5882 2.3529 +j10.58823.3333 +j10.0000 2.3529 +j10.5882 5.6863 j20.5882迭代運(yùn)行了3次,每次功率的不平衡量表3-1所示表3-1 迭代過(guò)程中各節(jié)點(diǎn)的功率不平衡量 表3-2 迭代過(guò)程中各節(jié)點(diǎn)電壓 P 3+jQ 3= 1.1188 +j0.5188下面計(jì)算網(wǎng)絡(luò)中的功率分布,根據(jù)(*0*.*.ij j i i i i ij i ij ij ij y V V y V V I V jQ P S -+=+=所有線路功率計(jì)算如下:S 12=(0.9632-j0.0423(0.9632+j0.04230+(0.9632+j0.0423-(0.9623+ j0.0457(3.33+j10=0.0355 - j 0