电力系统潮流上机程序孙英云.docx

附录 程序及结果#include pf.h#include complexusing namespace std;pf:pf(void) m_Line = NULL;m_Bus = NULL; bph=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 (m_Bus!=NULL) delete m_Bus;if (bph!=NULL) delete bph;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_str();/ 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.c_str();cout系统共有m_Bus_Num个节点endl;coutBussetw(6)Typesetw(6)Vsetw(7)Thetasetw(7)LoadPsetw(7)LoadQsetw(7)GenPsetw(7)GenQsetw(9)ShuntBsetw(9)ShuntGendl;/ 2.2 read each bus data here/ allocate memory for m_Busm_Bus = new Busm_Bus_Num;m_Bus_newIdx = new intm_Bus_Num;for (int i = 0; im_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.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(); /读入节点相角 strTemp = strLine.substr(33,6); m_Busi.theta= atof(strTemp.c_str(); /读入负载有功 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(); /读入发电机无功 strTemp = strLine.substr(67,8); m_Busi.GenQ= atof(strTemp.c_str(); /读入转移电导B strTemp = strLine.substr(114,8); m_Busi.ShuntB= atof(strTemp.c_str(); /读入转移电纳G strTemp = strLine.substr(106,8); m_Busi.ShuntG= atof(strTemp.c_str(); couti+1setw(7)m_Busi.Typesetw(7)m_Busi.Vsetw(7)m_Busi.thetasetw(7)m_Busi.LoadPsetw(7)m_Busi.LoadQsetw(9)m_Busi.GenPsetw(9)m_Busi.GenQsetw(9)m_Busi.ShuntBsetw(6)m_Busi.ShuntGendl; / 3. read Line Data here ;/ 3.1 read Line num;getline(fin,strLine);getline(fin,strLine);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_Line_Num = atoi(strTemp.c_str();/ 3.2 read each line data;/ allocate memory for m_Linem_Line = new Linem_Line_Num; coutendl系统共有m_Line_Num条线路endl;coutLinesetw(4)Rsetw(8)Xsetw(8)Bsetw(9)Ksetw(5)endlsetw(1);for (int i = 0; im_Line_Num; i+)getline(fin,strLine);/ 读取支路相关的节点NumINumJstrTemp = strLine.substr(0,4);m_Linei.NumI = atoi(strTemp.c_str();strTemp = strLine.substr(5,4);m_Linei.NumJ = atoi(strTemp.c_str();/读取线路电阻RstrTemp = strLine.substr(19,10);m_Linei.R = atof(strTemp.c_str();/读取线路电抗XstrTemp = strLine.substr(29,10);m_Linei.X = atof(strTemp.c_str();/读取线路电纳BstrTemp = strLine.substr(40,11);m_Linei.B = atof(strTemp.c_str();/读取变压器变比Km_Linei.K= 1;/读取线路类型Type strTemp = strLine.substr(18,1);m_Linei.Type= atoi(strTemp.c_str(); coutm_Linei.NumIm_Linei.NumJsetw(7)m_Linei.Rsetw(7) m_Linei.Xsetw(7)m_Linei.Bsetw(7)m_Linei.Kendl;/ 4. close the file;fin.close();elsecoutfile open fail!endl; return 0;int pf:initPFData()/ 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; im_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 = iSW + m_pq_Num + m_pv_Num;iSW+;break; std:coutendl共有iPQ个PQ节点endl; std:cout共有iPV个PV节点endl; std:cout共有iSW个sw节点endl;/ here we give the size of Jacobi matrix;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_Num,m_Bus_Num);m_Array_DeltaS.setSize(2*m_pq_Num+m_pv_Num,1);bph = new double2*m_pq_Num+m_pv_Num;return 0;int pf:getBusIdx(int busNum)/ return the index of bus from busNumint i,idx = -1;for (i = 0; im_Bus_Num; i+)if (m_Busi.Num = busNum)idx = i;return idx;int pf:makeYMatrix()/ TO DO for(int i=0;im_Line_Num;i+)int ii,jj; ii=m_Linei.NumI-1;jj=m_Linei.NumJ-1;double Z2=m_Linei.R*m_Linei.R+m_Linei.X*m_Linei.X;double YG=m_Linei.R/Z2;double YB=-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_Matrix_G.outputMatrixtoFile(g.txt);m_Matrix_B.outputMatrixtoFile(b.txt);coutendl节电导纳阵为endl;for(int i=0;im_Bus_Num;i+) for(int j=0;jm_Bus_Num;j+) coutm_Matrix_G(i,j)+jm_Matrix_B(i,j)setw(3); if (j=m_Bus_Num-1) coutendl; coutendl;return 0;int pf:makeJacobi()/ TODOdouble H1010=0;double N1010=0;double J1010=0;double L1010=0;int i,j; for(i=0;im_Bus_Num;i+)for(j=0;jm_Bus_Num;j+)if(i=j) for(int k=0;km_Bus_Num;k+)Hii+=-m_Busi.V*m_Busk.V*m_Matrix_G(i,k)*sin(m_Busi.theta-m_Busk.theta)+m_Busi.V*m_Busk.V*m_Matrix_B(i,k)*cos(m_Busi.theta-m_Busk.theta);Jii+=m_Busi.V*m_Busk.V*m_Matrix_G(i,k)*cos(m_Busi.theta-m_Busk.theta)+m_Busi.V*m_Busk.V*m_Matrix_B(i,k)*sin(m_Busi.theta-m_Busk.theta);Nii+=m_Busi.V*m_Busk.V*m_Matrix_G(i,k)*cos(m_Busi.theta-m_Busk.theta)+m_Busi.V*m_Busk.V*m_Matrix_B(i,k)*sin(m_Busi.theta-m_Busk.theta);Lii+=m_Busi.V*m_Busk.V*m_Matrix_G(i,k)*sin(m_Busi.theta-m_Busk.theta)-m_Busi.V*m_Busk.V*m_Matrix_B(i,k)*cos(m_Busi.theta-m_Busk.theta);Hii=Hii+m_Busi.V*m_Busi.V*m_Matrix_G(i,i)*sin(m_Busi.theta-m_Busi.theta)-m_Busi.V*m_Busi.V*m_Matrix_B(i,i)*cos(m_Busi.theta-m_Busi.theta);Jii=Jii-m_Busi.V*m_Busi.V*m_Matrix_G(i,i)*cos(m_Busi.theta-m_Busi.theta)-m_Busi.V*m_Busi.V*m_Matrix_B(i,i)*sin(m_Busi.theta-m_Busi.theta);Nii=Nii-m_Busi.V*m_Busi.V*m_Matrix_G(i,i)*cos(m_Busi.theta-m_Busi.theta)-m_Busi.V*m_Busi.V*m_Matrix_B(i,i)*sin(m_Busi.theta-m_Busi.theta)+2*m_Busi.V*m_Busi.V*m_Matrix_G(i,j);Lii=Lii-m_Busi.V*m_Busi.V*m_Matrix_G(i,i)*sin(m_Busi.theta-m_Busi.theta)+m_Busi.V*m_Busi.V*m_Matrix_B(i,i)*cos(m_Busi.theta-m_Busi.theta)-2*m_Busi.V*m_Busi.V*m_Matrix_B(i,j);else Hij=m_Busi.V*m_Busj.V*m_Matrix_G(i,j)*sin(m_Busi.theta-m_Busj.theta)-m_Busi.V*m_Busj.V*m_Matrix_B(i,j)*cos(m_Busi.theta-m_Busj.theta);Jij=-m_Busi.V*m_Busj.V*m_Matrix_G(i,j)*cos(m_Busi.theta-m_Busj.theta)-m_Busi.V*m_Busj.V*m_Matrix_B(i,j)*sin(m_Busi.theta-m_Busj.theta); Nij=m_Busi.V*m_Busj.V*m_Matrix_G(i,j)*cos(m_Busi.theta-m_Busj.theta)+m_Busi.V*m_Busj.V*m_Matrix_B(i,j)*sin(m_Busi.theta-m_Busj.theta);Lij=m_Busi.V*m_Busj.V*m_Matrix_G(i,j)*sin(m_Busi.theta-m_Busj.theta)-m_Busi.V*m_Busj.V*m_Matrix_B(i,j)*cos(m_Busi.theta-m_Busj.theta);for(i=0;im_pq_Num;i+)/上左for(j=0;jm_pq_Num;j+)m_Jacobi(2*i,2*j)=Hi+m_pv_Num+1j+1+m_pv_Num;m_Jacobi(2*i,2*j+1)=Ni+1+m_pv_Numj+1+m_pv_Num;for(i=0;im_pq_Num;i+)/上左for(j=0;jm_pq_Num;j+)m_Jacobi(2*i+1,2*j)=Ji+m_pv_Num+1j+1+m_pv_Num;m_Jacobi(2*i+1,2*j+1)=Li+1+m_pv_Numj+1+m_pv_Num;for(i=0;im_pq_Num;i+)/上右for(j=2*m_pq_Num;j2*m_pq_Num+m_pv_Num;j+)m_Jacobi(2*i,j)=Hi+m_pv_Num+1j+1-2*m_pq_Num;for(i=0;im_pq_Num;i+)/上右for(j=2*m_pq_Num;j2*m_pq_Num+m_pv_Num;j+)m_Jacobi(2*i+1,j)=Ji+m_pv_Num+1j+1-2*m_pq_Num;for(i=0;im_pv_Num;i+)/下左for(j=0;jm_pq_Num;j+)m_Jacobi(i+2*m_pq_Num,2*j)=Hi+1j+1+m_pv_Num;m_Jacobi(i+2*m_pq_Num,2*j+1)=Ni+1j+1+m_pv_Num;for(i=0;im_pv_Num;i+)/下右for(j=2*m_pq_Num;j2*m_pq_Num+m_pv_Num;j+)m_Jacobi(i+2*m_pq_Num,j)=Hi+m_pv_Num-1j+1-2*m_pq_Num;return 0;int pf:solveLF()/ TODOint i;/1. initializefor ( i = 0; im_Bus_Num; i+)if(m_Busi.Type=1)m_Busi.V = 1;m_Busi.theta = 0;if(m_Busi.Type=2)m_Busi.theta=0; /2. iterateint maxIter = 20;double *bph = new double2*m_pq_Num+m_pv_Num;double epsilon =0.000001;for(int t=0;tmaxIter;t+)calcDeltaS(bph);double error =bph0;for(int k=0;k2*m_pq_Num+m_pv_Num;k+)if(errorbphk)error= bphk;if(errorepsilon)cout一共迭代t次endlendl;break;cout第t+1次雅可比矩阵为makeJacobi(); m_Jacobi.outputMatrix(); m_Jacobi.outputMatrixtoFile(Jacobi.txt); m_Jacobi.solve(2*m_pq_Num+m_pv_Num,bph);for(int k=0;km_Bus_Num;k+)if(m_Busk.Type=1)m_Busk.theta+= bph2*(k-m_pv_Num-1);m_Busk.V +=bph2*(k-m_pv_Num)-1*m_Busk.V;if(m_Busk.Type=2)m_Busk.theta +=bphk+2*m_pq_Num-1; / 2.4