潮流计算代码c

合集下载

电力系统潮流计算C语言程序及说明

程序的稳定性分析
程序在不同计算机上的运行结果是否一致。
程序运行过程中，输入数据的变化对输出结果的影响程度。
程序在长时间运行过程中，输出结果是否保持稳定。
程序在处理异常情况时的表现和稳定性。
程序的扩展性分析
代码可读性：C语言程序应具备良好的可读性，方便后续维护和修改算法效率：C语言程序应采用高效的算法，提高计算速度内存占用：C语言程序应合理利用内存，避免内存泄漏和不必要的内存占用扩展性：C语言程序应具备良好的扩展性，方便添加新功能和优化性能
THANK YOU
汇报人：XX
程序的异常处理说明
异常类型：输入错误、计算错误、内存不足等
异常处理方式：使用try-catch 语句进行异常捕获和处理
异常处理流程：当异常发生时，程序会输出错误信息并终止运行
异常处理结果：确保程序在遇到异常时能够正确处理并给出相应的提示信息
C语言程序应用示例
示例程序的输入数据格式
添加标题
添加标题
添加标题Βιβλιοθήκη 输入输出函数：用于数据的输入和输出
函数：可重复使用的代码块，具有特定的功能
C语言程序中电力系统模型的建立
定义节点和支路：根据电力系统网络结构，定义节点和支路，为潮流计算做准备。
建立数学模型：根据电力系统的物理特性和元件参数，建立数学模型，包括节点电压、支路电流等。
实际运行时间测试
程序的内存占用性能分析
内存占用情况：分析程序运行过程中内存的占用情况，包括堆内存和栈内存。
内存泄漏检测：检查程序是否存在内存泄漏，即程序运行结束后未正确释放的内存。
内存使用优化：根据内存占用情况，优化程序中的数据结构或算法，降低内存占用。

电力系统潮流计算完整c语言程序(含网损计算的最终版)

