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

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

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

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

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

6、錯誤，得到結果。在形成雅克比矩陣的過程中，存在雅克比矩陣的坐標與節點的坐標對應和轉化的問題。附錄（設計流程圖、程序、表格、數據等）一、流程圖 輸入原始數據啟動設非平衡節點電壓初值形成節點導納矩陣 令迭代次數k=0對PQ節點計算對PV節點計算 令節點號i=1雅客比矩陣是否形成，i>n? 計算雅客比矩陣各元素增加迭代次數k=k+1增加節點號i=i+1解修正方程，由及雅客比矩陣用高斯法求各節點的電壓增量計算節點的新電壓求出迭代是否收斂停止計算平衡節點的功率及線路功率二、程序/*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();/讀入節點個數cout<<"系統共有"<<m_Bus_Num<<"個節點"<<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（為節點數據定義數組）m_Bus = new Busm_Bus_Num; /建立結構體保存節點數據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();/讀入節點電壓VstrTemp = strLine.substr(33,6);m_Busi.theta = atof(strTemp.c_str();/讀入節

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

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;/輸出新的節點編號/ 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()/做導納矩陣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區元素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區的元素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區元素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區元素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&