电力系统潮流计算程序实现C语言版(动态节点+稀疏技术)
2015-07-24 12:02
531 查看
概述
电力系统潮流计算程序采用的是牛顿-拉夫逊法(直角坐标)潮流计算程序的系统节点数是可以动态变化的,根据节点大小分配储存空间
为了减少储存空间和加快计算速度,根据电力系统节点导纳矩阵稀疏的特点,计算程序运用了稀疏储存技术,详细讲解请见参考目录。
关于本程序的计算原理和流程图均可参考后面参考目录中的图书
源数据格式说明
源数据有功功率、无功功率、电压、阻抗、感抗、对地充电容纳均以标幺值表示。数据文件必须命名为data.txt且与潮流计算程序放置于同一个文件中。数据文件data.txt包含两类参数:节点参数和线路参数。节点数据块与线路数据块之间用数字0作为间隔,即在节点数据块结束后,另起一行输入0,然后再在后面按格式要求录入线路参数。
1. 节点参数
节点参数包括:系统节点数Node、节点功率(有功P、无功Q),节点电压(有效值V、相角Delta) ,参数组织格式:
节点数Node 节点数Node写在参数文件的开头,如:4表明为四节点系统。
功率和电压P/Q/V/Delta
首先给出节点参数示例: 2 3 0.5 0 1.10 0 第一列数字2表明该节点的类型为2-PV节点;第二列数字3表明该行数据为节点3的参数;后面三列依次为P、Q、V的给定值,给定值为0,表明该项参数未知;第六列为相角δ,非平衡节点的δ即为PQ迭代的初始相角值,平衡节点的即为给定的相角值。 节点功率为各支路输入功率之和,且规定功率流入节点为正,流出为负。故而负荷功率为负值,发电机功率为正值。
节点类型的判断 根据节点的给定参数可以将节点分为三种类型:
1)PQ节点:给定P、Q初始值的节点,用数字代码1表示;
2)PV节点:给定P、V初始值的节点,用数字代码2表示;
3)平衡节点:没有给定P、Q初始值,仅给定V的初始值,用数字代码3表示;
注意:数据文件中,节点顺序为PQ、PV、平衡节点
2.线路参数
线路参数包括:线路阻抗(R,X)、变压器变比k和对地充电(接地支路)容纳b。每一条线路包括节点号在内共有6个参数,6个参数缺一不可。
线路阻抗Z
首先给出线路参数示例: 1 4 0.12 0.5 1 0.01920 前两列代表线路两端的节点编号,亦即线路编号14(即41);第三四列代表线路阻抗Z=0.12+j0.5;第五列为变压器变比k;最后一列为线路的对地充电电容的一半,即B/2。 k不为0。当k=1则表明该支路为普通支路;否则该支路为变压器支路。普通支路没有任何特殊要求,但对于变压器支路,有以下注意事项。
变压器支路
首先给出变压器线路示例: 3 1 0 0.3 0.909091 0
1)变比k为转换为变压器支路的标准等值电路(如下图)后的变比。
2)变压器线路的编号31有特定含义:3对应节点p,1对应节点q。即变压器支路的编号对节点顺序有要求,p节点编号在前,q节点在后。
3)q节点为全系统参考电压侧。
程序智能化读取线路数据,用户无需刻意对线路参数进行归类、排序,可随意输入,只要每行数据都各自符合格式要求即可。
输出数据说明
1. 节点参数表节点参数表实际上只是对源数据的重现,用于对源数据的整理以及检验程序是否成功读取原始数据。 其中,节点类型列中,BS代表平衡节点、PQ代表PQ节点、PV代表PV节点。 参数表后紧随的是对各类节点数目的统计结果。 该表属于对原始数据的基本统计处理。
2. 节点导纳矩阵
节点导纳矩阵以上三角方式按照节点编号顺序输出,括号内部逗号前为导纳Gij、逗号后为容纳Bij。下三角部分可由上三角部分对称得到。
3. 潮流计算结果(节点)
该部分计算结果包含PQ迭代次数k和节点参数计算结果。迭代次数k在65535次以内均可正常显示。一般PQ迭代不会达到这个数值,否则可以认为PQ迭代是发散的。
节点参数计算结果表,展示了PQ迭代收敛时的计算结果,包括非平衡节点的电压相角δ、PQ节点的电压V、PV节点的无功功率Q和平衡节点的功率。
4. 潮流计算结果(线路)(该部分暂未编写,请等待后续完善)
线路功率计算结果表,展示了PQ迭代收敛时的线路功率计算结果S,包括线路有功功率P和无功功率Q。
参考文献
陈亚明.电力系统计算程序.北京:水利电力出版社,1995张伯明,陈寿孙,严正.高等电力网络分析(第2版).北京:清华大学出版社,2007
源代码
/* FUNCTION : POWER FLOW */ /* WRITTEN BY : ZAHGN CHUN */ /* LAST EDITED : 2014-12-10 */ #include <stdio.h> #include <math.h> #include<stdlib.h> #include<string.h> #include<conio.h> #include<iostream> #include<vector> using namespace std; /*** 子函数声明 ***/ void in_data(); /* 读取数据 */ void Y_Form(); /* 生成导纳矩阵 */ void chuzhi(); /* 给节点电压赋初值 */ void B_Form(); /* 生成雅克比矩阵和求解修正方程*/ void xinzhi(); /*计算新值*/ void PrtNode(); /*输出潮流计算结果*/ void PrtY(); /*打印Y矩阵*/ /*** 全局变量声明 ***/ int Node; /*网络节点数*/ int Num; /*网络支路数*/ int i,j,k,k1,k2,k5,l,m,ig,i1,i2,i3,i5,t,s; /*通用下标值*/ int J1,J0,J3; float tmp; /*临时数据暂存*/ FILE *in,*out; /*输入、输出文件指针*/ int J5; //雅克比矩阵第i行非零元素的技术单元 int k0=1; //导纳矩阵非零互导纳元素的指针 int FL=0; //导纳矩阵非零互导纳元素的个数 float a1,b1,r1,x1,a,r,x,b; //寄存单元 float gij,bij,gi,bi,gj,bj; float A3; //第i节点电流的实部 float B3; //第i节点电流的虚部 float P2,Q2,P3,Q3,C3,D3,C5,D5,A5,B5,C6; float PG,QG,PL,QL; //全网功率 float E,F,PI; int N0; //平衡节点的节点号 float U0; //存平衡节点的给定电压值 int IQ; //发电机台数 int IP; //负荷个数 int N1; //PV节点个数 float UP=1;; //PQ节点的电压初值 int it; //迭代次数 float am=0; //节点功率误差的最大绝对值 int i0; //节点功率误差绝对值最大的节点号 //声明变长数组 vector<int> IZA(Num); //存支路状态数 vector<int> IZ1(Num); //存支路一端的节点号 vector<int> IZ2(Num); //存支路另一端的节点号 vector<float> Z1(Num); //存支路正序电阻 vector<float> Z2(Num); //存支路正序电抗 vector<float> Z3(Num); //存支路正序电纳或非标准变比 vector<int> IWGA(Node); //存发电机状态数 vector<int> IWG(Node); //存发电机节点号 vector<float> W1(Node); //存发电机有功功率 vector<float> W2(Node); //存发电机无功功率 vector<int> ILP(Node); //存负荷状态数 vector<int> ILD(Node); //存负荷节点号 vector<float> WL1(Node); //存负荷有功功率 vector<float> WL2(Node); //存负荷无功功率 vector<int> IPV(5); //存PV节点的节点号 vector<float> PV(5); //存PV节点的给定电压值 vector<float> YZ1(Node*Node); //不规则互导纳元素的实部 vector<float> YZ2(Node*Node); //不规则互导纳元素的虚部 vector<int> IY1(Node*Node); //不规则互导纳元素的行号 vector<int> IY2(Node*Node); //不规则互导纳元素的列号 vector<float> D1(Node); //导纳矩阵对角元素的实部 vector<float> D2(Node); //导纳矩阵对角元素的虚部 vector<float> Y1(Node*Node); //导纳矩阵非对角非零元素的实部 vector<float> Y2(Node*Node); //导纳矩阵非对角非零元素的虚部 vector<int> IY(Node*Node); //导纳矩阵非对角非零元素的列号 vector<int> IN(Node); //导纳矩阵每行非对角非零元素的个数 vector<float> U1(Node); //节点电压的实部 vector<float> U2(Node); //节点电压的虚部 vector<float> AK1(Node*Node); //存雅克比矩阵一行的非零元素 vector<float> AK2(Node*Node); vector<float> AK3(Node*Node); vector<float> AK4(Node*Node); vector<int> JK(Node*Node); //存雅克比矩阵一行的非零元素的列号 vector<float> AJ1(Node*Node); //存雅克比矩阵消去规格化后上三角非对角非零元素 vector<float> AJ2(Node*Node); vector<float> AJ3(Node*Node); vector<float> AJ4(Node*Node); vector<int> IJ(Node*Node); //存雅克比矩阵消去规格化后上三角非对角非零元素的列号 vector<int> JF(Node); //存雅克比矩阵消去规格化后上三角每行非对角非零元素的个数 vector<float> CK1(Node); //先存节点无功功率误差ΔQi(或ΔUi的平方),后存节点电压修正量的实部 vector<float> CK2(Node); //先存节点有功功率误差ΔPi,后存节点电压修正量的虚部 vector<float> PD(Node); //节点负荷的有功功率 vector<float> QD(Node); //节点负荷的无功功率 vector<float> PF(Node); //节点给定的发电有功功率 vector<float> QF(Node); //节点给定的发电无功功率 vector<float> GP(Node); //节点实际的发电有功功率 vector<float> GQ(Node); //节点实际的发电无功功率 vector<float> Q(Node); //节点的无功功率误差值 vector<int> IVI(Node); //PV节点的标致 /****************** 主函数 ******************/ /**↓****↓****↓** 主函数 **↓****↓****↓**/ int main(void) { if((in=fopen("Data.txt","r"))==NULL) printf("wrong\n"); if((out=fopen("out.txt","w"))==NULL) printf("wrong\n"); in_data(); /* 读取数据 */ for(i=0;i<Node*Node;i++) { Y1[i]=0; Y2[i]=0; IY[i]=0; YZ1[i]=0; YZ2[i]=0;IY1[i]=0; AK1[i]=AK2[i]=AK3[i]=AK4[i]=0; AJ1[i]=AJ2[i]=AJ3[i]=AJ4[i]=0; JK[i]=JF[i]=IJ[i]=0; } for(i=0;i<Node;i++) { D1[i]=D2[i]=0; } Y_Form(); PrtY(); chuzhi(); for(it=1;it<20;it++) { B_Form(); fprintf(out,"第%d迭代节点功率误差绝对值最大的节点号:%d,误差的最大绝对值:%f\n",it,i0,am); for(i=0;i<Node-1;i++) fprintf(out,"e[%d]=%f f[%d]=%f\n",i,CK1[i],i,CK2[i]); printf("IT=%d\n",it); if(am<0.0001) break; xinzhi(); } PrtNode(); /*输出潮流计算结果*/ fclose(out); system("cmd /c start out.txt"); return(0); } /**↑****↑****↑** 主函数 **↑****↑****↑**/ /**************** 计算新值 ****************/ void xinzhi() { for(i=Node-2;i>=0;i--) //解修正方程回代运算 //zhuyi!!!!!!!!!!!!!!!! { for(ig=0;ig<JF[i];ig++) { k--; m=IJ[k]; CK1[i]=CK1[i]-AJ1[k]*CK1[m]-AJ2[k]*CK2[m]; CK2[i]=CK2[i]-AJ3[k]*CK1[m]-AJ4[k]*CK2[m]; } } for(i=0;i<Node;i++) //计算节点新的电压值 { U1[i]=U1[i]+CK1[i]; U2[i]=U2[i]+CK2[i]; } for(i=0;i<Node-1;i++) printf("e[%d]=%f f[%d]=%f\n",i,U1[i],i,U2[i]); } /******** 生成雅克比矩阵和不平衡量 *********/ void B_Form() { am=0; k0=0; for(i=0;i<Node;i++) { a=D1[i]*U1[i]; b=D2[i]*U1[i]; r=D1[i]*U2[i]; x=D2[i]*U2[i]; A3=a-x; B3=r+b; J5=1; for(ig=0;ig<IN[i];ig++) { j=IY[k0]; A3=A3+Y1[k0]*U1[j]-Y2[k0]*U2[j]; B3=B3+Y1[k0]*U2[j]+Y2[k0]*U1[j]; if(i!=N0&&j!=N0) { JK[J5]=IY[k0]; AK1[J5]=Y1[k0]*U2[i]-Y2[k0]*U1[i]; AK3[J5]=Y1[k0]*U1[i]+Y2[k0]*U2[i]; AK2[J5]=-AK3[J5]; AK4[J5]=AK1[J5]; J5++; } k0++; } if(i==N0) { GQ[i]=A3*U2[i]-B3*U1[i]; GP[i]=A3*U1[i]+B3*U2[i]; CK1[i]=0; CK2[i]=0; JF[i]=0; } else { P2=PD[i]+(A3*U1[i]+B3*U2[i]); Q2=QD[i]+(A3*U2[i]-B3*U1[i]); GP[i]=P2; GQ[i]=Q2; P2=PF[i]-P2; Q2=QF[i]-Q2; P3=P2; Q3=Q2; if(IVI[i]==1) if(Q[i]>=0||it-1<5) Q3=0; AK3[0]=A3+a+x; AK4[0]=B3-b+r; JK[0]=i; CK2[i]=P2; if(IVI[i]==1) { Q[i]=Q2; if(it>=5) { if(Q2>=0) { for(j=1;j<J5;j++) AK1[j]=AK2[j]=0; AK1[0]=2*U1[i]; AK2[0]=2*U2[i]; Q2=0; CK1[i]=Q2; } else { AK1[0]=r-b-B3; AK2[0]=A3-a-x; CK1[i]=Q2; } } else { for(j=1;j<J5;j++) AK1[j]=AK2[j]=0; AK1[0]=2*U1[i]; AK2[0]=2*U2[i]; Q2=0; CK1[i]=Q2; } } else { AK1[0]=r-b-B3; AK2[0]=A3-a-x; CK1[i]=Q2; } C5=fabs(P3); D5=fabs(Q3); if(C5>am) { i0=i;am=C5; } if(D5>am) { i0=i;am=D5; } k=0; //消去运算 m=0; for(i1=0;i1<i;i1++) { for(i3=1;i3<J5;i3++) if(JK[i3]==i1) break; if(i3<J5) { // printf("xiaoq\n"); for(ig=0;ig<JF[i1];ig++) { for(i2=0;i2<J5;i2++) if(JK[i2]==IJ[k]) break; if(i2<J5) { AK1[i2]=AK1[i2]-AK1[i3]*AJ1[k]-AK2[i3]*AJ3[k]; AK2[i2]=AK2[i2]-AK1[i3]*AJ2[k]-AK2[i3]*AJ4[k]; AK3[i2]=AK3[i2]-AK3[i3]*AJ1[k]-AK4[i3]*AJ3[k]; AK4[i2]=AK4[i2]-AK3[i3]*AJ2[k]-AK4[i3]*AJ4[k]; } else { AK1[J5]=-AK1[i3]*AJ1[k]-AK2[i3]*AJ3[k]; AK2[J5]=-AK1[i3]*AJ2[k]-AK2[i3]*AJ4[k]; AK3[J5]=-AK3[i3]*AJ1[k]-AK4[i3]*AJ3[k]; AK4[J5]=-AK3[i3]*AJ2[k]-AK4[i3]*AJ4[k]; JK[J5]=IJ[k]; J5++; } k++; } CK1[i]=CK1[i]-AK1[i3]*CK1[i1]-AK2[i3]*CK2[i1]; CK2[i]=CK2[i]-AK3[i3]*CK1[i1]-AK4[i3]*CK2[i1]; } else { k=k+JF[i1]; } } k5=k; //规格化 a=AK1[0]*AK4[0]-AK2[0]*AK3[0]; if(a==0) { CK1[i]=0; CK2[i]=0; JF[i]=0; } else { A3=AK4[0]/a; B3=-AK2[0]/a; C3=-AK3[0]/a; D3=AK1[0]/a; for(j=1;j<J5;j++) if(JK[j]>i) { // printf("guife\n"); IJ[k]=JK[j]; AJ1[k]=A3*AK1[j]+B3*AK3[j]; AJ2[k]=A3*AK2[j]+B3*AK4[j]; AJ3[k]=C3*AK1[j]+D3*AK3[j]; AJ4[k]=C3*AK2[j]+D3*AK4[j]; k++; } A5=CK1[i]; B5=CK2[i]; CK1[i]=A3*A5+B3*B5; CK2[i]=C3*A5+D3*B5; JF[i]=k-k5; } } } } /**************** 子函数:读数据 ****************/ void in_data() { fscanf(in,"%d %d %d %d %d",&Node,&Num,&IQ,&IP,&N1); /*读取节点数Node、支路数、发电机个数、负荷个数、PV节点个数*/ for(i=0;i<Num;i++) { Z1[i]=Z2[i]=Z3[i]=0; } for(i=0;i<IQ;i++) { IWGA[i]=IWG[i]=W1[i]=W2[i]=0; } for(i=0;i<Num;i++) fscanf(in,"%d %d %d %f %f %f",&IZA[i],&IZ1[i],&IZ2[i],&Z1[i],&Z2[i],&Z3[i]); for(i=0;i<IQ;i++) fscanf(in,"%d %d %f %f",&IWGA[i],&IWG[i],&W1[i],&W2[i]); for(i=0;i<IP;i++) { fscanf(in,"%d %d %f %f",&ILP[i],&ILD[i],&WL1[i],&WL2[i]); } for(i=0;i<N1;i++) fscanf(in,"%d %f",&IPV[i],&PV[i]); fscanf(in,"%d %f",&N0,&U0); } /****************给节点电压赋初值****************/ void chuzhi() { for(i=0;i<Node;i++) { PD[i]=0; QD[i]=0; PF[i]=0; QF[i]=0; IVI[i]=0; } for(i=0;i<IP;i++) { if(ILP[i]!=0) { j=ILD[i]; PD[j]=WL1[i]; QD[j]=WL2[i]; } } for(i=0;i<IQ;i++) { if(IWGA[i]!=0) { j=IWG[i]; PF[j]=W1[i]; QF[j]=W2[i]; } } for(i=0;i<N1;i++) { j=IPV[i]; IVI[j]=1; } for(i=0;i<Node;i++) { U1[i]=UP; U2[i]=0; } for(i=0;i<N1;i++) { j=IPV[i]; U1[j]=PV[i]; } U1[N0]=U0; } /************形成节点导纳矩阵Y ************/ void Y_Form() { for(i=0;i<Node;i++) { D1[i]=0; D2[i]=0; } l=0; for(k=0;k<Num;k++) { ig=IZA[k]; i=IZ1[k]; j=IZ2[k]; r=Z1[k]; x=Z2[k]; b=Z3[k]; a=r*r+x*x; if(a!=0) { gij=r/a;bij=-x/a; if(ig==1) { gi=gij; gj=gij; bi=bij+b; bj=bi; D1[i]=D1[i]+gi; D2[i]=D2[i]+bi; D1[j]=D1[j]+gj; D2[j]=D2[j]+bj; YZ1[l]=-gij; YZ2[l]=-bij; IY1[l]=i; IY2[l]=j; l++; YZ1[l]=-gij; YZ2[l]=-bij; IY1[l]=j; IY2[l]=i; l++; } else { if(ig==2||ig==3) { gj=gij; bj=bij; gij=gij/b; bij=bij/b; gi=gij/b; bi=bij/b; D1[i]=D1[i]+gi; D2[i]=D2[i]+bi; D1[j]=D1[j]+gj; D2[j]=D2[j]+bj; YZ1[l]=-gij; YZ2[l]=-bij; IY1[l]=i; IY2[l]=j; l++; YZ1[l]=-gij; YZ2[l]=-bij; IY1[l]=j; IY2[l]=i; l++; } else { D1[i]=D1[i]+gij; D2[i]=D2[i]+bij; } } } } j=0; k0=0; m=0; //出口标志 for(i=0;i<Node;i++) { J1=0; for(k=0;k<l;k++) { if(IY1[k]==i) { J3=IY2[k]; for(k1=0;k1<J1;k1++) { k2=k0+k1; if(J3==IY[k2]) { break; m=1; } } if(m==1) { Y1[k2]=Y1[k2]+YZ1[k]; Y2[k2]=Y2[k2]+YZ2[k]; } else { Y1[j]=YZ1[k]; Y2[j]=YZ2[k]; IY[j]=J3; J1++; j++; } } } IN[i]=J1; k0=k0+J1; } FL=k0; } /*******************打印导纳矩阵**************/ void PrtY() { fprintf(out,"非对角线元素:\n"); for(i=0;i<FL;i++) { fprintf(out,"%10.6f,%10.6f,,,,,,%d\n",Y1[i],Y2[i],IY[i]); } fprintf(out,"\n"); fprintf(out,"对角线元素:\n"); for(i=0;i<Node;i++) { fprintf(out,"%10.6f,%10.6f\n",D1[i],D2[i]); } fprintf(out,"\n"); fprintf(out,"每行非对角线非零元素个数:\n"); for(i=0;i<Node;i++) { fprintf(out,"%d\n",IN[i]); } fprintf(out,"\n"); } /***************输出潮流计算结果****************/ void PrtNode() { GP[N0]=GP[N0]+PD[N0]; //输出节点信息 GQ[N0]=GQ[N0]+QD[N0]; PG=0; QG=0; PL=0; QL=0; PI=180/3.14159; fprintf(out,"\n"); for(i=0;i<Node;i++) { E=U1[i]; F=U2[i]; a=sqrt(E*E+F*F); b=PI*atan2(F,E); PG=PG+GP[i]; QG=QG+GQ[i]; PL=PL+PD[i]; QL=QL+QD[i]; fprintf(out,"第%d节点:\n幅值:%f,相角:%f\n",i+1,a,b); } // fprintf(out,"全网发电有功功率:%f,无功功率:%f,负荷有功功率:%f,负荷无功功率:%f\n",PG,QG,PL,QL); //输出支路信息 float ploss=0,qloss=0; //全网有功、无功功率损耗 float qb=0; //全网输电线充电功率 }