for(j=0;j<Bus_Num;j++)
{
ia[i]=ia[i]+gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j];
ib[i]=ib[i]+gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j];
}
}
for(i=0,n=1;i<Bus_Num-1;i++,n++)
{
gDelta_PQ[2*i]=gDelta_P[i];
gDelta_PQ[2*i+1]=gDelta_Q[i];
}
if((fp=fopen("C:\\Documents and Settings\\Zorro\\桌面\\1\\data\\unbalance.txt","w"))==NULL)
if(gBus[n].Type==1)
gDelta_Q[i]=gDelta_Q[i]-gf[n]*(gY_G[n][j]*ge[j]-gY_B[n][j]*gf[j])+ge[n]*(gY_G[n][j]*gf[j]+gY_B[n][j]*ge[j]);
}
}
for(i=0;i<Bus_Num-1;i++)
{
gY_G[i][j]=0.0;
gY_B[i][j]=0.0;
}
for(l=0;l<Line_Num;l++)
{
i=gLine[l].No_I-1;
j=gLine[l].No_J-1;
r=gLine[l].R;
x=gLine[l].X;

潮流计算源程序

潮流计算源程序/***********************************************************文件名称：潮流计算程序.cpp作者：版本：v1.0说明：修改记录：***********************************************************/#include<stdio.h>#include<math.h>#include<stdlib.h>/******************macro definition************************/#define M 100 /*矩阵阶数*/#define N 100 /*迭代次数*//*******************global variables***********************/int i,j,k,a,b,c; /*循环控制变量*/int t,l,g; /*中间变量*/int n, /* 节点数*/m, /* 支路数*/pq, /* PQ节点数*/pv, /* PV节点数*/duidi; /*对地支路数*/double eps; /* 精度*/double aa[M],bb[M],cc[M],dd[M],max,temp, rr,tt; /*中间变量*/double mo,c1,d1,c2,d2,fushuqiujiao; /*复数运算函数的返回值*/ double G[M][M],B[M][M],Y[M][M]; /*节点导纳矩阵中的实部、虚部及其模方值*/double ykb[M][M],D[M],dU[M]; /*雅克比矩阵、不平衡量矩阵*/struct jd /* 节点结构体*/ { int num,s; /* num为节点号，s为节点类型*/double p,q,S,e,f,U,zkj,dp,dq,du,de,df; /*节点有功、无功功率，功率模值，电压纵、横分量，电压模值,阻抗角，牛顿--拉夫逊中功率不平衡量、电压不平衡量*/} jd[M];struct zhl /* 支路结构体*/{ int numb; /*numb为支路号*/ int p1,p2; /*支路的两个节点*/double r,x; /*支路的电阻与电抗*/ } zhl[M];FILE *fp1,*fp2;/***********************************************************函数名称：void ReadData()函数功能：读取数据入口参数：出口参数：备注：***********************************************************/void ReadData(){int h,number;fp1=fopen("input.txt","r");if(fp1==NULL) /*如果输入为空则跳出*/{printf(" can not open file !\n");exit(0);}fscanf(fp1,"%d,%d,%d,%d,%d,%lf\n",&n,&m,&pq,&pv,&duidi,&eps); /*输入节点数,支路数,PQ节点数,PV节点数对地支路数和精度*/j=1;k=pq+1;for(i=1;i<=n;i++) /*输入节点编号、类型、输入功率和电压初值*/{fscanf(fp1,"%d,%d",&number,&h);if(h==1) /*类型h=1是PQ节点*/{fscanf(fp1,",%lf,%lf,%lf,%lf\n",&jd[j].p,&jd[j].q,&jd[j].e,&jd[j].f);jd[j].num=number;jd[j].s=h;j++;}if(h==2) /*类型h=2是pv节点*/{fscanf(fp1,",%lf,%lf\n",&jd[k].p,&jd[k].U);jd[k].num=number;jd[k].s=h;jd[k].q=0;k++;}if(h==3) /*类型h=3是平衡节点*/{fscanf(fp1,",%lf,%lf\n",&jd[n].e,&jd[n].f);jd[n].num=number;jd[n].s=h;}}for(i=1;i<=m;i++) /*输入支路阻抗*/fscanf(fp1,"%d,%d,%d,%lf,%lf\n",&zhl[i].numb,&zhl[i].p1,&zhl[i].p2,&zhl[i].r,&zhl[i].x);fclose(fp1);if((fp2=fopen("output.txt","w"))==NULL){printf(" can not open file!\n");exit(0);}fprintf(fp2,"\n 潮流上机实习华北电力大学电气化0807 段春姝200801000705\n");fprintf(fp2,"\n ********** 原始数据*********\n");fprintf(fp2,"====================================================================== ========\n");fprintf(fp2," 节点数:%d 支路数:%d PQ节点数:%d PV节点数:%d 对地支路数:%d 精度:%f\n",n,m,pq,pv,duidi,eps);fprintf(fp2," ------------------------------------------------------------------------------\n");for(i=1;i<=pq;i++)fprintf(fp2," PQ节点: 节点%d P[%d]=%f Q[%d]=%f\n",jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].q);for(i=pq+1;i<=pq+pv;i++)fprintf(fp2," PV节点: 节点%d P[%d]=%f U[%d]=%f 初值Q[%d]=%f\n",jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].U,jd[i].num,jd[i].q);fprintf(fp2," 平衡节点: 节点%d e[%d]=%f f[%d]=%f\n",jd[n].num,jd[n].num,jd[n].e,jd[n].num,jd[n].f);fprintf(fp2," -------------------------------------------------------------------------------\n");for(i=1;i<=m;i++)fprintf(fp2," 支路%d: 相关节点:%d,%d R=%f X=%f\n",i,zhl[i].p1,zhl[i].p2,zhl[i].r,zhl[i].x);fprintf(fp2,"====================================================================== ========\n");}/***********************************************************函数名称：double fushuqiumo(double a0,double b0)函数功能：复数运算函数入口参数：double a0,double b0出口参数：double备注：***********************************************************/double fushuqiumo(double a0,double b0) /*复数求模值函数*/ { mo=sqrt(a0*a0+b0*b0);return mo;}double fushuqiuji(double a1,double b1,double a2,double b2) /*复数求积函数*/{ c1=a1*a2-b1*b2;d1=a1*b2+a2*b1;return c1;return d1;}double fushuqiushang(double a3,double b3,double a4,double b4) /*复数求商函数*/{ c2=(a3*a4+b3*b4)/(a4*a4+b4*b4);d2=(a4*b3-a3*b4)/(a4*a4+b4*b4);return c2;return d2;}double fushuqiufushufujiao(double a5,double b5) /*复数求幅角函数*/{ if(a5>0)fushuqiujiao=atan(b5/a5);if(a5<0&&b5>=0)fushuqiujiao=3.1415926+atan(b5/a5);if(a5<0&&b5<0)fushuqiujiao=atan(b5/a5)-3.1415926;if(a5==0&&b5>0)fushuqiujiao=3.1415926/2;if(a5==0&&b5<0)fushuqiujiao=-3.1415926/2;if(a5==0&&b5==0)fushuqiujiao=0;return fushuqiujiao;}/***********************************************************函数名称：void Form_Y()函数功能：计算节点导纳矩阵入口参数：出口参数：备注：***********************************************************/void Form_Y(){for(i=1;i<=n;i++) /*节点导纳矩阵元素附初值*/ for(j=1;j<=n;j++)G[i][j]=B[i][j]=0;/*节点导纳矩阵的主对角线上的元素*/for(i=1;i<=m;i++) /*将对地导纳加入相应的主对角线元素中*/{if((zhl[i].p1==0)||(zhl[i].p2==0)){fushuqiumo(zhl[i].r,zhl[i].x);if(mo==0){G[zhl[i].p1][zhl[i].p2]=B[zhl[i].p1][zhl[i].p2]=G[zhl[i].p2][zhl[i].p1]=B[zhl[i].p2][zhl[i].p1] =0;continue;}fushuqiushang(1,0,zhl[i].r,zhl[i].x);if(zhl[i].p1==0){G[zhl[i].p2][zhl[i].p2]+=c2;B[zhl[i].p2][zhl[i].p2]+=d2;}if(zhl[i].p2==0){G[zhl[i].p1][zhl[i].p1]+=c2;B[zhl[i].p1][zhl[i].p1]+=d2;}}}for(i=1;i<=n;i++) /*非对地导纳加入相应的主对角线元素中*/for(j=1;j<=m;j++)if(zhl[j].p1==i||zhl[j].p2==i){fushuqiumo(zhl[j].r,zhl[j].x);if(mo==0) continue;fushuqiushang(1,0,zhl[j].r,zhl[j].x);G[i][i]+=c2;B[i][i]+=d2;}for(k=1;k<=m;k++) /*节点导纳矩阵非主对角线上的元素*/{i=zhl[k].p1;j=zhl[k].p2;fushuqiumo(zhl[k].r,zhl[k].x);if(mo==0) continue;fushuqiushang(1,0,zhl[k].r,zhl[k].x);G[i][j]-=c2;B[i][j]-=d2;G[j][i]=G[i][j];B[j][i]=B[i][j];}/*--------------------------输出节点导纳矩阵------------------------------*/fprintf(fp2,"\n\n ********* 潮流计算过程*********\n"); fprintf(fp2,"\n====================================================================== ========\n");fprintf(fp2,"\n 节点导纳矩阵为:");for(i=1;i<=n;i++){fprintf(fp2,"\n ");for(k=1;k<=n;k++){fprintf(fp2,"%f",G[i][k]);if(B[i][k]>=0){fprintf(fp2,"+j");fprintf(fp2,"%f ",B[i][k]);}else{B[i][k]=-B[i][k];fprintf(fp2,"-j");fprintf(fp2,"%f ",B[i][k]);B[i][k]=-B[i][k];}}}fprintf(fp2,"\n ------------------------------------------------------------------------------\n");}/***********************************************************函数名称：void fuchuzhi()函数功能：赋初值入口参数：出口参数：备注：***********************************************************/void fuchuzhi(){for(i=1;i<=pq;i++) /* e、f赋初值*/{jd[i].e=1; jd[i].f=0;}for(i=pq+1;i<n;i++){jd[i].e=jd[i].U; jd[i].f=0;}}/***********************************************************函数名称：void Calculate_Unbalanced_Value()函数功能：计算不平衡功率和电压入口参数：出口参数：备注：***********************************************************/void Calculate_Unbalanced_Value(){for(i=1;i<=n;i++){if(jd[i].s==1) /*计算PQ节点不平衡量*/{t=jd[i].num;cc[t]=dd[t]=0;for(j=1;j<=n;j++){for(a=1;a<=n;a++){if(jd[a].num==j)break;}fushuqiuji(G[t][j],-B[t][j],jd[a].e,-jd[a].f);cc[t]+=c1;dd[t]+=d1;}fushuqiuji(jd[i].e,jd[i].f,cc[t],dd[t]);jd[i].dp=jd[i].p-c1;jd[i].dq=jd[i].q-d1;}if(jd[i].s==2) /*计算PV节点不平衡量*/ {t=jd[i].num;cc[t]=dd[t]=0;for(j=1;j<=n;j++){for(a=1;a<=n;a++){if(jd[a].num==j)break;}fushuqiuji(G[t][j],-B[t][j],jd[a].e,-jd[a].f);cc[t]+=c1;dd[t]+=d1;}fushuqiuji(jd[i].e,jd[i].f,cc[t],dd[t]);jd[i].dp=jd[i].p-c1;jd[i].q=d1;jd[i].du=jd[i].U*jd[i].U-(jd[i].e*jd[i].e+jd[i].f*jd[i].f);}}for(i=1;i<=pq;i++) /*形成不平衡量矩阵D[M]*/ {D[2*i-1]=jd[i].dp;D[2*i]=jd[i].dq;}for(i=pq+1;i<=g;i++){D[2*i-1]=jd[i].dp;D[2*i]=jd[i].du;}fprintf(fp2,"\n 不平衡量为:");for(i=1;i<=pq;i++){t=jd[i].num;fprintf(fp2,"\n dp[%d]=%f",i,D[2*t-1]);fprintf(fp2,"\n dq[%d]=%f",i,D[2*t]);}for(i=pq+1;i<=g;i++){t=jd[i].num;fprintf(fp2,"\n dp[%d]=%f",i,D[2*t-1]);fprintf(fp2,"\n du[%d]=%f",i,D[2*t]);}}/***********************************************************函数名称：void Form_Jacobi_Matric()函数功能：形成雅克比矩阵入口参数：出口参数：备注：***********************************************************/void Form_Jacobi_Matric(){for(i=1;i<=pq;i++) /*形成pq节点子阵*/for(j=1;j<n;j++){int i1=jd[i].num;int j1=jd[j].num;double ei=jd[i].e;double ej=jd[j].e;double fi=jd[i].f;double fj=jd[j].f;if(i!=j) /*求i!=j时的H、N、J、L*/{ykb[2*i-1][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi; /* H */ykb[2*i-1][2*j]=G[i1][j1]*ei+B[i1][j1]*fi; /* N */ykb[2*i][2*j-1]=-G[i1][j1]*ei-B[i1][j1]*fi; /* J */ykb[2*i][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /* L */}else /*求i=j时的H、N、J、L*/{aa[i]=0;bb[i]=0;for(k=1;k<=n;k++){int k1=jd[k].num;fushuqiuji(G[i1][k1],B[i1][k1],jd[k].e,jd[k].f);aa[i]=aa[i]+c1;bb[i]=bb[i]+d1;}ykb[2*i-1][2*j-1]=-B[i1][i1]*ei+G[i1][i1]*fi+bb[i]; /*H*/ykb[2*i-1][2*j]=G[i1][i1]*ei+B[i1][i1]*fi+aa[i]; /*N*/ykb[2*i][2*j-1]=-G[i1][i1]*ei-B[i1][i1]*fi+aa[i]; /*J*/ykb[2*i][2*j]=-B[i1][i1]*ei+G[i1][i1]*fi-bb[i]; /*L*/}}for(i=pq+1;i<=g;i++) /*形成pv节点子阵*/ for(j=1;j<n;j++){int i1=jd[i].num;int j1=jd[j].num;double ei=jd[i].e;double ej=jd[j].e;double fi=jd[i].f;double fj=jd[j].f;if(i!=j) /*求i!=j时的H、N*/{ykb[2*i-1][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi; /*H*/ykb[2*i-1][2*j]=G[i1][j1]*ei+B[i1][j1]*fi; /*N*/ykb[2*i][2*j-1]=ykb[2*i][2*j]=0; /*R、S*/}else /*求i=j时的H、N、R、S*/{aa[i]=0;bb[i]=0;for(k=1;k<=n;k++){int k1=jd[k].num;fushuqiuji(G[i1][k1],B[i1][k1],jd[k].e,jd[k].f);aa[i]=aa[i]+c1;bb[i]=bb[i]+d1;}ykb[2*i-1][2*j-1]=-B[i1][i1]*ei+G[i1][i1]*fi+bb[i]; /*H*/ykb[2*i-1][2*j]=G[i1][i1]*ei+B[i1][i1]*fi+aa[i]; /*N*/ykb[2*i][2*j-1]=2*fi; /*R*/ykb[2*i][2*j]=2*ei; /*S*/ }}/*------------------------------输出雅克比矩阵--------------------------------*/fprintf(fp2,"\n\n 雅克比矩阵为: ");for(i=1;i<=2*g;i++){fprintf(fp2,"\n");for(j=1;j<=2*g;j++){fprintf(fp2," %f",ykb[i][j]);}}}/***********************************************************函数名称：void Solve_Equations()函数功能：求解修正方程组入口参数：出口参数：备注：用列主元消取法***********************************************************/void Solve_Equations(){for(i=1;i<=2*g;i++) /*把函数参量矩阵编入修正方程*/ ykb[i][2*g+1]=D[i];k=1;do{rr=ykb[k][k];l=k;i=k+1;do{if(fabs(ykb[i][k])>fabs(rr)) /*列选主元*/{rr=ykb[i][k];l=i;}i++;} while(i<=2*g);if(l!=k){for(j=k;j<=2*g+1;j++){tt=ykb[l][j];ykb[l][j]=ykb[k][j];ykb[k][j]=t;}}for(j=k+1;j<=2*g+1;j++) /*消元*/ykb[k][j]/=ykb[k][k];for(i=k+1;i<=2*g;i++)for(j=k+1;j<=2*g+1;j++)ykb[i][j]-=ykb[i][k]*ykb[k][j];k++;} while(k<=2*g);if(k!=1){for(i=2*g;i>=1;i--) /*回代*/{tt=0;for(j=i+1;j<=2*g;j++)tt+=ykb[i][j]*D[j];D[i]=ykb[i][2*g+1]-tt;}}for(i=1;i<=g;i++){jd[i].e+=D[2*i];jd[i].f+=D[2*i-1];}max=fabs(D[1]);for(i=1;i<=2*(pq+pv);i++)if(fabs(D[i])>max)max=fabs(D[i]);}/***********************************************************函数名称：void Niudun_Lafuxun()函数功能：牛顿--拉夫逊入口参数：出口参数：备注:***********************************************************/void Niudun_Lafuxun(){int z=1;fprintf(fp2,"\n --------------牛顿--拉夫逊迭代次数及结果:--------------\n");g=pq+pv;do{ max=1;if((z<N)&&(max>=eps)){fprintf(fp2,"\n 迭代次数: %d\n",z);}Calculate_Unbalanced_Value();Form_Jacobi_Matric();Solve_Equations();fprintf(fp2,"\n\n 输出df,de: ");for(c=1;c<=n;c++){for(a=1;a<=n;a++){if(jd[a].num==c)break;}fprintf(fp2,"\n");fprintf(fp2," 节点为%2d df=%8.5f de=%8.5f",c,D[2*a-1],D[2*a]);}fprintf(fp2,"\n\n 输出迭代过程中的电压值: ");for(c=1;c<=n;c++){for(a=1;a<=n;a++){if(jd[a].num==c)break;}fprintf(fp2,"\n U[%d]=%f",c,jd[a].e);if(jd[a].f>=0)fprintf(fp2,"+j%f",jd[a].f);elsefprintf(fp2,"-j%f",-jd[a].f);}fprintf(fp2,"\n\n ------------------------------------------------------------------------------");z++;} while((z<N)&&(max>=eps));}/***********************************************************函数名称：void Powerflow_Result()函数功能：输出潮流计算结果入口参数：出口参数：备注:***********************************************************/void Powerflow_Result(){int n1=jd[n].num;fprintf(fp2,"\n\n====================================================================== ========\n\n");fprintf(fp2," ******潮流计算结果******");fprintf(fp2,"\n\n====================================================================== ========\n\n");fprintf(fp2,"\n 各节点电压值: ");for(c=1;c<=n;c++){for(a=1;a<=n;a++){if(jd[a].num==c)break;}fprintf(fp2,"\n U[%d]=%f",c,jd[a].e);if(jd[a].f>=0)fprintf(fp2,"+j%f",jd[a].f);elsefprintf(fp2,"-j%f",-jd[a].f);}rr=tt=0;for(i=1;i<=n;i++){int i1=jd[i].num;fushuqiuji(G[n1][i1],-B[n1][i1],jd[i].e,-jd[i].f);rr+=c1;tt+=d1;}fushuqiuji(jd[n].e,jd[n].f,rr,tt);fprintf(fp2,"\n\n 各节点注入功率:\n");for(i=1;i<=pq;i++){int i1=jd[i].num;fprintf(fp2," PQ节点: 节点%d S[%d]=%f", i1,i1,jd[i].p);if(jd[i].q>=0)fprintf(fp2,"+j%f\n",jd[i].q);elsefprintf(fp2,"-j%f\n",-jd[i].q);}for(i=pq+1;i<=g;i++){int i1=jd[i].num;fprintf(fp2," PV节点: 节点%d S[%d]=%f", i1,i1,jd[i].p);if(jd[i].q>=0)fprintf(fp2,"+j%f\n",jd[i].q);elsefprintf(fp2,"-j%f\n",-jd[i].q);}fprintf(fp2," 平衡节点: 节点%d",jd[n].num,jd[n].num);fprintf(fp2," S[%d]=%f",n1,c1);if(d1>=0)fprintf(fp2,"+j%f",d1);elsefprintf(fp2,"-j%f",-d1);fprintf(fp2,"\n\n 线路功率:\n");rr=tt=0;for(i=1;i<=m;i++){int i1=zhl[i].p1;int j1=zhl[i].p2;aa[i]=bb[i]=0;for(a=1;a<=n;a++){if(jd[a].num==i1)break;}for(b=1;b<=n;b++){if(jd[b].num==j1)break;}fushuqiumo(zhl[i].r,zhl[i].x);if(mo==0)continue;fushuqiushang(1,0,zhl[i].r,zhl[i].x);fushuqiuji(jd[a].e-jd[b].e,-jd[a].f+jd[b].f,c2,-d2);fushuqiuji(jd[a].e,jd[a].f,c1,d1);fprintf(fp2," 支路%d: 线路%d--%d 的功率为: %f",i,i1,j1,c1);if(d1>=0)fprintf(fp2,"+j%f\n",d1);elsefprintf(fp2,"-j%f\n",-d1);aa[i]+=c1;bb[i]+=d1;fushuqiuji(jd[b].e-jd[a].e,-jd[b].f+jd[a].f,c2,-d2);fushuqiuji(jd[b].e,jd[b].f,c1,d1);fprintf(fp2," 支路%d: 线路%d--%d 的功率为: %f",i,j1,i1,c1);if(d1>=0)fprintf(fp2,"+j%f\n",d1);elsefprintf(fp2,"-j%f\n",-d1);aa[i]+=c1;bb[i]+=d1;rr+=aa[i];tt+=bb[i];}fprintf(fp2,"\n\n 线路损耗功率:\n");for(i=1;i<=m;i++){int i1=zhl[i].p1;int j1=zhl[i].p2;fprintf(fp2," 支路%d: 线路%d--%d 损耗的功率为: %f",i,i1,j1,aa[i]);if(bb[i]>=0)fprintf(fp2,"+j%f\n",bb[i]);elsefprintf(fp2,"-j%f\n",-bb[i]);}fprintf(fp2,"\n\n 网络总损耗功率为: %f",rr);if(tt>=0)fprintf(fp2,"+j%f\n",tt);elsefprintf(fp2,"-j%f\n",-tt);fprintf(fp2,"\n====================================================================== ========\n");fprintf(fp2,"\n\n ********* 潮流计算结束*********");}/***********************************************************函数名称：void main()函数功能：主函数入口参数：出口参数：备注:***********************************************************/void main(){ReadData(); /*读取数据*/Form_Y(); /*形成节点导纳矩阵*/fuchuzhi(); /*赋初值*/Niudun_Lafuxun(); /*牛顿--拉夫逊迭代*/Powerflow_Result();printf("!!******计算结果存放于output.txt文件中******!!\n");}潮流计算输入文件3,2,2,0,0,0.000101,1,-0.8,0.6,1,02,1,-0.8,0.6,1,03,3,1,01,1,3,0.049,0.22,2,3,0.01,0.02给定题目题号：1183个节点两条支路无对地支路节点个数支路条数PQ节点数PV节点数对地支路数精度要求（PQ节点）节点编号节点类型（1）输入有功输入无功电压实部初值电压虚部初值（PV节点）节点编号节点类型（2）输入有功电压实部初值（平衡节点）节点编号节点类型（3）输入有功电压实部初值（依次列写）支路号支路所连节点（两个）支路电阻支路电抗（依次列写）潮流计算输出程序文件潮流上机实习华北电力大学电气化0807 段春姝200801000705********** 原始数据*********====================================================================== ========节点数:3 支路数:2 PQ节点数:2 PV节点数:0 对地支路数:0 精度:0.000100------------------------------------------------------------------------------PQ节点: 节点1 P[1]=-0.800000 Q[1]=0.600000PQ节点: 节点2 P[2]=-0.800000 Q[2]=0.600000平衡节点: 节点3 e[3]=1.000000 f[3]=0.000000-------------------------------------------------------------------------------支路1: 相关节点:1,3 R=0.049000 X=0.200000支路2: 相关节点:2,3 R=0.010000 X=0.020000====================================================================== ========********* 潮流计算过程*********====================================================================== ========节点导纳矩阵为:1.155633-j4.716870 0.000000+j0.000000 -1.155633+j4.7168700.000000+j0.000000 20.000000-j40.000000 -20.000000+j40.000000-1.155633+j4.716870 -20.000000+j40.000000 21.155633-j44.716870--------------------------------------------------------------------------------------------牛顿--拉夫逊迭代次数及结果:--------------迭代次数: 1不平衡量为:dp[1]=-0.800000dq[1]=0.600000dp[2]=-0.800000dq[2]=0.600000雅克比矩阵为:4.716870 1.155633 0.000000 0.000000-1.155633 4.716870 0.000000 0.0000000.000000 0.000000 40.000000 20.0000000.000000 0.000000 -20.000000 40.000000输出df,de:节点为 1 df=-0.18940 de= 0.08080节点为 2 df=-0.02200 de= 0.00400节点为 3 df= 0.00000 de= 0.00000输出迭代过程中的电压值:U[1]=1.080800-j0.189400U[2]=1.004000-j0.022000U[3]=1.000000+j0.000000------------------------------------------------------------------------------ 迭代次数: 2不平衡量为:dp[1]=-0.049000dq[1]=-0.200000dp[2]=-0.010000dq[2]=-0.020000雅克比矩阵为:4.279116 1.342383 0.000000 0.000000-2.942383 5.479116 0.000000 0.0000000.000000 0.000000 39.120000 20.1600000.000000 0.000000 -21.760000 40.320000输出df,de:节点为 1 df=-0.00000 de=-0.03650节点为 2 df=-0.00000 de=-0.00050节点为 3 df= 0.00000 de= 0.00000输出迭代过程中的电压值:U[1]=1.044298-j0.189400U[2]=1.003504-j0.022000U[3]=1.000000+j0.000000------------------------------------------------------------------------------ 迭代次数: 3不平衡量为:dp[1]=-0.001540dq[1]=-0.006285dp[2]=-0.000005dq[2]=-0.000010雅克比矩阵为:4.279116 1.258017 0.000000 0.000000-2.942383 5.134763 0.000000 0.0000000.000000 0.000000 39.120000 20.1401590.000000 0.000000 -21.760000 40.280317输出df,de:节点为 1 df= 0.00000 de=-0.00122节点为 2 df= 0.00000 de=-0.00000节点为 3 df= 0.00000 de= 0.00000输出迭代过程中的电压值:U[1]=1.043074-j0.189400U[2]=1.003504-j0.022000U[3]=1.000000+j0.000000------------------------------------------------------------------------------ 迭代次数: 4不平衡量为:dp[1]=-0.000002dq[1]=-0.000007dp[2]=-0.000000dq[2]=-0.000000雅克比矩阵为:4.279116 1.255188 0.000000 0.000000-2.942383 5.123217 0.000000 0.0000000.000000 0.000000 39.120000 20.1401490.000000 0.000000 -21.760000 40.280298输出df,de:节点为 1 df=-0.00000 de=-0.00000节点为 2 df=-0.00000 de=-0.00000节点为 3 df= 0.00000 de= 0.00000输出迭代过程中的电压值:U[1]=1.043072-j0.189400U[2]=1.003504-j0.022000U[3]=1.000000+j0.000000------------------------------------------------------------------------------====================================================================== ========******潮流计算结果******====================================================================== ========各节点电压值:U[1]=1.043072-j0.189400U[2]=1.003504-j0.022000U[3]=1.000000+j0.000000各节点注入功率:PQ节点: 节点1 S[1]=-0.800000+j0.600000PQ节点: 节点2 S[2]=-0.800000+j0.600000平衡节点: 节点3 S[3]=1.653525-j1.002193线路功率:支路1: 线路1--3 的功率为: -0.800000+j0.600000支路1: 线路3--1 的功率为: 0.843599-j0.422044支路2: 线路2--3 的功率为: -0.800000+j0.600000支路2: 线路3--2 的功率为: 0.809926-j0.580149线路损耗功率:支路1: 线路1--3 损耗的功率为: 0.043599+j0.177956支路2: 线路2--3 损耗的功率为: 0.009926+j0.019851网络总损耗功率为: 0.053525+j0.197807====================================================================== ========********* 潮流计算结束*********。

C语言潮流计算-牛顿-拉夫逊法(直角坐标)

( k 1)fi(k ) f i
(k )
计算平衡节点功率 S s 和线路功率
~
停止
1 / 20
程序代码如下：
#include<stdio.h> #include<stdlib.h> #include<math.h> #include<string.h> #include<conio.h> struct linetype // 线路参数 { int jiedian[2]; // 若为变压器，则左端为“低”压侧（换算成Π型等值电路）,变压器的阻抗在“低”压侧 double R,X,K,B0; }line[30]; struct Nodetype // 节点功率 { int lei; // PQ 定义 1，PV 定义 2，平衡节点定义 3 int jie; // 节点编号 double P,Q; double y1,y2; // 初始电压 }poin[30]; int point,road; // 节点数 point 支路数 road int p1,p2; // PQ PV 节点数目 //************************************************* 自定义函数 *************************************************************** void chargescreen() // 调节屏幕 { int mode; printf("\t 请选择界面模式： ①. 106*45 ②. 134*45\n\t>>"); a: scanf("%d",&mode); if(mode!=1 && mode!=2) { printf("\n\t 错误，请重新输入...\n\t>>"); goto a; } printf("\n\t"); system("pause"); if(mode==1) system("mode con:cols=106 lines=45"); // 调整屏幕大小 else system("mode con:cols=134 lines=45"); } void pqpv() // 统计 PQ、PV 节点数目

电力系统分析潮流计算C语言编程-pq分解法2

void factor();/*求因子表*/
void solve(float **B,float *X,int N);/*解方程组*/
void PrtNode();/*打印输出节点参数*/
void ErrorMsg(int Flag);/*错误信息*/
int Node;/*节点数*/
int num;/*原始节点序号*/
kp=0;
for(i=0;i<NP;i++)
{
dPi=dP+i;
Yi=*(Y+i)-i;
Dltai=*(Dlta+i);
*dPi=0;
for(j=0;j<Node;j++)
{
temp=Dltai-*(Dlta+j);
if(i>j)*dPi+=*(V+j)*(Pji);
tP=*(V+j)*(Pij);
tP=*(V+i)*Yij.G-tP;
tP*=*(V+i);
tQ=*(V+j)*(Qij);
tQ-=*(V+i)*(Yij.B-Yij.B0);
tQ*=*(V+i);
}
fprintf(out,"S[%d,%d]=(%10.6f,%10.6f)\n",k+1,m+1,-tP,-tQ)
*(num+i)=k;
fscanf(in,"%d",&k);
}
if(NQ+j!=Node)ErrorMsg(4);
fprintf(out,"【节点参数表】\n");

C语言计算潮流程序

节点数：3 支路数：3 计算精度：0.00010支路1：0.0300+j0.09001┠—————□—————┨2支路2：0.0200+j0.09002┠—————□—————┨3支路3：0.0300+j0.09003┠—————□—————┨1节点1：PQ节点，S（1）=-0.5000-j0.2000节点2：PQ节点，S（2）=-0.6000-j0.2500节点3：平衡节点，U（3）=1.0000∠0.0000n=5;nl=5;isb=1;pr=0.00001;B1=[1 2 0.03i 0 1.05 0;2 3 0.08+0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1];B2=[0 0 1.05 1.05 0 1;0 3.7+1.3i 1.05 0 0 2;0 2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0 0 2;5 0 1.05 1.05 0 3];X=[1 0;2 0;3 0;4 0;5 0];na=3;Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n); for i=1:nif X(i,2)~=0;p=X(i,1);Y(p,p)=1./X(i,2);endendfor i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));YI(p,q)=YI(p,q)-1./B1(i,3);Y(q,p)=Y(p,q);YI(q,p)=YI(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;YI(q,q)=YI(q,q)+1./B1(i,3);Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;YI(p,p)=YI(p,p)+1./B1(i,3);endG=real(Y);B=imag(YI);BI=imag(Y);for i=1:nS(i)=B2(i,1)-B2(i,2);BI(i,i)=BI(i,i)+B2(i,5);endP=real(S);Q=imag(S);for i=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfor i=1:nif B2(i,6)==2V(i)=sqrt(e(i)^2+f(i)^2);O(i)=atan(f(i)./e(i));endendfor i=2:nif i==nB(i,i)=1./B(i,i);else IC1=i+1;for j1=IC1:nB(i,j1)=B(i,j1)./B(i,i);endB(i,i)=1./B(i,i);for k=i+1:nfor j1=i+1:nB(k,j1)=B(k,j1)-B(k,i)*B(i,j1);endendendendp=0;q=0;for i=1:nif B2(i,6)==2p=p+1;k=0;for j1=1:nif B2(j1,6)==2k=k+1;A(p,k)=BI(i,j1);endendendendfor i=1:naif i==naA(i,i)=1./A(i,i);else k=i+1;for j1=k:naA(i,j1)=A(i,j1)./A(i,i);endA(i,i)=1./A(i,i);for k=i+1:nafor j1=i+1:naA(k,j1)=A(k,j1)-A(k,i)*A(i,j1);endendendendICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;while ICT2~=0|ICT3~=0ICT2=0;ICT3=0;for i=1:nif i~=isbC(i)=0;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k)));endDP1(i)=P(i)-V(i)*C(i);DP(i)=DP1(i)./V(i);DET=abs(DP1(i));if DET>=prICT2=ICT2+1;endendendNp(K)=ICT2;if ICT2~=0for i=2:nDP(i)=B(i,i)*DP(i);if i~=nIC1=i+1;for k=IC1:nDP(k)=DP(k)-B(k,i)*DP(i);endelsefor LZ=3:iL=i+3-LZ;IC4=L-1;for MZ=2:IC4I=IC4+2-MZ;DP(I)=DP(I)-B(I,L)*DP(L);endendendendfor i=2:nO(i)=O(i)-DP(i);endkq=1;L=0;for i=1:nif B2(i,6)==2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i));if DET>=prICT3=ICT3+1;endendendelse kp=0;if kq~=0;L=0;for i=1:nif B2(i,6)==2C(i)=0;L=L+1;for k=1:nC(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));endDQ1(i)=Q(i)-V(i)*C(i);DQ(L)=DQ1(i)./V(i);DET=abs(DQ1(i));endendendendNq(K)=ICT3;if ICT3~=0L=0;for i=1:naDQ(i)=A(i,i)*DQ(i);if i==nafor LZ=2:iL=i+2-LZ;IC4=L-1;for MZ=1:IC4I=IC4+1-MZ;DQ(I)=DQ(I)-A(I,L)*DQ(L);endendelseIC1=i+1;for k=IC1:naDQ(k)=DQ(k)-A(k,i)*DQ(i);endendendL=0;for i=1:nif B2(i,6)==2L=L+1;V(i)=V(i)-DQ(L);endendkp=1;K=K+1;elsekq=0;if kp~=0K=K+1;endendfor i=1:nDy(K-1,i)=V(i);endenddisp('迭代次数')disp(K);disp('每次没有达到精度要求的有功功率个数为'); disp(Np);disp('每次没有达到精度要求的无功功率个数为'); disp(Nq);for k=1:nE(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;O(k)=O(k)*180./pi;enddisp('各节点的电压标么值E为');disp(E);disp('各节点的电压V大小');disp(V);disp('各节点的电压相角O');disp(O);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率为');disp(S);disp('各条支路的首端功率为');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1. /(B1(i,3)*B1(i,5))));disp(Si(p,q));enddisp('各条支路的末端功率为');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1 ./(B1(i,3)*B1(i,5))));disp(Sj(q,p));enddisp('各条支路的功率损耗为');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);disp(DS(i));endfor i=1:KCs(i)=i;for j=1:nDy(K,j)=Dy(K-1,j);endenddisp('每次迭代后各节点的电压值如图所示'); plot(Cs,Dy)xlabel('迭代次数')ylabel('电压')title('电压迭代次数曲线');题号：2节点数：4 支路数：4 计算精度：0.00010支路1：0.0200+j0.08001┠—————□—————┨3支路2：0.0400+j0.12001┠—————□—————┨4支路3：0.0500+j0.14002┠—————□—————┨4支路4：0.0400+j0.12003┠—————□—————┨4节点1：PQ节点，S(1)=-0.6000-j0.2500节点2：PQ节点，S(2)=-0.8000-j0.3500节点3：PV节点，P(3)=0.4000 V(3)=0.9500节点4：平衡节点，U(4)=1.0000∠0.0000#include "stdio.h"#include "math.h"#define M 20 /*节点数、支路数极限值*/#define N 5 /*迭代次数极限值*/int n,m,dd=0,pq=0,pv=0,tt,qy;float eps; /*节点数、支路数、对地支路数、pq节点数、pv节点数、控制打印、互联网数、精度*/struct jiedian /*节点数据结构体*/{int s; /*节点类型：1-PQ节点；2-PV节点；3-平衡节点*/ float p,q,e,f,v; /*节点的有功、无功、电压实部、虚部、电压辐值*/}jiedian[M]; /*如引用节点1的无功，则为jiedian[1].q */struct zhilu /*支路数据结构体*/{int p1,p2,s; /*支路两端节点号，支路类型：1-普通支路；2-变压器支路；3-对地支路*/float r,x,b,kt; /*支路的电阻、电抗、导纳、变压器非标准变比*/}zhilu[M];struct hulianwang /*互联网结构体*/{int num,pv; /*互联区域号,指定的pv节点号*/float p,eps2; /*规定的有功功率及其允许误差*/int count; /*每个互联区域包括的节点数*/int a[M]; /*包括的节点*/} hulianwang[M];static float G[M][M],B[M][M],G1[M][M],B1[M][M]; /*节点导纳阵矩阵*/static float ykb[2*M][2*M]; /*节点导纳阵矩阵*/static float yinzi[2*M][2*M]; /*因子表*/static float P[M][M],Q[M][M]; /*潮流计算结果*/FILE *fp1,*fp2; /*文件输入、输出指针*/void input() /*从文件in.txt 输入线路基本参数、节点数据、支路数据*/{int i,j,h; /*循环变量、节点（支路）类型*//**********************************打开电网参数文件**********************/fp1=fopen("in.txt","r");if (fp1==NULL){printf("Can not open file in.txt !\n");exit(0);}/***********读取节点数,支路数,互联网数,PQ节点数PV节点数和精度***********/fscanf(fp1,"%d,%d,%d,%f\n",&n,&m,&qy,&eps);/************************读取节点信息*************************************/for(i=1;i<=n;i++){fscanf(fp1,"%d",&h);if(h==1) /*类型h=1是PQ节点*/{pq++;fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].q,&jiedian[i].e,&jiedian[i].f);jiedian[i].s=1;jiedian[i].v=sqrt(jiedian[i].e*jiedian[i].e+jiedian[i].f*jiedian[i].f);}else if(h==2) /*类型h=2是pv节点*/{pv++;fscanf(fp1,",%f,%f,%f,%f\n",&jiedian[i].p,&jiedian[i].v,&jiedian[i].e,&jiedian[i].f);jiedian[i].s=2;jiedian[i].q=0;}else /*类型h=3是平衡节点*/{jiedian[i].p=0;jiedian[i].q=0;jiedian[i].e=1;jiedian[i].f=0;jiedian[i].v=1;jiedian[i].s=3;}}/*******************************读取支路信息*************************************/for(i=1;i<=m;i++){fscanf(fp1,"%d",&h);if(h==1) /*类型h=1是普通支路*/{fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,& zhilu[i].b);zhilu[i].kt=1;zhilu[i].s=1;}if(h==2) /*类型h=2是变压器支路*/{fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,& zhilu[i].kt);zhilu[i].s=2;}if(h==3) /*类型h=3是接地支路*/{fscanf(fp1,",%d,%d,%f,%f,%f\n",&zhilu[i].p1,&zhilu[i].p2,&zhilu[i].r,&zhilu[i].x,& zhilu[i].b);zhilu[i].kt=1;zhilu[i].s=3;dd++;}}/*******************************读互联网信息*************************************/if(qy!=0)for(i=1;i<=qy;i++) /*输入互联网状况*/{fscanf(fp1,"%d,%d,%f,%f,%d",&hulianwang[i].num,&hulianwang[i].pv,&hulian wang[i].p,&hulianwang[i].eps2,&hulianwang[i].count);for(j=1;j<=hulianwang[i].count;j++)fscanf(fp1,",%d",&(hulianwang[i].a[j]));}fclose(fp1);/**********************************打开输出结果文件********************************/fp2=fopen("out.txt","w");if(fp2==NULL){printf("Can not open file!\n");exit(0);}fprintf(fp2,"\n**************************** 原始数据如下***********************************\n\n ");fprintf(fp2," 节点数:%2d 支路数:%2d 对地支路数:%2d PQ节点数:%2d PV节点数:%2d 精度1:%f \n",n,m,dd,pq,pv,eps);fprintf(fp2,"\n------------------------------------------------------------------------------\n\n ");for(i=1;i<=pq;i++)fprintf(fp2," 节点%2d PQ节点P[%d]=%fQ[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].q);for(i=pq+1;i<=pq+pv;i++)fprintf(fp2," 节点%2d PV节点P[%d]=%fV[%d]=%f\n",i,i,jiedian[i].p,i,jiedian[i].v);fprintf(fp2," 节点%2d 平衡节点\n",i);fprintf(fp2,"\n-------------------------------------------------------------------------------\n\n ");for(i=1;i<=m;i++){if(zhilu[i].s==1)fprintf(fp2," 支路%2d：普通支路相关节点:%2d,%2d R=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].b);else if(zhilu[i].s==2)fprintf(fp2," 支路%2d：变压器支路相关节点:%2d,%2dR=%f X=%f Kt=%f\n",i,zhilu[i].p1,zhilu[i].p2,zhilu[i].r,zhilu[i].x,zhilu[i].kt);elsefprintf(fp2," 支路%2d：对地支路相关节点:%2dR=%f X=%f B=%f\n",i,zhilu[i].p1,zhilu[i].r,zhilu[i].x,zhilu[i].b);}for(i=1;i<=qy;i++){fprintf(fp2," 互联网区域:%2d 指定pv节点:%2d 规定有功功率:%8.5f 允许误差:%8.5f 包括的节点:",hulianwang[i].num,hulianwang[i].pv,hulianwang[i].p,hulianwang[i].eps2);for (j=1;j<=hulianwang[i].count;j++)fprintf(fp2,"%2d ",hulianwang[i].a[j]);fprintf(fp2,"\n");}}void youhua() /*利用节点数据、支路数据形成新编号*/ {int i,k,j,jd[N]; /*jd[N]记录节点*/struct jiedian tem; /*中间临时节点所连支路数*/for(i=1;i<=n;i++) /*对节点连接支路数赋0*/{jd[i]=0;}for(i=1;i<=n;i++) /*考虑每个节点*/{for(j=1;j<=m;j++) /*考虑每条支路*/{if((zhilu[j].p1==i||zhilu[j].p2==i)&&(zhilu[j].s!=3)) /*如果是这条支路的节点且不是接地支路*/jd[i]++;}}for(i=1;i<pq;i++) /*对pq节点按所连支路的个数进行重新编号*/{for(j=i+1;j<=pq;j++) /*从小到大排序算法*/{if(jd[i]>jd[j]){tem=jiedian[i];jiedian[i]=jiedian[j];jiedian[j]=tem;for(k=1;k<=m;k++) /*更新支路信息*/{if(zhilu[k].p1==j)zhilu[k].p1=i;else if(zhilu[k].p2==j)zhilu[k].p2=i;else if(zhilu[k].p1==i)zhilu[k].p1=j;else if(zhilu[k].p2==i)zhilu[k].p2=j;}}}}for(i=pq+1;i<pq+pv;i++) /*对pv节点按所连支路的个数进行重新编号*/{for(j=i+1;j<=pq+pv;j++) /*从小到大排序算法*/{if(jd[i]>jd[j]){tem=jiedian[i];jiedian[i]=jiedian[j];jiedian[j]=tem;for(k=1;k<=m;k++) /*更新支路信息*/{if(zhilu[k].p1==j)zhilu[k].p1=i;else if(zhilu[k].p2==j)zhilu[k].p2=i;else if(zhilu[k].p1==i)zhilu[k].p1=j;else if(zhilu[k].p2==i)zhilu[k].p2=j;}}}}fprintf(fp2,"\n**************************** 节点优化结果如下*********************************\n\n ");for(i=1;i<=n;i++){if(jiedian[i].s==1)fprintf(fp2," 节点%2d PQ节点P[%d]=%fQ[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].q);if(jiedian[i].s==2)fprintf(fp2," 节点%2d PV节点P[%d]=%fV[%d]=%f \n",i,i,jiedian[i].p,i,jiedian[i].v);if(jiedian[i].s==3)fprintf(fp2," 节点%2d 平衡节点\n",i);}}void form_y() /*利用支路数据形成Y，注意对地支路*/ {int i,j,k;float S;for(i=1;i<=n;i++)for(j=1;j<=n;j++)G[i][j]=B[i][j]=0;for(i=1;i<=m;i++) /*节点导纳矩阵的主对角线上的导纳*/ for(j=1;j<=n;j++)if(zhilu[i].p1==j||zhilu[i].p2==j){S=zhilu[i].r*zhilu[i].r+zhilu[i].x*zhilu[i].x;if(S==0) continue;G[j][j]+=zhilu[i].r/S;B[j][j]+=-zhilu[i].x/S;if(zhilu[i].s==1) /*如果是普通支路*/{B[j][j]+=zhilu[i].b/2;}else if(zhilu[i].s==2) /*如果是普通变压器支路*/{if(zhilu[i].p1==j){G[j][j]+=(zhilu[i].r/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt);B[j][j]+=(-zhilu[i].x/S*(1-zhilu[i].kt))/(zhilu[i].kt*zhilu[i].kt);}else if(zhilu[i].p2==j){G[j][j]+=(zhilu[i].r/S*(zhilu[i].kt-1))/zhilu[i].kt;B[j][j]+=(-zhilu[i].x/S*(zhilu[i].kt-1))/zhilu[i].kt;}}else if(zhilu[i].s==3) /*如果是对地支路*/{B[j][j]+=zhilu[i].b;}}for(k=1;k<=m;k++) /*节点导纳矩阵非主对角线上的导纳*/ {i=zhilu[k].p1;j=zhilu[k].p2;S=zhilu[k].r*zhilu[k].r+zhilu[k].x*zhilu[k].x;if(S==0) continue;G[i][j]+=-zhilu[k].r/S;B[i][j]+=zhilu[k].x/S;if(zhilu[k].kt!=1.0){G[i][j]/=zhilu[k].kt;B[i][j]/=zhilu[k].kt;}G[j][i]=G[i][j];B[j][i]=B[i][j];}fprintf(fp2,"\n\n****************************** 节点导纳矩阵为**********************************\n");for(i=1;i<=n;i++){fprintf(fp2,"\n ");for(j=1;j<=n;j++)fprintf(fp2,"%8.5f+j%8.5f ",G[i][j],B[i][j]);}}void form_j() /*利用节点数据和Y形成J*/{float ei,fi,a=0,b=0;int i1,j1,k1,i,j,k;for(i=1;i<=2*(pq+pv)+1;i++)for(j=1;j<=2*(pq+pv)+1;j++)ykb[i][j]=0;for(i=1;i<=pq;i++)for(j=1;j<n;j++){ i1=i;j1=j;ei=jiedian[i].e;fi=jiedian[i].f;if(i!=j) /*求i!=j时的H、N、J、L*/ { ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /* H */ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /* N */ykb[2*i][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi; /* J */ykb[2*i][2*j]=-G[i1][j1]*ei-B[i1][j1]*fi; /* L */}else /*求i=j时的H、N、J、K*/ { a=0;b=0;for(k=1;k<=n;k++)if(k!=i){ k1=k;a=a+G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;b=b+G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;}ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/ykb[2*i][2*j]=-2*B[i1][i1]*jiedian[i].f+a; /*L*/ykb[2*i-1][2*j]= 2*G[i1][i1]*jiedian[i].f+b; /*N*/ykb[2*i][2*j-1]=-2*B[i1][i1]*jiedian[i].e-b; /*J*/}}for(i=pq+1;i<=pq+pv;i++) /*形成pv节点子阵*/for(j=1;j<n;j++){ i1=i;j1=j;ei=jiedian[i].e;fi=jiedian[i].f;if(i!=j) /*求i!=j时的H、N*/ {ykb[2*i-1][2*j-1]=G[i1][j1]*ei+B[i1][j1]*fi; /*H*/ykb[2*i-1][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi; /*N*/}else /*求i=j时的H、N、R、S*/{ a=0;b=0;for(k=1;k<=n;k++)if(k!=i){ k1=k;a+=G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;b+=G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;}ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].e+a; /*H*/ykb[2*i-1][2*i]=2*G[i1][i1]*jiedian[i].f+b; /*N*/ykb[2*i][2*j-1]=2*ei; /*R*/ykb[2*i][2*j]=2*fi; /*S*/}}fprintf(fp2,"\n\n\n******************************* 雅可比矩阵*************************************** \n");for(i=1;i<=2*(pq+pv);i++){fprintf(fp2,"\n ");for(j=1;j<=2*(pq+pv);j++)fprintf(fp2,"%8.5f ",ykb[i][j]);}}void form_yz() /*利用J形成因子表*/{int i,j,k;float a[2*M][2*M],L[2*M][2*M],R[2*M][2*M];float x,y,z;for(i=1;i<=2*(pq+pv);i++)for(j=1;j<=2*(pq+pv);j++){ a[i][j]=0;R[i][j]=0; L[i][j]=0;}for(i=1;i<=2*(pq+pv);i++)for(j=1;j<=2*(pq+pv);j++)a[i][j]=ykb[i][j];for(i=1;i<=2*(pq+pv);i++){for(j=1;j<=2*(pq+pv);j++){if(i>j){R[i][j]=0;x=0;for(k=1;k<=j-1;k++){x=x+L[i][k]*R[k][j];}L[i][j]=a[i][j]-x;}if(i==j){R[i][j]=1;y=0;for(k=1;k<=i-1;k++){y=y+L[i][k]*R[k][i];}L[i][i]=a[i][j]-y;}if(i<j){L[i][j]=0;z=0;for(k=1;k<=i-1;k++){z=z+L[i][k]*R[k][j];}R[i][j]=(a[i][j]-z)/L[i][i];}}}for(i=1;i<=2*(pq+pv);i++){for(j=1;j<=2*(pq+pv);j++){if(i<j)yinzi[i][j]=R[i][j];if(i==j)yinzi[i][j]=1/L[i][j];if(i>j)yinzi[i][j]=L[i][j];}}fprintf(fp2,"\n\n\n******************************** 因子表***************************************** \n");for(i=1;i<=2*(pq+pv);i++){fprintf(fp2,"\n ");for(j=1;j<=2*(pq+pv);j++)fprintf(fp2,"%8.5f ",yinzi[i][j]);}}。

电力系统潮流计算C语言程序及说明

实验目的根据所给的电力系统，编制潮流计算程序，通过计算机进行调试，最后完成一个切实可行的电力系统计算应用程序.通过自己设计电力系统计算程序使同学对电力系统分析有进一步理解，同时加强计算机实际应用能力的训练。

程序计算原理1、概述应用计算机进行电力系统计算，首先要掌握电力系统相应计算的数学模型；其次是运用合理的计算方法；第三则是选择合适的计算机语言编制计算程序。

建立电力系统计算的相关数学模型，就是建立用于描述电力系统相应计算的有关参数间的相互关系的数学方程式。

该数学模型的建立往往要突出问题的主要方面，即考虑影响问题的主要因素，而忽略一些次要因素，使数学模型既能正确地反映实际问题，又使计算不过于复杂。

运用合理的计算方法，就是要求所选用的计算方法能快速准确地得出正确结果，同时还应要求在解算过程中占用内存少，以利提高计算机的解题规模。

选择合适的语言编写程序，就是首先确定用什么计算机语言来编制程序；其次是作出计算的流程图；第三根据流程图用选择的语言编写计算程序.然后上机调试，直到语法上无错误。

本程序采用C语言进行编程。

所编制的程序难免存在逻辑错误,因此先用一个已知结果的系统作为例题进行计算.用程序计算的结果和已知结果相比较,如果结果相差甚远就要逐步分析程序的计算步骤，查出问题的出处；如果结果比较接近,则逐步分析误差来源;直到结果正确为止。

2、电力系统潮流计算的程序算法潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角）、网络中的功率分布及功率损耗等.目前计算机潮流计算的方法主要有牛顿—拉夫逊算法和PQ分解法。

牛顿-拉夫逊算法是数学上求解非线形方程组的有效方法,具有较好的收敛性,曾经是潮流计算中应用比较普遍的方法。

PQ快速分解法是从牛顿—拉夫逊算法演变而来的，是将纯数学的牛顿—拉夫逊算法与电力系统具体特点相结合并进行简化与改进而得出的。

PQ快速分解法比牛顿—拉夫逊算法大大提高了计算速度和节省了内存，故而本程序以PQ快速分解法进行潮流计算.1）形成节点导纳矩阵（1)自导纳的形成对节点i其自导纳Y ii是节点i以外的所有节点都接地时节点i对地的总导纳。

