matlab實(shí)驗(yàn) 電力系統(tǒng)潮流計(jì)算

《matlab實(shí)驗(yàn) 電力系統(tǒng)潮流計(jì)算》由會(huì)員分享，可在線閱讀，更多相關(guān)《matlab實(shí)驗(yàn) 電力系統(tǒng)潮流計(jì)算（9頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

1、實(shí)驗(yàn)一電力系統(tǒng)潮流計(jì)算 一、一元非線性方程求解 例1-1試求非線性方程f(x)=0的解。 解：(1)取一個(gè)合理的初值x(o)作為方程f(x)=0的解，如果正好f(x(0))二0，則方程的解x*=x(0)。否則做下一步。 (2)取x(0)+Ax(o)為第一次修正值。Ax(0)充分小，將f(x(0)+Ax⑼)在x(0)附近展開 成泰勒級(jí)數(shù)，并且將的高次項(xiàng)略去，取其線性部分，得到 f(x(0)+Ax⑼)af(x(0))+f'(x(0))Ax(0)=0(1-1) 上式表明，在x(0)處把非線性方程f(x)=0線性化，變成求x(0)附近修正量Ax(0)的線性方 程，這個(gè)方程也稱為修正方程

2、式。從而可求得 f(x(0)) Ax(0)二一(1-2) f(x(0)) 所以，可以確定第一次修正值x⑴=x(O)+Ax(0)。若f(x⑴)二0,則x*=x⑴。 (3) 若f(x⑴)豐0，則用步驟(2)闡述的方法由x(1)確定出第二次修正值x⑵。如此迭代下去，在第(k+1)次迭代時(shí)，x(k+1)應(yīng)為 f(x(k)) x(k+1)=x(k)+Ax(k)=x(k)—(1-3) 八x(k)) 式中k為迭代次數(shù)。 如果f(x(k+l))|<￡(￡是預(yù)設(shè)的一個(gè)小的正數(shù)，如￡=10-5)，則方程的解x*=x(k+1),迭 代停止。 例1-2應(yīng)用牛頓—拉夫遜法求解非線性方程 f(x

3、)=x3一2x2+x一12=0 解：設(shè)初始近似解x(0)=2.0，首先根據(jù)(1-1)計(jì)算f(x(0)) f(x(o))=-10 然后計(jì)算f'（x（0）） 八x（0））二5 根據(jù)（1-2）式計(jì)算Ax（o） Ax（o） f（x（0）） 八x（0）） -10 "I" -9- 再根據(jù)（1-3）式計(jì)算Ax（1） x（1）=x⑼+Ax⑼=2+2=4 重復(fù)以上計(jì)算直到|f（x（k+1））|<10-5，得到的計(jì)算過程量和結(jié)果見表1-1。 表1-1 k 0 1 2 3 4 5 x(k) f(x(k)) f'(x(k)) Ax(k) 2.0

4、 -10 5 2 4.0 24 33 -0.7273 3.2727 4.9046 20.0413 -0.2447 3.0280 0.4536 16.3944 -0.0277 3.003 0.0054 16.0047 -3.3747X10-4 3.000 7.9727X10-7 16 -4.9829X10-8 非線性方程的解x*=3.0000。 、二元非線性方程組求解 例1-3應(yīng)用牛頓—拉夫遜法求非線性方程組的近似解 f(x,x)=3x2 1121 f(x,x)=2x2 2121 +2x2+xx 212 +x2+2

5、xx 212 -x-11.04=0 2 +x-11.31=0 1 解：令X=[x,x]t，F(xiàn)(X)=[f(X),f(X)]t，迭代次數(shù)為k。 1212 F(X)的Jacobi矩陣為 x+4x—1 12 2x+2x 12 6x+x J(X)=12 4x+2x+1 12 設(shè)初始近似解為X（0）=[1.0,2.0]T，X迭代精度取0.0001。計(jì)算過程量和結(jié)果見表1-2。 表1-2 k 0 1 2 3 X(K) 1 X(K) 2 彳(x(K),x2K)) F2(X(k),x(k)) 1.0 2.0 -0.0400 -0.3

6、100 1.0933 1.9117 0.0335 0.0087 1.0999 1.9001 0.3210X10-3 0.0680X10-3 1.1000 1.9000 0.3518X10-7 0.0740X10-7 則x=1.1000,x=1.9000。本例題中經(jīng)過3次迭代就得到了原方程的精確解。當(dāng)然, 12 這是一個(gè)特例。一般情況下只能得到近似解。從例題也能看出牛頓—拉夫遜法的收斂速度是比較快的。 下面給出用MATLAB5.3語(yǔ)言寫的源程序。 clear x(1)=1.0;x(2)=2.0; k=0;precision=1; k,x whil

