華北電力大學(xué)潮流計算上機(jī)實(shí)驗(yàn)

1、 電力系統(tǒng)潮流上機(jī)課程設(shè)計報告 院 系：電氣與電子工程學(xué)院 班 級： 電氣1006 學(xué) 號 學(xué)生姓名： 指導(dǎo)教師： 孫英云 設(shè)計周數(shù)： 兩周 成 績： 日期： 2013 年 1 月 10 日一、 目的與要求培養(yǎng)學(xué)生的電力系統(tǒng)潮流計算機(jī)編程能力，掌握計算機(jī)潮流計算的相關(guān)知識二、 主要內(nèi)容1 手算.要求應(yīng)用牛頓-拉夫遜法或P-Q分解法手算求解，要求精度為0.001MW。節(jié)點(diǎn)1為平衡節(jié)點(diǎn)，電壓，節(jié)點(diǎn)2為PQ節(jié)點(diǎn)，負(fù)荷功率，節(jié)點(diǎn)3是PV節(jié)點(diǎn)，兩條支路分別為，對地支路。2 計算機(jī)計算(詳見附錄)3 思考題3.1潮流計算的方法有哪些？各有何特點(diǎn)？答：潮流計算分為手算和機(jī)算兩大類，其中機(jī)算又有高斯-賽德爾迭

2、代法、牛頓-拉夫遜迭代法、P-Q分解法等算法。迭特點(diǎn)：手算求解潮流一般只用在簡單的網(wǎng)絡(luò)中，其計算量大，對于多節(jié)點(diǎn)的網(wǎng)絡(luò)用手算一般難以解決問題，但通過手算可以加深物理概念的理解，還可以在運(yùn)用計算機(jī)計算前以手算求取某些原始數(shù)據(jù)。高斯-賽德爾迭代法：算法簡單，對初值的要求不高，但需要迭代的次數(shù)多，收斂的速度慢，在早期的潮流計算程序中應(yīng)用很多，之后逐漸被牛頓-拉夫遜迭代法所取代，但仍可作為計算程序前幾次迭代的算法，以彌補(bǔ)后者對初值要求高的缺點(diǎn)。牛頓-拉夫遜迭代法：是常用的解非線性方程組的方法，也是當(dāng)前廣泛采用的計算潮流的方法，其收斂速度快，幾次迭代就可以得到最終的結(jié)果。但其缺點(diǎn)是要求初值的選擇得比較接

3、近它們的精確值，否則迭代過程可能不收斂。P-Q分解法潮流計算：派生于以極坐標(biāo)表示時的牛頓-拉夫遜法，其根據(jù)電力系統(tǒng)的特點(diǎn)，對后者的修正方程做了簡化，P-Q分解法的系數(shù)矩陣B和B”代替了牛拉法中的雅可比矩陣J，階數(shù)降低，其中的元素在迭代過程中不發(fā)生變化，而且元素對稱，這些都大大提高了運(yùn)算速度，而且精確度幾乎不受影響。P-Q分解法的收斂特性接近于直線，而牛頓-拉夫遜的收斂速度要比P-Q分解法快，但是由于牛頓-拉夫遜每次迭代都要形成雅客比矩陣，所以一次代的時間比P-Q分解法長。3.2如果交給你一個任務(wù)，請你用已有的潮流計算軟件計算北京城市電網(wǎng)的潮流，你應(yīng)該做哪些工作？（收集哪些數(shù)據(jù)，如何整理，計算結(jié)

4、果如何分析答： （1）所需要收集的數(shù)據(jù)：A. 電網(wǎng)中所有節(jié)點(diǎn)的數(shù)據(jù)：a各節(jié)點(diǎn)的類型,包括平衡節(jié)點(diǎn)、PV 節(jié)點(diǎn)、PQ 節(jié)點(diǎn)；b. 對于平衡節(jié)點(diǎn)要了解節(jié)點(diǎn)的電壓大小相位、及節(jié)點(diǎn)所能提供的最大最小有功無功功率；c. PV節(jié)點(diǎn)要知道節(jié)點(diǎn)電壓大小注入有功功率及節(jié)點(diǎn)所能提供的最大和最小無功功.率；d. PQ節(jié)點(diǎn)要知道節(jié)點(diǎn)的注入有功和無功功率B電網(wǎng)中所有支路的數(shù)據(jù)： a.各支路類型，即是否含有變壓器；b.各支路的電阻、電感、電納；c.各變壓器的變比。（2）數(shù)據(jù)整理：將上述數(shù)據(jù)資料進(jìn)行分類整理，并為每個節(jié)點(diǎn)及支路編上編號。將整理的結(jié)果寫成本 實(shí)驗(yàn)中所要求的格式（原始數(shù)據(jù)的 txt 文檔） ，再用本實(shí)驗(yàn)所編制的