电力系统潮流计算的C语言实现

//////////////////////////////////////////////////////////////////////// PQ分解法潮流////文件输入格式：节点总数n(包括联络节点),支路数zls ////节点数(发电机和负荷)nb,接地电抗数mdk,迭代精度eps // //考虑负荷静特性标志kk2(0考虑),平衡节点号，优化标志(0不优化) ////最大迭代次数it1,支路左右节点号izl[],jzl[],支路电阻zr[],电抗zx[] ////支路容纳zyk[],节点号nob[]及标志nobt[](0-PQ,-1-PV) ////发电机和负荷有功、无功pg[],qg[],pl[],ql[] ////电压v0[](pv节点输入实际值，PQ节点任输入一值) // //电抗节点号idk[],电抗值dkk[] ////////////////////////////////////////////////////////////////////////#include "math.h"#include "stdio.h"#define NS 2000 //最大节点数#define NS2 NS * 2#define NS4 1000 //NS4、NS必须大于2*zls。

#define ZS 3000 //最大支路数#define ZS2 ZS * 2#define DKS 200 //最大电抗器数#define N2 ZS * 4#define N3 ZS * 8 + NS * 4FILE *fp1, *fp2;char inname[12], outname[12];// fp1输入数据文件指针fp2输出文件指针// inname[]输入数据文件名outname[]输出数据文件名int n, zls, nb, mdk, mpj, bnsopton, it1, dsd, kk2, nzls;// 节点总数n(包括联络节点) 支路数(回路数)zls 节点数nb(发电机和负荷) // 接地电抗数mdk 精度eps 平衡节点号mpj// 节点优化(标志)bnsopton(=0节点不优化，!=0节点优化)// 最大迭代次数it1 最低电压或最大功率误差节点号dsd// 负荷静特性标志(=0考虑负荷静特性)// 支路数(双回线算一条支路)int izl[ZS], jzl[ZS], idk[DKS], yds[NS], ydz[NS], iy[ZS2];// izl[],jzl[],idk[]:分别存放左、右节点号和电抗器节点号。

