潮流计算代码c

《电力系统潮流上机》课程设计报告院系：电气与电子工程学院班级：电气1405学号： **********学生姓名：指导教师：***设计周数：两周成绩：日期：2017年7月5日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文1．掌握计算机潮流计算的原理：a)复习电力系统分析基础中潮流的计算机算法一章，重点掌握节点分类、潮流算法介绍b)详细阅读牛拉法部分，掌握潮流方程（极坐标、直角坐标）的写法，掌握雅可比矩阵的公式及排列顺序和潮流方程、变量顺序的关系，掌握迭代法收敛条件及迭代法的基本原理c)设计程序框图，划分功能模块、并对每个模块的输入输出量进行细化。

2．编写计算机潮流计算程序a)学习了解IEEE标准格式数据，学习掌握C/C++读取数据的方法b)设计计算机数据存储母线、支路数据的结构，并将所读取的数据存放于所设计的结构当中c)学习节点排序、节点导纳阵计算方法，编写节点导纳阵生成模块d)编写潮流方程不平衡量计算模块e)编写雅可比矩阵生成子模块f)利用给定的pfMatrix类，编写修正量计算模块g)实现潮流计算主程序，并利用IEEE标准节点数据进行校验，要求能够输出计算结果、支路潮流等必要信息3．思考题1．潮流计算的方法有哪些？各有何特点？答：潮流计算分为简单电力网络的手算和复杂电力网络的机算两大类，其中机算又有高斯-赛德尔法、牛顿-拉夫逊法和P-Q分解法。

各方法特点如下所示：手算求解潮流一般只用于简单的网络中，计算量大，对于多节点的网络用手算一般难以解决问题。

但是通过手算可以对物理概念的理解，还可以在运用计算机计算前由手算的形式求取某些原始数据。

2．如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）答：①.所需要收集的数据：A.电网中所有节点的数据：a．各节点的类型,包括平衡节点、PV 节点、PQ 节点b. 对于平衡节点要了解节点的电压大小相位、及节点所能提供的最大最小有功无功功率c. PV节点要知道节点电压大小注入有功功率及节点所能提供的最大和最小无功功.率d. PQ节点要知道节点的注入有功和无功功率B．电网中所有支路的数据：a.各支路类型，即是否含有变压器b.各支路的电阻、电感、电纳c.各变压器的变比。

②．数据整理：将上述数据资料进行分类整理，并为每个节点及支路编上编号。

将整理的结果写成本实验中所要求的格式（原始数据的txt文档），再用本实验所编制的程序进行求解，得到各节点电压、相位，各线路传输功率、损耗，平衡节点注入功率等数值。

③.计算结果分析：考虑PQ节点的电压是否过高或过低；分析PV节点的电压幅值是否正常及无功功率是否超出范围；分析平衡节点有功、无功功率是否在节点所能提供的范围之内；分析给定之路的功率，看是否超出线路的最大传输容量；分析整个系统的网损是否达到标准。

3．设计中遇到的问题和解决的办法。

c++好久没用，有些生疏。

经过复习与百度，渐渐回忆起来。

潮流计算机解法已经遗忘，经过复习查书，很快熟悉起来。

对老师的思路不是很理解，经过与同学一起探讨，得到了正确答案。

三、课程设计总结或结论2016下半年学历电力系统潮流计算，当时并没有编程实践，就背了背矩阵公式。

现在真让我们上手实践，感觉还是略有难度，很有挑战性，毕竟平时没多少机会接触程序。

通过这两周的摸索与交流，最终完成了潮流的编程计算。

由于是在老师的工作基础上进行补充与改造，所以要读懂老师的代码。

我觉得老师的注释还是太少，而且还是英文（虽然英语也能看懂，但还是觉得中文环境用中文好）。

在对节点数据的处理上，我们对老师的思路并不能感到理解，因此在后面雅克比矩阵生成与不平衡量计算模块绕了些许弯路，我最后还是没采用老师的办法。

除了算法的设计外，最恼人的当属开发工具了，机房是vs2010，而我电脑上是vc++6.0与vs2015，一开始用vc写，然后出了一个迷之bug，换到了vs2010才解决。

