Matlab及C語言在潮流計(jì)算運(yùn)用

1、Matlab 及 C 語言在潮流計(jì)算中的運(yùn)用(三峽大學(xué)電氣信息學(xué)院 20071096 班)摘要：潮流計(jì)算是研究電力系統(tǒng)穩(wěn)態(tài)運(yùn)行情況的一種基本電氣計(jì)算，常規(guī) 潮流計(jì)算的任務(wù)是根據(jù)給定的運(yùn)行條件和網(wǎng)路結(jié)構(gòu)確定整個(gè)系統(tǒng)的運(yùn)行狀 態(tài)，如各母線上的電壓 (幅值及相角) 、網(wǎng)絡(luò)中的功率分布以及功率損耗等。 潮流計(jì)算的結(jié)果是電力系統(tǒng)穩(wěn)定計(jì)算和故障分析的基礎(chǔ)。關(guān)鍵詞：電力系統(tǒng)分析潮流計(jì)算 牛頓-拉夫遜法 C 語言 Matlab一，潮流計(jì)算算法原理：牛頓拉夫遜法的基本原理 牛頓 -拉夫遜法是一種求解非線性方程的數(shù)值解法，由于便于編寫程序用計(jì) 算機(jī)求解，應(yīng)用較廣。下面以一元非線性代數(shù)方程的求解為例，來說明牛 頓-