C语言进行潮流计算

C语言进行潮流计算C语言是一种通用的编程语言，被广泛用于开发各种应用程序。

虽然它的历史可以追溯到上个世纪70年代，但C语言在当今仍然非常流行。

它被广泛用于系统编程、嵌入式系统开发、游戏开发和科学计算等领域。

本文将重点介绍C语言的潮流计算。

潮流计算是一种用于分析电力系统稳态和动态行为的工具和技术。

它可以用于预测电力系统中潮流（电压和电流）的分布和变化。

潮流计算的目标是找到系统中各个节点的电压和相应的功率流向。

这对于电力系统的规划和运行至关重要。

C语言可以通过编写相应的程序来实现潮流计算。

在这个程序中，需要定义一组节点和支路，以及相应的导纳矩阵（节点之间的电导和电纳）。

导纳矩阵描述了电力系统中节点之间的电压和电流关系。

一旦定义了这些参数，就可以使用C语言编写算法来解决潮流计算问题。

潮流计算通常使用迭代的方法来求解。

在每一次迭代中，系统的节点电压和功率分布都会进行更新，直到达到稳定状态。

C语言提供了丰富的控制结构和数据类型，使得编写这样的迭代算法变得相对容易。

例如，可以使用for循环来实现迭代过程，同时使用浮点数类型来存储节点电压和功率分布。

潮流计算的核心是求解节点电压的方程组。

这个方程组通常是一个非线性方程组，需要使用牛顿-拉夫逊方法或高斯-赛德尔方法等迭代算法进行求解。

C语言提供了一些数值计算库，如GSL（GNU科学库），可以用于实现这些迭代算法。

这些库提供了各种数值计算函数，如求解非线性方程组、矩阵运算和插值等，可以大大简化潮流计算程序的编写。

除了求解节点电压方程组，潮流计算还需要考虑各种系统参数和条件。

例如，支路的参数（如电阻和电感）、发电机的容量和调节策略、负荷的需求和变化等。

C语言可以通过定义和存储这些参数，并在计算过程中引用它们来实现对这些条件的考虑。

此外，潮流计算还需要处理一些特殊情况，如无源系统和有环路系统。

无源系统是指没有外部电压源的电力系统，只由发电机和负荷组成。

有环路系统是指电力系统中存在环路，这会导致方程组无解。

电力系统通用潮流计算C语言程序

＃include <iostream〉#include 〈fstream>＃include<iomanip>#include〈math。

h〉using namespace std；//节点号类型负荷有功负荷无功母线数据（类型1＝PV节点，2＝PQ节点，3＝平衡节点）struct BUS｛int busno；int type；float Pd;float Qd;｝;//发电机数据节点号有功发电电压幅值struct Generator{int busno;float Pg;float Vg；｝;//支路信息节点I 节点J R X B/2 kstruct Line{int busi;int busj;float R;float X；float B;float k；｝；//deltaP deltaQ deltaV^2//void fun1（double YG[][50],double YB［][50]，double e［］，double f［］,int type[］，int N，double W［]，double P［],double Q［］，double V[］）｛double dP=0，dQ=0,dV=0;int i,j;for（i=0；i<N-1；i++)｛double A=0,B=0；for(j=0；j〈N；j++）{A+=YG［i]［j]＊e[j］-YB［i］[j］*f［j];B+=YG［i］［j］＊f[j]+YB［i][j]＊e[j］;}dV=V[i]＊V[i］-e［i]*e［i］—f［i］*f［i];dP=P[i]—e［i］＊A—f［i］＊B；W［2＊i]=dP;dQ=Q[i］-f［i]*A+e［i］＊B;if(type[i］==1)W[2＊i+1］=dQ;else W[2＊i+1]=dV；｝｝//Jacobi矩阵//void Jacobi（double YG［]［50］，double YB[］[50]，double e[50］,double f［50],int type［50]，int N ，double Ja［100］[101]）｛int i，j；for（i=0;i<N;i++)｛for(j=0；j〈N；j++){if(i！=j）{if（type［i]==1){Ja[2*i][2＊j]=—(YG[i]［j］＊e［i］+YB［i][j］＊f［i]）;Ja［2＊i]［2＊j+1]=YB［i]［j]＊e［i]-YG［i］［j］*f[i］;Ja［2＊i+1]［2*j]=Ja［2＊i］［2＊j+1］;Ja［2*i+1］[2*j+1]=—Ja［2*i］[2＊j]；}else {Ja［2*i］[2*j］=—YG［i］［j］*e[i］+YB[i][j］*f［i］;Ja［2＊i][2*j+1]=YB[i］[j］*e[i］-YG［i]［j]＊f［i]；Ja[2＊i+1]［2*j+1]=Ja[2＊i+1］［2＊j]=0；}｝else {double a[50］={0｝，b[50]=｛0};for（int k=0；k<N；k++){a[i]+=(YG[i］［k]*e[k]—YB［i]［k］*f[k］）;b[i]+=（YG[i][k]*f[k]+YB[i］[k]*e［k］）；Ja［2＊i]［2*j］=—a［i]-YG［i]［i］*e［i］-YB[i]［i]＊f［i］；Ja[2*i]［2*j+1］=-b[i］+YB[i]［i］*e[i]-YG［i］［i]＊f[i]；if(type［i]==1）{Ja[2*i+1][2＊j］=b［i]+YB[i］［i]＊e［i］—YG［i］[i]＊f[i]；Ja[2＊i+1][2*j+1]=-a［i]+YG［i]［i]＊e［i]+YB［i］[i］＊f[i];}else ｛Ja[2*i+1］［2*j］=—2＊e[i]；Ja［2*i+1]［2*j+1］=-2＊f[i];}}｝｝}}//高斯消元法解方程组函数//void gauss（double a[］［101］，int n）｛int i，j,k；double c;for（k=0；k〈n-1；k++) ｛c=a[k］[k］；for（j=k；j〈=n；j++) a［k][j］/=c；for（i=k+1；i〈n;i++）｛c=a［i][k］；for(j=k；j<=n；j++）a[i][j］—=c*a［k］[j];｝}a［n-1]［n]/=a[n-1][n-1];for（k=n-2;k>=0;k——）for(j=k+1；j<n；j++) a［k]［n]-=a[k］［j]＊a［j]［n］;}void main（){ifstream fin；int N=0，GS=0，LD=0,ZLs=0; //节点数发电机数负荷数支路数// BUS *B；Generator *G；Line ＊L;//从文本中读入原始数据到数组中//fin。

C语言-_潮流计算实现

C 语言程序设计潮流计算学院自动化学院专业班级学号姓名联系方式本程序潮流计算部分采用牛顿拉夫逊极坐标法进行计算，求解一次多元方程采用高斯列主元分解法进行求解。

根据工程实际，在存储文件实时记录产生文件时间。

此外本程序特意增加了文件查看功能，方便文件的查看。