7、eprecision>0.0001 f1=3*x(1)入2+2*x(2)入2+x(1)*x(2)-x(2)-11.04; f2=2*x(1)入2+x(2)入2+2*x(1)*x(2)+x(1)-11.31; f=[f1f2]' k=k+1; k J=[6*x(1)+x(2)x(1)+4*x(2)-1 4*x(1)+2*x(2)+12*x(1)+2*x(2)]; xx=-J\f; x(1)=x(1)+xx(1); x(2)=x(2)+xx(2); x precision=max(abs(xx)); end 說明：(1)MATLAB是目前國(guó)際上最流行的科學(xué)與工程計(jì)算的軟

8、件工具，它具有強(qiáng)大的數(shù)值計(jì)算和圖形功能。有關(guān)MATLAB的內(nèi)容請(qǐng)讀者參考專門介紹它的書籍。 (2)程序中語(yǔ)句xx=-J\f的功能相當(dāng)于xx=-inv(J)*f，即矩陣xx等于J的逆矩陣的負(fù)數(shù)左乘矩陣f。但是前者比后者的運(yùn)算速度快得多。 三、電力系統(tǒng)潮流計(jì)算例題 例1-4網(wǎng)絡(luò)接線如圖7-2所示，各支路導(dǎo)納均以標(biāo)幺值標(biāo)于圖1-1中。節(jié)點(diǎn)注入功率分 別為:S=0.20+j0.20，S=-0.45—j0.15，S=-0.40—j0.05，S=-0.60—j0.10，其中節(jié) 1234 點(diǎn)1連接的實(shí)際上相當(dāng)于給定功率的發(fā)電廠。設(shè)節(jié)點(diǎn)5電壓保持定值，V5=1.06。試運(yùn)用以 極坐標(biāo)表示的牛頓—

9、拉夫遜法計(jì)算該系統(tǒng)的潮流分布。計(jì)算精度要求個(gè)節(jié)點(diǎn)電壓修正量不大 于10-5。 「1.25-j3.75 io-j3oq 5 t3 51j-5 2.5-j7.5 4 圖1-1以導(dǎo)納表示的等值電路 解：在該系統(tǒng)中，節(jié)點(diǎn)5為平衡節(jié)點(diǎn)，電壓保持定值，V5=1.06。其余4個(gè)節(jié)點(diǎn)都是PQ節(jié) 點(diǎn)，給定的輸入功率分別為：S=0.20+j0.20，S=-0.45-j0.15，S=-0.40-j0.05， 123 S4=-0.60-j0.10。 圖1-2牛頓—拉夫遜法潮流計(jì)算原理框圖 1．給出潮流計(jì)算的基本步驟。 （1）形成節(jié)點(diǎn)導(dǎo)納矩陣 (2)根據(jù)文獻(xiàn)

10、［1］的公式(2-60)和(2-61)式計(jì)算個(gè)節(jié)點(diǎn)功率的不平衡量 (3)根據(jù)文獻(xiàn)［1］的公式(2-64)和(2-65)計(jì)算雅可比矩陣中各元素 (4) 解修正方程式求各節(jié)點(diǎn)電壓 (5) 計(jì)算平衡節(jié)點(diǎn)出功率和線路功率。 2．給出計(jì)算原理框圖如圖1-2。 3．再給出牛頓—拉夫遜法潮流計(jì)算的源程序。 %ThefollowingProgramforloadflowcalculationisbasedonMATLAB5.3.clear G(1,1)=10.834;B(1,1)=-32.500;G(1,2)=-1.667;B(1,2)=5.000;G(1,3)=-1.667; B(1,3)=

11、5.000;G(1,4)=-2.500;B(1,4)=7.5000;G(1,5)=-5.000;B(1,5)=15.000; G(2,1)=-1.667;B(2,1)=5.000;G(2,2)=12.917;B(2,2)=-38.750;G(2,3)=-10.000; B(2,3)=30.0000;G(2,4)=0;B(2,4)=0;G(2,5)=-1.250;B(2,5)=3.750; B(3,1)=5.000;G(3,2)=-10.000;B(3,2)=30.000;G(3,3)=12.917;B(3,3)=-38.750; G(3,4)=-1.250;B(3,4)=3.750;G

