电网1202刘妍君1121600218.doc

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

院 系：电气与电子工程学院
班 级： 电网1202
学 号： 1121600218
学生姓名： 刘妍君
指导教师： 毛安家
设计周数： 两周

成 绩：

日期：2014年12月30日
课程 课程设计报告
一、课程设计的目的与要求
培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识

二、设计正文（详细内容见附录）
1． 手算：
要求应用牛顿-拉夫逊法或P-Q分解法手算求解，要求迭代两次。
节点1为平衡节点，电压，节点2为PQ节点，负荷功率，节点3是PV节点，，两条支路分别为，，对地支路。手算结果见附录。
根据题意得出节点导纳矩阵为：

1）第一次迭代：
①不平衡量为：

②雅克比矩阵为：

③由修正方程

解得：

④节点电压新值为：

2）第二次迭代：
①不平衡量为：

②雅克比矩阵：

③由修正方程可得：

④迭代后的新值为：

2． 编写潮流计算程序，计算简单系统.
简单系统如下图所示，支路数据如下：
支路14，27，39为变压器支路，参数为
，

其余支路为线路支路，参数为
,

节点数据如下：

，
,,
牛顿-拉夫逊法的求解过程为：

（程序代码以及结果见附录）

3．思考题
3.1潮流计算的方法有哪些？各有何特点？
潮流计算的方法主要有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法和P-Q分解法。 它们各自的特点如下： 
（1）高斯-赛德尔迭代法分为以节点导纳矩阵为基础的高斯-赛德尔迭代法和以以节点阻抗矩阵为基础的高斯-赛德尔迭代法的原理比较简单，要求的数字计算机的内存量也比较小，但收敛性差，当系统规模变大时，迭代次数急剧上升，往往出现 迭代不收敛的情况；而阻抗法改善了电力系统潮流计算导纳法德收敛性问题，在当时获得了广泛的应用，但是，阻抗法的主要缺点是占用计算机的内存很大，每次迭代的计算量很大。当系统不断扩大时，这些缺点就更加突出。 
（2）牛顿-拉夫逊法是数学中求解非线性方程式的典型方法，有较好的收敛性。只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流计算程序的计算效率，牛顿法在收敛性、 内存要求、计算速度方面都超过了阻抗法，成为知道目前仍被广泛采用的方法。 
（3）P-Q分解法潮流计算派生于以极坐标表示时的牛顿-拉夫逊法，它根据电力系统的特点，抓住主要矛盾，对纯数学的牛顿法进行了改造。与牛顿法相比，P-Q分解法的修正方程的系数矩阵B’和B”分别是（n-1）和（m-1）的方阵，替代了原有的（n+m-2）阶系数矩阵J；B’、B”在迭代过程中可以保持不变且为对称的系数矩阵，提高了计算速度，降低了对存储容量的要求。P-Q分解法在计算速度方面有显著地提高，迅速得到了推广。 
3.2如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）
在进行北京城市电网的潮流计算之前需要了解北京城市电网中所有的节点支路的相关数据，并对节点和支路分类。处理PQ节点时要了解节点的注入有功和无功功率；PV节点要了解节点电压大小，注入有功功率及节点所能提供的最大和最小无功功率；对于平衡节点要了解节点的电压大小、相位及节点所能提供的最大和最小有功无功功率，此外还需的节点数据有：支路的电阻电抗电纳支路变压器的变比及线路所能传输的最大容量等。 
3.3设计中遇到的问题和解决的办法。
（2）实验中需要使用文件的输入输出，学习了简单的文件读入读出函数及文件指针的使用方法。
（3）关于导纳矩阵的形成，实验中有两种类型的线路，通过循环读入每条线路的信息并修改导纳矩阵，在读入输电线线路是由于疏忽讲电纳取了倒数，以致导纳矩阵gY_B中部分数值很大，幸好及时发现问题并改正。其次是变压器线路，变压器线路的等效模型参考《电力系统分析基础》书中第89页图4-3，程序要求读入数据中变压器线路的起点与终点的比值为1：k。
（7）由于程序设计初期设置变量大多将一个复数的实部虚部分开保存，在计算线路功率是要将公式展开，对应计算实部虚部，但是如果程序定义了复数类并且重载运算符，会使线路功率计算更加便捷。