程序代码：#include<stdio.h>#include<malloc.h>#include<stdlib.h>#include<math.h>#include<string.h>#include<time.h>#include<conio.h>void dataprepare(void);void initial(void);void Yc(void);void showy(void);double detpqc(void);void showdetav(void);void Jrc(void);void showdetpq(void);void Gauss(void);void showj(void);void showsolution(void);char save2file(void);char option(void);int chose(void);int list(char filename[],int i,char dstn[]);struct PQV{char name[5];double vb;double p;double q;double v;double ag;}*pqv;struct Y{double G;double B;}**y;struct LZ{char name1[5];double vb1;char name2[5];double vb2;double r;double x;double b;}*lz;struct TZ{char name1[5];double vb1;char name2[5];double vb2;double x;double k1;double k2;}*tz;char fsource[20];int pqnum=0,pvnum=0,lznum=0,tznum=0;int temp1=0,temp2=0,temp3=0,temp4=0;int diedai=0;int numw=1,numofw=0;double *detpq,*detav,*Pi,*Qi,error,wucha,**J; /* detpq节点功率的误差量，detav修正量，Pi 节点的有功功率*/char c,ch,dstname[20],filename[20];int main(){opt:c=option();if(c=='1')goto opt1;else if(c=='2'){list("1.bin",1,filename);system("mode con:cols=100 lines=30 & color 07");while(c!='0'){system("cls");printf(" ################################计算结果文件查看################################\n");printf(" [8]上一个[2]下一个[5]查看[-]删除文件【0】退出\n\n");list("1.bin",0,filename);chose();goto opt;}else if(c=='0')exit(0);else goto opt;opt1:dataprepare();init:initial();Yc();printf("输入允许误差量：");scanf("%lf",&wucha);printf("选择模式：a.详情；b.简单：");redo:fflush(stdin);c=getchar();if(c!='a'&&c!='b'&&c!='A'&&c!='B'){printf("输入错误！重新输入：");goto redo;}if(c!='B'&&c!='b')showy();error=detpqc();showdetpq();while(error>wucha&&error<100){if(c!='B'&&c!='b'){printf("\t\t是否继续迭代<空格：是R/r:重赋初值E/e:退出潮流计算>？\n");ifdie: fflush(stdin);ch=getch();switch(ch){case ' ':break;case 'R':;case 'r':goto init;break;case 'E':;case 'e':goto opt;break;default: goto ifdie;}}diedai++;Jrc();if(c!='B'&&c!='b'){printf("\n>>>>>>>>>>>>>>>>>>>第%d 次迭代<<<<<<<<<<<<<<<<<<<<<\n",diedai);showj();Gauss();if(c!='B'&&c!='b')showdetav();for(temp1=0;temp1<pqnum+pvnum;temp1++) /* 修正电压*/pqv[temp1].ag=pqv[temp1].ag-detav[temp1];for(temp1=0;temp1<pqnum;temp1++)pqv[temp1].v=pqv[temp1].v-detav[temp1+pqnum+pvnum];error=detpqc();if(c!='B'&&c!='b')showdetpq();}if(error>=300||error<0){printf("\t\t迭代不收敛<R/r:重赋初值E/e:退出潮流计算>？\n");what: fflush(stdin);ch=getch();switch(ch){case 'R':;case 'r':goto init;break;case 'E':;case 'e':goto opt;break;default: goto what;}}showsolution();ch=save2file();if(ch=='1')goto opt;elsegoto init;}void dataprepare(void){FILE *fp;char ch,type[4];pqnum=0,pvnum=0,lznum=0,tznum=0,temp1=0,temp2=0,temp3=0,temp4=0;printf("\n 请输入潮流计算的数据文件\n --");getf:fflush(stdin);gets(fsource);while((fp=fopen(fsource,"r"))==NULL) /* 打开文件*/ {printf("\t\t不能打开文件！！\n\t\t重新输入请键入“Y”，其他任意键退出:");fflush(stdin);ch=getchar();if(ch=='Y'||ch=='y'){printf("\n\t\t潮流计算的数据文件：");goto getf;}elseexit(0);}system("mode con:cols=130 lines=50 & color 06");while(!feof(fp)) /* 计算各种节点数目*/{ch=fgetc(fp);if(ch==10){if(fgets(type,3,fp)){if(strcmp(type,"pq")==0)pqnum++;if(strcmp(type,"pv")==0)pvnum++;if(strcmp(type,"lz")==0)lznum++;if(strcmp(type,"tz")==0)tznum++;}}}pqv=(struct PQV*)malloc((pqnum+pvnum+1)*sizeof(struct PQV)); /* 根据节点数目开辟各节点储存空间*/lz=(struct LZ*)malloc((lznum)*sizeof(struct LZ));tz=(struct TZ*)malloc((tznum)*sizeof(struct TZ));rewind(fp);while(!feof(fp)) /* 读取和输入各节点数据*/{ch=fgetc(fp);if(ch==10){if(!fgets(type,3,fp))break;if(strcmp(type,"pq")==0){fscanf(fp," %s %lf %lf %lf",pqv[temp1].name,&pqv[temp1].vb,&pqv[temp1].p,&pqv[temp1]. q);temp1++;}if(strcmp(type,"pv")==0){fscanf(fp," %s %lf %lf\t%lf",pqv[temp2+pqnum].name,&pqv[temp2+pqnum].vb,&pqv[temp2 +pqnum].p,&pqv[temp2+pqnum].v);temp2++;}if(strcmp(type,"ph")==0){fscanf(fp," %s %lf\t\t%lf %lf",pqv[pqnum+pvnum].name,&pqv[pqnum+pvnum].vb,&pqv[pqn um+pvnum].v,&(pqv[pqnum+pvnum].ag));pqv[pqnum+pvnum].ag=pqv[pqnum+pvnum].ag/180*3.1415926;}if(strcmp(type,"lz")==0){fscanf(fp," %s %lf %s %lf %lf %lf %lf",lz[temp3].name1,&lz[temp3].vb1,lz[temp3].name2,&lz[ temp3].vb2,&lz[temp3].r,&lz[temp3].x,&lz[temp3].b);temp3++;}if(strcmp(type,"tz")==0){fscanf(fp," %s %lf %s %lf %lf %lf %lf",tz[temp4].name1,&tz[temp4].vb1,tz[temp4].name2,&tz [temp4].vb2,&tz[temp4].x,&tz[temp4].k1,&tz[temp4].k2);temp4++;}}}fclose(fp);puts("******************************************************************************* ********************************************");puts(" 电力系统潮流计算");puts("******************************************************************************* ********************************************");printf("从文件中获得的数据如下：\n");printf("节点类型名称额定电压有功标幺无功标幺电压标幺相角弧度\n");for(temp1=0;temp1<=pqnum+pvnum;temp1++){if(temp1<pqnum)printf("PQ %s %9.3lf %9.3lf %9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[temp1].p,pqv[temp1].q);else if(temp1<pqnum+pvnum)printf("PV %s %9.3lf %9.3lf %9.3lf\n",pqv[temp1].name,pqv[temp1].vb,pqv[te mp1].p,pqv[temp1].v);elseprintf(" 平衡%s %9.3lf %9.3lf %9.3lf\n",pqv[temp1].name,pqv[temp 1].vb,pqv[temp1].v,pqv[temp1].ag);}printf("\n");printf(" 线路节点1 额定电压节点2 额定电压电阻电抗对地导纳\n");for(temp1=0;temp1<lznum;temp1++)printf("输电线路%s %9.3lf %s %9.3lf %9.3lf %9.3lf %9.5lf\n",lz[temp1].name1,lz[temp1].vb1,lz[te mp1].name2,lz[temp1].vb2,lz[temp1].r,lz[temp1].x,lz[temp1].b);printf("\n");printf(" 线路节点1 额定电压节点2 额定电压电抗值一次比二次\n");for(temp1=0;temp1<tznum;temp1++)printf(" 变压器%s %9.3lf %s %9.3lf %9.3lf %9.3lf ：%9.3lf\n",tz[temp1].name1,tz[temp1].vb1, tz[temp1].name2,tz[temp1].vb2,tz[temp1].x,tz[temp1].k1,tz[temp1].k2);}void initial(void){char c;re_insert:for(temp1=0;temp1<pqnum+pvnum;temp1++){if(temp1<pqnum){printf("请输入第%d个节点%s(pq)的电压幅值：",temp1+1,pqv[temp1].name);fflush(stdin);scanf("%lf",&pqv[temp1].v);printf("请输入第%d个节点%s(pq)的电压相角：",temp1+1,pqv[temp1].name);fflush(stdin);scanf("%lf",&pqv[temp1].ag);pqv[temp1].ag=pqv[temp1].ag/180*3.1415926;}if(temp1>=pqnum){printf("请输入第%d个节点%s(pv)的电压相角：",temp1+1,pqv[temp1].name);fflush(stdin);scanf("%lf",&pqv[temp1].ag);}}printf("\t\t\t按ESC重新赋值，其他任意键继续！\n");c=getch();if(c==27)goto re_insert;}void Yc(void){int i, j,temp,flag=0;y=(struct Y **)malloc((pqnum+pvnum+1)*sizeof(struct Y *));for(temp1=0;temp1<pqnum+pvnum+1;temp1++)y[temp1]=(struct Y *)malloc((pqnum+pvnum+1)*sizeof(struct Y));for(i=0;i<=pqnum+pvnum;i++){for(j=0;j<=(pqnum+pvnum);j++){if(i!=j){flag=0;for(temp=0;temp<lznum;temp++) /* 当两节点间为线路时*/if((strcmp(pqv[i].name,lz[temp].name1)==0&&strcmp(pqv[j].name,lz[temp].name2)==0)||(s trcmp(pqv[i].name,lz[temp].name2)==0&&strcmp(pqv[j].name,lz[temp].name1)==0)){ flag++;y[i][j].G=-lz[temp].r/(pow(lz[temp].r,2)+pow(lz[temp].x,2));y[i][j].B=lz[temp].x/(pow(lz[temp].r,2)+pow(lz[temp].x,2));}for(temp=0;temp<tznum;temp++) /* 当两节点间为变压器*/if((strcmp(pqv[i].name,tz[temp].name1)==0&&strcmp(pqv[j].name,tz[temp].name2)==0)||( strcmp(pqv[i].name,tz[temp].name2)==0&&strcmp(pqv[j].name,tz[temp].name1)==0)){ flag++;y[i][j].G=0;y[i][j].B=1/tz[temp].x*tz[temp].k1/tz[temp].k2;}if(flag==0){y[i][j].G=0;y[i][j].B=0;}if(flag>1)printf("\n两节点间出现多条支路，本程序暂不支持!!\n");}}}for(i=0;i<=pqnum+pvnum;i++){y[i][i].G=0,y[i][i].B=0;for(temp=0;temp<lznum;temp++) /* 当与另一节点间为线路时*/if((strcmp(pqv[i].name,lz[temp].name1)==0)||(strcmp(pqv[i].name,lz[temp].name2)==0)) {y[i][i].G+=lz[temp].r/(pow(lz[temp].r,2)+pow(lz[temp].x,2));y[i][i].B+=-lz[temp].x/(pow(lz[temp].r,2)+pow(lz[temp].x,2))+lz[temp].b/2;}for(temp=0;temp<tznum;temp++) /* 当与另一节点间为变压器*/{if(strcmp(pqv[i].name,tz[temp].name1)==0){y[i][i].B+=-(1/tz[temp].x*tz[temp].k1/tz[temp].k2+1/tz[temp].x*(tz[temp].k2-tz[temp].k1)/tz[ temp].k2);}if(strcmp(pqv[i].name,tz[temp].name2)==0){y[i][i].B+=-(1/tz[temp].x*tz[temp].k1/tz[temp].k2+1/tz[temp].x*(tz[temp].k1-tz[temp].k2)*tz[ temp].k1/pow(tz[temp].k2,2));}}}}void showy(void){printf("\n");printf("导纳矩阵Y:\n");for(temp1=0;temp1<=pqnum+pvnum;temp1++){for(temp2=0;temp2<=pqnum+pvnum;temp2++)printf("%9.6lf+j%9.6lf ",y[temp1][temp2].G,y[temp1][temp2].B);printf("\n");}printf("\n");}double detpqc(void){double max=0;detpq=(double *)malloc((2*pqnum+pvnum)*sizeof(double)); /* 为detpq误差量*/ detav=(double *)malloc((2*pqnum+pvnum)*sizeof(double)); /* 修正量*/Pi=(double *)malloc((pqnum+pvnum)*sizeof(double)); /* */Qi=(double *)malloc((pqnum+pvnum)*sizeof(double));for(temp1=0;temp1<2*(pqnum+pvnum);temp1++){if(temp1<pqnum+pvnum) /* 求P（pq,pv)误差量pqnum+pvnum个*/{Pi[temp1]=0;for(temp2=0;temp2<=(pqnum+pvnum);temp2++)Pi[temp1]=Pi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*cos(pqv[temp1].ag-pqv[temp2].ag )+y[temp1][temp2].B*sin(pqv[temp1].ag-pqv[temp2].ag));Pi[temp1]=pqv[temp1].v*Pi[temp1];detpq[temp1]=pqv[temp1].p-Pi[temp1];if(fabs(detpq[temp1])>max)max=fabs(detpq[temp1]);}if(temp1>=pqnum+pvnum) /* 求Q（pq）误差量pqnum个*/{Qi[temp1-pqnum-pvnum]=0;for(temp2=0;temp2<=(pqnum+pvnum);temp2++)Qi[temp1-pqnum-pvnum]=Qi[temp1-pqnum-pvnum]+pqv[temp2].v*(y[temp1-pqnum-pvnu m][temp2].G*sin(pqv[temp1-pqnum-pvnum].ag-pqv[temp2].ag)-y[temp1-pqnum-pvnum][temp2 ].B*cos(pqv[temp1-pqnum-pvnum].ag-pqv[temp2].ag));Qi[temp1-pqnum-pvnum]=pqv[temp1-(pqnum+pvnum)].v*Qi[temp1-pqnum-pvnum];detpq[temp1]=pqv[temp1-(pqnum+pvnum)].q-Qi[temp1-pqnum-pvnum];if(fabs(detpq[temp1])>max&&temp1!=2*pqnum+pvnum)max=fabs(detpq[temp1]);}}return max;}void showdetav(void){printf("解得:\n");for(temp1=0;temp1<2*pqnum+pvnum;temp1++)if(temp1<pqnum)printf("Δδ%d：%lf ",temp1,detav[temp1]);elseprintf("ΔV %d: %lf ",temp1-pqnum,detav[temp1]);printf("\n");}void Jrc(void){int i,j;double agij;J=(double **)malloc((2*pqnum+pvnum)*sizeof(double *)); /* 为雅可比矩阵分配内存空间*/for(temp1=0;temp1<2*pqnum+pvnum;temp1++)J[temp1]=(double *)malloc((2*pqnum+pvnum)*sizeof(double));for(i=0;i<pqnum+pvnum;i++){for(j=0;j<pqnum+pvnum;j++) /* Hij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i][j]=-pqv[i].v*pqv[j].v*(y[i][j].G*sin(agij)-y[i][j].B*cos(agij));}elseJ[i][j]=pqv[i].v*pqv[i].v*y[i][i].B+Qi[i];}for(j=0;j<pqnum;j++) /* Nij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i][j+pqnum+pvnum]=-pqv[i].v*(y[i][j].G*cos(agij)+y[i][j].B*sin(agij));}elseJ[i][j+pqnum+pvnum]=-pqv[i].v*y[i][i].G-Pi[i];}}for(i=0;i<pqnum;i++){for(j=0;j<pqnum+pvnum;j++) /* Kij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i+pqnum+pvnum][j]=pqv[i].v*pqv[j].v*(y[i][j].G*cos(agij)+y[i][j].B*sin(agij));}elseJ[i+pqnum+pvnum][j]=pqv[i].v*pqv[i].v*y[i][i].G-Pi[i];}for(j=0;j<pqnum;j++) /* Lij */{if(i!=j){agij=pqv[i].ag-pqv[j].ag;J[i+pqnum+pvnum][j+pqnum+pvnum]=-pqv[i].v*(y[i][j].G*sin(agij)-y[i][j].B*cos(agij));}elseJ[i+pqnum+pvnum][j+pqnum+pvnum]=pqv[i].v*y[i][i].B-Qi[i];}}}void showj(void){printf("\n");printf("雅可比矩阵J:\n");for(temp1=0;temp1<2*pqnum+pvnum;temp1++){for(temp2=0;temp2<2*pqnum+pvnum;temp2++)printf("%9.6lf ",J[temp1][temp2]);printf("\n");}printf("\n");}void Gauss(void){int xiao,big,i,j,n=2*pqnum+pvnum;double biggest,sta,E;double conj=0,conpq=0;for(xiao=0;xiao<2*pqnum+pvnum-1;xiao++){big=xiao;biggest=fabs(J[xiao][xiao]);for(i=xiao+1;i<2*pvnum+pqnum;i++) /* 找出列主元*/ if(fabs(J[i][xiao])>biggest){big=i;biggest=fabs(J[i][xiao]);}for(i=0;i<2*pqnum+pvnum;i++) /* 交换两行*/{conj=J[xiao][i];J[xiao][i]=J[big][i];J[big][i]=conj;}conpq=detpq[xiao];detpq[xiao]=detpq[big];detpq[big]=conpq;for(i=xiao+1;i<n;i++) /* 消元*/{sta=J[i][xiao];for(j=0;j<n;j++)J[i][j]=J[i][j]-sta*J[xiao][j]/J[xiao][xiao];detpq[i]=detpq[i]-sta*detpq[xiao]/J[xiao][xiao];}}detav[n-1]=detpq[n-1]/J[n-1][n-1]; /* 回代*/for(i=n-2;i>=0;i--){ E=0;for(j=i+1;j<2*pqnum+pvnum;j++)E=E+J[i][j]*detav[j];detav[i]=(detpq[i]-E)/J[i][i];}}void showdetpq(void){printf("\n");printf("功率不平衡量:\n");for(temp1=0;temp1<2*pqnum+pvnum;temp1++){if(temp1<pqnum+pvnum)printf("ΔP%d: %lf ",temp1+1,detpq[temp1]);elsepr intf("ΔQ%d: %lf ",temp1-pqnum,detpq[temp1]);}printf("\n最大功率不平衡量为（绝对值）：%lf\n",error);printf("\n");}void showsolution(void){printf("\n>>>>>>>>>>>>>>>>>>> 迭代结束<<<<<<<<<<<<<<<<<<<<<\n");printf("\t\t总共进行%d次迭代:\n",diedai);for(temp1=0;temp1<pqnum+pvnum;temp1++)printf("节点%d<%s>电压：幅值%lf 相角%lf\n",temp1+1,pqv[temp1].name,pqv[temp1].v,pqv[temp1].ag);for(temp1=0;temp1<(pqnum+pvnum+1);temp1++){Pi[temp1]=0;for(temp2=0;temp2<=(pqnum+pvnum);temp2++){Pi[temp1]=Pi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*cos(pqv[temp1].ag-pqv[temp2].ag )+y[temp1][temp2].B*sin(pqv[temp1].ag-pqv[temp2].ag));Qi[temp1]=Qi[temp1]+pqv[temp2].v*(y[temp1][temp2].G*sin(pqv[temp1].ag-pqv[temp2].a g)-y[temp1][temp2].B*cos(pqv[temp1].ag-pqv[temp2].ag));}Pi[temp1]=Pi[temp1]*pqv[temp1].v;Qi[temp1]=Qi[temp1]*pqv[temp1].v;printf("节点%d<%s>功率: %8.6lf+j%8.6lf\n",temp1+1,pqv[temp1].name,Pi[temp1],Qi[temp1]);}printf("\n");}char save2file(void){FILE *fp;int i;time_t t;char filename[20],c;printf("\t\t是否保存文件<Y/y>：");fflush(stdin);c=getchar();if(c!='Y'&&c!='y')goto next;printf("输入要储存数据的文件名(默认后缀.data)：");fflush(stdin);scanf("%s",filename);strcat(filename,".data");if((fp=fopen(filename,"w"))==NULL){printf("无法执行存储操作！！");exit(0);}time(&t);fprintf(fp,"\t\t\t潮流计算结果\n保存时间：%s\n对应数据文件：%s\n",ctime(&t),fsource);fprintf(fp,"<标幺值>\n");for(i=0;i<pqnum;i++)fprintf(fp," PQ节点%s，电压幅值：%10.6lf 相角:%10.6lf度有功功率：%10.6lf 无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3.14159265337,Pi[i],Qi[i]);for(i=pqnum;i<pqnum+pvnum;i++)fprintf(fp," PV节点%s，电压幅值：%10.6lf 相角:%10.6lf度有功功率：%10.6lf 无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3.14159265337,Pi[i],Qi[i]);fprintf(fp,"平衡节点%s，电压幅值：%10.6lf 相角:%10.6lf度有功功率：%10.6lf 无功功率：%10.6lf\n",pqv[i].name,pqv[i].v,pqv[i].ag*180/3.14159265337,Pi[i],Qi[i]);fclose(fp);printf("\t\t\t文件保存成功\n");next:puts("\t1.返回主界面 2.重新计算此文件数据:");nex:fflush(stdin);c=getch();if(c!='1'&&c!='2')goto nex;return c;}char option(void){char c;system("mode con:cols=75 lines=25 &color 16");cls:puts(" ***************************************************************");puts(" 电力系统潮流计算");puts("***************************************************************\n\n");printf(" 选择您的操作：\n\n\t\t1.潮流计算 2.查看数据文件0.退出程序\n\t\t\t\t");fflush(stdin);c=getch();return c;}int list(char filename[],int i,char dstn[]){FILE *fp;char ch,c;int shu=0,j=0,k=0;char name[50];system("dir /x *.data>1.bin 2>nul ");if((fp=fopen("1.bin","r"))==NULL){printf("读取列表有误！");system("pause");return 0;}while((ch=fgetc(fp))!=EOF){if(ch=='2'&&c==10){shu++;if(shu!=i){fseek(fp,47L,1);fgets(name,50L,fp);}if(shu==i){fseek(fp,35L,1);fgets(dstn,13,fp);if(dstn[0]==' '){fseek(fp,1L,1);fgets(dstn,13,fp);dstn[strlen(dstn)-1]=0;}return 0;}if(i==0){printf("\t\t");if(numw==shu)printf("\b\b->");printf(" %3d . %s",shu,name);numofw=shu;}ch=10;}c=ch;}fclose(fp);return 0;}int chose(void){char del[30]={"del "};FILE *fp;fflush(stdin);c=getch();if(c=='2'&&numw<numofw)numw++;if(c=='8'&&numw>1)numw--;if(c=='-'){list("1.bin",numw,filename);strcat(del,filename);strcat(del," 2>null");system(del);}if(c=='5') /* 5 is play key */{list("1.bin",numw,filename);if((fp=fopen(filename,"r"))==NULL){printf("读取列表有误！");system("pause");return 0;}while((ch=fgetc(fp))!=EOF)putchar(ch);system("pause & del /f /Q 2>nul");fclose(fp);}return 0;}运行效果：启动主界面：1.键入1，开始潮流计算，输入潮流计算数据文件：程序将自动读取数据，并进入如图界面，列出数据：此时需为潮流计算设定初值，如下图：若需重新给定初值则键入按ESC键，否则按任意键继续：按任意键给定允许误差，并选择显示模式：a.详情可列出计算过程的各种数据，如导纳矩阵，雅可比矩阵，误差量等。

