潮流计算代码c++

课程设计报告

《电力系统潮流上机》课程设计报告

院 系：电气与电子工程学院 班 级： 电气1405 学 号： 1141180505 学生姓名：

指导教师： 孙英云 设计周数： 两周

成 绩：

日期：2017年7月5日

课程设计报告

一、课程设计的目的与要求

培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识

二、 设计正文

1． 掌握计算机潮流计算的原理：

a) 复习电力系统分析基础中潮流的计算机算法一章，重点掌握节点分类、潮流算法介

绍 b) 详细阅读牛拉法部分，掌握潮流方程（极坐标、直角坐标）的写法，掌握雅可比矩

阵的公式及排列顺序和潮流方程、变量顺序的关系，掌握迭代法收敛条件及迭代法的基本原理 c) 设计程序框图，划分功能模块、并对每个模块的输入输出量进行细化。 2． 编写计算机潮流计算程序

a) 学习了解IEEE标准格式数据，学习掌握C/C++读取数据的方法

b) 设计计算机数据存储母线、支路数据的结构，并将所读取的数据存放于所设计的结

构当中 c) 学习节点排序、节点导纳阵计算方法，编写节点导纳阵生成模块 d) 编写潮流方程不平衡量计算模块 e) 编写雅可比矩阵生成子模块 f)

利用给定的pfMatrix类，编写修正量计算模块

g) 实现潮流计算主程序，并利用IEEE标准节点数据进行校验，要求能够输出计算结

果、支路潮流等必要信息 3． 思考题

1．潮流计算的方法有哪些？各有何特点？

答：潮流计算分为简单电力网络的手算和复杂电力网络的机算两大类，其中机算又有高斯-赛德尔法、牛顿-拉夫逊法和P-Q分解法。 各方法特点如下所示：

手算求解潮流一般只用于简单的网络中，计算量大，对于多节点的网络用手算一般难以解决问题。但是通过手算可以对物理概念的理解，还可以在运用计算机计算前由手算的形式求取某些原始数据。

课程设计报告

方法 初值要求 迭代次数 收敛速度 精度 应用 早期应用多，现在较少 高斯-赛德尔法 不高 多 慢 牛顿-拉夫逊法 高 少 较快 三者一样 广泛应用 应用较多 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 #include #include #include \

using namespace std;