2、拉夫遜法的基本思想。設(shè)欲求解的非線性代數(shù)方程為f(x)=o設(shè)方程的真實(shí)解為x*，貝U必有f(x*)=O。用牛頓-拉夫遜法求方程真實(shí)解x*的步驟如下：首先選取余割合適的初始估值x 作為方程f(x)=O的解，若恰巧有f(x )=0 ,則方程的真實(shí)解即為x*= x。若f(x)千0則做下一步。取x仁x +為第一次的修正估值，則f(x1)=f(x + x)其中 X為初始估值的增量，即 x =x1設(shè)函數(shù)f(x)具有任意階導(dǎo)數(shù)，即可將上式在X。的鄰域展開為泰勒級(jí)數(shù)，即：f(x1)=f(x + x )=f(x )+f(x ) x +f(x )( x )2/2+若所取的I X。足夠小，則含( x 的項(xiàng)及其余的一

3、切高階項(xiàng)均可略去，并使 其等于零，即：f(x1)f(x )+f(x ) x =0 故得 x =(x )/f(x ) 從而x1= x -f(x )/f(x ) 可見，只要f(x )產(chǎn)即可根據(jù)上式求出第一次的修正估值x1,若恰巧有f(x1)=0，則方程的真實(shí)解即為 x*=x1。若f(x1)藝，則用上述方法由 X1再確定 第二次的修正估值 x2。如此反復(fù)疊代下去，直到求得真實(shí)解x*為止。二，節(jié)點(diǎn)電壓用極坐標(biāo)牛頓拉夫遜法潮流計(jì)算 節(jié)點(diǎn)的功率方程寫成PGi PLiUiU j (Gij cos ijBij sin ij )jii 1,2,3, ,nQGiQLiUiU j (Gij sin ijBij co

4、s ij )ji式中， ij i j 是兩節(jié)點(diǎn)電壓的相角差。把節(jié)點(diǎn)功率表示為節(jié)點(diǎn)電壓的幅值和相角的函數(shù)。 在有 n 個(gè)節(jié)點(diǎn)的系統(tǒng) 中，假定第1m號(hào)節(jié)點(diǎn)為節(jié)點(diǎn)，第 m+1n-1號(hào)節(jié)點(diǎn)為節(jié)點(diǎn)，第 n號(hào)節(jié)點(diǎn) 為平衡節(jié)點(diǎn)。Un和B n是給定的，PV節(jié)點(diǎn)的電壓幅值 Um+1Un-1也是 給定的。因此，只剩下 n-1個(gè)節(jié)點(diǎn)的電壓相角B1，9 2，,& n-1和m個(gè)節(jié)點(diǎn)的電壓幅值 U1，U2，Um是未知量。實(shí)際上，對(duì)于每一個(gè) PQ 節(jié)點(diǎn)或每一個(gè) PV 節(jié)點(diǎn)都可以列寫一個(gè)有功功 率不平衡量方程式，一共包含了 n-1+m 個(gè)方程式，正好同未知數(shù)的數(shù)目相 同，而比直角坐標(biāo)形式的方程少了 n-1-m 個(gè)。對(duì)于方程式

5、可以寫出修正方程式如下：H ij ijU iU j (G ij sinB ij ijcosij)N ij ijU iU j (Gij cosBij ijsinij )J1j UJJij ( j i)M ijN ij ijLijH ij ijH ii iiU i U j (Gij sin ij jiBij cos ij)NiiiijiU i U j (Gij cos ij jiBij sinij)2U i 2GiiM ii iijiUiU j(Gij cos ijBij sinij)LiiiijiU iU j (Gij sin ijji jiBij cosij)2Ui2BiiH是(n-1) x

6、(n-1)階方陣，其元素為；N是(n-1) x m階矩陣，其元素為；K是m x (n-1)階矩陣，其元素為;L是m x m階方陣，其元素為。在這里把節(jié)點(diǎn)不平衡功率對(duì)節(jié)點(diǎn)電壓幅值的偏導(dǎo)數(shù)都乘以該節(jié)點(diǎn)電壓， 相應(yīng)地把節(jié)點(diǎn)電壓的修正量都除以該節(jié)點(diǎn)的電壓幅值，這樣，雅可比矩陣 的表達(dá)式就具有比較整齊的形式。三， matlab 及 c 語言運(yùn)算過程的具體實(shí)例 例：網(wǎng)絡(luò)接線如圖，各支路阻抗和各節(jié)點(diǎn)功率均以標(biāo)幺值標(biāo)于圖中，其 中節(jié)點(diǎn) 2連接的實(shí)際是發(fā)額定功率的發(fā)電廠， 設(shè)節(jié)點(diǎn) 1 的電壓保持為 1.06， 用牛頓 -拉夫遜法計(jì)算系統(tǒng)中的潮流分布。四，程序運(yùn)行結(jié)果：Matlab程序運(yùn)行結(jié)果如下:Y =3.75

7、00 -11.2500i-2.5000 + 7.5000i0-1.2500 + 3.7500i02.5000 + 7.5000i10.8333 -32.5000i-1.6667 + 5.0000i-1.6667 + 5.0000i-5.0000 +15.0000i0-1.6667 + 5.0000i12.9167 -38.7500i-10.0000 +30.0000i-1.2500 + 3.7500i1.2500 + 3.7500i-1.6667 + 5.0000i -10.0000 +30.0000i12.9167 -38.7500i00-5.0000 +15.0000i-1.2500 +

8、3.7500i06.2500 -18.7500icount =4U =1.05551.05721.03191.03311.0600Angle =0.04170.0055-0.0261-0.02360Error_Angle_U =1.0e-008 *-0.5980-0.03660.13340.1379-0.8207-0.4688-0.5584-0.6216S =0.6000 - 0.1000i0.2000 + 0.2000i-0.4500 - 0.1500i-0.4000 - 0.0500i0.0677 + 0.1531iLine_S =00.2905 - 0.1333i00.2991 - 0.

9、0163i0-0.2956 + 0.1229i0.2185 + 0.0764i0.2028 + 0.0746i0.0764 - 0.0760i0-0.2109 - 0.0812i0-0.0918 - 0.0132i-0.1389 - 0.0789i-0.2931 - 0.0031i-0.1960 - 0.0786i0.0919 + 0.0130i0Line=0-0.0771 + 0.0758i0.1448 + 0.0773i0.2697 + 0.1377i00.2824 + 0.0273i-0.2790 - 0.1178i00.2071 - 0.0711i0.1922 - 0.0695i0.0

10、719 + 0.0723i-0.2836 + 0.0097i-0.2022 + 0.0840i0-0.0886 + 0.0151i-0.1326 + 0.0799i-0.1879 + 0.0806i0.0886 - 0.0146i0-0.0727 - 0.0715i0.1366 - 0.0729i0C語言運(yùn)行結(jié)果如下：節(jié)點(diǎn)導(dǎo)納矩陣：3.750000+j-11.250000 -2.500000+j7.500000 -0.000000+j-0.000000 -1.250000+j3.750000 -0.000000+j-0.000000 -2.500000+j7.500000 10.833333+

11、j-32.500000 -1.666667+j5.000000 -1.666667+j5.000000 -5.000000+j15.000000 -0.000000+j-0.000000-1.666667+j5.00000012.916667+j-38.750000-10.000000+j30.000000-1.250000+j3.750000-1.250000+j3.750000-1.666667+j5.000000-10.000000+j30.00000012.916667+j-38.750000-0.000000+j-0.000000-0.000000+j-0.000000 -5.000

12、000+j15.000000 -1.250000+j3.750000 -0.000000+j-0.000000 6.250000+j-18.750000 迭代次數(shù)：count=4節(jié)點(diǎn)相角分別為：9 =0.041710 0.005480 -0.026114 -0.023594 0.000000節(jié)點(diǎn)電壓分別為:U=1.055525 1.057178 1.031931 1.033072 1.060000節(jié)點(diǎn)功率為:S=0.600000+j-0.100000 0.200000+j0.200000 -0.450000+j-0.150000 -0.400000+j-0.050000 0.067700+j0

13、.153099請(qǐng)按任意鍵繼續(xù).經(jīng)過比較matlab和c語言運(yùn)行的結(jié)果是一樣的！最終的潮流分布如下圖：0.0771 - 0.0758i0.0764 - 0.0760i0.2185 + 0.0764i0.0677 + 0.1531i0.45+j0.15-0.1389 - 0.0789i3n-0.0918 - 0.0132i -0.0918 - 0.0132i0.4+j0.05-0.1448 - 0.0773i0.2109 + 0.0812i-0.092 - 0.013i0.1960 + 0.0786i0.2931 + 0.0031i0.2991 - 0.0163i0.2028 + 0.0746i0

14、.2+j0.23956 - 0 1229i0.2905 - 0.1333i1-0.06+0.1i五，小結(jié)在計(jì)算機(jī)潮流計(jì)算中， 使用牛頓 -拉夫遜法進(jìn)行潮流計(jì)算， 非?？旖?。 在程序 編寫過程中，對(duì)公式的由來應(yīng)該非常清楚 ，不然在編寫過程中易出錯(cuò)，比如 for 循環(huán)的次數(shù)，一些細(xì)節(jié)問題都應(yīng)該反復(fù)琢磨。 在程序調(diào)試過程中， 應(yīng)該分步調(diào)試， 不要一開始就運(yùn)行完整的程序。 對(duì)程序中出現(xiàn)的問題應(yīng)該從運(yùn)行結(jié)果中分析問題 的所在，有什么問題解決什么問題，抓住問題的要害。在潮流計(jì)算中運(yùn)用 matlab 比運(yùn)用c語言要方便好多，matlab里有大量的數(shù)學(xué)函數(shù)直接調(diào)用。用c語言寫的 話每一個(gè)函數(shù)都要自己編寫，寫起

15、來很痛苦,其中求逆矩陣就用了一百多行。最后想說的是，想學(xué)沒有什么不可以。六，程序代碼Matlab 程序： %Writer:ChenZhou%Time:2009年 11 月 9日 -2009 年 11 月 15日 %編譯環(huán)境： matlab R2007bclc%說明：為了使節(jié)點(diǎn)按照先PQ,再PV節(jié)點(diǎn)，最后平衡節(jié)點(diǎn)的次序編號(hào)，以便與公式對(duì)照，節(jié)點(diǎn)1與節(jié)點(diǎn) 5 對(duì)調(diào)。%節(jié)點(diǎn)阻抗矩陣 Z=0,0.04+0.12i,0,0.08+0.24i,0;0.04+0.12i,0,0.06+0.18i,0.06+0.18i,0.0 2+0.06i;0,0.06+0.18i,0,0.01+0.03i,0.08+0.

16、24i;0.08+0.24i,0.06+0.18i,0. 01+0.03i,0,0;0,0.02+0.06i,0.08+0.24i,0,0;%求互導(dǎo)納for i=1:5for j=1:5if Z(i,j)=0y(i,j)=0;elsey(i,j)=1/Z(i,j);endendend%求導(dǎo)納for i=1:5for j=1:5if i=jY(i,j)=-y(i,j);elsefor k=1:5Y(i,j)=sum(y(i,:);endendend end %導(dǎo)納矩陣 G=real(Y);B=imag(Y);%賦初值U=1,1,1,1,1.06;Angle=0,0,0,0,0;%節(jié)點(diǎn)功率S=0.

17、6-0.1i,0.2+0.2i,-0.45-0.15i,-0.4-0.05i,0;P=real(S);Q=imag(S);Error_Power=ones(8,1);count=0;% 迭代次數(shù)while max(Error_Power)1e-5for i=1:4H(i,i)=0;N(i,i)=0;M(i,i)=0;L(i,i)=0;Error_P(i)=0;Error_Q(i)=0; endfor i=1:4for j=1:5Error_P(i)=Error_P(i)-U(i)*U(j)*(G(i,j)*cos(Angle(i)-Angle(j)+B(i,j) *sin(Angle(i)-A

18、ngle(j);Error_Q(i)=Error_Q(i)-U(i)*U(j)*(G(i,j)*sin(Angle(i)-Angle(j)-B(i,j) *cos(Angle(i)-Angle(j);endError_P(i)=Error_P(i)+P(i); Error_Q(i)=Error_Q(i)+Q(i);endError_Power=Error_P,Error_Q;%求雅可比矩陣%當(dāng) i=j時(shí)H,N, M丄如下for i=1:4for j=1:4if i=jH(i,j)=-U(i)*U(j)*(G(i,j)*sin(Angle(i)-Angle(j)-B(i,j)*cos(Angle

19、(i)- Angle(j);N(i,j)=-U(i)*U(j)*(G(i,j)*cos(Angle(i)-Angle(j)+B(i,j)*sin(Angle(i)- Angle(j);L(i,j)=-U(i)*U(j)*(G(i,j)*sin(Angle(i)-Angle(j)-B(i,j)*cos(Angle(i)- Angle(j);M(i,j)=U(i)*U(j)*(G(i,j)*cos(Angle(i)-Angle(j)+B(i,j)*sin(Angle(i)-A ngle(j);endendend%當(dāng)i=j 時(shí)H, N, M, L如下for i=1:4for j=1:5if i=jH

20、(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(Angle(i)-Angle(j)-B(i,j)*cos(Ang le(i)-Angle(j);N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(Angle(i)-Angle(j)+B(i,j)*sin(Ang le(i)-Angle(j);M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(Angle(i)-Angle(j)+B(i,j)*sin(Ang le(i)-Angle(j);L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(Angle(i)-Angl

21、e(j)-B(i,j)*cos(Ang le(i)-Angle(j);endendN(i,i)=N(i,i)-2*(U(i)F2*G(i,i);L(i,i)=L(i,i)+2*(U(i)F2*B(i,i);endJ=H,N;M,L;Error_Angle_U=-(inv(J)*Error_Power);Error_Angle=zeros(1,4);Error_U=zeros(1,4);for i=1:4Error_Angle(i)=Error_Angle_U(i);Error_U(i)=Error_Angle_U(i+4)*U(i); endfor i=1:4Angle(i)=Angle(i)

22、+Error_Angle(i); U(i)=U(i)+Error_U(i);endcount=count+1;end%count%平衡節(jié)點(diǎn)的注入功率for j=1:5P(5)=U(5)*U(j)*(G(5,j)*cos(Angle(5)-Angle(j)+B(5,j)*sin(Angle(5)-Ang le(j)+P(5);Q(5)=U(5)*U(j)*(G(5,j)*sin(Angle(5)-Angle(j)-B(5,j)*cos(Angle(5)-Ang le(j)+Q(5);endS(5)=P(5)+Q(5)*sqrt(-1);% 節(jié)點(diǎn)電壓用復(fù)數(shù)表示for i=1:5Complex_U(

23、i)=U(i)*exp(Angle(i)*sqrt(-1);end% 元件端功率for i=1:5for j=1:5Line_S(i,j)=conj(Complex_U(i)*(conj(Complex_U(i)-conj(Complex_U(j)*conj(y(i,j);endend% 元件端電流for i=1:5for j=1:5Line_I(i,j)=conj(Line_S(i,j)/conj(Complex_U(i);endend Y,count,U,Angle,Error_Angle_U,S,Line_S,Line_IC語言程序：此程序一共三個(gè)文件，分別為 InverseMatrix

24、.h ， Functions.h ， ChaoLiuJiSuan.c 主函數(shù)：/Writer:ChenZhou/Time:2009 年11月9日2009 年 11月15日/ 編譯環(huán)境： visual stdio 2008/說明：為了使節(jié)點(diǎn)按照先PQ再PV節(jié)點(diǎn)，最后平衡節(jié)點(diǎn)的次序編號(hào)，以便與公 式對(duì)照，節(jié)點(diǎn)與節(jié)點(diǎn)對(duì)調(diào)。#include #include #include InverseMatrix.h#include Functions.hvoid main()unsigned int i,j,count=0;Z_to_Y();printf( 節(jié)點(diǎn)導(dǎo)納矩陣： n );for (i=0;i5;i+

25、 )for (j=0;j0.00001)|(max_array(Error_Q)0.00001) One_Times();count+;printf( 迭代次數(shù)： count=%d ,count+1);printf( n );printf(節(jié)點(diǎn)相角分別為：0 =);for (i=0;i5;i+)printf( %f ,Anglei);printf( n );printf( 節(jié)點(diǎn)電壓分別為 ：U=);for (i=0;i5;i+)printf( %f ,Ui);Calculate_S();printf( n 節(jié)點(diǎn)功率為 ：S=n );for (i=0;i5;i+ )printf( %f+j%f

26、,S_Pi,S_Qi);printf( n );求逆矩陣的頭文件( InverseMatrix.h 文件 )int max=8;/= 雅可比矩陣 = double J88=0;int switchnum28;void switchrow( int r, int k) /* 將矩陣中的 R,K行交換 */int i;double temp;for (i=0;imax;i+ )temp=Jki;Jki=Jri;Jri=temp;void switchcol( int r, int k) /* 將矩陣中的 R,K行交換 */ int i;double temp;for (i=0;imax;i+ )t

27、emp=Jik;Jik=Jir;Jir=temp;int Inverse_Matrix()int i,j,r,k;double temp;for (k=0;kmax ;k+ )temp=0;r=k;for (j=k;j=max; j+)/* 找主元 */if (tempfabs(Jrj)temp=fabs(Jrj);r=j;if (0=temp) /* 非奇異矩陣 , 退出 */return (0);if (r!=k)switchrow(r,k);switchnum0k=r;switchnum1k=k;elseswitchnum0k=0;switchnum1k=0;Jkk=1/Jkk;for

28、(j=0;jmax;j+ ) if (j!=k)Jkj=-Jkj*Jkk;for (i=0;imax ;i+ )for (j=0;jmax;j+ )if (i!=k & j!=k)Jij=Jij+Jik*Jkj;for (i=0;i=0 ;k-)if (switchnum0k!=0) switchcol(switchnum0k,switchnum1k); 求解過程中的一些功能函數(shù)( Functions.h 文件)/= 計(jì)算導(dǎo)納矩陣/ 阻抗的實(shí)部 doubleZ_Re55=0,0.04,0,0.08,0,0.04,0,0.06,0.06,0.02,0,0.06,0,0.01,0.08,0.08,

29、0.06,0.01,0,0,0,0.02,0.08,0,0;/ 阻抗的虛部doubleZ_Im55=0,0.12,0,0.24,0,0.12,0,0.18,0.18,0.06,0,0.18,0,0.03,0.24,0.24,0.18,0.03,0,0,0,0.06,0.24,0,0;/ 導(dǎo)納矩陣double Y_Re55=0;double Y_Im55=0;/= 節(jié)點(diǎn)電壓 = double U5=1,1,1,1,1.06;/= 節(jié)點(diǎn)電壓相角 = double Angle5=0;/= 節(jié)點(diǎn)電壓相角補(bǔ)償 = double Error_Angle5=0;/= 功率不平衡量 = double Erro

30、r_P5=0;double Error_Q5=0;/= 節(jié)點(diǎn)注入的有功和無功功率 =double S_P5=0;double S_Q5=0;void Z_to_Y() /z=m+nj y_re=m/(m*m+n*n) y_im=-n/(m*m+n*n)junsigned char i=0,j=0,k;double y_temp_Re55=0;double y_temp_Im55=0;for (i=0;i5;i+)for (j=0;j5;j+)if (Z_Reij=0)&(Z_Imij=0)y_temp_Reij=0;y_temp_Imij=0;elsey_temp_Reij=Z_Reij/(Z

31、_Reij*Z_Reij+Z_Imij*Z_Imij);y_temp_Imij=-Z_Imij/(Z_Reij*Z_Reij+Z_Imij*Z_Imij);/= 導(dǎo)納矩陣 = for (i=0;i5;i+)for (j=0;j5;j+)if (i=j)for (k=0;k5;k+)Y_Reij=Y_Reij+y_temp_Reik; Y_Imij=Y_Imij+y_temp_Imik;elseY_Reij=-y_temp_Reij; Y_Imij=-y_temp_Imij;/= 計(jì)算功率不平衡量 = void Error_Power()unsigned char i,j;double Sum_

32、P5=0;double Sum_Q5=0;double Start_P5=0.6,0.2,-0.45,-0.4,0;double Start_Q5=-0.1,0.2,-0.15,-0.05,0;for (i=0;i5;i+)for (j=0;j5;j+)Sum_Pi=Sum_Pi+Ui*Uj*(Y_Reij*cos(Anglei-Anglej)+Y_Imij*sin(Anglei-Anglej);Sum_Qi=Sum_Qi+Ui*Uj*(Y_Reij*sin(Anglei-Anglej)-Y_Imij*cos( Anglei-Anglej);for (i=0;i5;i+)Error_Pi=St

33、art_Pi-Sum_Pi;Error_Qi=Start_Qi-Sum_Qi;/= 形成修正方程式 = void Calculate_J()unsigned char i=0,j=0,k=0;double Sum=0;for (i=0;i4;i+) /=H 參數(shù) =for (j=0;j4;j+)if (i=j)for (k=0;k5;k+)if (i!=k)Sum=Sum+Ui*Uk*(Y_Reik*sin(Anglei-Anglek)-Y_Imik*cos(Anglei-A nglek);elseSum=Sum+0;Jij=Sum;Sum=0;else 宀0_6U)U_S*N=_E_IA+(M26U 宀0_6UV)S8*M2 MIA)*Mn* 曰 nEnsuens(丘一)七(+X9VX0艾)0(H上)七(+&vohdoj懸( + +寸 V 一0上)OJ宀宀宀S