C语言进行潮流计算

电力系统课程设计C语言潮流计算学院：电气工程班级：电092 班学号：0912002020 学生姓名：闵凯2013.3.7电力系统的潮流计算是对电力系统分析的最基本步骤也是最重要的步骤，是指在一定的系统结构和运行条件下，确定系统运行状态的计算，也即是对各母线（节点）电压，各元件（支路）传输电线或功率的计算。

通过计算出的节点电压和功率分布用以检查系统各元件是否过负荷，各点电压是否合理，以及功率损耗等。

即使对于一个简单的电力系统，潮流计算也不是一件简单就可以完成的事，其运算量很大，因此如果对于一个大的、复杂的电网来说的话，由于其节点多，分支杂，其计算量可想而知，人工对其计算也更是难上加难了。

特别是在现实生活中，遇到一个电力系统不会像我们期望的那样可以知道它的首端电压和首端功率或者是末端电压和末端功率，而是只知道它的首端电压和末端功率，更是使计算变的头疼万分。

为了使计算变的简单，我们就可以利用计算机，用C语言编程来实现牛顿-拉夫逊（Newton-Raphson）迭代法，最终实现对电力系统潮流的计算。

用牛顿-拉夫逊迭代法进行电力系统潮流计算的相关概念1.节点导纳矩阵如图所示的电力网络，将节点i和j的电压用U i和U j表示,它们之间的支路导纳表示为y j,那么有基尔霍夫电流定律可知注入接点I的电流I i（设流入节点的电流为正）等于离开节点I的电流之和，因此有y u上式也可以写为n-zj =0y j u=Y HI i:I = YU-y j =Yj则可将（1-2）改写为：Y j U j I=1,2,…,n.(1-1)(1-2)(1-3)(1-4)' y iJ(U i-u J)I i 八I ij其中Y 为节点导纳矩阵，也称为稀疏的对称矩阵，它是 nM 阶方阵。

对角元 Yu 称为自导纳，它等于与该节点I 直接相连的所有支路导纳总和；非对角元 Y j ( i 并)称为互导纳或转移导纳，它等于连结节点 I , j 支路导纳的负数，且有 Y j =Y ji ，当节点I , j 之间没有支路直接相连时，Y j =Y ji =O 。

潮流计算C语言编程

// flow.cpp : Defines the entry point for the console application.//#include "stdafx.h"#include <iostream.h>#include <fstream.h>#include "NEquation.h"#include "math.h"#include "config.h"void test(){NEquation ob1;ob1.SetSize(2);ob1.Data(0,0)=1;ob1.Data(0,1)=2;ob1.Data(1,0)=2;ob1.Data(1,1)=1;ob1.V alue(0)=5;ob1.V alue(1)=4;ob1.Run();printf("x1=%f\n",ob1.V alue(0));printf("x2=%f\n",ob1.V alue(1));}void GetData()//Read the data{FILE *fp;int i;if((fp=fopen("F:\\1061180523\\flow\\data\\data.txt","r"))==NULL){printf("Can not open the file named 'data.txt' \n");return;}for(i=0;i<=Bus_Num-1;i++){fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&gBus[i].No,&gBus[i].V oltage,&gBus[i].Phase, &gBus[i].GenP,&gBus[i].GenQ,&gBus[i].LoadP,&gBus[i].LoadQ,&gBus[i].Type);}for(i=0;i<=Line_Num-1;i++){fscanf(fp,"%d,%d,%d,%f,%f,%f,%f",&gLine[i].No,&gLine[i].No_I,&gLine[i].No_J, &gLine[i].R,&gLine[i].X,&gLine[i].B,&gLine[i].k);}fclose(fp);}void GetYMatrix(){int i,j;FILE *fp;// calculate the Y matrixint bus1,bus2;float r,x,d,g,b;for(i=0;i<=Bus_Num-1;i++) //导纳矩阵赋初值为零{for(j=0;j<=Bus_Num-1;j++){gY_G[i][j]=0;gY_B[i][j]=0;}}for(i=0; i<=Line_Num-1; i++){bus1=gLine[i].No_I-1; //线路一侧连接的节点号减1bus2=gLine[i].No_J-1; //线路另一侧连接的节点号减1r=gLine[i].R;x=gLine[i].X;d=r*r+x*x;g=r/d;b=-x/d;if(gLine[i].k==0) //没有变压器的情况{gY_G[bus1][bus1]=gY_G[bus1][bus1]+g;//gY_G[bus1][bus1]=ggY_G[bus2][bus2]=gY_G[bus2][bus2]+g;//gY_G[bus2][bus2]=ggY_G[bus1][bus2]=gY_G[bus1][bus2]-g;//gY_G[bus1][bus2]=-g 矩阵是对称的gY_G[bus2][bus1]=gY_G[bus2][bus1]-g;//gY_G[bus2][bus1]=-ggY_B[bus1][bus1]=gY_B[bus1][bus1]+b+gLine[i].B;//gY_B[bus1][bus1]=b+gLine[i].BgY_B[bus2][bus2]=gY_B[bus2][bus2]+b+gLine[i].B;//gY_B[bus2][bus2]=b+gLine[i].BgY_B[bus1][bus2]=gY_B[bus1][bus2]-b;//gY_B[bus1][bus2]=-bgY_B[bus2][bus1]=gY_B[bus2][bus1]-b;//gY_B[bus2][bus1]=-b}else //有变压器的情况{gY_G[bus1][bus1]=gY_G[bus1][bus1]+g/gLine[i].k+g*(1-gLine[i].k)/gLine[i].k/gLine[i].k;gY_G[bus2][bus2]=gY_G[bus2][bus2]+g/gLine[i].k+g*(gLine[i].k-1)/gLine[i].k;gY_G[bus1][bus2]=gY_G[bus1][bus2]-g/gLine[i].k;gY_G[bus2][bus1]=gY_G[bus2][bus1]-g/gLine[i].k;gY_B[bus1][bus1]=gY_B[bus1][bus1]+b/gLine[i].k+b*(1-gLine[i].k)/gLine[i].k/gLine[i].k;gY_B[bus2][bus2]=gY_B[bus2][bus2]+b/gLine[i].k+b*(gLine[i].k-1)/gLine[i].k;gY_B[bus1][bus2]=gY_B[bus1][bus2]-b/gLine[i].k;gY_B[bus2][bus1]=gY_B[bus2][bus1]-b/gLine[i].k;}}// output the Y matrixif((fp=fopen("F:\\1061180523\\flow\\data\\ymatrix.txt","w"))==NULL){printf("Can not open the file named 'ymatrix.txt' \n");return ;}fprintf(fp,"---Y Matrix---\n");for(i=0;i<=Bus_Num-1;i++){for(j=0;j<=Bus_Num-1;j++){fprintf(fp,"Y(%d,%d)=(%10.5f,%10.5f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);}}fclose(fp);}void SetInitial(){int i;for(i=0;i<=Bus_Num-1;i++){if(gBus[i].Type==3){gf[i]=gBus[i].V oltage*sin(gBus[i].Phase);ge[i]=gBus[i].V oltage*cos(gBus[i].Phase);}else{gf[i]=0;ge[i]=1;}}}void GetUnbalance(){int i,j;FILE *fp;for(i=0;i<=Bus_Num-2;i++){gDelta_P[i]=gBus[i+1].GenP-gBus[i+1].LoadP;if(gBus[i+1].Type==2) //PV节点gDelta_Q[i]=gBus[i+1].V oltage*gBus[i+1].V oltage-(ge[i+1]*ge[i+1]+gf[i+1]*gf[i+1]);elsegDelta_Q[i]=gBus[i+1].GenQ-gBus[i+1].LoadQ;for(j=0;j<=Bus_Num-1;j++){gDelta_P[i]=gDelta_P[i]-ge[i+1]*(gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j])-gf[i+1]*(gY_G[i +1][j]*gf[j]+gY_B[i+1][j]*ge[j]);if(gBus[i+1].Type==1) //PQ节点gDelta_Q[i]=gDelta_Q[i]-gf[i+1]*(gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j])+ge[i+1]*(gY_G[i+1][ j]*gf[j]+gY_B[i+1][j]*ge[j]);}}for(i=0;i<=Bus_Num-2;i++) //合并{gDelta_PQ[2*i]=gDelta_P[i];gDelta_PQ[2*i+1]=gDelta_Q[i];}if((fp=fopen("F:\\1061180523\\flow\\data\\unbalance.txt","w"))==NULL){printf("Can not open the file named 'unbalance.txt' \n");return ;}fprintf(fp,"---Unbalance---\n");for(i=0;i<=2*Bus_Num-3;i++){fprintf(fp,"Unbalance[%d]=%10.5f\n",i+1,gDelta_PQ[i]);}fclose(fp);}void GetJaccobi(){int i,j;float ga[Bus_Num-1],gb[Bus_Num-1];FILE *fp;for(i=0;i<=Bus_Num-2;i++) //计算注入电流{ga[i]=0;gb[i]=0;for(j=0;j<=Bus_Num-1;j++){ga[i]=ga[i]+gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j];gb[i]=gb[i]+gY_G[i+1][j]*gf[j]+gY_B[i+1][j]*ge[j];}}for(i=0;i<=Bus_Num-2;i++){for(j=0;j<=Bus_Num-2;j++){if(i!=j){gJaccobi[2*i][2*j]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1];gJaccobi[2*i][2*j+1]=gY_G[i+1][j+1]*ge[i+1]+gY_B[i+1][j+1]*gf[i+1];if(gBus[i+1].Type==2) //PV节点{gJaccobi[2*i+1][2*j]=0;gJaccobi[2*i+1][2*j+1]=0;}else //PQ{gJaccobi[2*i+1][2*j]=-gJaccobi[2*i][2*j+1];gJaccobi[2*i+1][2*j+1]=gJaccobi[2*i][2*j];}}else{gJaccobi[2*i][2*j]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1]+gb[i];gJaccobi[2*i][2*j+1]=gY_G[i+1][j+1]*ge[i+1]+gY_B[i+1][j+1]*gf[i+1]+ga[i];if(gBus[i+1].Type==2) //PV节点{gJaccobi[2*i+1][2*j]=2*gf[i+1];gJaccobi[2*i+1][2*j+1]=2*ge[i+1];}else //PQ{gJaccobi[2*i+1][2*j]=-gY_G[i+1][j+1]*ge[i+1]-gY_B[i+1][j+1]*gf[i+1]+ga[i];gJaccobi[2*i+1][2*j+1]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1]-gb[i];}}}}if((fp=fopen("F:\\1061180523\\flow\\data\\jaccobi.txt","w"))==NULL){printf("Can not open the file named 'jaccobi.txt' \n");return ;}fprintf(fp,"---Jaccobi Matrix---\n");for(i=0;i<=2*Bus_Num-3;i++){for(j=0;j<=2*Bus_Num-3;j++){fprintf(fp,"jaccobi(%d,%d)=%10.5f\n",i+1,j+1,gJaccobi[i][j]);}}fclose(fp);}void GetRevised(){int i,j;FILE *fp;NEquation ob1; //解矩阵方程ob1.SetSize(2*(Bus_Num-1));for(i=0;i<=2*Bus_Num-3;i++)for(j=0;j<=2*Bus_Num-3;j++)ob1.Data(i,j)=gJaccobi[i][j];for(i=0;i<=2*Bus_Num-3;i++)ob1.V alue(i)=gDelta_PQ[i];ob1.Run();for(i=0;i<=Bus_Num-2;i++){gDelta_f[i]=ob1.V alue(2*i);gDelta_e[i]=ob1.V alue(2*i+1);gDelta_fe[2*i]=gDelta_f[i];gDelta_fe[2*i+1]=gDelta_e[i];}if((fp=fopen("F:\\1061180523\\flow\\data\\revised.txt","w"))==NULL) {printf("Can not open the file named 'revised.txt' \n");return ;}fprintf(fp,"---Revised---\n");for(i=0;i<=2*Bus_Num-3;i++){fprintf(fp,"revised[%d]=%10.5f\n",i+1,gDelta_fe[i]);}fclose(fp);}void GetNew V alue(){int i;FILE *fp;for(i=0;i<=Bus_Num-2;i++){gf[i+1]=gf[i+1]+gDelta_f[i];ge[i+1]=ge[i+1]+gDelta_e[i];}if((fp=fopen("F:\\1061180523\\flow\\data\\newvalue.txt","w"))==NULL) {printf("Can not open the file named 'newvalue.txt' \n");return ;}fprintf(fp,"---New V alue---\n");for(i=0;i<=Bus_Num-2;i++){fprintf(fp,"f(%d)=%10.5f,e(%d)=%10.5f\n",i+1,gf[i+1],i+1,ge[i+1]);}fclose(fp);}int main(int argc, char* argv[]){int i,Count_Num;float maxV alue;test();GetData();GetYMatrix();SetInitial();for(Count_Num=0;Count_Num<=100;Count_Num++){GetUnbalance();GetJaccobi();GetRevised();GetNewV alue();maxV alue=fabs(gDelta_fe[0]);for(i=1;i<=2*(Bus_Num-1)-1;i++){if(maxV alue<fabs(gDelta_fe[i])){maxV alue=fabs(gDelta_fe[i]);}}if(maxV alue<Precision){break;}}printf("%d\n",Count_Num);for(i=0;i<=Bus_Num-1;i++){printf("%10.5f\n",sqrt(ge[i]*ge[i]+gf[i]*gf[i]));}return 0;}。