三、课程设计总结或结论

四、参考文献
1. 《电力系统计算：电子数字计算机的应用》，西安交通大学等合编。北京：水利电力出版社；
2. 《现代电力系统分析》，王锡凡主编，科学出版社；
3. 《电力系统稳态分析》，陈珩，中国电力出版社，1995年，第三版；
4. 《算法竞赛入门经典训练指南》，刘汝佳清华大学出版社，2012年，第一版；
附录（设计流程图、程序、表格、数据等）
1、程序设计流程图
开始
读入数据
形成导纳矩阵
设置初值，设cas=1
计算不平衡量，雅克比矩阵
高斯消元求修正值
更新新值
存在△e，△f大于精度或者cas<规定值
Cas++
计算不平衡量，雅克比矩阵
高斯消元求修正值
更新新值
计算各节点电压值及各线路功率，网损
将结果输出到文件中
结束程序
否
是

2、计算机潮流程序
（1）头文件 ”DeclearVar.h”
#ifndef DECLEARVAR_H_INCLUDED
#define DECLEARVAR_H_INCLUDED
const int N=100;
const int NLINE=100;
const double Limit=1e-4;
int n_node,n_line,cas;

//平衡节点 type=0,e,f;PQ节点 type=1，PQ;PV节点 type=2,p v;
struct tynode
{
int id,type;
double noda,nodb;
}gBus[N];

//变压器type=1,输电线type=2;
struct tyline
{
int id,type,start,ends;
double x,r,bk;
}gLine[NLINE];
double gY_G[N][N],gY_B[N][N];
double gDelta_P[N],gDelta_Q[N],gDelta_PQ[2*N];
double gJaccobi[2*N][2*N];
double gDelta_f[N],gDelta_e[N],gDelta_fe[2*N];
double gf[N],ge[N];
struct tynodeans
{
int id,type;
double e,f,v,g,p,q;
}nodeans[N];
struct fushu
{
double p,q;
}S[N][N];
#endif // DECLEARVAR_H_INCLUDED
（2）头文件 ”DeclearFunc.h”
#ifndef DECLEARFUNC_H_INCLUDED
#define DECLEARFUNC_H_INCLUDED
void GetData();
void GetYMatrix();
void SetInitial();
void GetUnbalance();
void GetJaccobi();
void GetRevised();
void GetNewValue();
void CalFlow();
void GetAns();
void OutPut();
#endif // DECLEARFUNC_H_INCLUDED
（3）主函数”main.cpp”
#include<cstdio>
#include<cmath>
#include "DeclearVar.h"
#include "DeclearFunc.h"
using namespace std;
#define DataPath "data.txt"
/******************/
//复数由实部虚部转化为幅值p和相角q；
fushu Getvg(fushu a)
{
fushu tmpf;
tmpf.p=sqrt(a.p*a.p+a.q*a.q);
tmpf.q=atan(a.q/a.p);
return tmpf;
}

/******************/
void GetData()
{
int i;
FILE *fp;
fp=fopen(DataPath,"r");
if(fp==NULL)
printf("data.txt open error!\n");
else
{
fscanf(fp,"%d %d",&n_node,&n_line);
for(i=0;i<n_node;i++)
fscanf(fp,"%d%d%lf%lf",&gBus[i].id,&gBus[i].type,&gBus[i].noda,&gBus[i].nodb);
for(i=0;i<n_line;i++)
fscanf(fp,"%d%d%d%d%lf%lf%lf",&gLine[i].id,&gLine[i].type,&gLine[i].start,&gLine[i].ends,&gLine[i].r,&gLine[i].x,&gLine[i].bk);
}
/*
//DataRead Judge
printf("%d %d\n",n_node,n_line);
for(i=0;i<n_node;i++)
printf("%d %d %lf %lf\n",gBus[i].id,gBus[i].type,gBus[i].noda,gBus[i].nodb);
for(i=0;i<n_line;i++)
printf("%d %d %d %d %lf %lf %lf\n",gLine[i].id,gLine[i].type,gLine[i].start,gLine[i].ends,gLine[i].r,gLine[i].x,gLine[i].bk);
*/
printf("文件读入搞定^-^\n");
return ;
}