void main() { pf A;

A.readDataFromFile(\ A.initPFData(); A.makeYMatrix(); A.makeJacobi(); A.solveLF(); A.outputResult(); system(\}

Pf.cpp

#include \

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(\ pos_begin = pos_begin + size_t(10); pos_end = strLine.find(\

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_Bus m_Bus = new Bus[m_Bus_Num];

m_Bus_newIdx = new int[m_Bus_Num];

for (int i = 0; i

getline(fin,strLine);

strTemp = strLine.substr(0,4);

// read bus num

m_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; //判断条件YY

strTemp = 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 Voltage

strTemp = strLine.substr(27,6); m_Bus[i].V = atof(strTemp.c_str()); //read bus angle

strTemp = strLine.substr(33,7);

m_Bus[i].theta = atof(strTemp.c_str())/180*3.1415926; //read bus Load P

strTemp = strLine.substr(40,9);

m_Bus[i].LoadP = atof(strTemp.c_str())/m_SBase; //read bus Load Q

strTemp = strLine.substr(49,10);

m_Bus[i].LoadQ = atof(strTemp.c_str())/m_SBase; //read bus Gen P

strTemp = strLine.substr(59,8);

m_Bus[i].GenP = atof(strTemp.c_str())/m_SBase; //read bus Gen Q

strTemp = strLine.substr(67,8);

m_Bus[i].GenQ = atof(strTemp.c_str())/m_SBase; //read bus Shunt conductance G strTemp = strLine.substr(106,8);

//read bus Shunt susceptance B strTemp = strLine.substr(114,8);

课程设计报告

}

// 3. read Line Data here ; // 3.1 read Line num; getline(fin,strLine); getline(fin,strLine);

pos_begin = strLine.find(\ pos_begin = pos_begin + size_t(10); pos_end = strLine.find(\

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

getline(fin,strLine); //read Tap bus number

strTemp = strLine.substr(0,4);

m_Line[i].NumI = atoi(strTemp.c_str()); //read Z bus number

strTemp = strLine.substr(5,4);

m_Line[i].NumJ = atoi(strTemp.c_str()); //read line type

strTemp = strLine.substr(18,1);

m_Line[i].Type = atoi(strTemp.c_str()); //read line resistance R

strTemp = strLine.substr(19,10); m_Line[i].R = atof(strTemp.c_str()); //read line reactance X

strTemp = strLine.substr(29,11); m_Line[i].X = atof(strTemp.c_str()); //read line charging B

strTemp = strLine.substr(40,10); m_Line[i].B = atof(strTemp.c_str()); //read transformer ratio

strTemp = strLine.substr(76,6);

m_Line[i].K = atof(strTemp.c_str()); }

// 4. close the file; fin.close(); } else

cout<<\

课程设计报告

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

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

// 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 busNum int i,idx = -1;

for (i = 0; i

if (m_Bus[i].Num == busNum) idx = i; }

return idx; }

int pf::makeYMatrix() {

// TO DO int i;

float Line_G; float Line_B;

for(i = 0;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(\

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(\实部 m_Matrix_B.outputMatrixtoFile(\虚部 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

H=0.0; J=0.0; L=0.0; N=0.0; if(i==j) {

for(int k=0;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].Num))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx(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].Num))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx(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].Num),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].theta)-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_Matrix_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_Matrix_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].Num))*sin(m_Bus[i].theta-m_Bus[j].theta)-m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx(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].Num))*cos(m_Bus[i].theta-m_Bus[j].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx(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_Matrix_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_Matrix_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 elements

m_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;

课程设计报告

H;

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 elements

m_Jacobi(2*m_Bus_newIdx[i],m_Bus_newIdx[j]+m_pq_Num) =

m_Jacobi(2*m_Bus_newIdx[i]+1,m_Bus_newIdx[j]+m_pq_Num) = J;

break;

case 3: // SW bus break; }

break;

case 2: // PV bus;

switch(m_Bus[j].Type) {

case 1: // PQ bus

// H N two elements

m_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 element

m_Jacobi(m_Bus_newIdx[i]+m_pq_Num,m_Bus_newIdx[j]+m_pq_Num) = H;

break;

case 3: // SW bus break; }

break;

case 3: // SW bus; break; } } }

ofstream fout(\ for(int i=0;i

课程设计报告

for(int j=0;j

fout<

fout<

fout.close();

return 0; }

int pf::solveLF() {

// TODO int i;

int equNum = 2*m_pq_Num + m_pv_Num; double * bph = new double[equNum]; // 1. initialize

for ( i = 0; i

switch (m_Bus[i].Type) {

case 1: //PQ node {

m_Bus[i].V = 1; m_Bus[i].theta = 0; break; }

case 2: //PV node m_Bus[i].theta = 0; break;

case 3: //SW node break; } }

// 2. iterate int maxIter = 20; int p,k;

for (i = 0; i

// 2.1 calDeltaS

if(calcDeltaS(bph)==1) //判断是否收敛 {

课程设计报告

cout<<\一共迭代\次收敛\ break; }

// 2.2 calJacobi;

cout<<\第\次\雅可比矩阵为\ this->makeJacobi(); m_Jacobi.outputMatrix();

m_Jacobi.outputMatrixtoFile(\ // 2.3 output deltaf to screen

m_Jacobi.solve(equNum,bph);

printf(\第%d次修正方程为\\n\

for (k = 0; k<2*m_pq_Num + m_pv_Num;k++) {

printf(\ }

printf(\

printf(\第%d次不平衡量为\\n\

ofstream fout (\int p;

for(p=0;p

// 2.4 update variables 更新电压 相角

/*

for(int p=0;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

课程设计报告

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 node

m_Bus[p].theta=m_Bus[p].theta+bph[m_pq_Num+m_Bus_newIdx[p]];

break;

case 3: // SW node break; } } }

// 3. output the result. if(calcDeltaS(bph)!=1)

printf(\ 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

switch (m_Bus[i].Type) {

case 1:

// PQ node;

for(k=0; 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].Nu

课程设计报告

m))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx(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].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));

}

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

// 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].Num))*cos(m_Bus[i].theta-m_Bus[k].theta)+m_Matrix_B(getBusIdx(m_Bus[i].Num),getBusIdx(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(\ if(maximum

printf(\

if(maximum<1e-6)

课程设计报告

answer=1;

ofstream fout (\ int p;

for(p=0;p<2*m_pq_Num + m_pv_Num;p++) fout<

return answer; }

int pf::outputResult() {

ofstream fout(\ int p,q;

for ( p = 0 ; p

fout<

fout<

fout.close();

int j;

for(j=0;j

cout<<\节点\的电压为\相角为\ return 0;

}

PfMatrix.cpp

#include \

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 (\ m_value = new double*[rows]; m_idx = new int[m_row_num]; for(i = 0 ; i

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

double& pfMatrix::operator() (int row, int col) {

if (row >= m_row_num || col >= m_col_num) {

#ifdef _DEBUG_MATRIX char strTemp[100]; sprintf(strTemp,\= ] col = ], Matrix bounds\\n\

DumpInfo(strTemp); #endif

throw (\ }

return m_value[row][col]; }

double pfMatrix::operator() (int row, int col) const {

subscript out of

课程设计报告

if (row >= m_row_num || col >= m_col_num) {

#ifdef _DEBUG_MATRIX char strTemp[100]; sprintf(strTemp,\= ] col = ], Matrix subscript out of bounds\

DumpInfo(strTemp); #endif

throw (\ }

return m_value[row][col]; }

void pfMatrix::outputMatrix() { int i,j; for ( i = 0; i

printf(\ } printf(\ }

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 big) big=temp; if (big < TINY) throw(\ vv[i]=1.0/big; } for (k=0;k

课程设计报告

big=0.0; for (i=k;i big) { big=temp; imax=i; } } if (k != imax) { for (j=0;j

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_MATRIX char strTemp[100]; sprintf(strTemp,\matrix is ]*%d,but the vector is ], so it can't be solved!\ DumpInfo(strTemp); #endif throw(\ }

课程设计报告

for (i=0;i=0;i--) { sum=b[i]; for (j=i+1;j

void pfMatrix::setSize(int rows, int cols) { int i,j;

if (rows <= 0 || cols <= 0)

throw (\ m_row_num = rows; m_col_num = cols;

m_value = new double*[rows]; m_idx = new int[m_row_num]; for(i = 0 ; i

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<

课程设计报告

for ( i = 0 ; i

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>value[j]; m_value[i][j]=value[j]; } } delete[] value; }

void pfMatrix::Tokenize(const string& str, vector& tokens,

const string& delimiters = \{

// Skip delimiters at beginning.

string::size_type lastPos = str.find_first_not_of(delimiters, 0); // Find first \

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 \

lastPos = str.find_first_not_of(delimiters, pos);

课程设计报告

// Find next \

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

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

for(i=0;i

课程设计报告

{ for(j = 0; j

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

课程设计报告

} } //计算平衡节点功率及PV节点的无功功率 for(i = 0; i< m_Bus_Num; i++) { if(m_Bus[i].Type==3) for(j = 0; 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[i].V*m_Bus[j].V*(m_Matrix_G(i,j)*sin(m_Bus[i].theta-m_Bus[j].theta)-m_Matrix_B(i,j)*cos(m_Bus[i].theta-m_Bus[j].theta)); } k++; } } //输出支路功率 ofstream fout1(\ fout1 << setprecision(6) << endl; for ( p = 0 ; p

//输出各支路网损和系统总网损 ofstream fout2(\ fout2 << setprecision(6) <

课程设计报告

for( q = p+1; q

fout2<<\总网损：\ //输出PV节点的无功功率 ofstream fout3(\ fout3 << setprecision(6) << endl; for ( p = 0 ; p

//输出平衡节点的有功和无功功率 ofstream fout4(\ fout4 << setprecision(6) << endl; fout4<

课程设计报告

三、结果

第一次雅克比矩阵：

随后两次雅克比矩阵

课程设计报告

收敛后的电压幅值相角

迭代次数：3

各节点电压幅值和相角：（节点类型3为平衡节点，2为PV节点，3为PQ节点） 幅值 相角 节点类型3 1.04 0 节点类型2 1.025 9.28471 节点类型2 1.025 4.66712 节点类型1 1.02579 -2.21791

课程设计报告

节点类型1 0.995631 -3.99083 节点类型1 1.01265 -3.68927 节点类型1 1.02577 3.72159 节点类型1 1.01588 0.727907 节点类型1 1.03235 1.96772

各支路功率：(S1-6代表节点1流向节点6的功率） S0-3支路：0.71641+j0.270459 S1-6支路： 1.63+j0.0665362 S2-8支路： 0.85+j-0.108597 S3-0支路：-0.71641+j-0.160839 S3-4支路：0.409373+j0.39992 S3-5支路：0.307037+j0.171819 S4-3支路：-0.406798+j-0.196051 S4-6支路：-0.843202+j0.142127 S5-3支路：-0.305373+j0.0724752 S5-8支路：-0.594628+j0.205889 S6-1支路：-1.63+j0.280126 S6-4支路：0.866202+j0.265523 S6-7支路：0.763799+j0.25876 S7-6支路：-0.759046+j0.0606611 S7-8支路：-0.240954+j-0.0442948 S8-2支路：-0.85+j0.233748

S8-5支路：0.608166+j0.0942158 S8-7支路：0.241834+j0.226761 各支路网损与总的网损：（0-3支路代表0节点到3节点之间支路上损耗的功率） 0-3支路： 0+j0.10962 1-6支路： 0+j0.346662 2-8支路： 0+j0.125151 3-4支路：0.00257514+j0.20387 3-5支路：0.00166407+j0.244294 4-6支路：0.0229997+j0.407651 5-8支路：0.0135385+j0.300105 6-7支路：0.00475284+j0.319421 7-8支路：0.000879965+j0.182466 总网损：0.0464102+j2.23924 平衡节点功率： 0+j0

PV节点的无功功率：

PV节点功率: 0.0665362+j PV节点功率: -0.108597+j

课程设计报告

节点类型1 0.995631 -3.99083 节点类型1 1.01265 -3.68927 节点类型1 1.02577 3.72159 节点类型1 1.01588 0.727907 节点类型1 1.03235 1.96772

各支路功率：(S1-6代表节点1流向节点6的功率） S0-3支路：0.71641+j0.270459 S1-6支路： 1.63+j0.0665362 S2-8支路： 0.85+j-0.108597 S3-0支路：-0.71641+j-0.160839 S3-4支路：0.409373+j0.39992 S3-5支路：0.307037+j0.171819 S4-3支路：-0.406798+j-0.196051 S4-6支路：-0.843202+j0.142127 S5-3支路：-0.305373+j0.0724752 S5-8支路：-0.594628+j0.205889 S6-1支路：-1.63+j0.280126 S6-4支路：0.866202+j0.265523 S6-7支路：0.763799+j0.25876 S7-6支路：-0.759046+j0.0606611 S7-8支路：-0.240954+j-0.0442948 S8-2支路：-0.85+j0.233748

S8-5支路：0.608166+j0.0942158 S8-7支路：0.241834+j0.226761 各支路网损与总的网损：（0-3支路代表0节点到3节点之间支路上损耗的功率） 0-3支路： 0+j0.10962 1-6支路： 0+j0.346662 2-8支路： 0+j0.125151 3-4支路：0.00257514+j0.20387 3-5支路：0.00166407+j0.244294 4-6支路：0.0229997+j0.407651 5-8支路：0.0135385+j0.300105 6-7支路：0.00475284+j0.319421 7-8支路：0.000879965+j0.182466 总网损：0.0464102+j2.23924 平衡节点功率： 0+j0

PV节点的无功功率：

PV节点功率: 0.0665362+j PV节点功率: -0.108597+j

本文来源：https://www.bwwdw.com/article/uthr.html