电力系统潮流计算C程序精编版

电力系统潮流计算注：这是一个基于N－R法的潮流计算通用程序，仅提供了子程序，需要做些处理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享！！！/******************************************************************* * 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是 ** Newton_Raphson法.** 程序中所用的变量说明如下:** N:网络节点总数. M:网络的PQ节点数. ** L:网络的支路总数. N0:雅可比矩阵的行数. ** N1:N0+1K:打印开关.K=1,则打印;否则,不打印.** K1:子程序PLSC中判断输入电压的形式.K1=1,则为极座标形式.否则** 为直角坐标形式.** D:有功及无功功率误差的最大值. * * G(I,J):Ybus的电导元素(实部).** B(I,J):Ybus的电纳元素(虚部).** G1(I) :第I支路的串联电导. B1(I):第I支路的串联电纳. ** C1(I) :第I支路的pie型对称接地电纳. ** C(I,J):第I节点J支路不对称接地电纳. ** CO(I) :第I节点的接地电纳.** S1(I) :第I节点的起始节点号. E1(I):第I节点的终止节点号. ** P(I) :第I节点的注入有功功率. Q(I):第I节点的注入无功功率.** P0(I) :第I节点有功功率误差. Q0(I):第I节点无功功率误差. ** V0(I) :第I节点(PV节点)的电压误差(平方误差). ** V(I) :第I节点的电压误差幅值. * * E(I) :第I节点的电压的实部. F(I):第I节点的电压的虚部. ** JM(I,J):Jacoby矩阵的第I行J列元素. ** A(I,J):修正方程的增广矩阵,三角化矩阵的第I行J列元素,运算结 ** 束后A矩阵的最后一列存放修正的解. ** P1(I) :第I支路由S1(I)节点注入的有功功率. ** Q1(I) :第I支路由S1(I)节点注入的无功功率. ** P2(I) :第I支路由E1(I)节点注入的有功功率. ** Q2(I) :第I支路由E1(I)节点注入的无功功率. ** P3(I) :第I支路的有功功率损耗. * * Q3(I) :第I支路的无功功率损耗. * * ANGLE(I):第I节点电压的角度.********************************************************************/ #include <math.h>#include <stdio.h>#define f1(i) (i-1)/* 把习惯的一阶矩阵的下标转化为C语言数组下标*/#define f2(i,j,n) ((i-1)*(n)+j-1)/* 把习惯的二阶矩阵的下标转化为C语言数组下标*//***************************************************** 本子程序根据所给的支路导纳及有关信息,形成结点 ** 导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B *****************************************************/void ybus(int n,int l,int m,float *g,float *b,float *g1,float *b1,float *c1,\float *c,float *co,int k,int *s1,int *e1) {extern FILE *file4;FILE *fp;int i,j,io,i0;int pos1,pos2;int st,en;if(file4==NULL){fp=stdout;}else{fp=file4; /* 输出到文件 */}/* 初始化矩阵G,B */for(i=1;i<=n;i++){for(j=1;j<=n;j++){pos2=f2(i,j,n);g[pos2]=0;b[pos2]=0;}}/* 计算支路导纳 */for(i=1;i<=l;i++){/* 计算对角元 */pos1=f1(i);st=s1[pos1];en=e1[pos1];pos2=f2(st,st,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];pos2=f2(en,en,n);g[pos2]+=g1[pos1];b[pos2]+=b1[pos1]+c1[pos1];/* 计算非对角元 */pos2=f2(st,en,n);g[pos2]-=g1[pos1];b[pos2]-=b1[pos1];g[f2(en,st,n)]=g[f2(st,en,n)];b[f2(en,st,n)]=b[f2(st,en,n)];}/* 计算接地支路导纳 */for(i=1;i<=n;i++){/* 对称部分*/b[f2(i,i,n)]+=co[f1(i)];/* 非对称部分*/for(j=1;j<=l;j++){b[f2(i,i,n)]+=c[f2(i,j,l)];}}if(k!=1){return; /* 如果K不为 1,则返回;否则,打印导纳矩阵 */ }fprintf(fp,"\n BUS ADMITTANCE MATRIXY(BUS):");fprintf(fp,"\n ******************* ARRAY G ********************"); for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",g[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n ******************* ARRAY B ********************");for(io=1;io<=n;io+=5){i0=(io+4)>n?n:(io+4);fprintf(fp,"\n");for(j=io;j<=i0;j++){fprintf(fp,"%13d",j);}for(i=1;i<=n;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%13.6f",b[f2(i,j,n)]);}}fprintf(fp,"\n");}fprintf(fp,"\n************************************************");}/******************************************** 本子程序根据所给的功率及电压等数据 ** 求出功率及电压误差量,并返回最大有功功率 ** 以用于与给定误差比较.如打印参数K=1,则输 ** 出P0,Q0(对PQ结点),V0(对PV结点). ********************************************/void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,int m,\int n,float *e,float *f,int k,float *g,float *b,float *dd){extern FILE *file4;FILE *fp;int i,j,l;int pos1,pos2;float a1,a2,d1,d;if(file4==NULL){fp=stdout; /* 输出到屏幕 */}else{fp=file4; /* 输出到文件*/}l=n-1;if(k==1){fprintf(fp,"\n CHANGE OFP0,V**2,P0(I),Q0(I),V0(I) ");fprintf(fp,"\n I P0(I)Q0(I)");}for(i=1;i<=l;i++){a1=0;a2=0;pos1=f1(i);for(j=1;j<=n;j++){/* a1,a2对应课本p185式(4-86)中括号内的式子 */pos2=f2(i,j,n);a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)];a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)];}/* 计算式(4-86)(4-87)中的deltaPi */p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2;if(i <= m){ /* 计算PQ结点中的deltaQi */q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2;}else{ /* 计算PV结点中的deltaVi平方*/v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1];}/* 输出结果 */if(k==1){if(i<m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);}else if(i==m){fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]);fprintf(fp,"\n I P0(I)V0(I)");}else{fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],v0[pos1]);}}}/* 找到deltaP和deltaQ中的最小者，作为收敛指标, 存在dd中 */d=0;for(i=1;i<=l;i++){pos1=f1(i);d1=p0[pos1]>0 ? p0[pos1] : -p0[pos1];if(d<d1){d=d1;}if(i<=m){d1=q0[pos1]>0?q0[pos1]:-q0[pos1];if(d<d1){d=d1;}}}(*dd)=d;}/**************************************************** 本子程序根据节点导纳及电压求Jacoby矩阵,用于求** 电压修正量,如打印参数K=1,则输出Jacoby矩阵. ** 对应于课本P186式(4-89)(4-90) ****************************************************/void jmcc(int m,int n,int n0,float *e,float *f,float *g,float *b,float *jm,int k){extern FILE *file4;FILE *fp;int i,j,i1,io,i0,ns;int pos1,pos2;if(file4==NULL){fp=stdout;}else{fp=file4;}/* 初始化矩阵jm */for(i=1;i<=n0;i++){for(j=1;j<=n0;j++){jm[f2(i,j,n0)]=0;}}ns=n-1; /* 去掉一个平衡结点 *//* 计算式(4-89)(4-90) */for(i=1;i<=ns;i++){/* 计算式(4-90) */for(i1=1;i1<=n;i1++){/* pos1是式(4-90)中的j */pos1=f1(i1);/* pos2是式(4-90)中的ij */pos2=f2(i,i1,n);if(i<=m) /* i是PQ结点 */{/* 计算式(4-90)中的Jii等式右侧第一部分 */jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1 ];/* 计算式(4-90)中的Lii等式右侧第一部分 */jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];}/* 计算式(4-90)中的Hii等式右侧第一部分 */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1];/* 计算式(4-90)中的Nii等式右侧第一部分 */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1];}/* pos2是式(4-90)中的ii */pos2=f2(i,i,n);/* pos1是式(4-90)中的i */pos1=f1(i);if(i<=m) /* i是PQ结点 */{/* 计算式(4-90)中的Jii */jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1] ;/* 计算式(4-90)中的Lii */jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[po s1];}/* 计算式(4-90)中的Hii */jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1];/* 计算式(4-90)中的Jii */jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1];if(i>m) /* PV结点 */{/* 计算式(4-90)中的Rii */jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1];/* 计算式(4-90)中的Sii */jm[f2(2*i-1,2*i,n0)]=-2*f[pos1];}/* 计算式(4-89) */for(j=1;j<=ns;j++){if(j!=i){/* pos1是式(4-89)中的i */pos1=f1(i);/* pos2是式(4-89)中的ij */pos2=f2(i,j,n);/* 计算式(4-89)中的Nij */jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1];/* 计算式(4-89)中的Hij */jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1];if(i<=m) /* i是PQ结点 */{/* 计算式(4-89)中的Lij (=-Hij) */jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)];/* 计算式(4-89)中的Jij (=Nij) */jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)];}else /* i是PV结点 */{/* 计算式(4-89)中的Rij (=0) */jm[f2(2*i-1,2*j-1,n0)]=0;/* 计算式(4-89)中的Sij (=0) */jm[f2(2*i-1,2*j,n0)]=0;}}}}if(k!=1){return;}/* 输出Jacoby矩阵 */fprintf(fp,"\n J MATRIX(C)"); for(io=1;io<=n0;io+=5){i1=(io+4)>n0?n0:(io+4);fprintf(fp,"\n");for(j=io;j<=i1;j++){fprintf(fp,"%10d",j);}for(i=1;i<=n0;i++){fprintf(fp,"\n%2d",i);for(j=io;j<=i1;j++){fprintf(fp,"%12.6f",jm[f2(i,j,n0)]);}}}fprintf(fp,"\n");}/*********************************************** 本子程序用选列主元素的高斯消元法求解组 ** 性方程组求各结点电压修正量,如打印参数K=1,则** 输出增广矩阵变换中的上三角及电压修正量.如果** 无唯一解,则给出信息,并停止程序运行. ***********************************************/void sevc ( float a[], int n0, int k, int n1){extern FILE *file4;FILE *fp;int i,j,l,n2,n3,n4,i0,io,j1,i1;float t0,t,c;if(file4==NULL) fp=stdout;else fp=file4;for(i=1;i<=n0;i++){l=i;for(j=i;j<=n0;j++){if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) ){l=j; /* 找到这行中的最大元 */}}if(l!=i){ /* 行列交换 */for (j=i;j<=n1;j++){t=a[f2(i,j,n1)];a[f2(i,j,n1)]=a[f2(l,j,n1)];a[f2(l,j,n1)]=t;}}if (fabs(a[f2(i,i,n1)]-0)<1e-10){ /* 对角元近似于0, 无解 */printf("\nNo Solution\n");exit (1);}t0=a[f2(i,i,n1)];for(j=i;j<=n1;j++){/* 除对角元 */a[f2(i,j,n1)]/=t0;}if(i==n0){ /* 最后一行，不用消元 */continue;}/* 消元 */j1=i+1;for(i1=j1;i1<=n0;i1++){c=a[f2(i1,i,n1)];for(j=i;j<=n1;j++){a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c;}}}if(k==1){ /* 输出上三角矩阵 */fprintf(fp,"\nTrianglar Angmentex Matrix ");for(io=1;io<=n1;io+=5){i0=(io+4)>n1?n1:(io+4);fprintf(fp,"\n");fprintf(fp," ");for(i=io;i<=i0;i++){fprintf(fp,"%12d",i);}for(i=1;i<=n0;i++){fprintf(fp,"\n");fprintf(fp,"%2d",i);for(j=io;j<=i0;j++){fprintf(fp,"%15.6f", a[f2(i,j,n1)]);}}}}/* 回代求方程解 */n2=n1-2;for(i=1;i<=n2;i++){n3=n1-i;for(i1=n3;i1<=n0;i1++){n4=n0-i;a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)];}}if(k!=1){return;}/* 输出电压修正值 */fprintf(fp,"\nVoltage correction E(i), F(i) :");for(io=1;io<=n0;io+=4){i1=(io+1)/2;i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2);fprintf(fp,"\n");for(j=i1;j<=i0;j++){fprintf(fp,"%16d%16d",j,j);}i1 = 2*i0;fprintf(fp,"\n");for(i=io;i<=i1;i++){fprintf(fp,"%15.6f", a[f2(i,n1,n1)]);}}}#define Pi 3.1415927/180void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\int e1[],int s1[],float g1[],floatb1[],float c1[],float c[],\float co[],float p1[],float q1[],floatp2[],float q2[],float p3[],\float q3[],float p[],float q[],floatv[],float angle[],int k1){extern FILE *file4;FILE *fp;float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2;int i,i1,j,m1,ns,pos1,pos2,km;ns=n-1;if(file4==NULL){fp=stdout;}else{fp=file4;}fprintf(fp,"\nTHE RESULT ARE:");if(k1==1){for(i=0;i<n;i++){angle[i]*=Pi;e[i]=v[i]*cos(angle[i]);f[i]=v[i]*sin(angle[i]);}}t1=0.0;t2=0.0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(n,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(n);p[pos1]=t1*e[pos1];q[pos1]=-t2*e[pos1];m1=m+1;for(i1=m1;i1<=ns;i1++){t1=0;t2=0;for(i=1;i<=n;i++){pos1=f1(i);pos2=f2(i1,i,n);t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1];t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1];}pos1=f1(i1);q[pos1]=f[pos1]*t1-e[pos1]*t2;}for(i=0;i<n; i++){cm=co[i];if(cm!=0){q[i]-=(e[i]*e[i]+f[i]*f[i])*cm;}}fprintf(fp,"\nBUS DATA");fprintf(fp,"\nBUS VOLTAGE ANGLE(DEG S.) BUS P BUS Q");for(i=0;i<n;i++){v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);x=e[i];y=f[i];z=y/x;angle[i]=atan(z);angle[i]/=Pi;fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[ i],p[i],q[i]);}fprintf(fp,"\n LINE FLOW ");for(i=1;i<=l;i++){pos1=f1(i);st=s1[pos1];en=e1[pos1];x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)];x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)];y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)];y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)];p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1];q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1];p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1];q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1];for(j=1;j<=n;j++){cm=c[f2(j,i,l)];if(cm!=0.0){km=1;if(en==j){km=2;}if(km==1){q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}else{q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm;}}}p3[pos1]=p1[pos1]+p2[pos1] ;q3[pos1]=q1[pos1]+q2[pos1] ;fprintf(fp,"\n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d%11d%1 3.6e%13.6e",\i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1 ],\e1[pos1],s1[pos1],p2[pos1],q2[pos1]);}}。

电力系统潮流计算程序

电力系统潮流计算c语言程序，两行，大家可以看看，仔细研究，然后在这个基础上修改。