void GetYMatrix()
{
int i,j,tst,ten;
double tmp,tg,tb,k;
for(i=0;i<n_node;i++)
for(j=0;j<n_node;j++)
{
gY_G[i][j]=0;
gY_B[i][j]=0;
}
for(i=0;i<n_line;i++)
{
tmp=gLine[i].x*gLine[i].x+gLine[i].r*gLine[i].r;
tg=gLine[i].r/tmp;
tb=-gLine[i].x/tmp;
//printf("r=%.3lf x=%.3lf tmp=%.3lf tg=%.3lf tb=%3lf\n",gLine[i].r,gLine[i].x,tmp,tg,tb);
tst=gLine[i].start-1;
ten=gLine[i].ends-1;
if(gLine[i].type==1)
{
k=gLine[i].bk;
gY_G[tst][tst]+=(tg);
gY_B[tst][tst]+=(tb);

gY_G[ten][ten]+=(tg/k/k);
gY_B[ten][ten]+=(tb/k/k);

gY_G[tst][ten]-=(tg/k);
gY_B[tst][ten]-=(tb/k);

gY_G[ten][tst]-=(tg/k);
gY_B[ten][tst]-=(tb/k);
}
else if(gLine[i].type==2)
{
k=gLine[i].bk;
gY_G[tst][tst]+=(tg);
gY_B[tst][tst]+=(tb+k);

gY_G[ten][ten]+=(tg);
gY_B[ten][ten]+=(tb+k);

gY_G[tst][ten]-=(tg);
gY_B[tst][ten]-=(tb);

gY_G[ten][tst]-=(tg);
gY_B[ten][tst]-=(tb);
}
}
FILE *fp;
fp=fopen("YMatrix.txt","w");
if(fp==NULL)
printf("YMatrix.txt文件打开失败！\n");
for(i=0;i<n_node;i++)
{
for(j=0;j<n_node;j++)
fprintf(fp,"%.3lf %.3lf\t",gY_G[i][j],gY_B[i][j]);
fprintf(fp,"\n");
}
printf("导纳矩阵搞定^-^\n");
}

void SetInitial()
{
int i;
for(i=0;i<n_node;i++)
{
if(gBus[i].type==0) //平衡节点
{
ge[i]=gBus[i].noda;
gf[i]=gBus[i].nodb;
}
else if(gBus[i].type==1) //PQ节点
{
ge[i]=1;
gf[i]=0.0;
}
else if(gBus[i].type==2) //PV节点
{
ge[i]=1;//gBus[i].nodb;
gf[i]=0.0;
}
else
printf("输入的节点类型出错！\n");
}
/*
for(i=0;i<n_node;i++) printf("%.3lf ",ge[i]);printf("\n");
for(i=0;i<n_node;i++) printf("%.3lf ",gf[i]);printf("\n");
*/
printf("初始化数据完成^-^\n");
}

void GetUnbalance()
{
int i,j;
double tmp;
for(i=1;i<n_node;i++)
{
if(gBus[i].type==2)
{
tmp=0;
for(j=0;j<n_node;j++)
tmp+=(ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]));
gDelta_P[i]=gBus[i].noda-tmp;
gDelta_Q[i]=gBus[i].nodb*gBus[i].nodb-ge[i]*ge[i]-gf[i]*gf[i];
}
else if(gBus[i].type==1)
{

tmp=0;
for(j=0;j<n_node;j++)
tmp+=ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
gDelta_P[i]=gBus[i].noda-tmp;
tmp=0;
for(j=0;j<n_node;j++)
tmp+=gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
gDelta_Q[i]=gBus[i].nodb-tmp;
}
}
for(i=1;i<n_node;i++) //gDelta_PQ从2开始
{
gDelta_PQ[i*2]=gDelta_P[i];
gDelta_PQ[i*2+1]=gDelta_Q[i];
//printf("%d %.5lf %d %.5lf\n",i*2,gDelta_PQ[i*2],i*2+1,gDelta_PQ[i*2+1]);
}
FILE *fp;
fp=fopen("Unbalance.txt","w");
if(fp==NULL)
printf("Unbalance.txt打开失败\n");
fprintf(fp,"节点类型1为pv节点，2为pq节点\n");
for(i=2;i<n_node*2;i+=2)
fprintf(fp,"节点%d类型%d,Delta_p=%lf Delta_q=%lf\n",i/2+1,gBus[i/2].type,gDelta_PQ[i],gDelta_PQ[i+1]);
}