5、程序進(jìn)行求解，得到各節(jié)點(diǎn)電壓、相位，各線路傳輸功率、損耗，平衡節(jié)點(diǎn)注入功率等數(shù)值。（3）計算結(jié)果分析：可以大致了解整個電網(wǎng)的運(yùn)行狀況并按照電力系統(tǒng)知識依次進(jìn)行所需的調(diào)整。如調(diào)節(jié)各節(jié)點(diǎn)電壓的大小以防止超出規(guī)定限度，調(diào)整各線路的有功無功分布狀況分析整個系統(tǒng)的網(wǎng)損是否達(dá)到標(biāo)準(zhǔn)等。3.3設(shè)計中遇到的問題和解決的辦法。1）設(shè)計中遇到是對C+部分指令的不熟悉，通過老師講解、看書和向他人詢問減小了編程語言上的障礙；2）在設(shè)計中，由于一開始知識脈絡(luò)以及思維脈絡(luò)都不是很清晰，導(dǎo)致在程序設(shè)計過程中的各種誤操作，例如在寫不平衡量的時候?qū)Χx數(shù)組的下腳標(biāo)標(biāo)注不清。最好經(jīng)過老師的教導(dǎo)和我們自己的仔細(xì)排查，終于排除掉各類

6、錯誤，得到結(jié)果。在形成雅克比矩陣的過程中，存在雅克比矩陣的坐標(biāo)與節(jié)點(diǎn)的坐標(biāo)對應(yīng)和轉(zhuǎn)化的問題。附錄（設(shè)計流程圖、程序、表格、數(shù)據(jù)等）一、流程圖 輸入原始數(shù)據(jù)啟動設(shè)非平衡節(jié)點(diǎn)電壓初值形成節(jié)點(diǎn)導(dǎo)納矩陣 令迭代次數(shù)k=0對PQ節(jié)點(diǎn)計算對PV節(jié)點(diǎn)計算 令節(jié)點(diǎn)號i=1雅客比矩陣是否形成，i>n? 計算雅客比矩陣各元素增加迭代次數(shù)k=k+1增加節(jié)點(diǎn)號i=i+1解修正方程，由及雅客比矩陣用高斯法求各節(jié)點(diǎn)的電壓增量計算節(jié)點(diǎn)的新電壓求出迭代是否收斂停止計算平衡節(jié)點(diǎn)的功率及線路功率二、程序/*Purpose: for undergraduate power flow calculation courseExa

7、mple2Author: Sun YingyunDate: 9-Jun-2010Version: 1.0.0This class is to help the students to solve pf problems.*/#include "pf.h"using namespace std;pf:pf(void)m_Line = NULL;m_Bus = NULL;m_Bus_newIdx = NULL;m_pv_Num = 0;m_sw_Num = 0;m_pq_Num = 0;pf:pf(void)if (m_Line!=NULL) delete m_Line;if