12、(3,5)=0;B(3,5)=0;G(3,1)=-1.667; B(4,1)=7.500;G(4,2)=0;B(4,2)=0;G(4,3)=-1.250;B(4,3)=3.750;G(4,4)=3.750;B(4,4)=-11.250;G(4,5)=0;B(4,5)=0;G(4,1)=-2.500; G(5,1)=-5.000;B(5,1)=15.000;G(5,2)=-1.250;B(5,2)=3.750;G(5,3)=0;B(5,3)=0;G(5,4)=0;B(5,4)=0;G(5,5)=6.250;B(5,5)=-18.750; Y=G+j*B; delt(1)=0;delt(2

13、)=0;delt(3)=0;delt(4)=0;u(1)=1.0;u(2)=1.0;u(3)=1.0;u(4)=1.0; p(1)=0.20;q(1)=0.20;p(2)=-0.45;q(2)=-0.15; p(3)=-0.40;q(3)=-0.05;p(4)=-0.60;q(4)=-0.10;k=0;precision=1; N1=4;%theN1istheamountofthePQbus whileprecision>0.00001 delt(5)=0;u(5)=1.06; form=1:N1 forn=1:N1+1 pt(n)=u(m)*u(n)*(G(m,n)*cos(

14、delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n))); qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); end pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt); end form=1:N1 forn=1:N1+1 h0(n)= u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); n0(n)=-u(m)*u(n)*(G(m,n)

15、*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n))); j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n))); L0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); end H(m,m)=sum(hO)-u(m)入2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m))); N(m,m)

16、=sum(n0)-2*u(m)入2*G(m,m)+u(m)入2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m))); J(m,m)=sum(jO)+u(m)入2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m))); L(m,m)=sum(L0)+2*u(m)A2*B(m,m)+u(m)A2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m))); end form=1:N1 JJ(2*m-1,2*m

17、-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m); end form=1:N1 forn=1:N1 ifm==n else H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n))); J(m,n)= u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n))); N(m,n)=-J(m,n);L(m,n)=H(m,

18、n); JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n); end end end form=1:N1 PP(2*m-1)=pp(m);PP(2*m)=qq(m); enduu=-inv(JJ)*PP';precision=max(abs(uu)); forn=1:N1 delt(n)=delt(n)+uu(2*n-1); u(n)=u(n)+uu(2*n); end k=k+1; end k-1,delt',u' %thefollowingp

19、rogramisusedtocalculatetheS5andSmn forn=1:N1+1 U(n)=u(n)*(cos(delt(n))+j*sin(delt(n))); endform=1:N1+1 I(m)=Y(5,m)*U(m); end S5=U(5)*sum(conj(I)) form=1:N1+1 forn=1:N1+1S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n)); end end S說明：(1)編寫程序時(shí)把平衡節(jié)點(diǎn)標(biāo)為最大號(hào)，如在本例中標(biāo)為5號(hào)。 (2)因?yàn)樵贛ATLAB中i和j是作為虛數(shù)單位，所以

20、表示節(jié)點(diǎn)導(dǎo)納矩陣的行號(hào)和列號(hào) 的變量用m和n。 4．潮流計(jì)算結(jié)果 表1-3迭代過程中各節(jié)點(diǎn)的電壓 k V 5 V 5 V 5 V 5 4 4 0 1.0 0 1.0 0 1.0 0 1.0 0 1 2 3 4 1.0430 -0.0473 1.0154 -0.0863 1.0141 -0.0922 1.0093 -0.1076 1.0368 -0.0461 1.0089 -0.0839 0.0074 -0.0896 1.0017 -0.01044

21、 1.0365 -0.0461 1.0088 -0.0839 0.0073 -0.0896 1.0016 -0.1044 1.0365 -0.0461 1.0087 -0.0839 1.0073 -0.0896 1.0016 -0.1044 表1-4各線路功率Smn 1 2 3 4 1 0 0.2469+j0.0815 0.2793+j0.0806 0.5489+j0.1333 2 -0.2431-j0.0701 0 0.1891-j0.0121 0 3 -0.2746-j0.0664 -0.1887+j0.0132 0 0.0633+j0.0033 4 -0.5370-j0.0977 0 -0.0630-j0.0023 0 5 0.8895+j0.1387 0.4087+j0.1058 0 0 5 O.8751-jO.O954 0.3960-j0.0677 0 0 0 參考文獻(xiàn) [1]邱曉燕劉天琪編著電力系統(tǒng)分析的計(jì)算機(jī)算法中國(guó)電力出版社2009年8月第一版