void GetJaccobi()
{
int i,j,k;
double tmp;
for(i=1;i<n_node;i++)
{
if(gBus[i].type==1)
{
for(j=1;j<n_node;j++)
{
if(i!=j)
{
gJaccobi[i*2][j*2]=gY_B[i][j]*ge[i]-gY_G[i][j]*gf[i];
gJaccobi[i*2][j*2+1]=-gY_G[i][j]*ge[i]-gY_B[i][j]*gf[i];
gJaccobi[i*2+1][j*2]=gY_B[i][j]*gf[i]+gY_G[i][j]*ge[i];
gJaccobi[i*2+1][j*2+1]=-gY_G[i][j]*gf[i]+gY_B[i][j]*ge[i];
}
else
{
tmp=0;
for(k=0;k<n_node;k++)
tmp+=(gY_G[i][k]*gf[k]+gY_B[i][k]*ge[k]);
gJaccobi[i*2][j*2]=-tmp-gY_G[i][i]*gf[i]+gY_B[i][i]*ge[i];
tmp=0;
for(k=0;k<n_node;k++)
tmp+=(gY_G[i][k]*ge[k]-gY_B[i][k]*gf[k]);
gJaccobi[i*2][j*2+1]=-tmp-gY_G[i][i]*ge[i]-gY_B[i][i]*gf[i];
tmp=0;
for(k=0;k<n_node;k++)
tmp+=(gY_G[i][k]*ge[k]-gY_B[i][k]*gf[k]);
gJaccobi[i*2+1][j*2]=-tmp+gY_G[i][i]*ge[i]+gY_B[i][i]*gf[i];
tmp=0;
for(k=0;k<n_node;k++)
tmp+=(gY_G[i][k]*gf[k]+gY_B[i][k]*ge[k]);
gJaccobi[i*2+1][j*2+1]=+tmp-gY_G[i][i]*gf[i]+gY_B[i][i]*ge[i];
}
}
}
else if(gBus[i].type==2)
{
for(j=1;j<n_node;j++)
{
if(i!=j)
{
gJaccobi[i*2][j*2]=gY_B[i][j]*ge[i]-gY_G[i][j]*gf[i];
gJaccobi[i*2][j*2+1]=-gY_G[i][j]*ge[i]-gY_B[i][j]*gf[i];
gJaccobi[i*2+1][j*2]=0.0;
gJaccobi[i*2+1][j*2+1]=0.0;
}
else
{
tmp=0;
for(k=0;k<n_node;k++)
tmp+=(gY_G[i][k]*gf[k]+gY_B[i][k]*ge[k]);
gJaccobi[i*2][j*2]=-tmp-gY_G[i][i]*gf[i]+gY_B[i][i]*ge[i];
tmp=0;
for(k=0;k<n_node;k++)
tmp+=(gY_G[i][k]*ge[k]-gY_B[i][k]*gf[k]);
gJaccobi[i*2][j*2+1]=-tmp-gY_G[i][i]*ge[i]-gY_B[i][i]*gf[i];
gJaccobi[i*2+1][j*2]=-2*gf[i];
gJaccobi[i*2+1][j*2+1]=-2*ge[i];
}
}
}
}
FILE *fp1;
fp1=fopen("Jaccobi.txt","w");
if(fp1==NULL)
printf("Jaccobi.txt文件打开失败！\n");
else
{
for(i=2;i<2*n_node;i++)
{
for(j=2;j<2*n_node;j++)
fprintf(fp1,"%8.3lf",gJaccobi[i][j]);
fprintf(fp1,"\n");
}
}
//printf("Jaccobi矩阵搞定^-^\n");
}
//高斯消元，参考刘汝佳算法竞赛入门经典训练指南154页
void GetRevised()
{
int i,j,k,r,n=n_node*2-2;
double tmp,f,A[N*2-2][N*2-1];//存放增广矩阵
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
A[i][j]=gJaccobi[i+2][j+2];
A[i][n]=gDelta_PQ[i+2];
}
for(i=0;i<n;i++)
{
r=i;
for(j=i+1;j<n;j++)
if(fabs(A[j][i])>fabs(A[r][i]))
r=j;
if(r!=i)
{
for(j=0;j<=n;j++)
{
tmp=A[r][j];
A[r][j]=A[i][j];
A[i][j]=tmp;
}
}
//i+1~n行与第i行消元
for(k=i+1;k<n;k++)
{
f=A[k][i]/A[i][i];
for(j=i;j<=n;j++)
A[k][j]-=f*A[i][j];
}
}
//回代过程
for(i=n-1;i>=0;i--)
{
for(j=i+1;j<n;j++)
A[i][n]-=A[j][n]*A[i][j];
A[i][n]/=A[i][i];
}
for(i=0;i<n;i++)
gDelta_fe[i+2]=A[i][n];
FILE *fp;
fp=fopen("Delta_fe.txt","w");
if(fp==NULL)
printf("Delata_fe.txt打开失败\n");
for(i=2;i<n+2;i+=2)
{
fprintf(fp,"Delta_f%d=%.9lf\t Delta_e%d=%.9lf\n",i/2+1,gDelta_fe[i],i/2+1,gDelta_fe[i+1]);
}
}