谢谢#include "stdafx.h"#include <iostream>#include <fstream>#include <process.h>#include"Complex.h"#include"wanjing.h"#include"gauss.h"using namespace std;int _tmain(int argc, _TCHAR* argv[]){int i; //i作为整个程序的循环变量int N=Bus::ScanfBusNo(); //输入节点个数int L=Line::ScanflineNo(); //输入支路个数if((L&&N)==0){return 0;} //如果找不到两个文件中的任意一个，退出Line *line=new Line[L]; //动态分配支路结构体Line::ScanfLineData(line); //输入支路参数Line::PrintfLineData(line,L); //输出支路参数Bus *bus=new Bus[N]; //动态分配结点结构体for(int i=0;i<N;i++){bus[i].Sdelta.real=0;bus[i].Sdelta.image=0;}Bus::ScanfBusData(bus); //输入节点参数Bus::PrintfBusData(bus,N); //输出结点参数Complex **X; X=new Complex *[N]; for(i=0;i<N;i++) {X[i]=new Complex[N];}//动态分配二维数组作节点导纳阵Bus::JisuanNodeDnz(X,line,bus,L,N); //计算节点导纳矩阵Bus::PrintfNodeDnz(X,N); //输出节点导纳矩阵int NN=(N-1)*2;double **JacAug; JacAug=new double *[NN]; for(i=0;i<NN;i++) {JacAug[i]=new double[NN+1];}//为雅可比增广动态分配一个二维数组double *x;x=new double[NN];int count=1;LOOP:Bus::JisuanNodeI(X,bus,N); //计算节点注入电流Bus::JisuanNodeScal(X,bus,N); //计算节点功率Bus::JisuanNodeScal(X,bus,N);//计算节点功率Bus::JisuanNodeSdelta(bus,N); //计算节点功率差值Bus::PrintfNodeScal(X,bus,N);//输出节点功率差值int icon=wehcon1(bus,N);//whether converbence看迭代是否结束if(icon==1){cout<<"icon="<<icon<<"进行第"<<count<<"次迭代"<<endl;Bus::JisuanJacAug(JacAug,X,bus,N); //计算雅可比增广矩阵// Bus::PrintfJacAug(JacAug,N);gauss::gauss_slove(JacAug,x,NN);//解方程组求出电压差值Bus::ReviseNodeV(bus,x,N);//修正节点电压// Bus::PrintfNodeV(bus,N);count++;goto LOOP;}else{for(i=0;i<L;i++)//计算输电线路上的传输功率{int statemp,endtemp;Complex aa,bb,cc,dd,B;B.real=0;B.image=-line[i].B;statemp=line[i].start;endtemp=line[i].end;aa=Complex::productComplex (Complex::getconj(bus[statemp-1].V), B);bb=Complex::subComplex (Complex::getconj(bus[statemp-1].V), Complex::getconj(bus[endtemp-1].V));cc=Complex::productComplex (bb , Complex::getconj(line[i].Y));dd=Complex::CaddC(aa,cc);line[i].stoe=Complex::productComplex(b us[statemp-1].V,dd);aa=Complex::productComplex (Complex::getconj(bus[endtemp-1].V), B);bb=Complex::subComplex (Complex::getconj(bus[endtemp-1].V), Complex::getconj(bus[statemp-1].V));cc=Complex::productComplex (bb , Complex::getconj(line[i].Y));dd=Complex::CaddC(aa,cc);line[i].etos=Complex::productComplex(b us[endtemp-1].V,dd);}cout<<"icon="<<icon<<",迭代结束。

潮流计算C++程序

潮流计算C++程序#include <fstream> #include<iomanip> #include<math、h> using namespace std; struct node{//节点类 int i;//节点编号 double U,P,Q,delta;//额定电压计算负荷电压相角 }; struct line{//线路类连接父节点子节点 nodef_node,s_node;//父节点子节点 double R,X,B;//线路参数R X B/2 double P_in,Q_in,P_out,Q_out,d_P,d_Q,D_U,d_U;//线路输入输出功率以及线路消耗功率 void Set_node(node nod1,node nod2){ f_node=nod1; s_node=nod2; } }; void fun1(line&lin){//由后往前递推功率 double p=lin、P_out; doubleq=lin、Q_out; double u=lin、s_node、U; lin、d_P=(p*p+q*q)/u/u*lin、R; lin、d_Q=(p*p+q*q)/u/u*lin、X; lin、P_in=lin、d_P+lin、P_out; lin、Q_in=lin、d_Q+lin、Q_out; }; void fun2(line &lin){//由前往后推电压 doublep=lin、P_in; double q=lin、Q_in; double u=lin、f_node、U; lin、D_U=(p*lin、R+q*lin、X)/u; lin、d_U=(p*lin、X-q*lin、R)/u; lin、s_node、U=sqrt(pow(lin、f_node、U-lin、D_U,2)+pow(lin、d_U,2));//子节点电压 lin、s_node、delta=lin、f_node、delta-atan(lin、d_U/(lin、f_node、U-lin、D_U)); }; void fun3(line &lin){//由前往后推电压不计横向分量 double p=lin、P_in; double q=lin、Q_in; doubleu=lin、f_node、U; lin、D_U=(p*lin、R+q*lin、X)/u; lin、d_U=(p*lin、X-q*lin、R)/u; lin、s_node、U=lin、f_node、U-lin、D_U,2;//子节点电压 lin、s_node、delta=lin、f_node、delta-atan(lin、d_U/(lin、f_node、U-lin、D_U)); }; void main(){ int num_l; int num_n;//支路数节点数 ifstream fin; fin、open("E:\\data、txt"); fin>>num_n>>num_l;//输入节点数支路数 ofstream fout; fout、open("E:\\databak、txt"); node *nod; nod=new node[num_n];//节点数目 line *lin;lin=new line[num_l];//线路数目 nod[0]、delta=0; double*u;//节点额定电压 u=new double[num_n]; for(inti=0;i<num_n;i++){ fin>>u[i]; }; double *p;//节点有功功率p=new double[num_n]; for(inti=0;i<num_n;i++){ fin>>p[i]; }; double *q;//节点无功功率q=new double[num_n]; for(inti=0;i<num_n;i++){ fin>>q[i]; }; for(inti=0;i<num_n;i++){//设定节点标号参数 nod[i]、i=i;nod[i]、P=p[i]; nod[i]、Q=q[i]; nod[i]、U=u[i]; }; double *r;//线路电阻 r=new double[num_l]; for(inti=0;i<num_l;i++){ fin>>r[i]; }; double *x;//线路电抗x=new double[num_l]; for(inti=0;i<num_l;i++){ fin>>x[i]; }; double *b;//线路电纳b=new double[num_l]; for(inti=0;i<num_l;i++){ fin>>b[i]; }; for(inti=0;i<num_l;i++){//设定线路参数 lin[i]、R=r[i]; lin[i]、X=x[i]; lin[i]、B=b[i]; }; for(int i=0;i<num_l;i++){//确定线路节点关系 fin>>lin[i]、f_node、i; fin>>lin[i]、s_node、i; }; for(int i=0;i<num_l;i++){//计算节点运算负荷nod[lin[i]、f_node、i]、Q-=lin[i]、B*nod[lin[i]、f_node、i]、U*nod[lin[i]、f_node、i]、U; nod[lin[i]、s_node、i]、Q-=lin[i]、B*nod[lin[i]、s_node、i]、U*nod[lin[i]、s_node、i]、U; }; double *P_c,*Q_c;//保存运算负荷数据P_c=new double[num_n]; Q_c=new double[num_n]; for(inttc=0;tc<num_n;tc++){ P_c[tc]=nod[tc]、P;Q_c[tc]=nod[tc]、Q; }; for(int i=0;i<num_l;i++){//设定线路节点数据 lin[i]、Set_node(nod[lin[i]、f_node、i],nod[lin[i]、s_node、i]); }; for(int re=0;re<3;re++){//迭代运算开始 fout<<"第"<<re+1<<"次迭代"<<endl; doubleto_P(0),to_dP(0); for(inti=1;i<num_n;i++){ to_P+=nod[i]、P; }; for(inttc=0;tc<num_n;tc++){//重置运算负荷 nod[tc]、P=P_c[tc]; nod[tc]、Q=Q_c[tc]; }; for(int ts=0;ts<num_l;ts++){//置各线路初始输出功率为子节点运算负荷 lin[ts]、P_out=lin[ts]、s_node、P; lin[ts]、Q_out=lin[ts]、s_node、Q; }; for(int i=0;i<num_l;i++){//设定线路节点数据 lin[i]、Set_node(nod[lin[i]、f_node、i],nod[lin[i]、s_node、i]); }; int j=num_l-1;//反向求各支路功率损耗和功率分布for(j=num_l-1;j>0;j--){ fun1(lin[j]); nod[lin[j]、f_node、i]、P+=lin[j]、P_in; nod[lin[j]、f_node、i]、Q+=lin[j]、Q_in; for(int i=0;i<num_l;i++){//设定线路节点数据 lin[i]、Set_node(nod[lin[i]、f_node、i],nod[lin[i]、s_node、i]); }; for(int ts=0;ts<num_l;ts++){ lin[ts]、P_out=lin[ts]、s_node、P; lin[ts]、Q_out=lin[ts]、s_node、Q; }; }; fun1(lin[j]); for(inti=0;i<num_l;i++){//设定线路节点数据 lin[i]、Set_node(nod[lin[i]、f_node、i],nod[lin[i]、s_node、i]); }; int t=0;//求线路各点电压 for(t=0;t<num_l-1;t++){ fun2(lin[t]); nod[lin[t]、s_node、i]、U=lin[t]、s_node、U; nod[lin[t]、s_node、i]、delta=lin[t]、s_node、delta; for(int i=0;i<num_l;i++){//设定线路节点数据lin[i]、Set_node(nod[lin[i]、f_node、i],nod[lin[i]、s_node、i]); }; }; fun2(lin[t]); nod[lin[t]、s_node、i]、U=lin[t]、s_node、U; nod[lin[t]、s_node、i]、delta=lin[t]、s_node、delta; fout<<"支路信息:"<<endl;for(int i=0;i<num_l;i++){//输出线路信息 fout<<"支路"<<lin[i]、f_node、i+1<<"-"<<lin[i]、s_node、i+1<<":"<<endl; fout<<"始端功率:"<<lin[i]、P_in<<"+j"<<lin[i]、Q_in<<endl; fout<<"末端功率:"<<lin[i]、P_out<<"+j"<<lin[i]、Q_out<<endl; fout<<"功率损耗:"<<lin[i]、d_P<<"+j"<<lin[i]、d_Q<<endl; fout<<"电压损耗"<<lin[i]、f_node、U-lin[i]、s_node、U<<endl; };fout<<"节点信息:"<<endl; for(int i=0;i<num_n;i++){//输出节点信息 fout<<"节点"<<i+1<<endl; fout<<"电压:"<<nod[i]、U<<"相角:"<<nod[i]、delta*180/3、14<<endl; }; double *lu;//求最低电压及最低电压点lu=new double[num_n]; int *lua; lua=new int[num_n];for(int i=0;i<num_n;i++){ lu[i]=nod[i]、U; lua[i]=i; };for(int i=0;i<num_n-1;i++){ if(lu[i]<lu[i+1]){ double st; int a; st=lu[i]; a=lua[i]; lu[i]=lu[i+1]; lua[i]=lua[i+1]; lu[i+1]=st; lua[i+1]=a; }; }; for(inti=0;i<num_l;i++){ to_dP+=lin[i]、d_P; }; fout<<"全网信息:"<<endl; fout<<"总电源有功:"<<lin[0]、P_in<<endl;fout<<"总负荷有功:"<<to_P<<endl; fout<<"总有功损耗:"<<to_dP<<endl; fout<<"网损率:"<<to_dP/(to_P+lin[0]、P_in)<<endl; fout<<"最低电压:"<<lu[num_n-1]<<"最低电压点:"<<lua[num_n-1]+1<<endl; delete[]lu; delete[]lua; }; delete[]nod; delete[]lin; delete[]u; delete[]p; delete[]q; delete[]r; delete[]x; delete[]b; }附：3-4的data文件32113100 0 0、1720 01、7158、51、2220、520、2 0、 0 0112输出为：第1次迭代支路信息: 支路1-2: 始端功率:1、02165+j0、末端功率:1、00434+j0、功率损耗:0、+j0、电压损耗0、支路2-3: 始端功率:0、5034+j0、3068 末端功率:0、5+j0、3 功率损耗:0、0034+j0、0068 电压损耗0、109 支路2-4: 始端功率:0、+j0、末端功率:0、2+j0、15 功率损耗:0、+j0、电压损耗0、节点信息: 节点1 电压:10、5相角:0 节点2 电压:10、2259相角:-0、节点3 电压:10、1169相角:-1、25281 节点4 电压:10、152相角:-1、072 全网信息: 总电源有功:1、02165 总负荷有功:1 总有功损耗:0、网损率:0、最低电压:10、1169最低电压点:3第2次迭代支路信息: 支路1-2: 始端功率:1、02078+j0、69156 末端功率:1、00423+j0、功率损耗:0、+j0、电压损耗0、支路2-3: 始端功率:0、+j0、末端功率:0、5+j0、3 功率损耗:0、+j0、电压损耗0、支路2-4: 始端功率:0、20211+j0、末端功率:0、2+j0、15 功率损耗:0、+j0、电压损耗0、节点信息: 节点1 电压:10、5相角:0 节点2 电压:10、2264相角:-0、86489 节点3 电压:10、1175相角:-1、25273 节点4 电压:10、1525相角:-1、07194 全网信息: 总电源有功:1、02078 总负荷有功:1、70434 总有功损耗:0、网损率:0、最低电压:10、1175最低电压点:3第3次迭代支路信息: 支路1-2: 始端功率:1、02078+j0、末端功率:1、00423+j0、功率损耗:0、+j0、电压损耗0、支路2-3: 始端功率:0、+j0、末端功率:0、5+j0、3 功率损耗:0、+j0、电压损耗0、支路2-4: 始端功率:0、20211+j0、末端功率:0、2+j0、15 功率损耗:0、+j0、电压损耗0、节点信息: 节点1 电压:10、5相角:0 节点2 电压:10、2264相角:-0、86489 节点3 电压:10、1175相角:-1、25273 节点4 电压:10、1525相角:-1、07194 全网信息: 总电源有功:1、02078 总负荷有功:1、70423 总有功损耗:0、网损率:0、最低电压:10、1175最低电压点:3例3-2的data：4310、5010 0 0、3 0、5 0、2 0 0、2 0、3 0、151、211、52、423 0 0 0 011213输出：第1次迭代支路信息:支路1-2:始端功率:1、02165+j0、末端功率:1、00434+j0、功率损耗:0、+j0、电压损耗0、支路2-3:始端功率:0、5034+j0、3068末端功率:0、5+j0、3功率损耗:0、0034+j0、0068电压损耗0、支路2-4:始端功率:0、+j0、末端功率:0、2+j0、15功率损耗:0、+j0、电压损耗0、节点信息:节点1电压:10、5相角:0节点2电压:10、2248相角:-0、节点3电压:10、1155相角:-1、2529节点4电压:10、1507相角:-1、07205全网信息:总电源有功:1、02165总负荷有功:1总有功损耗:0、网损率:0、最低电压:10、1155最低电压点:3第2次迭代支路信息:支路1-2:始端功率:1、02078+j0、69157末端功率:1、00423+j0、功率损耗:0、+j0、电压损耗0、支路2-3:始端功率:0、+j0、末端功率:0、5+j0、3功率损耗:0、+j0、电压损耗0、支路2-4:始端功率:0、20211+j0、15182末端功率:0、2+j0、15功率损耗:0、+j0、电压损耗0、节点信息:节点1电压:10、5相角:0节点2电压:10、2253相角:-0、86489节点3电压:10、1161相角:-1、25282节点4电压:10、1513相角:-1、07199全网信息:总电源有功:1、02078总负荷有功:1、70434总有功损耗:0、网损率:0、最低电压:10、1161最低电压点:3第3次迭代支路信息:支路1-2:始端功率:1、02078+j0、末端功率:1、00423+j0、功率损耗:0、+j0、电压损耗0、支路2-3:始端功率:0、+j0、末端功率:0、5+j0、3功率损耗:0、+j0、电压损耗0、支路2-4:始端功率:0、20211+j0、15182末端功率:0、2+j0、15功率损耗:0、+j0、电压损耗0、节点信息:节点1电压:10、5相角:0节点2电压:10、2253相角:-0、86489节点3电压:10、1161相角:-1、25282节点4电压:10、1513相角:-1、07199全网信息:总电源有功:1、02078总负荷有功:1、70423总有功损耗:0、网损率:0、最低电压:10、1161最低电压点:3。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

《电力系统潮流上机》课程设计报告院系：电气与电子工程学院班级：电气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++){。