但我电脑装上vs2010却因为2015的存在无法运行，vs2015也无法运行我在2010下写好的程序。

不想卸掉花老长时间才装上的巨大2015，为此浪费了许多时间，很令人生气。

能够亲手实践电力系统潮流的计算机计算，还是学到了很多知识，对潮流计算那一部分知识又有了更深的印象。

四、参考文献1.《电力系统稳态分析》，陈珩，中国电力出版社，2007年，第三版；2.《高等电力网络分析》，张伯明，陈寿孙，严正，清华大学出版社，2007年，第二版3.《电力系统计算：电子数字计算机的应用》，西安交通大学等合编。

北京：水利电力出版社；4.《现代电力系统分析》，王锡凡主编，科学出版社；二、程序Example.cpp#include <string>#include <iostream>#include <fstream>#include "pf.h"using namespace std;void main(){pf A;A.readDataFromFile("009ieee.dat");A.initPFData();A.makeYMatrix();A.makeJacobi();A.solveLF();A.outputResult();system("pause");}Pf.cpp#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 (m_Bus!=NULL) delete [] m_Bus;if (m_Bus_newIdx!=NULL) delete[] m_Bus_newIdx;}int pf::readDataFromFile(string fileName){int i;string strLine,strTemp;ifstream fin;// open file to 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());// 2.2 read each bus data here// allocate memory for m_Busm_Bus = new Bus[m_Bus_Num];m_Bus_newIdx = new int[m_Bus_Num];for (int i = 0; i<m_Bus_Num; i++){getline(fin,strLine);strTemp = strLine.substr(0,4);// read bus numm_Bus[i].Num = atoi(strTemp.c_str());// read bus Name;strTemp = strLine.substr(5,7);m_Bus[i].Name = strTemp;// read bus type PQ: Type = 1; PV: Type = 2; swing: Type = 3;//判断条件YYstrTemp = strLine.substr(24,2);if (atoi(strTemp.c_str())<=1){m_Bus[i].Type = 1;m_pq_Num ++;}else if(atoi(strTemp.c_str())==2){m_Bus[i].Type = 2;m_pv_Num ++;}else if(atoi(strTemp.c_str())==3){m_Bus[i].Type = 3;m_sw_Num ++;}//read bus VoltagestrTemp = strLine.substr(27,6);m_Bus[i].V = atof(strTemp.c_str());//read bus anglestrTemp = strLine.substr(33,7);m_Bus[i].theta = atof(strTemp.c_str())/180*3.1415926;//read bus Load PstrTemp = strLine.substr(40,9);m_Bus[i].LoadP = atof(strTemp.c_str())/m_SBase;//read bus Load QstrTemp = strLine.substr(49,10);m_Bus[i].LoadQ = atof(strTemp.c_str())/m_SBase;//read bus Gen PstrTemp = strLine.substr(59,8);m_Bus[i].GenP = atof(strTemp.c_str())/m_SBase;//read bus Gen QstrTemp = strLine.substr(67,8);m_Bus[i].GenQ = atof(strTemp.c_str())/m_SBase;//read bus Shunt conductance GstrTemp = strLine.substr(106,8);//read bus Shunt susceptance BstrTemp = strLine.substr(114,8);}// 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;m_Line = new Line[m_Line_Num];for (i = 0;i <m_Line_Num;i++){getline(fin,strLine);//read Tap bus numberstrTemp = strLine.substr(0,4);m_Line[i].NumI = atoi(strTemp.c_str());//read Z bus numberstrTemp = strLine.substr(5,4);m_Line[i].NumJ = atoi(strTemp.c_str());//read line typestrTemp = strLine.substr(18,1);m_Line[i].Type = atoi(strTemp.c_str());//read line resistance RstrTemp = strLine.substr(19,10);m_Line[i].R = atof(strTemp.c_str());//read line reactance XstrTemp = strLine.substr(29,11);m_Line[i].X = atof(strTemp.c_str());//read line charging BstrTemp = strLine.substr(40,10);m_Line[i].B = atof(strTemp.c_str());//read transformer ratiostrTemp = strLine.substr(76,6);m_Line[i].K = atof(strTemp.c_str());}// 4. close the file;fin.close();}elsecout<<"file 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; i<m_Bus_Num; i++)//按PQPVSW排序{switch (m_Bus[i].Type){case 1:m_Bus_newIdx[i] = iPQ;iPQ++;break;case 2:m_Bus_newIdx[i] = iPV + m_pq_Num;iPV++;break;case 3:m_Bus_newIdx[i] = iSW + m_pq_Num + m_pv_Num;iSW++;break;}}for (i =0; i<m_Bus_Num; i++)cout<<m_Bus_newIdx[i]<<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_G.setSize(m_Bus_Num,m_Bus_Num);m_Matrix_B.setSize(m_Bus_Num,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_Bus[i].Num == busNum)idx = i;}return idx;}int pf::makeYMatrix(){// TO DOint i;float Line_G;float Line_B;for(i = 0;i <m_Line_Num;i++){Line_G = m_Line[i].R/(m_Line[i].R*m_Line[i].R+m_Line[i].X*m_Line[i].X);Line_B = -m_Line[i].X/(m_Line[i].R*m_Line[i].R+m_Line[i].X*m_Line[i].X);if(m_Line[i].Type != 2){//m_Matrix_G.DumpInfo("here");m_Matrix_G(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) = m_Matrix_G(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) + Line_G;m_Matrix_G(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumJ)) += -Line_G;m_Matrix_G(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumI)) += -Line_G;m_Matrix_G(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) = m_Matrix_G(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) + Line_G;m_Matrix_B(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) = m_Matrix_B(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) + Line_B + m_Line[i].B/2;m_Matrix_B(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumJ)) += -Line_B;m_Matrix_B(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumI)) += -Line_B;m_Matrix_B(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) = m_Matrix_B(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) + Line_B + m_Line[i].B/2;}else{m_Matrix_G(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) = m_Matrix_G(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) + Line_G/(m_Line[i].K*m_Line[i].K);m_Matrix_G(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumJ)) += -Line_G/m_Line[i].K;m_Matrix_G(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumI)) += -Line_G/m_Line[i].K;m_Matrix_G(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) = m_Matrix_G(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) + Line_G;m_Matrix_B(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) = m_Matrix_B(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumI)) + Line_B/(m_Line[i].K*m_Line[i].K);m_Matrix_B(getBusIdx(m_Line[i].NumI),getBusIdx(m_Line[i].NumJ)) += -Line_B/m_Line[i].K;m_Matrix_B(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumI)) += -Line_B/m_Line[i].K;m_Matrix_B(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) = m_Matrix_B(getBusIdx(m_Line[i].NumJ),getBusIdx(m_Line[i].NumJ)) + Line_B;}}m_Matrix_G.outputMatrixtoFile("G.txt");//实部m_Matrix_B.outputMatrixtoFile("B.txt");//虚部return 0;}int pf::makeJacobi(){int equNum = 2*m_pq_Num+m_pv_Num;int i,j,k;double H,J,N,L;// init Jacobi matrix;for(i = 0;i < equNum;i++){for(j = 0;j < equNum;j++){m_Jacobi(i,j) = 0.0;}}for(i = 0; i< m_Bus_Num; i++){for(j = 0; j<m_Bus_Num; j++){H=0.0;J=0.0;L=0.0;N=0.0;if(i==j){for(int k=0;k<m_Bus_Num;k++){if(i!=k){//H+=-m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m _Bus[k].Num))*sin(m_Bus[i].theta-m_Bus[k].theta)-m_Matrix_B(getBusIdx(m_Bus[i]. Num),getBusIdx(m_Bus[k].Num))*cos(m_Bus[i].theta-m_Bus[k].theta));//J+=m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[k].Nu m))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBu sIdx(m_Bus[k].Num))*sin(m_Bus[i].theta-m_Bus[k].theta));//N+=m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[k].Nu m))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBu sIdx(m_Bus[k].Num))*sin(m_Bus[i].theta-m_Bus[k].theta));//L+=m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_ Bus[k].Num))*sin(m_Bus[i].theta-m_Bus[k].theta)-m_Matrix_B(getBusIdx(m_Bus[i].N um),getBusIdx(m_Bus[k].Num))*cos(m_Bus[i].theta-m_Bus[k].theta));H+=-m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(i,k)*sin(m_Bus[i].theta-m_Bus[k].thet a)-m_Matrix_B(i,k)*cos(m_Bus[i].theta-m_Bus[k].theta));J+=m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(i,k)*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Mat rix_B(i,k)*sin(m_Bus[i].theta-m_Bus[k].theta));N+=m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(i,k)*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Mat rix_B(i,k)*sin(m_Bus[i].theta-m_Bus[k].theta));L+=m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(i,k)*sin(m_Bus[i].theta-m_Bus[k].theta )-m_Matrix_B(i,k)*cos(m_Bus[i].theta-m_Bus[k].theta));}}//N+=2*m_Bus[i].V*m_Bus[i].V*m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m _Bus[i].Num));//L+=-2*m_Bus[i].V*m_Bus[i].V*m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx( m_Bus[i].Num));N+=2*m_Bus[i].V*m_Bus[i].V*m_Matrix_G(i,i);L+=-2*m_Bus[i].V*m_Bus[i].V*m_Matrix_B(i,i);}else{//H = m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[j].Nu m))*sin(m_Bus[i].theta-m_Bus[j].theta)-m_Matrix_B(getBusIdx(m_Bus[i].Num),getBu sIdx(m_Bus[j].Num))*cos(m_Bus[i].theta-m_Bus[j].theta));//J=-m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[j].N um))*cos(m_Bus[i].theta-m_Bus[j].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getB usIdx(m_Bus[j].Num))*sin(m_Bus[i].theta-m_Bus[j].theta));//N=-J;//L=H;H = m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*sin(m_Bus[i].theta-m_Bus[j].theta)-m_Mat rix_B(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta));J=-m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta)+m_Ma trix_B(i,j)*sin(m_Bus[i].theta-m_Bus[j].theta));N=-J;L=H;}//雅可比矩阵int a,b;switch(m_Bus[i].Type){case 1: // PQ bus;switch(m_Bus[j].Type){case 1: // PQ bus// H N J L four elementsm_Jacobi(2*m_Bus_newIdx[i],2*m_Bus_newIdx[j]) = H;m_Jacobi(2*m_Bus_newIdx[i],2*m_Bus_newIdx[j]+1) = N;m_Jacobi(2*m_Bus_newIdx[i]+1,2*m_Bus_newIdx[j]) = J;m_Jacobi(2*m_Bus_newIdx[i]+1,2*m_Bus_newIdx[j]+1) = L;break;case 2: // PV bus// H J two elementsm_Jacobi(2*m_Bus_newIdx[i],m_Bus_newIdx[j]+m_pq_Num) = H;m_Jacobi(2*m_Bus_newIdx[i]+1,m_Bus_newIdx[j]+m_pq_Num) = J;break;case 3: // SW busbreak;}break;case 2: // PV bus;switch(m_Bus[j].Type){case 1: // PQ bus// H N two elementsm_Jacobi(m_Bus_newIdx[i]+m_pq_Num,2*m_Bus_newIdx[j]) = H;m_Jacobi(m_Bus_newIdx[i]+m_pq_Num,2*m_Bus_newIdx[j]+1) = N;break;case 2: // PV bus// H, one elementm_Jacobi(m_Bus_newIdx[i]+m_pq_Num,m_Bus_newIdx[j]+m_pq_Num) = H;break;case 3: // SW busbreak;}break;case 3: // SW bus;break;}}}ofstream fout("J.txt");for(int i=0;i<equNum;i++){for(int j=0;j<equNum;j++)fout<<setw(12)<<m_Jacobi(i,j)<<" ";fout<<endl;}fout.close();return 0;}int pf::solveLF(){// TODOint i;int equNum = 2*m_pq_Num + m_pv_Num;double * bph = new double[equNum];// 1. initializefor ( i = 0; i<m_Bus_Num; i++){switch (m_Bus[i].Type){case 1: //PQ node{m_Bus[i].V = 1;m_Bus[i].theta = 0;break;}case 2: //PV nodem_Bus[i].theta = 0;break;case 3: //SW nodebreak;}}// 2. iterateint maxIter = 20;int p,k;for (i = 0; i<maxIter; i++){// 2.1 calDeltaSif(calcDeltaS(bph)==1) //判断是否收敛{cout<<"一共迭代"<<i+1<<"次收敛"<<endl;break;}// 2.2 calJacobi;cout<<"第"<<i+1<<"次"<<"雅可比矩阵为"<<endl;this->makeJacobi();m_Jacobi.outputMatrix();m_Jacobi.outputMatrixtoFile("Jacobi.txt");// 2.3 output deltaf to screenm_Jacobi.solve(equNum,bph);printf("第%d次修正方程为\n",i+1);for (k = 0; k<2*m_pq_Num + m_pv_Num;k++){printf("%7.7f\n",bph[k]);}printf("\n");printf("第%d次不平衡量为\n",i+1);ofstream fout ("bph.txt");int p;for(p=0;p<equNum;p++)fout<<bph[p]<<"\n";fout<<endl;fout.close();// 2.4 update variables 更新电压相角/*for(int p=0;p<m_Bus_Num;p++){if(m_Bus[p].Type==1){m_Bus[p].theta+= bph[2*(p-m_pv_Num-1)];m_Bus[p].V +=bph[2*(p-m_pv_Num)-1]*m_Bus[p].V;}if(m_Bus[p].Type==2){m_Bus[p].theta +=bph[p+2*m_pq_Num-1];}}*/for (p = 0; p<m_Bus_Num; p++){switch (m_Bus[p].Type){case 1: //PQ node{m_Bus[p].theta=m_Bus[p].theta+bph[2*m_Bus_newIdx[p]];m_Bus[p].V=m_Bus[p].V+bph[2*m_Bus_newIdx[p]+1]*m_Bus[p].V;break;}case 2: //PV nodem_Bus[p].theta=m_Bus[p].theta+bph[m_pq_Num+m_Bus_newIdx[p]];break;case 3: // SW nodebreak;}}}// 3. output the result.if(calcDeltaS(bph)!=1)printf("output unsolved\n");return 0;}int pf::calcDeltaS(double* bph){int i,j,k;int answer = -1;double maximum;for (i = 0; i<2*m_pq_Num + m_pv_Num; i++)bph[i] = 0;for (i = 0; i<m_Bus_Num; i++){switch (m_Bus[i].Type){case 1:// PQ node;for(k=0; k<m_Bus_Num; k++){// active power unblance;bph[2*m_Bus_newIdx[i]] = bph[2*m_Bus_newIdx[i]] + m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[k].Num))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBu sIdx(m_Bus[k].Num))*sin(m_Bus[i].theta-m_Bus[k].theta));// reactive power unblance;bph[2*m_Bus_newIdx[i]+1] = bph[2*m_Bus_newIdx[i]+1] + m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[k].Nu m))*sin(m_Bus[i].theta-m_Bus[k].theta)-m_Matrix_B(getBusIdx(m_Bus[i].Num),getBu sIdx(m_Bus[k].Num))*cos(m_Bus[i].theta-m_Bus[k].theta));}bph[2*m_Bus_newIdx[i]] = m_Bus[i].GenP-m_Bus[i].LoadP-bph[2*m_Bus_newIdx[i]];bph[2*m_Bus_newIdx[i]+1] = m_Bus[i].GenQ-m_Bus[i].LoadQ-bph[2*m_Bus_newIdx[i]+1];break;case 2:// PV node;for(k=0; k<m_Bus_Num; k++){// active power unblance;bph[m_Bus_newIdx[i]+m_pq_Num] = bph[m_Bus_newIdx[i]+m_pq_Num] +m_Bus[i].V*m_Bus[k].V*(m_Matrix_G(getBusIdx(m_Bus[i].Num),getBusIdx(m_Bus[k].Nu m))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBu sIdx(m_Bus[k].Num))*sin(m_Bus[i].theta-m_Bus[k].theta));}bph[m_Bus_newIdx[i]+m_pq_Num] = m_Bus[i].GenP-m_Bus[i].LoadP-bph[m_Bus_newIdx[i]+m_pq_Num];break;case 3:// SW node;break;}}maximum=bph[0];for (i = 0; i<2*m_pq_Num + m_pv_Num;i++){printf("%7.9f\n",bph[i]);if(maximum<fabs(bph[i]))maximum=fabs(bph[i]);}printf("\n");if(maximum<1e-6)answer=1;ofstream fout ("bph1.txt");int p;for(p=0;p<2*m_pq_Num + m_pv_Num;p++)fout<<bph[p]<<"\n";fout<<endl;fout.close();return answer;}int pf::outputResult(){ofstream fout("Result.txt");int p,q;for ( p = 0 ; p<m_Bus_Num; p++){fout<<m_Bus[p].V<<"\t";fout<<m_Bus[p].theta*180/3.14<<"\t";fout<<endl;}fout.close();int j;for(j=0;j<m_Bus_Num;j++)cout<<"节点"<<j+1<<"的电压为"<<m_Bus[j].V<<"\t"<<"相角为"<<m_Bus[j].theta*180/3.14<<endl;return 0;}PfMatrix.cpp#include "pfMatrix.h"pfMatrix::pfMatrix()//data_ <--initialized below (after the 'if/throw' statement){m_row_num = 0;m_col_num = 0;m_value = NULL;m_idx = NULL;}pfMatrix::pfMatrix(int rows, int cols) : m_row_num (rows), m_col_num (cols)//data_ <--initialized below (after the 'if/throw' statement){int i,j;if (rows <= 0 || cols <= 0)throw ("Matrix constructor has wrong size");m_value = new double*[rows];m_idx = new int[m_row_num];for(i = 0 ; i<m_row_num; i++){m_value[i] = new double[m_col_num];m_idx[i] = 0;}for(i = 0; i < m_row_num; i++)for(j = 0; j < m_col_num; j++)m_value[i][j] = 0.0;}pfMatrix::~pfMatrix(){for (int i = 0; i<m_row_num; i++)delete[] m_value[i];delete[] m_value;delete[] m_idx;}double& pfMatrix::operator() (int row, int col){if (row >= m_row_num || col >= m_col_num){#ifdef _DEBUG_MATRIXchar strTemp[100];sprintf(strTemp,"row = %5d col = %5d, Matrix subscript out of bounds\n",row,col);DumpInfo(strTemp);#endifthrow ("const Matrix subscript out of bounds");}return m_value[row][col];}double pfMatrix::operator() (int row, int col) const{if (row >= m_row_num || col >= m_col_num){#ifdef _DEBUG_MATRIXchar strTemp[100];sprintf(strTemp,"row = %5d col = %5d, Matrix subscript out of bounds",row,col);DumpInfo(strTemp);#endifthrow ("const Matrix subscript out of bounds");}return m_value[row][col];}void pfMatrix::outputMatrix(){int i,j;for ( i = 0; i <m_row_num; i++){for ( j = 0; j <m_col_num;j++){printf("%7.3f\t",m_value[i][j]);}printf("\n"); //here is the end of line}printf("\n");}void pfMatrix::LUdcmp(){const double TINY=1.0e-40;int i,imax,j,k;double big,temp;double *vv = new double[m_row_num];double d=1.0;for (i=0;i<m_row_num;i++) {big=0.0;for (j=0;j<m_row_num;j++)if ((temp=abs(m_value[i][j])) > big) big=temp;if (big < TINY) throw("pfMatrix::singular matrix");vv[i]=1.0/big;}for (k=0;k<m_row_num;k++) {big=0.0;for (i=k;i<m_row_num;i++) {temp=vv[i]*abs(m_value[i][k]);if (temp > big) {big=temp;imax=i;}}if (k != imax) {for (j=0;j<m_row_num;j++) {temp=m_value[imax][j];m_value[imax][j]=m_value[k][j];m_value[k][j]=temp;}d = -d;vv[imax]=vv[k];}m_idx[k]=imax;if (m_value[k][k] == 0.0) m_value[k][k]=TINY;for (i=k+1;i<m_row_num;i++) {temp=m_value[i][k] /= m_value[k][k];for (j=k+1;j<m_row_num;j++)m_value[i][j] -= temp*m_value[k][j];}}delete [] vv;}void pfMatrix::solve(int n, double * b){LUdcmp();int i,ii=0,ip,j;double sum;if (m_col_num != n || m_row_num!= n){#ifdef _DEBUG_MATRIXchar strTemp[100];sprintf(strTemp,"the matrix is %5d*%d,but the vector is %5d, so it can't be solved!",m_col_num,m_row_num,n);DumpInfo(strTemp);#endifthrow("pfMatrix::solve bad sizes");}for (i=0;i<n;i++) {ip=m_idx[i];sum=b[ip];b[ip]=b[i];if (ii != 0)for (j=ii-1;j<i;j++) sum -= m_value[i][j]*b[j];else if (sum != 0.0)ii=i+1;b[i]=sum;}for (i=n-1;i>=0;i--) {sum=b[i];for (j=i+1;j<n;j++) sum -= m_value[i][j]*b[j];b[i]=sum/m_value[i][i];}}void pfMatrix::setSize(int rows, int cols){int i,j;if (rows <= 0 || cols <= 0)throw ("Matrix constructor has wrong size");m_row_num = rows;m_col_num = cols;m_value = new double*[rows];m_idx = new int[m_row_num];for(i = 0 ; i<m_row_num; i++){m_value[i] = new double[m_col_num];m_idx[i] = 0;}for(i = 0; i < m_row_num; i++)for(j = 0; j < m_col_num; j++)m_value[i][j] = 0.0;}void pfMatrix::outputMatrixtoFile(string fileName){ofstream fout(fileName.c_str());int i,j;fout<<m_row_num<<"\t"<<m_col_num<<endl;for ( i = 0 ; i<m_row_num; i++){for( j = 0; j<m_col_num; j++)fout<<m_value[i][j]<<"\t";fout<<endl;}fout.close();}void pfMatrix::readMatrixFromFile(string fileName){ifstream ob(fileName.c_str());int i,j;double *value;ob>>m_row_num>>m_col_num;value = new double[m_col_num];setSize(m_row_num,m_col_num);for(i=0;i<m_row_num;i++){for(j = 0; j<m_col_num;j++){ob>>value[j];m_value[i][j]=value[j];}}delete[] value;}void pfMatrix::Tokenize(const string& str,vector<string>& tokens,const string& delimiters = " "){// Skip delimiters at beginning.string::size_type lastPos = str.find_first_not_of(delimiters, 0); // Find first "non-delimiter".string::size_type pos = str.find_first_of(delimiters, lastPos);while (string::npos != pos || string::npos != lastPos){// Found a token, add it to the vector.tokens.push_back(str.substr(lastPos, pos - lastPos));// Skip delimiters. Note the "not_of"lastPos = str.find_first_not_of(delimiters, pos);// Find next "non-delimiter"pos = str.find_first_of(delimiters, lastPos);}}/////////////////////////////////////////////////////////double **branchresistive;//支路功率有功double **branchreactive;//支路功率无功double **deltaresistive;//支路网损有功double **deltareactive;//支路网损无功功率double *PVreactive;//PV节点无功功率//分配内存空间branchresistive=(double **)malloc(m_Bus_Num*sizeof(double));branchreactive=(double **)malloc(m_Bus_Num*sizeof(double));deltaresistive=(double **)malloc(m_Bus_Num*sizeof(double));deltareactive=(double **)malloc(m_Bus_Num*sizeof(double));PVreactive=(double *)malloc(m_pv_Num*sizeof(double));for(i=0;i<m_Bus_Num;i++){branchresistive[i]=(double *)malloc(m_Bus_Num*sizeof(double)); branchreactive[i]=(double *)malloc(m_Bus_Num*sizeof(double)); deltaresistive[i]=(double *)malloc(m_Bus_Num*sizeof(double)); deltareactive[i]=(double *)malloc(m_Bus_Num*sizeof(double)); }//初始化数据for(i=0;i<m_pv_Num;i++){PVreactive[i]=0.0;}for(i=0;i<m_Bus_Num;i++)for(j=0;j<m_Bus_Num;j++){branchresistive[i][j]=0.0;branchreactive[i][j]=0.0;deltaresistive[i][j]=0.0;deltareactive[i][j]=0.0;}//计算各支路功率和网损for(i = 0; i< m_Bus_Num; i++){for(j = 0; j<m_Bus_Num; j++){if(i!=j){branchresistive[i][j]=-m_Bus[i].V*m_Bus[i].V*m_Matrix_G(i,j)+m_Bus[i].V*m_B us[j].V*(m_Matrix_G(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta)+m_Matrix_B(i,j)*sin (m_Bus[i].theta-m_Bus[j].theta));branchreactive[i][j]=m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*sin(m_Bus[i].th eta-m_Bus[j].theta)-m_Matrix_B(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta))+m_Bus[i ].V*m_Bus[i].V*(m_Matrix_B(i,j));}}}for(i = 0; i< m_Bus_Num; i++){for(j = i+1; j<m_Bus_Num; j++){deltaresistive[i][j]=branchresistive[i][j]+branchresistive[j][i];deltareactive[i][j]=branchreactive[i][j]+branchreactive[j][i];if(deltareactive[i][j]!=0)//对于彼此之间有支路的两个节点，应该计入对地导纳m_Line[i].B/2的功率损耗，所以重新计算相应的支路无功{branchreactive[i][j]=m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*sin(m_Bus[i].th eta-m_Bus[j].theta)-m_Matrix_B(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta))+m_Bus[i ].V*m_Bus[i].V*(m_Matrix_B(i,j)+m_Line[i].B/2);branchreactive[j][i]=m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(j,i)*sin(m_Bus[j].theta-m_Bus[i].theta)-m_Matrix_B(j,i)*cos(m_Bus[j].theta-m_Bus[i].theta))+m_Bus[j].V* m_Bus[j].V*(m_Matrix_B(j,i)+m_Line[j].B/2);}}}for(i = 0; i< m_Bus_Num; i++){for(j = i+1; j<m_Bus_Num; j++){deltareactive[i][j]=branchreactive[i][j]+branchreactive[j][i];zongwangsunresistive=zongwangsunresistive+deltaresistive[i][j];zongwangsunreactive=zongwangsunreactive+deltareactive[i][j];}}//计算平衡节点功率及PV节点的无功功率for(i = 0; i< m_Bus_Num; i++){if(m_Bus[i].Type==3)for(j = 0; j<m_Bus_Num; j++){SWresistive=m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*cos(m_Bus[j].theta)+m_Matrix _B(i,j)*sin(m_Bus[j].theta));SWreactive=-m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*sin(m_Bus[j].theta)+m_Matrix _B(i,j)*cos(m_Bus[j].theta));}if(m_Bus[i].Type==2){for(j=0; j<m_Bus_Num; j++){PVreactive[k]=PVreactive[k]+m_Bus[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*sin(m_Bus[i].theta-m_Bus[j].theta)-m_Mat rix_B(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta));}k++;}}//输出支路功率ofstream fout1("branchpower.txt");fout1 << setprecision(6) << endl;for ( p = 0 ; p<m_Bus_Num; p++){for( q = 0; q<m_Bus_Num; q++)if((branchreactive[p][q]!=0)||(branchresistive[p][q]!=0))fout1<<"S"<<p<<"-"<<q<<"支路："<<setw(5)<<branchresistive[p][q]<<"+j"<<branchreactive[p][q]<<"\n";fout1<<endl;}fout1.close();//输出各支路网损和系统总网损ofstream fout2("wangsun.txt");fout2 << setprecision(6) <<endl;for ( p = 0 ; p<m_Bus_Num; p++){。