8、(m_Bus!=NULL) delete m_Bus;if (m_Bus_newIdx!=NULL) delete m_Bus_newIdx;int pf:readDataFromFile(string fileName)string strLine,strTemp;ifstream fin;/ open file for read;fin.open(fileName.c_str();if(!fin.fail()/ 1. read SBase;getline(fin,strLine);strTemp.assign(strLine,31,6);m_SBase = atof(strTemp.c_s

9、tr();/ 2. read Bus Data here ;/ 2.1 read Bus num;getline(fin,strLine);size_t pos_begin, pos_end;pos_begin = strLine.find("FOLLOWS");pos_begin = pos_begin + size_t(10);pos_end = strLine.find("ITEM");strTemp = strLine.substr(pos_begin, pos_end - pos_begin);m_Bus_Num = atoi(strTemp.

10、c_str();/讀入節(jié)點(diǎn)個數(shù)cout<<"系統(tǒng)共有"<<m_Bus_Num<<"個節(jié)點(diǎn)"<<endl;cout<<"Type"<<setw(6)<<"V"<<setw(7)<<"Theta"<<setw(7)<<"LoadP"<<setw(7)<<"LoadQ"<<setw(7)<&

11、lt;"GenP"<<setw(7)<<"GenQ"<<setw(9)<<"ShuntB"<<setw(9)<<"ShuntG"<<endl;/ 2.2 read each bus data here/ allocate memory for m_Bus（為節(jié)點(diǎn)數(shù)據(jù)定義數(shù)組）m_Bus = new Busm_Bus_Num; /建立結(jié)構(gòu)體保存節(jié)點(diǎn)數(shù)據(jù)m_Bus_newIdx = new intm_Bus_Num; for (int

12、i = 0; i<m_Bus_Num; i+)getline(fin,strLine);strTemp = strLine.substr(0,4);/ read bus numm_Busi.Num = atoi(strTemp.c_str();/ read bus Name;strTemp = strLine.substr(5,7);m_Busi.Name = strTemp;/ read bus type PQ: Type = 1; PV: Type = 2; swing: Type = 3;strTemp = strLine.substr(24,2);if (atoi(strTemp

13、.c_str()<=1)m_Busi.Type = 1;m_pq_Num +;else if(atoi(strTemp.c_str()=2)m_Busi.Type = 2;m_pv_Num +;else if(atoi(strTemp.c_str()=3)m_Busi.Type = 3;m_sw_Num +;strTemp = strLine.substr(27,6);m_Busi.V = atof(strTemp.c_str();/讀入節(jié)點(diǎn)電壓VstrTemp = strLine.substr(33,6);m_Busi.theta = atof(strTemp.c_str();/讀入節(jié)

14、點(diǎn)相角strTemp = strLine.substr(39,9);m_Busi.LoadP = atof(strTemp.c_str();/讀入節(jié)點(diǎn)有功負(fù)荷strTemp = strLine.substr(49,10);m_Busi.LoadQ = atof(strTemp.c_str();/讀入節(jié)點(diǎn)無功負(fù)荷strTemp = strLine.substr(59,8);m_Busi.GenP = atof(strTemp.c_str();/讀入節(jié)點(diǎn)電壓VstrTemp = strLine.substr(67,8);m_Busi.GenQ= atof(strTemp.c_str();/讀入節(jié)點(diǎn)

15、電壓VstrTemp = strLine.substr(114,8);m_Busi.ShuntB = atof(strTemp.c_str();/input into ShuntB in systemstrTemp = strLine.substr(106,8);m_Busi.ShuntG = atof(strTemp.c_str();/input into ShuntG in systemcout<<m_Busi.Type<<setw(7)<<m_Busi.V<<setw(7)<<m_Busi.theta<<setw(7

16、)<<m_Busi.LoadP<<setw(7)<<m_Busi.LoadQ<<setw(9)<<m_Busi.GenP<<setw(9)<<m_Busi.GenQ<<setw(9)<<m_Busi.ShuntB<<setw(6)<<m_Busi.ShuntG<<endl;/ 3. read Line Data here ;/ 3.1 read Line num;getline(fin,strLine);getline(fin,strLine);pos

17、_begin = strLine.find("FOLLOWS");pos_begin = pos_begin + size_t(10);pos_end = strLine.find("ITEM");strTemp = strLine.substr(pos_begin, pos_end - pos_begin);m_Line_Num = atoi(strTemp.c_str();m_Line = new Linem_Line_Num;for (int i = 0; i<m_Line_Num; i+)getline(fin,strLine);strTe

18、mp = strLine.substr(0,4);m_Linei.NumI = atoi(strTemp.c_str();strTemp = strLine.substr(5,4);m_Linei.NumJ = atoi(strTemp.c_str();/ read linestrTemp = strLine.substr(19,10);m_Linei.R= atof(strTemp.c_str();/ read bus R;strTemp = strLine.substr(29,10);m_Linei.X= atof(strTemp.c_str();/ read bus X;strTemp

19、= strLine.substr(40,11);m_Linei.B= atof(strTemp.c_str();/ read bus B;m_Linei.K=1; / read bus K;strTemp = strLine.substr(18,1);m_Linei.Type= atoi(strTemp.c_str();/ read bus Type/ 4. close the file;fin.close();elsecout<<"file open fail!"<<endl;return 0; /重新編號pq、pv、swint pf:initPF

20、Data()/ according to Page 132 of ref book 3,/ reindex the bus num ase the sequence PQ PV SW;int iPQ,iPV,iSW;iPQ = 0;iPV = 0;iSW = 0;int i;for( i = 0; i<m_Bus_Num; i+)switch (m_Busi.Type)case 1:m_Bus_newIdxi = iPQ;iPQ+;break;case 2:m_Bus_newIdxi = iPV + m_pq_Num;iPV+;break;case 3:m_Bus_newIdxi = i

21、SW + m_pq_Num + m_pv_Num;iSW+;break;for (i =0; i<m_Bus_Num; i+)cout<<m_Bus_newIdxi<<endl;/輸出新的節(jié)點(diǎn)編號/ here we give the size of Jacobi matrix;/分別給出jacobi矩陣和BG矩陣的空間m_Jacobi.setSize(2*m_pq_Num+m_pv_Num,2*m_pq_Num+m_pv_Num);m_Matrix_B.setSize( m_Bus_Num,m_Bus_Num);m_Matrix_G.setSize(m_Bus_N

22、um,m_Bus_Num);return 0;/？int pf:getBusIdx(int busNum)/ return the index of bus from busNumint i,idx = -1;for (i = 0; i<m_Bus_Num; i+)if (m_Busi.Num = busNum)idx = i;return idx;/做jacobi矩陣int pf:makeYMatrix()/做導(dǎo)納矩陣int ii,jj;double YG,YB;double Z2; for(int i=0;i<m_Line_Num;i+)ii=m_Linei.NumI-1;jj

23、=m_Linei.NumJ-1;Z2= m_Linei.R*m_Linei.R+ m_Linei.X*m_Linei.X;YG=m_Linei.R/Z2;YB=-1*m_Linei.X/Z2;m_Matrix_G(ii,ii)+=YG;m_Matrix_B(ii,ii)+=YB+m_Linei.B/2;m_Matrix_G(jj,jj)+=YG;m_Matrix_B(jj,jj)+=YB+m_Linei.B/2;m_Matrix_G(ii,jj)-=YG;m_Matrix_B(ii,jj)-=YB;m_Matrix_G(jj,ii)-=YG;m_Matrix_B(jj,ii)-=YB;m_Ma

24、trix_B.outputMatrix();m_Matrix_G.outputMatrix();return 0;/做jacobi矩陣int pf:makeJacobi()/ TODOint i,j,k=0,l=0,n=0,h; struct Bus a9;int b9;for(i=0;i<m_Bus_Num;i+)if(m_Busi.Type=1)ak=m_Busi;bk=i;k+; if(m_Busi.Type=2)al+m_pq_Num=m_Busi;bl+m_pq_Num=i;l+; if(m_Busi.Type=3)an+m_pq_Num+m_pv_Num=m_Busi;bn+

25、m_pq_Num+m_pv_Num=i;n+;/計算jacobi矩陣中A區(qū)元素for(i=0;i<k;i+)for(j=0;j<k;j+)if(i=j)for(h=0;h<m_Bus_Num;h+)if(h!=i)m_Jacobi(2*i,2*i)+=ai.V*ah.V*(m_Matrix_G(bi,bh)*sin(ai.theta*3.14/180-ah.theta*3.14/180)-m_Matrix_B(bi,bh)*cos(ai.theta*3.14/180-ah.theta*3.14/180); m_Jacobi(2*i,2*i+1)-=ai.V*ah.V*(m_M

26、atrix_G(bi,bh)*cos(ai.theta*3.14/180-ah.theta*3.14/180)+m_Matrix_B(bi,bh)*sin(ai.theta*3.14/180-ah.theta*3.14/180); m_Jacobi(2*i+1,2*i)-=ai.V*ah.V*(m_Matrix_G(bi,bh)*cos(ai.theta*3.14/180-ah.theta*3.14/180)+m_Matrix_B(bi,bh)*sin(ai.theta*3.14/180-ah.theta*3.14/180); m_Jacobi(2*i+1,2*i+1)-=ai.V*ah.V*

27、(m_Matrix_G(bi,bh)*sin(ai.theta*3.14/180-ah.theta*3.14/180)-m_Matrix_B(bi,bh)*cos(ai.theta*3.14/180-ah.theta*3.14/180); elsecontinue;m_Jacobi(2*i,2*i+1)-=2*ai.V*ai.V*m_Matrix_G(bi,bi);m_Jacobi(2*i+1,2*i+1)+=2*ai.V*ai.V*m_Matrix_B(bi,bi);elsem_Jacobi(2*i,2*j)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*sin(ai.thet

28、a*3.14/180-aj.theta*3.14/180)-m_Matrix_B(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/180);m_Jacobi(2*i,2*j+1)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/180)+m_Matrix_B(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180); m_Jacobi(2*i+1,2*j)=ai.V*aj.V*(m_Matrix_G(bi,bj)*cos(ai.thet

29、a*3.14/180-aj.theta*3.14/180)+m_Matrix_B(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180);m_Jacobi(2*i+1,2*j+1)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180)-m_Matrix_B(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/180);/計算jacobi矩陣中B區(qū)的元素for(i=k;i<k+l;i+)for(j=k;j<k+l;j+)if

30、(i=j)for(h=0;h<m_Bus_Num;h+)if(i!=h)m_Jacobi(i+k,i+k)+=ai.V*ah.V*(m_Matrix_G(bi,bh)*sin(ai.theta*3.14/180-ah.theta*3.14/180)-m_Matrix_B(bi,bh)*cos(ai.theta*3.14/180-ah.theta*3.14/180);elsecontinue;elsem_Jacobi(k+i,k+j)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180)-m_Matri

31、x_B(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/180);/計算jacobi矩陣中C區(qū)元素for(i=0;i<k;i+)for(j=k;j<k+l;j+)m_Jacobi(2*i,k+j)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180)-m_Matrix_B(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/180);m_Jacobi(2*i+1,k+j)=ai.V*aj.V*(m_Matrix_G(bi,b

32、j)*cos(ai.theta*3.14/180-aj.theta*3.14/180)+m_Matrix_B(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180);/計算jacobi矩陣中D區(qū)元素for(i=k;i<k+l;i+)for(j=0;j<k;j+)m_Jacobi(k+i,2*j)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180)-m_Matrix_B(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/

33、180);m_Jacobi(k+i,2*j+1)=-ai.V*aj.V*(m_Matrix_G(bi,bj)*cos(ai.theta*3.14/180-aj.theta*3.14/180)+m_Matrix_B(bi,bj)*sin(ai.theta*3.14/180-aj.theta*3.14/180);m_Jacobi.outputMatrix();return 0;int pf:solveLF()/ TODOint i;/ 1. initializefor ( i = 0; i<m_Bus_Num; i+)if(m_Busi.Type=1)m_Busi.V = 1;m_Busi.

34、theta = 0;if(m_Busi.Type=2)m_Busi.theta = 0;/ 2. iteratedouble *bph=new double2*m_pq_Num+m_pv_Num;int maxIter = 20;double epsilon=0.000001;for(int iter=0;iter<maxIter;iter+)int k=0,l=0;/ 2.1 calDeltaScalcDeltaS(bph);double error=0;for(int k=0;k<2*m_pq_Num+m_pv_Num;k+)error+=bphk*bphk;/cout<

35、<"誤差"<<error<<endl;if(error<epsilon)/cout<<iter<<endl;break;this->makeJacobi();/ 2.2 calJacobi;m_Jacobi.solve(2*m_pq_Num+m_pv_Num,bph);for ( i = 0; i<2*m_pq_Num+m_pv_Num; i+)for ( int j = 0; j<2*m_pq_Num+m_pv_Num;j+)m_Jacobi(i,j)=0.0;/ 2.4 update vari

36、ablesfor ( i = 0; i<m_Bus_Num; i+)if(m_Busi.Type=1) m_Busi.V -=bph2*k+1*m_Busi.V;m_Busi.theta -= bph2*k*180/3.14;k+; if(m_Busi.Type=2)m_Busi.theta -= bph2*m_pq_Num+l*180/3.14;l+; outputResult();/cout<<iter<<endl;/ 3. output the result.outputResult();return 0;int pf:calcDeltaS(double*

37、bph)/ TODO/ calculate the active and reactive power equation and store the result to the double array bphdouble deltaP,deltaQ;double Ui,Uj;double thetai,thetaj;int k=0,l=0;for(int i=0;i<m_Bus_Num;i+)Ui=m_Busi.V;thetai=m_Busi.theta*3.14/180;deltaP=(m_Busi.GenP-m_Busi.LoadP)/100;deltaQ=(m_Busi.GenQ

38、-m_Busi.LoadQ)/100;for(int j=0;j<m_Bus_Num;j+) Uj=m_Busj.V; thetaj=m_Busj.theta*3.14/180; if(m_Busi.Type=1) deltaP-=Ui*Uj*(m_Matrix_G(i,j)*cos(thetai-thetaj)+m_Matrix_B(i,j)*sin(thetai-thetaj); deltaQ-=Ui*Uj*(m_Matrix_G(i,j)*sin(thetai-thetaj)-m_Matrix_B(i,j)*cos(thetai-thetaj); if(m_Busi.Type=2)

39、 deltaP-=Ui*Uj*(m_Matrix_G(i,j)*cos(thetai-thetaj)+m_Matrix_B(i,j)*sin(thetai-thetaj); if(m_Busi.Type=1)bphk=deltaP;bph+k=deltaQ;/cout<<deltaP<<endl;/cout<<deltaQ<<endl;k+;if(m_Busi.Type=2)bphl+2*m_pq_Num=deltaP;/cout<<deltaP<<endl;l+; for(k=0;k<2*m_pq_Num+m_pv_Num;k+)cout<<bphk&