void GetNewValue()
{
int i;
for(i=1;i<n_node;i++)
{
gf[i]-=gDelta_fe[i*2];
ge[i]-=gDelta_fe[i*2+1];
}
}

int is_ok() //返回值为1需要继续迭代，返回值为0则不需要继续迭代
{
int i;
double tmp=0;
for(i=0;i<n_node*2;i++)
{
if(fabs(gDelta_fe[i])>tmp)
tmp=fabs(gDelta_fe[i]);
}
if(fabs(tmp)>Limit)
return 1;
return 0;
}

void GetAns()
{
int i,j;
fushu tmpf,a;
double tmp,tmq;
for(i=0;i<n_node;i++)
{
nodeans[i].id=gBus[i].id;
nodeans[i].type=gBus[i].type;
if(nodeans[i].type==0)
{
nodeans[i].e=gBus[i].noda;
nodeans[i].f=gBus[i].nodb;
a.p=gBus[i].noda;
a.q=gBus[i].nodb;
tmpf=Getvg(a);
nodeans[i].v=tmpf.p;
nodeans[i].g=tmpf.q;
tmp=0;
tmq=0;
for(j=0;j<n_node;j++)
{
tmp+=ge[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])+gf[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
tmq+=gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
}
nodeans[i].p=tmp;
nodeans[i].q=tmq;
}
else if(nodeans[i].type==1)
{
nodeans[i].e=ge[i];
nodeans[i].f=gf[i];
nodeans[i].p=gBus[i].noda;
nodeans[i].q=gBus[i].nodb;
a.p=ge[i];
a.q=gf[i];
tmpf=Getvg(a);
nodeans[i].v=tmpf.p;
nodeans[i].g=tmpf.q;
}
else if(nodeans[i].type==2)
{
nodeans[i].e=ge[i];
nodeans[i].f=gf[i];
a.p=ge[i];
a.q=gf[i];
tmpf=Getvg(a);
nodeans[i].v=tmpf.p;
nodeans[i].g=tmpf.q;
nodeans[i].p=gBus[i].noda;
tmq=0;
for(j=0;j<n_node;j++)
tmq+=gf[i]*(gY_G[i][j]*ge[j]-gY_B[i][j]*gf[j])-ge[i]*(gY_G[i][j]*gf[j]+gY_B[i][j]*ge[j]);
nodeans[i].q=tmq;
}
}
}

void CalFlow()
{
int i,j,k;
double tmpg,tmpb,tmp,x,r,kb;
double B[N],G[N];//对地导纳统计
for(i=0;i<n_node;i++)
B[i]=G[i]=0;
for(i=0;i<n_line;i++)
{
if(gLine[i].type==2)
{
B[gLine[i].start-1]+=gLine[i].bk;
B[gLine[i].ends-1]+=gLine[i].bk;