潮流计算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;

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语言程序

/*利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。

所有参数应归算至标幺值下。

/*可计算最大节点数为100，可计算PQ,PV,平衡节点*//*可计算非标准变比和平行支路*/#include<stdio.h>#include<math.h>#include<stdlib.h>#define M 100 /*最大矩阵阶数*/#define Nl 100 /*迭代次数*/int i,j,k,a,b,c; /*循环控制变量*/int t,l;double P,Q,H,J; /*中间变量*/int n, /*节点数*/m, /*支路数*/pq, /*PQ节点数*/pv; /*PV节点数*/double eps; /*迭代精度*/double aa[M],bb[M],cc[M],dd[M],max, rr,tt; /*中间变量*/double mo,c1,d1,c2,d2; /*复数运算函数的返回值*/double G[M][M],B[M][M],Y[M][M]; /*节点导纳矩阵中的实部、虚部及其模方值*/double ykb[M][M],D[M],d[M],dU[M]; /*雅克比矩阵、不平衡量矩阵*/struct jd /*节点结构体*/{ int num,ty; /* num为节点号，ty为节点类型*/ double p,q,S,U,zkj,dp,dq,du,dj; /*节点有功、无功功率，功率模值，电压模值,阻抗角牛顿--拉夫逊中功率不平衡量、电压修正量*/} jd[M];struct zl /*支路结构体*/{ int numb; /*numb为支路号*/int p1,p2; /*支路的两个节点*/double kx; /*非标准变比*/double r,x; /*支路的电阻与电抗*/} zl[M];FILE *fp1,*fp2;void data() /* 读取数据*/{int h,number;fp1=fopen("input.txt","r");fscanf(fp1,"%d,%d,%d,%d,%lf\n",&n,&m,&pq,&pv,&eps); /*输入节点数,支路数,PQ节点数,PV节点数和迭代精度*/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[i].p,&jd[i].q,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;}if(h==2) /*类型h=2是pv节点*/{fscanf(fp1,",%lf,%lf,%lf\n",&jd[i].p,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;jd[i].q=-1.567;}if(h==3) /*类型h=3是平衡节点*/ {fscanf(fp1,",%lf,%lf\n",&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;}}for(i=1;i<=m;i++) /*输入支路阻抗*/ fscanf(fp1,"%d,%lf,%d,%d,%lf,%lf\n",&zl[i].numb,&zl[i].kx,&zl[i].p1,&zl[i].p2,&zl[i].r,&zl[i].x);fclose(fp1);if((fp2=fopen("output.txt","w"))==NULL){printf(" can not open file!\n");exit(0);}fprintf(fp2," 电力系统潮流计算\n ");fprintf(fp2," ********** 原始数据*********\n");fprintf(fp2,"====================================================================== ==========\n");fprintf(fp2," 节点数:%d 支路数:%d PQ节点数:%d PV节点数:%d 精度:%f\n",n,m,pq,pv,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].U,jd[n].num,jd[n].zkj);fprintf(fp2," -------------------------------------------------------------------------------\n");for(i=1;i<=m;i++)fprintf(fp2," 支路%d: 相关节点:%d,%d 非标准变比:kx=%f R=%f X=%f \n",i,zl[i].p1,zl[i].p2,zl[i].kx,zl[i].r,zl[i].x);fprintf(fp2,"==============================================================================\n ");}/*------------------------------------复数运算函数--------------------------------------*/double mozhi(double a0,double b0) /*复数求模值函数*/{ mo=sqrt(a0*a0+b0*b0);return mo;}double ji(double a1,double b1,double a2,double b2) /*复数求积函数a1为电压模值，a2为阻抗角，a3为导纳实部，a4为导纳虚部*/{ a1=a1*cos(b1);b1=a1*tan(b1);c1=a1*a2-b1*b2;d1=a1*b2+a2*b1;return c1;return d1;}double shang(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;}/*--------------------------------计算节点导纳矩阵----------------------------------*/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<=n;i++) /*节点导纳矩阵的主对角线上的元素,非对地导纳加入相应的主对角线元素中(考虑非标准变比）*/for(j=1;j<=m;j++)if(zl[j].p1==i){if(zl[j].kx==1){mozhi(zl[j].r,zl[j].x);if(mo==0) continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2;B[i][i]+=d2;}else{mozhi(zl[j].r,zl[j].x);if(mo==0) continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2/zl[j].kx+c2*(1-zl[j].kx)/(zl[j].kx*zl[j].kx);B[i][i]+=d2/zl[j].kx+d2*(1-zl[j].kx)/(zl[j].kx*zl[j].kx);}}else if(zl[j].p2==i){if(zl[j].kx==1){mozhi(zl[j].r,zl[j].x);if(mo==0) continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2;B[i][i]+=d2;}else{mozhi(zl[j].r,zl[j].x);if(mo==0) continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2/zl[j].kx+c2*(zl[j].kx-1)/zl[j].kx;B[i][i]+=d2/zl[j].kx+d2*(zl[j].kx-1)/zl[j].kx;}}for(k=1;k<=m;k++) /*节点导纳矩阵非主对角线上（考虑非标准变比）的元素*/if(zl[k].kx==1){i=zl[k].p1;j=zl[k].p2;mozhi(zl[k].r,zl[k].x);if(mo==0) continue;shang(1,0,zl[k].r,zl[k].x);G[i][j]-=c2;B[i][j]-=d2;G[j][i]=G[i][j];B[j][i]=B[i][j];}else{i=zl[k].p1;j=zl[k].p2;mozhi(zl[k].r,zl[k].x);if(mo==0) continue;shang(1,0,zl[k].r,zl[k].x);G[i][j]-=c2/zl[k].kx;B[i][j]-=d2/zl[k].kx;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 Calculate_Unbalanced_Para(){for(i=1;i<=n;i++){if(jd[i].ty==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;}P=Q=0;P=jd[a].U*(G[t][j]*cos(jd[i].zkj-jd[a].zkj)+B[t][j]*sin(jd[i].zkj-jd[a].zkj));Q=jd[a].U*(G[t][j]*sin(jd[i].zkj-jd[a].zkj)-B[t][j]*cos(jd[i].zkj-jd[a].zkj));cc[t]+=P;dd[t]+=Q;}jd[i].dp=jd[i].p-jd[i].U*cc[t];jd[i].dq=jd[i].q-jd[i].U*dd[t];}if(jd[i].ty==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;}P=Q=0;P=jd[a].U*(G[t][j]*cos(jd[i].zkj-jd[a].zkj)+B[t][j]*sin(jd[i].zkj-jd[a].zkj));Q=jd[a].U*(G[t][j]*sin(jd[i].zkj-jd[a].zkj)-B[t][j]*cos(jd[i].zkj-jd[a].zkj));cc[t]+=P;dd[t]+=Q;}jd[i].dp=jd[i].p-jd[i].U*cc[t];jd[i].q=jd[i].U*dd[t];}}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<=n-1;i++){D[pq+i]=jd[i].dp;}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<=n-1;i++){t=jd[i].num;fprintf(fp2,"\n dp[%d]=%f",i,D[pq+t]);}}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 Ui=jd[i].U;double Uj=jd[j].U;double zi=jd[i].zkj;double zj=jd[j].zkj;if(j<=pq) /*求j<=pq时的H、N、J、L */{if(i!=j) /*求i!=j时的H、N、J、L*/{ykb[2*i-1][2*j-1]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ykb[2*i-1][2*j]=Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* N */ykb[2*i][2*j-1]=-ykb[2*i-1][2*j]; /* J */ykb[2*i][2*j]=ykb[2*i-1][2*j-1]; /* 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;H=J=0;H=jd[k].U*(G[i1][k1]*sin(jd[i].zkj-jd[k].zkj)-B[i1][k1]*cos(jd[i].zkj-jd[k].zkj));J=jd[k].U*(G[i1][k1]*cos(jd[i].zkj-jd[k].zkj)+B[i1][k1]*sin(jd[i].zkj-jd[k].zkj));aa[i]=aa[i]+H;bb[i]=bb[i]+J;}ykb[2*i-1][2*j-1]=-Ui*(aa[i]-Ui*(G[i1][i1]*sin(jd[i].zkj-jd[i].zkj)-B[i1][i1]*cos(jd[i].zkj-jd[i].zkj)));/*H*/ykb[2*i][2*j-1]=Ui*(bb[i]-Ui*(G[i1][i1]*cos(jd[i].zkj-jd[i].zkj)+B[i1][i1]*sin(jd[i].zkj-jd[i].zkj)));/*J*/ykb[2*i-1][2*j]=ykb[2*i][2*j-1]+2*Ui*Ui*G[i1][i1];/*N*/ykb[2*i][2*j]=-ykb[2*i-1][2*j-1]-2*Ui*Ui*B[i1][i1];/*L*/}}else /*求j>pq时的H、J */{ ykb[2*i-1][pq+j]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ykb[2*i][pq+j]=-Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* J */}}for(i=pq+1;i<=n-1;i++) /*形成pv节点子阵*/ for(j=1;j<n;j++){int i1=jd[i].num;int j1=jd[j].num;double Ui=jd[i].U;double Uj=jd[j].U;double zi=jd[i].zkj;double zj=jd[j].zkj;if(j<=pq) /*求j<=pq时的H、N */{ykb[pq+i][2*j-1]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ykb[pq+i][2*j]=Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* N */ }else /*求j>pq时的H*/{if(i!=j) /*求i!=j时的H*/ykb[pq+i][pq+j]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */else /*求i=j时的H*/{aa[i]=0;for(k=1;k<=n;k++){int k1=jd[k].num;H=0;H=jd[k].U*(G[i1][k1]*sin(jd[i].zkj-jd[k].zkj)-B[i1][k1]*cos(jd[i].zkj-jd[k].zkj));aa[i]=aa[i]+H;}ykb[pq+i][pq+j]=-Ui*(aa[i]-Ui*(G[i1][i1]*sin(jd[i].zkj-jd[i].zkj)-B[i1][i1]*cos(jd[i].zkj-jd[i].zkj))); /*H*/ }}}/*------------------------------输出雅克比矩阵--------------------------------*/fprintf(fp2,"\n\n 雅克比矩阵为: ");for(i=1;i<=(2*pq+pv);i++){fprintf(fp2,"\n");for(j=1;j<=2*pq+pv;j++){fprintf(fp2," %f",ykb[i][j]);}}}void Solve_Equations() /* 求解修正方程组(LU分解法)*/{double l[Nl][Nl]={0}; //定义L矩阵double u[Nl][Nl]={0}; //定义U矩阵double y[Nl]={0}; //定义数组Ydouble x[Nl]={0}; //定义数组Xdouble a[Nl][Nl]={0}; //定义系数矩阵double b[Nl]={0}; //定义右端项double sum=0;int i,j,k,s;int n;n=2*pq+pv;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=ykb[i+1][j+1];}}for(i=0; i<n; i++){b[i]=D[i+1];}for(i=0; i<n; i++) /*初始化矩阵l*/{for(j=0; j<n; j++){if(i==j) l[i][j] = 1;}}for(i=0;i<n;i++) /*开始LU分解*/{ u[0][i]=(float)(a[0][i]);} /*第一步：对矩阵U的首行进行计算*/ for(k=0;k<n-1;k++) /*第二步：逐步进行LU分解*/ { for(i=k+1;i<n;i++) /*对L的第k列进行计算*/{ for(s=0,sum=0;s<n;s++){ if(s!=k)sum+=l[i][s]*u[s][k];}l[i][k]=(float)((a[i][k]-sum)/u[k][k]);}for(j=k+1;j<n;j++) /*对U的第k+1行进行计算*/ { for(s=0,sum=0;s<n;s++){ if(s!=k+1)sum+=l[k+1][s]*u[s][j];}u[k+1][j]=(float)((a[k+1][j]-sum));}}y[0]=b[0] ; /*回代法计算数组Y*/for(i=1;i<n;i++){ for(j=0,sum=0;j<i;j++){ sum+=y[j]*l[i][j];}y[i]=(float)(b[i]-sum);}x[n-1]=(float)(y[n-1]/u[n-1][n-1]); /*回代法计算数组X*/for(i=n-2;i>=0;i--){ for(j=n-1,sum=0;j>i;j--){ sum+=x[j]*u[i][j];}x[i]=(float)((y[i]-sum)/u[i][i]);}for(i=1; i<=n; i++){d[i]=x[i-1];}max=fabs(d[1]); /*选出最大的修正量的值*/for(i=1;i<=n;i++)if((fabs(d[i]))>max)max=fabs(d[i]);}void Niudun_Lafuxun(){int z=1;fprintf(fp2,"\n --------------极坐标形式的牛顿-拉夫逊迭代法结果:--------------\n");do{ max=1;if((z<Nl)&&(max>=eps)){fprintf(fp2,"\n 迭代次数: %d\n",z); /*开始迭代计算*/ }Calculate_Unbalanced_Para();Form_Jacobi_Matric();Solve_Equations();for(i=1;i<=pq;i++){jd[i].zkj+=d[2*i-1];jd[i].U+=d[2*i]*jd[i].U;}for(i=pq+1;i<=n-1;i++){jd[i].zkj+=d[pq+i];}fprintf(fp2,"\n\n 输出dδ,dU: "); /*输出修正量的值*/for(c=1;c<=n;c++){for(a=1;a<=n;a++){if(jd[a].num==c)break;}fprintf(fp2,"\n");if(jd[a].ty==1)fprintf(fp2," 节点为%2d dδ=%8.5f dU=%8.5f",c,d[2*a-1],d[2*a]);if(jd[a].ty==2)fprintf(fp2," 节点为%2d dδ=%8.5f",c,d[pq+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].U);fprintf(fp2,"∠%f",jd[a].zkj);}fprintf(fp2,"\n\n ------------------------------------------------------------------------------");z++;} while((z<Nl)&&(max>=eps)); /*判断是否达到精度要求*/}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].U);fprintf(fp2,"∠%f",jd[a].zkj);}rr=tt=0; /*计算节点的注入功率*/for(i=1;i<=n;i++){int i1=jd[i].num;ji(jd[i].U,-jd[i].zkj,G[n1][i1],-B[n1][i1]);rr+=c1;tt+=d1;}ji(jd[n].U,jd[n].zkj,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); /*PQ节点注入功率*/ 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<=n-1;i++){int i1=jd[i].num;fprintf(fp2," PV节点: 节点%d S[%d]=%f", i1,i1,jd[i].p); /*PV节点注入功率*/ 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=zl[i].p1; /*计算线路功率*/int j1=zl[i].p2;aa[i]=bb[i]=P=Q=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;}mozhi(zl[i].r,zl[i].x);if(mo==0)continue;shang(1,0,zl[i].r,zl[i].x);ji(jd[a].U/zl[i].kx/zl[i].kx,-jd[a].zkj,c2,-d2); /*考虑非标准变比*/P+=c1;Q+=d1;ji(jd[b].U/zl[i].kx,-jd[b].zkj,c2,-d2);P-=c1;Q-=d1;ji(jd[a].U,jd[a].zkj,P,Q);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;P=Q=0;ji(jd[b].U,-jd[b].zkj,c2,-d2); /*考虑非标准变比*/ P+=c1;Q+=d1;ji(jd[a].U/zl[i].kx,-jd[a].zkj,c2,-d2);P-=c1;Q-=d1;ji(jd[b].U,jd[b].zkj,P,Q);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=zl[i].p1;int j1=zl[i].p2;fprintf(fp2," 线路%d损耗的功率为: %f",i,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(){printf("请仔细阅读附录里的输入文件使用说明以便于你正确输入!\n");data(); /*读取数据*/Form_Y(); /*形成节点导纳矩阵*/for(i=1;i<=pq;i++) /* U、zkj附初值*/{jd[i].U=1; jd[i].zkj=0;}for(i=pq+1;i<n;i++){jd[i].U=jd[i].U; jd[i].zkj=0;}Niudun_Lafuxun(); /*牛顿--拉夫逊迭代*/Powerflow_Result();printf("潮流计算结果存放于output.txt文件中!\n");}附录：输入文件按以下规则输入节点数，支路数，PQ节点数，PV节点数，精度节点编号，节点类型（1为PQ节点，2为PV节点，3为平衡节点）节点输入有功功率，无功功率(节点电压纵分量，横分量)支路编号，非标准变比（若没有则填1即可），相关节点，相关节点，支路电阻，支路电抗。

电力系统潮流计算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程序精编版

电力系统潮流计算注：这是一个基于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语言进行复杂电力系统的潮流计算(1)

用C 语言进行复杂电力系统的潮流计算杨汉明电力系统的潮流计算是对电力系统分析的最基本步骤也是最重要的步骤，是指在一定的系统结构和运行条件下，确定系统运行状态的计算，也即是对各母线（节点）电压，各元件（支路）传输电线或功率的计算。

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

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

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

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

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

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

C语言潮流计算研究报告

潮流计算研究报告一、潮流计算概述电力工业是国民经济的基础产业,为保证电力系统安全、稳定及可靠的运行,首先需计算电力系统网络潮流。

潮流计算指根据给定的电力系统运行条件及系统接线情况确定整个电力系统各部分的运行状态:各母线电压、各元件流过的功率及系统功率损耗等等。

在电力系统规划、设计及运行方式的研究中,都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性及经济性,所以潮流计算是研究电力系统的一种很重要和很基础的计算。

在数学上潮流计算问题是一组多元非线性方程式求解问题,其解法离不开迭代。

二、潮流计算算法框图三、节点分类本次设计采用直角坐标牛顿拉夫逊法，用C语言编制潮流计算程序。

牛顿法拉夫逊：把非线性方程式的求解过程变成反复对相应的线性方程式的求解过程，通常称为逐次线性化过程。

节点分类：1）平衡节点已知节点的电压幅值和相角，求有功和无功，一般设在调频厂（一个）。

2）PQ节点已知节点的有功、无功，求电压幅值和相角，一般设在变电所母线，按给定频率发电的电厂母线（多个）。

3）PV节点已知节点的有功和电压幅值，求无功和电压相角，一般设在有无功补偿的变电所母线，无功可调的发电厂母线（少个）。

四、基本公式1、节点注入电流2、节点功率3、不平衡量4、雅可比矩阵（1）i ≠ j时（2）i=j 时5、线路功率及损耗五、IEEE30拓扑结构IEEE30节点结构六、算法输入原始数据1、支路参数（支路数L=41）始端末端回线数线路类型线路电阻线路电抗线路电纳(B/2)/ 变比2 1 1 0 0.0192 0.0575 0.0528 03 1 1 0 0.0452 0.1852 0.0408 04 2 1 0 0.057 0.1737 0.0368 04 3 1 0 0.0132 0.0379 0.0084 05 2 1 0 0.0472 0.1983 0.0418 06 2 1 0 0.0581 0.1763 0.0374 06 4 1 0 0.0119 0.0414 0.009 07 5 1 0 0.046 0.116 0.0204 07 6 1 0 0.0267 0.082 0.017 08 6 1 0 0.012 0.042 0.009 06 9 1 1 0 0.208 0 0.9786 10 1 1 0 0.556 0 0.96911 9 1 0 0 0.208 0 010 9 1 0 0 0.11 0 04 12 1 1 0 0.256 0 0.93213 12 1 0 0 0.14 0 014 12 1 0 0.1231 0.2559 0 015 12 1 0 0.0662 0.1304 0 016 12 1 0 0.0945 0.1987 0 015 14 1 0 0.221 0.1997 0 017 16 1 0 0.0524 0.1923 0 018 15 1 0 0.1073 0.2185 0 019 18 1 0 0.0639 0.1292 0 020 19 1 0 0.034 0.068 0 020 10 1 0 0.0936 0.209 0 017 10 1 0 0.0324 0.0845 0 021 10 1 0 0.0348 0.0749 0 022 10 1 0 0.0727 0.1499 0 022 21 1 0 0.0116 0.0236 0 023 15 1 0 0.1 0.202 0 024 22 1 0 0.115 0.179 0 024 23 1 0 0.132 0.27 0 025 24 1 0 0.1885 0.3292 0 026 25 1 0 0.2544 0.38 0 027 25 1 0 0.1093 0.2087 0 028 27 1 1 0 0.396 0 0.96829 27 1 0 0.2198 0.4153 0 030 27 1 0 0.3202 0.6027 0 030 29 1 0 0.2399 0.4533 0 028 8 1 0 0.0636 0.2 0.0428 028 6 1 0 0.0169 0.0599 0.013 02、节点参数（节点数N=30）编号节点类型电压实部电压虚部负荷有功负荷无功发电有功发电无功PV节点的V 并联电容1 1 1.05 0 0 0 0 0 1.05 02 3 1.045 0 0.217 0.127 0.8 0 1.045 03 2 1 0 0.024 0.012 0 0 1 04 2 1 0 0.076 0.016 0 0 1 05 3 1.01 0 0.942 0.190 0.5 0 1.01 06 2 1 0 0 0 0 0 1 07 2 1 0 0.228 0.109 0 0 1 08 3 1.01 0 0.30 0.30 0.2 0 1.01 09 2 1 0 0 0 0 0 1 010 2 1 0 0.058 0.02 0 0 1 0.1911 3 1.05 0 0 0 0.2 0 1.05 012 2 1 0 0.112 0.075 0 0 1 013 3 1.05 0 0 0 0.2 0 1.05 014 2 1 0 0.062 0.016 0 0 1 015 2 1 0 0.082 0.025 0 0 1 016 2 1 0 0.035 0.018 0 0 1 017 2 1 0 0.09 0.058 0 0 1 018 2 1 0 0.032 0.009 0 0 1 019 2 1 0 0.095 0.034 0 0 1 020 2 1 0 0.022 0.007 0 0 1 021 2 1 0 0.175 0.112 0 0 1 022 2 1 0 0 0 0 0 1 023 2 1 0 0.032 0.016 0 0 1 024 2 1 0 0.087 0.067 0 0 1 0.04325 2 1 0 0 0 0 0 1 026 2 1 0 0.035 0.023 0 0 1 027 2 1 0 0 0 0 0 1 028 2 1 0 0 0 0 0 1 029 2 1 0 0.024 0.009 0 0 1 0302 1 0 0.106 0.019 0 0 1 0七、潮流程序#include<iostream.h>#include<fstream.h>#include <iomanip.h> //使用setw函数#include<math.h>#include "stdlib.h"#include<stdio.h>#include "math.h"#include <iostream>#include <stdio.h>#include <tchar.h>#include <fstream>#include <cmath>#include <process.h>#include <math.h>#include <malloc.h>using namespace std;void main(){FILE *fp;fp=fopen("result.txt","w");ifstream fin;fin.open("node information.txt");int node_num;fin>>node_num;int n = node_num*(node_num-1);int Jocabi_rank=2*(node_num-1);double * * Jocabi=new double * [Jocabi_rank];int *II=new int[n];int *JJ=new int[n];double *WR=new double[n];double *WX=new double[n];double *WXC=new double[n];double *WNT=new double[n];for(int i=0;i<n;i++){II[i]=0;JJ[i]=0;WR[i]=0;WX[i]=0;WXC[i]=0;WNT[i]=0;}double *GF=new double[n];for(int GF_i=0;GF_i<=n;GF_i++)GF[GF_i]=0;double *BF=new double[n];for(int BF_i=0;BF_i<=n;BF_i++)BF[BF_i]=0;double * * GF1=new double * [node_num+1];for (int GF1_i=0;GF1_i<node_num;GF1_i++){GF1[GF1_i]=new double[node_num+1];for (int GF1_j=0;GF1_j<node_num;GF1_j++){GF1[GF1_i][GF1_j]=0;}} double * * BF1=new double * [node_num+1];for (int BF1_i=0;BF1_i<=node_num;BF1_i++) {BF1[BF1_i]=new double[node_num+1];for (int BF1_j=0;BF1_j<node_num;BF1_j++) {BF1[BF1_i][BF1_j]=0;}} double *GD=new double[node_num];for(int GD_i=0;GD_i<=node_num;GD_i++){GD[GD_i]=0;}double *BD=new double[node_num];for(int BD_i=0;BD_i<=node_num;BD_i++){BD[BD_i]=0;}double * BC=new double[node_num];for(int BC_i=0;BC_i<=node_num;BC_i++){BC[BC_i]=0;}for(i=0;i<n;i++)fin>>II[i]>>JJ[i]>>WR[i]>>WX[i]>>WXC[i]>>WNT[i];//节点导纳矩阵int p,q;double K;for(i=0;i<n;i++){p=0;q=0;K=0.0;if(WNT[i]==0){GF[i]=WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]);BF[i]=(-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]);}if(WNT[i]>0){K=WNT[i];p=II[i];GF[i]=(WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]))/K;BF[i]=((-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]))/K;GD[p-1]=GD[p-1]+((K-1)/K)*(WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]));BD[p-1]=BD[p-1]+((K-1)/K)*(-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]);}if(WNT[i]<0){K=fabs(WNT[i]);q=JJ[i];GF[i]=(WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]))/K;BF[i]=((-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]))/K;GD[q-1]=GD[q-1]+((1-K)/(K*K))*WR[i]/(WR[i]*WR[i]+WX[i]*WX[i]);BD[q-1]=BD[q-1]+((1-K)/(K*K))*(-WX[i])/(WR[i]*WR[i]+WX[i]*WX[i]);}}for(int I=1;I<=node_num;I++){BC[I-1]=0.0;for(int i=0;i<n;i++){int J=0;if(II[i]==I){if(WXC[i]<1000)BC[I-1]=BC[I-1]+WXC[i];J=JJ[i];GF1[I-1][J-1]=-GF[i];BF1[I-1][J-1]=-BF[i];GF1[I-1][I-1]=GF1[I-1][I-1]+GF[i];BF1[I-1][I-1]=BF1[I-1][I-1]+BF[i];}}GF1[I-1][I-1]=GF1[I-1][I-1]+GD[I-1];BF1[I-1][I-1]=BF1[I-1][I-1]+BD[I-1]+BC[I-1];} //输出节点导纳矩阵fprintf(fp,"节点导纳矩阵是： \n");for(int m=0;m<node_num;m++){for(int n=0;n<node_num;n++){fprintf(fp,"%9.6f",GF1[m][n]);fprintf(fp,"+j");fprintf(fp,"%9.6f",BF1[m][n]);fprintf(fp," ");}fprintf(fp,"\n");}fprintf(fp,"\n");fprintf(fp,"\n");fprintf(fp,"\n");//$$$$$$$$$$$$$$$$$输出节点导纳矩阵完毕$$$$$$$$$$$$$$$$$$$$$//$$$$$$$$$$$$$$$$$$$求节点初始差值$$$$$$$$$$$$$$$$$$$$$$$$$ double *wv=new double [node_num];double *we=new double [node_num];double *wf=new double [node_num];double *wp_g=new double [node_num];double *wq_g=new double [node_num];double *wp_l=new double [node_num];double *wq_l=new double [node_num];double *wn=new double [node_num];for(int p_i=0;p_i<node_num;p_i++){wv[p_i]=0;we[p_i]=0;wf[p_i]=0;wp_g[p_i]=0;wq_g[p_i]=0;wp_l[p_i]=0;wq_l[p_i]=0;wn[p_i]=0;}double * dp=new double [node_num];for (int dp_i=0;dp_i<node_num;dp_i++)dp[dp_i]=0;double * dq=new double [node_num];for (int dq_i=0;dq_i<node_num;dq_i++)dq[dq_i]=0;double * dvv=new double [node_num];for (int dvv_i=0;dvv_i<node_num;dvv_i++)dvv[dvv_i]=0;double * sm=new double [node_num];double * sn=new double [node_num];int n_m=0 ;ifstream ffin;ffin.open("node power.txt");for(int pi=0;pi<node_num;pi++){ffin>>wp_g[pi]>>wq_g[pi]>>wp_l[pi]>>wq_l[pi]>>we[pi]>>wf[pi]>> wn[pi];}for(int wv_i=0;wv_i<node_num;wv_i++){wv[wv_i]=we[wv_i];}double * w_p=new double [node_num];double * w_q=new double [node_num];for (int w_i=0;w_i<node_num;w_i++){w_p[w_i]=0;w_q[w_i]=0;}//开始迭代，次数为iterative_frequencyfor (intiterative_frequency=0;iterative_frequency<12;iterative_frequency++) {fprintf(fp,"第[");fprintf(fp,"%2d",iterative_frequency);fprintf(fp,"]次迭代结果如下：");fprintf(fp,"\n");n_m=0;for(int Jocabi_i=0;Jocabi_i<Jocabi_rank;Jocabi_i++){Jocabi[Jocabi_i]=new double[Jocabi_rank];for(int Jocabi_j=0;Jocabi_j<Jocabi_rank;Jocabi_j++)Jocabi[Jocabi_i][Jocabi_j]=0;}for(i=0;i<node_num;i++){sm[i]=0;sn[i]=0;}for(i=0;i<node_num;i++){if(wn[i]==1){for(int j=0;j<node_num;j++){sm[i]=sm[i]+(GF1[i][j]*we[j]-BF1[i][j]*wf[j]);sn[i]=sn[i]+(GF1[i][j]*wf[j]+BF1[i][j]*we[j]);}dp[i]=wp_l[i]-we[i]*sm[i]-wf[i]*sn[i];dq[i]=wq_l[i]-wf[i]*sm[i]+we[i]*sn[i];n_m++;}if(wn[i]==2){for(int j=0;j<node_num;j++){sm[i]=sm[i]+(GF1[i][j]*we[j]-BF1[i][j]*wf[j]);sn[i]=sn[i]+(GF1[i][j]*wf[j]+BF1[i][j]*we[j]);}dp[i]=wp_g[i]-we[i]*sm[i]-wf[i]*sn[i];dvv[i]=wv[i]*wv[i]-(we[i]*we[i]+wf[i]*wf[i]);w_p[i]=we[i]*sm[i]+wf[i]*sn[i];w_q[i]=wf[i]*sm[i]-we[i]*sn[i];fprintf(fp,"第[" );fprintf(fp,"%2d",i+1 );fprintf(fp, "]个节点是PV节点，且其节点功率为：");fprintf(fp,"%9.6f",w_p[i]);fprintf(fp,"+j" );fprintf(fp,"%9.6f",w_q[i]);fprintf(fp,"\n" );} if(wn[i]==0){for(int j=0;j<node_num;j++){sm[i]=sm[i]+(we[j]*GF1[i][j]-wf[j]*BF1[i][j]);sn[i]=sn[i]+(wf[j]*GF1[i][j]+we[j]*BF1[i][j]);}w_p[i]=we[i]*sm[i]+wf[i]*sn[i];;w_q[i]=wf[i]*sm[i]-we[i]*sn[i];fprintf(fp,"第[" );fprintf(fp,"%2d",i+1 );fprintf(fp, "]个节点是平衡节点，其节点功率为：");fprintf(fp,"%9.6f",w_p[i]);fprintf(fp,"+j" );fprintf(fp,"%9.6f",w_q[i]);fprintf(fp,"\n" );} //$$$$$$$$$$$$$$$$$$$$求Jocabi矩阵$$$$$$$$$$$$$$$$$$$$$for(int j=0;j<node_num;j++){if(j==i&&wn[i]==1){Jocabi[2*i][2*j]=-sm[i]-GF1[i][i]*we[i]-BF1[i][i]*wf[i]; Jocabi[2*i][2*j+1]=-sn[i]+BF1[i][i]*we[i]-GF1[i][i]*wf[i]; Jocabi[2*i+1][2*j]=sn[i]+BF1[i][i]*we[i]-GF1[i][i]*wf[i]; Jocabi[2*i+1][2*j+1]=-sm[i]+GF1[i][i]*we[i]+BF1[i][i]*wf[i];} if(j==i&&wn[i]==2){Jocabi[2*i][2*j]=-sm[i]-GF1[i][i]*we[i]-BF1[i][i]*wf[i]; Jocabi[2*i][2*j+1]=-sn[i]+BF1[i][i]*we[i]-GF1[i][i]*wf[i]; Jocabi[2*i+1][2*j]=(-2)*we[i];Jocabi[2*i+1][2*j+1]=(-2)*wf[i];} if(j!=i&&wn[i]==1){Jocabi[2*i][2*j]=(-1)*(GF1[i][j]*we[i]+BF1[i][j]*wf[i]);Jocabi[2*i][2*j+1]=BF1[i][j]*we[j]-GF1[i][j]*wf[i];Jocabi[2*i+1][2*j]=BF1[i][j]*we[j]-GF1[i][j]*wf[i];Jocabi[2*i+1][2*j+1]=GF1[i][j]*we[i]+BF1[i][j]*wf[i];} if(j!=i&&wn[i]==2){Jocabi[2*i][2*j]=(-1)*(GF1[i][j]*we[i]+BF1[i][j]*wf[i]);Jocabi[2*i][2*j+1]=BF1[i][j]*we[j]-GF1[i][j]*wf[i];Jocabi[2*i+1][2*j]=0;Jocabi[2*i+1][2*j+1]=0;}}} //$$$$$$$$$$$$$$$$$$$高斯法解方程$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ double * Dalta_w=new double[Jocabi_rank];for(i=0;i<node_num;i++){if(wn[i]==1){Dalta_w[2*i]=dp[i];Dalta_w[2*i+1]=dq[i];}if(wn[i]==2){Dalta_w[2*i]=dp[i];Dalta_w[2*i+1]=dvv[i];}}double * M=new double[Jocabi_rank];double * X=new double[Jocabi_rank];for(i=0;i<Jocabi_rank;i++){X[i]=0;M[i]=0;}double t;int * js=new int [Jocabi_rank];for(i=0;i<Jocabi_rank;i++)js[i]=0;//交换Jocabi矩阵的行，使每列的对角线元素为所有行中最大的元素double * temp_swap_vector=new double[Jocabi_rank+1];for(inttemp_swap_vector_i=0;temp_swap_vector_i<Jocabi_rank;temp_swap_vector_ i++)temp_swap_vector[temp_swap_vector_i]=0;double temp_variable_for_compare=0;int temp_line=0;for(int diagonal_i=0;diagonal_i<Jocabi_rank;diagonal_i++){temp_variable_for_compare=Jocabi[diagonal_i][diagonal_i];for(intJocabi_line=diagonal_i;Jocabi_line<Jocabi_rank;Jocabi_line++) {if(temp_variable_for_compare<fabs(Jocabi[Jocabi_line][diagonal _i])){temp_variable_for_compare=fabs(Jocabi[Jocabi_line][diagonal_i] );temp_line=Jocabi_line;for(int Jocabi_j_j=0;Jocabi_j_j<Jocabi_rank;Jocabi_j_j++) {temp_swap_vector[Jocabi_j_j]=Jocabi[diagonal_i][Jocabi_j_j];Jocabi[diagonal_i][Jocabi_j_j]=Jocabi[temp_line][Jocabi_j_j];Jocabi[temp_line][Jocabi_j_j]=temp_swap_vector[Jocabi_j_j];}t=Dalta_w[diagonal_i];Dalta_w[diagonal_i]=Dalta_w[temp_line];Dalta_w[temp_line]=t;}}}//$$$$$$$$$$$$$$交换Jocabi矩阵的行完毕$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ double * * temp=new double * [Jocabi_rank];for(i=0;i<Jocabi_rank;i++){temp[i]=new double[Jocabi_rank];for(int j=0;j<Jocabi_rank;j++)temp[i][j]=0.0;}for(int k=0;k<Jocabi_rank;k++){if(Jocabi[k][k]==0) cout<<"fail"<<endl;else{temp[k][k]=Jocabi[k][k];for(int j=0;j<Jocabi_rank;j++)Jocabi[k][j]=Jocabi[k][j]/temp[k][k];Dalta_w[k]=Dalta_w[k]/temp[k][k];for(int ii=k+1;ii<Jocabi_rank;ii++){M[ii]=0.0;M[ii]=Jocabi[ii][k];for(int jj=k;jj<Jocabi_rank;jj++)Jocabi[ii][jj]=Jocabi[ii][jj]-M[ii]*Jocabi[k][jj];Dalta_w[ii]=Dalta_w[ii]-M[ii]*Dalta_w[k];}}} X[Jocabi_rank-1]=Dalta_w[Jocabi_rank-1]/(-Jocabi[Jocabi_rank-1][J ocabi_rank-1]);for(i=Jocabi_rank-2;i>=0;i--){t=0;for(k=i+1;k<Jocabi_rank;k++){t=t+(-Jocabi[i][k])*X[k];}X[i]=(Dalta_w[i]-t)/(-Jocabi[i][i]);}for(i=0;i<node_num;i++){if(wn[i]==1){we[i]=we[i]+X[2*i];wf[i]=wf[i]+X[2*i+1];fprintf(fp,"第[");fprintf(fp,"%2d",i+1);fprintf(fp,"]个节点是PQ节点其节点电压实部修正量是：");fprintf(fp,"%9.6f",X[2*i]);fprintf(fp," 其节点电压虚部修正量是：");fprintf(fp,"%9.6f",X[2*i+1]);fprintf(fp,"\n");} if(wn[i]==2){we[i]=we[i]+X[2*i];wf[i]=wf[i]+X[2*i+1];fprintf(fp,"第[");fprintf(fp,"%2d",i+1);fprintf(fp,"]个节点是PV节点其节点电压实部修正量是：");fprintf(fp,"%9.6f",X[2*i]);fprintf(fp," 其节点电压虚部修正量是：");fprintf(fp,"%9.6f",X[2*i+1]);fprintf(fp,"\n");}} fprintf(fp,"\n");fprintf(fp,"\n");//$$$$$$$$$$$$$$$$判断是否收敛并求最后的差值$$$$$$$$$$$$$$$$$$$$$$ double ch=0;double * * s_p=new double * [node_num+1];for(int s_p_i=1;s_p_i<=node_num+1;s_p_i++){s_p[s_p_i]=new double[node_num+1];for(int s_p_j=1;s_p_j<=node_num+1;s_p_j++)s_p[s_p_i][s_p_j]=0;}double * * s_q=new double * [node_num+1];for(int s_q_i=1;s_q_i<=node_num+1;s_q_i++){s_q[s_q_i]=new double[node_num+1];for(int s_q_j=1;s_q_j<=node_num+1;s_q_j++)s_q[s_q_i][s_q_j]=0;}for(int ci=1;ci<=node_num-1;ci++){if((fabs(Dalta_w[ci]))>(fabs(ch))){ch=fabs(Dalta_w[ci]);}}fprintf(fp,"判断是否迭代的最大值是:");fprintf(fp,"%9.6f",ch);fprintf(fp,"\n");fprintf(fp,"\n");if(ch<1e-005){double *UU=new double [node_num];for(int k=0;k<node_num;k++)UU[k]=0.0;for(k=0;k<node_num;k++)UU[k]=we[k]*we[k]+wf[k]*wf[k];double *G_ground=new double [n];double *B_ground=new double [n];for(int i=0;i<n;i++){G_ground[i]=0;B_ground[i]=0;}cout<<"满足迭代要求，计算完成，迭代次数是"<<iterative_frequency<<endl;fprintf(fp,"||||||||||||||||||||||||||||||||||||||||||||||||||||| ||||||||\n");fprintf(fp,"\n");fprintf(fp,"满足迭代要求，计算完成，迭代次数是");fprintf(fp,"%2d",iterative_frequency);fprintf(fp,"\n");for(i=0;i<node_num;i++){fprintf(fp,"节点");fprintf(fp,"[");fprintf(fp,"%2d",i+1);fprintf(fp,"]的电压是：");fprintf(fp,"%9.6f",we[i]);if(wf[i]>0){fprintf(fp," + j");fprintf(fp,"%9.6f",wf[i]);}if(wf[i]<0){fprintf(fp," - j");fprintf(fp,"%9.6f",fabs(wf[i]));}fprintf(fp," 幅值是：");fprintf(fp,"%9.6f",sqrt(we[i]*we[i]+wf[i]*wf[i]));fprintf(fp," 相角是：");fprintf(fp,"%9.6f",atan(wf[i]/we[i])*180/3.1415926);fprintf(fp,"度");fprintf(fp,"\n");}fprintf(fp,"线路功率是:\n");for(int j=0;j<n;j++){double K=0.0;int p=0;int q=0;p=II[j];q=JJ[j];if(p==0&&q==0)break;GF[j]=WR[j]/(WR[j]*WR[j]+WX[j]*WX[j]);BF[j]=(-WX[j])/(WR[j]*WR[j]+WX[j]*WX[j]);if(WNT[j]<0){K=fabs(WNT[j]);G_ground[j]=((1-K)/(K*K))*GF[j];B_ground[j]=((1-K)/(K*K))*BF[j];}if(WNT[j]>0){K=WNT[j];G_ground[j]=((K-1)/K)*GF[j];B_ground[j]=((K-1)/K)*BF[j];}if(WXC[j]<1000)B_ground[j]=B_ground[j]+WXC[j];s_p[p][q]=UU[p-1]*G_ground[j]+UU[p-1]*GF1[p-1][q-1]-(we[p-1]*we[q -1]+wf[p-1]*wf[q-1])*GF1[p-1][q-1]+(we[p-1]*wf[q-1]-we[q-1]*wf[p-1])*BF1[p-1][q-1];s_q[p][q]=-UU[p-1]*B_ground[j]-UU[p-1]*BF1[p-1][q-1]+(we[p-1]*we[ q-1]+wf[q-1]*wf[q-1])*BF1[p-1][q-1]+(we[p-1]*wf[q-1]-we[q-1]*wf[p-1])*GF1[p-1][q-1];fprintf(fp,"S[");fprintf(fp,"%2d",p);fprintf(fp,"][");fprintf(fp,"%2d",q);fprintf(fp,"]=");fprintf(fp,"%9.6f",s_p[p][q]);if(s_q[p][q]>0){fprintf(fp," + j");fprintf(fp,"%9.6f",s_q[p][q]);}if(s_q[p][q]<0){fprintf(fp," - j");fprintf(fp,"%9.6f",fabs(s_q[p][q]));}fprintf(fp,"\n");}break;}}//迭代结束八、IEEE30节点的仿真结果1、节点电压编号电压实部电压虚部电压幅值电压相角节点1 1.05 0i 1.05 0度节点2 1.04445 -0.0337697i 1.045 -1.85186度节点3 1.01843 -0.067472i 1.02066 -3.79037度节点4 1.01033 -0.0798977i 1.01348 -4.52159度节点5 1.0035 -0.114413i 1.01 -6.50445度节点6 1.00662 -0.0943613i 1.01104 -5.35528度节点7 0.996459 -0.111039i 1.00263 -6.35844度节点8 1.00511 -0.0993071i 1.01 -5.64266度节点9 1.03035 -0.122151i 1.03756 -6.76104度节点10 1.02229 -0.156054i 1.03413 -8.67927度节点11 1.04666 -0.0837098i 1.05 -4.57268度节点12 1.03608 -0.141734i 1.04573 -7.78959度节点13 1.0436 -0.115738i 1.05 -6.32836度节点14 1.01898 -0.156444i 1.03092 -8.72848度节点15 1.01403 -0.157951i 1.02626 -8.85357度节点16 1.02183 -0.151553i 1.033 -8.43636度节点17 1.01658 -0.157803i 1.02876 -8.82354度节点18 1.00279 -0.167848i 1.01674 -9.50214度节点19 0.999803 -0.170757i 1.01428 -9.69206度节点20 1.00449 -0.16805i 1.01845 -9.49755度节点21 1.00872 -0.162656i 1.02175 -9.16004度节点22 1.00931 -0.162655i 1.02234 -9.15476度节点23 1.00254 -0.16504i 1.01604 -9.34827度节点24 0.996609 -0.169594i 1.01094 -9.65755度节点25 0.995082 -0.169371i 1.00939 -9.65966度节点26 0.976247 -0.173651i 0.991571 -10.0861度节点27 1.00354 -0.166036i 1.01718 -9.3945度节点28 1.00168 -0.102188i 1.00688 -5.82495度节点29 0.980068 -0.184114i 0.997212 -10.6395度节点30 0.96576 -0.197073i 0.985662 -11.5334度2、线路传输功率首端—末端首末功率末首功率2-1 -57.9109+(8.36395)i 58.5228+(-12.325)i3-1 -39.536+(-7.37948)i 40.2261+(5.83313)i4-2 -30.8909+(-9.48018)i 31.4524+(7.29204)i4-3 -36.9557+(-6.53076)i 37.136+(6.17948)i5-2 -44.4784+(-7.6178)i 45.4077+(7.10769)i6-2 -38.4904+(-7.58429)i 39.3508+(6.2414)i6-4 -34.7811+(4.04816)i 34.9241+(-4.02213)i 7-5 -0.260435+(-7.29197)i 0.278436+(5.2715)i7-6 -22.5396+(-3.60803)i 22.6765+(2.30528)i8-6 -11.9244+(0.486521)i 11.9413+(-1.34664)i 6-9 12.3726+(-12.7417)i -12.3726+(13.3835)i6-10 10.9034+(-3.88366)i -10.9034+(4.61235)i 11-9 20+(6.66081)i -20+(-5.82246)i10-9 -32.6509+(-2.67651)i 32.6509+(3.78043)i4-12 23.6005+(-12.0947)i -23.6005+(13.8474)i 13-12 20+(3.45425)i -20+(-2.93116)i14-12 -7.91606+(-2.10553)i 7.99377+(2.26709)i15-12 -18.28+(-5.90486)i 18.5119+(6.36176)i16-12 -7.55825+(-2.98863)i 7.61675+(3.11163)i15-14 -1.7094+(-0.499524)i 1.71605+(0.505537)i 17-16 -4.04947+(-1.1564)i 4.05825+(1.18863)i18-15 -6.09617+(-1.40631)i 6.13679+(1.48904)i19-18 -2.89082+(-0.495509)i 2.89617+(0.506313)i 20-19 6.6264+(2.93893)i -6.60918+(-2.90448)i 20-10 -8.82641+(-3.63893)i 8.90866+(3.82259)i17-10 -4.95054+(-4.64359)i 4.96464+(4.68038)i21-10 -16.1704+(-9.32613)i 16.2865+(9.57613)i22-10 -7.88783+(-4.19568)i 7.94335+(4.31016)i22-21 1.33022+(1.87505)i -1.32963+(-1.87386)i 23-15 -5.6167+(-2.34288)i 5.65258+(2.41535)i24-22 -6.50437+(-2.23775)i 6.55761+(2.32062)i24-23 -2.40853+(-0.726163)i 2.4167+(0.742882)i25-24 -0.212582+(-0.351348)i 0.212894+(0.351893)i 26-25 -3.5+(-2.29999)i 3.54539+(2.36778)i27-25 3.34908+(2.04752)i -3.3328+(-2.01643)i 28-27 16.1025+(-2.11682)i -16.1025+(3.14712)i 29-27 -6.10422+(-1.50686)i 6.1916+(1.67196)i30-27 -6.92977+(-1.35732)i 7.09411+(1.66665)i30-29 -3.67023+(-0.542645)i 3.70422+(0.60687)i28-8 -1.92156+(-3.12492)i 1.92445+(-1.21858)i28-6 -14.7132+(-3.43133)i 14.7506+(2.24037)i 3、发电机注入的功率编号功率实部功率虚部节点1 98.7489 -6.4919i 节点2 80 41.7051i 节点5 50 16.6537i 节点8 20 29.2679i 节点11 20 6.66081i 节点13 20 3.45425i。

电力系统潮流计算程序

电力系统潮流计算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语言编写的PQ分解法解潮流程序

这里要提一下，本人使用的数据出自鞠平主编《电力工程》，PQ分解法可能由于求逆程序的精确度问题而比matlab迭代次数稍多，可将float改为dobule float一试，若有不当之处，还望包涵（另外，本人使用c语言编译器为tc3.0）//PQ分解法解潮流程序#include<stdio.h>#include<math.h>#define N 4 //节点数#define n_PQ 2 //PQ节点数#define n_PV 1 //PV节点数#define n_br 5 //串联支路数#define PI 3.1415void main(){void disp_matrix(float *disp_p,int disp_m,int disp_n); //矩阵显示函数float Us[2*N]={1.0,0,1.0,0,1.05,0,1.05,0}; //电压初值float Ps[N]={0,-0.5,0.2}; //有功初值float Qs[N]={0,-0.3}; //无功初值floatG[N][N],B[N][N],B1[N-1][N-1],B2[n_PQ][n_PQ],invB1[N-1][N-1],invB2[n_PQ][n_PQ];struct //阻抗参数{int nl; //左节点int nr; //右节点float R; //串联电阻值float X; //串联电抗值float Bl; //左节点并联电纳float Br; //右节点并联电纳}ydata[n_br]={{1,2,0,0.1880,-0.6815,0.6040},{1,3,0.1302,0.2479,0.0129,0.0129},{1,4,0.1736,0.3306,0.0172,0.0172},{3,4,0.2603,0.4959,0.0259,0.0259},{2,2,0,0.05,0,0}};float Z2; //Z^2=R^2+X^2 各串联阻抗值的平方float U_amp[N],U_ang[N]; //e,f存储电压的x轴分量和y轴分量,dfe存储电压修正值float dU_amp[N],dU_ang[N]; //电压幅值和相角修正值float mid1[N],mid2[N],dP[N-1],dQ[n_PQ]; //mid1、mid2存储计算雅克比行列式对角线元素的中间值,dS存储PQU的不平衡量float dPQU=1.0; //PQU不平衡量最大值int kk=0; //迭代次数int i,j,k;float t;//形成导纳矩阵-------------------------------------------------------------------------------------------------- for(i=0;i<N;i++)for(j=0;j<N;j++){G[i][j]=0;B[i][j]=0;}for(i=0;i<n_br;i++){if(ydata[i].nl!=ydata[i].nr){Z2=(ydata[i].R)*(ydata[i].R)+(ydata[i].X)*(ydata[i].X);//串联阻抗等效导纳值//非对角元素G[ydata[i].nl-1][ydata[i].nr-1]=(-ydata[i].R)/Z2;B[ydata[i].nl-1][ydata[i].nr-1]=ydata[i].X/Z2;G[ydata[i].nr-1][ydata[i].nl-1]=(-ydata[i].R)/Z2;B[ydata[i].nr-1][ydata[i].nl-1]=ydata[i].X/Z2;//对角元素G[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].R/Z2;G[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].R/Z2;B[ydata[i].nl-1][ydata[i].nl-1]+=(-ydata[i].X/Z2);B[ydata[i].nr-1][ydata[i].nr-1]+=(-ydata[i].X/Z2);//并联导纳等效导纳值B[ydata[i].nl-1][ydata[i].nl-1]+=ydata[i].Bl;B[ydata[i].nr-1][ydata[i].nr-1]+=ydata[i].Br;}else{G[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].R;B[ydata[i].nl-1][ydata[i].nr-1]+=ydata[i].X;}}printf("G=\n");disp_matrix(*G,N,N);printf("B=\n");disp_matrix(*B,N,N);//分离U,deltafor(i=0;i<N;i++){U_amp[i]=Us[2*i];U_ang[i]=Us[2*i+1];}=================================================================== ==================================while(dPQU>0.00001){for(i=0;i<N-1;i++){mid1[i]=0;mid2[i]=0;for(j=0;j<N;j++){mid1[i]=mid1[i]+G[i][j]*U_amp[j]*cos(U_ang[j])-B[i][j]*U_amp[j]*sin(U_ang[j]);mid2[i]=mid2[i]+G[i][j]*U_amp[j]*sin(U_ang[j])+B[i][j]*U_amp[j]*cos(U_ang[j]);}dP[i]=Ps[i]-(U_amp[i]*cos(U_ang[i])*mid1[i]+U_amp[i]*sin(U_ang[i])*mid2[i]);if(i<n_PQ)dQ[i]=Qs[i]-(U_amp[i]*sin(U_ang[i])*mid1[i]-U_amp[i]*cos(U_ang[i])*mid2[i]);}dPQU=0;for(i=0;i<N-1;i++){if(dP[i]<0&&dPQU<-dP[i])dPQU=-dP[i];else if(dP[i]>0&&dPQU<dP[i])dPQU=dP[i];}for(i=0;i<n_PQ;i++){if(dQ[i]<0&&dPQU<-dQ[i])dPQU=-dQ[i];else if(dQ[i]>0&&dPQU<dQ[i])dPQU=dQ[i];}if(dPQU>0.00001){kk++;//形成B',B''矩阵for(i=0;i<N-1;i++)for(j=0;j<N-1;j++){B1[i][j]=B[i][j];if(i<n_PQ&&j<n_PQ)B2[i][j]=B[i][j];//求P,delta修正值for(i=0;i<N-1;i++){dP[i]=dP[i]/U_amp[i];if(i<n_PQ)dQ[i]=dQ[i]/U_amp[i];}//B'求逆----------------------------------------------------------------------------------- for(i=0;i<N-1;i++)for(j=0;j<N-1;j++){if(i!=j)invB1[i][j]=0;elseinvB1[i][j]=1;}for(i=0;i<N-1;i++){for(j=0;j<N-1;j++){if(i!=j){t=B1[j][i]/B1[i][i];for(k=0;k<N-1;k++){B1[j][k]-=B1[i][k]*t;invB1[j][k]-=invB1[i][k]*t;}}}}for(i=0;i<N-1;i++)if(B1[i][i]!=1){t=B1[i][i];for(j=0;j<N-1;j++)invB1[i][j]=invB1[i][j]/t;}for(i=0;i<N-1;i++)dU_ang[i]=0;for(i=0;i<N-1;i++){for(j=0;j<N-1;j++)dU_ang[i]+=invB1[i][j]*dP[j]; dU_ang[i]=dU_ang[i]/U_amp[i]; }//Q,U_amp修正值//B''求逆for(i=0;i<n_PQ;i++)for(j=0;j<n_PQ;j++){if(i!=j)invB2[i][j]=0;elseinvB2[i][j]=1;}for(i=0;i<n_PQ;i++){for(j=0;j<n_PQ;j++){if(i!=j){t=B2[j][i]/B2[i][i];for(k=0;k<n_PQ;k++){B2[j][k]-=B2[i][k]*t;invB2[j][k]-=invB2[i][k]*t;}}}}for(i=0;i<n_PQ;i++)if(B2[i][i]!=1){t=B2[i][i];for(j=0;j<n_PQ;j++)invB2[i][j]=invB2[i][j]/t;}for(i=0;i<n_PQ;i++)dQ[i]=dQ[i]/U_amp[i];for(i=0;i<n_PQ;i++)dU_amp[i]=0;for(i=0;i<n_PQ;i++){for(j=0;j<n_PQ;j++)dU_amp[i]+=invB2[i][j]*dQ[j]; }for(i=0;i<N;i++){U_amp[i]-=dU_amp[i];U_ang[i]-=dU_ang[i];}}elsebreak;}//循环结束---------------------------------------------------------------------------------------------------------------//求平衡节点功率---------------------------------------------------------------------------------------------------mid1[N-1]=0;mid2[N-1]=0;for(j=0;j<N;j++){mid1[N-1]=mid1[N-1]+G[N-1][j]*U_amp[j]*cos(U_ang[j])-B[N-1][j]*U_amp[j]*sin(U_ang[j] );mid2[N-1]=mid2[N-1]+G[N-1][j]*U_amp[j]*sin(U_ang[j])+B[N-1][j]*U_amp[j]*cos(U_ang[j ]);}Ps[N-1]=U_amp[N-1]*cos(U_ang[N-1])*mid1[N-1]+U_amp[N-1]*sin(U_ang[N-1])*mid2[N-1];Qs[N-1]=U_amp[N-1]*sin(U_ang[N-1])*mid1[N-1]-U_amp[N-1]*cos(U_ang[N-1])*mid2[N-1 ];for(i=n_PQ;i<N-1;i++)Qs[i]=U_amp[i]*sin(U_ang[i])*mid1[i]-U_amp[i]*cos(U_ang[i])*mid2[i];//---------------------------------------------------------------------------------------------------------------------------// 显示输出结果printf("kk=%d\n",kk);printf("P=");for(i=0;i<N;i++)printf("%9.4f",Ps[i]);printf("\nQ=");for(i=0;i<N;i++)printf("%9.4f",Qs[i]);printf("\nU_amp=");for(i=0;i<N;i++)printf("%9.4f",U_amp[i]);printf("\nU_ang(deg)=");for(i=0;i<N;i++)printf("%9.4f",U_ang[i]*180/PI);printf("\n\n");}//================================================================== ======================================void disp_matrix(float *disp_p,int disp_m,int disp_n){int i,j;for(i=0;i<disp_m;i++){for(j=0;j<disp_n;j++)printf("%9.4f",*(disp_p+i*disp_m+j));printf("\n");}printf("\n");}。

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

程序为计算书3-4的过程程序可以解决开式单直网络和树状网络的计算。

树状网络计算时要自己先设定好支路的起始节点和终止节点标号以及计算顺序源代码：#include <iostream>#include <fstream>#include<iomanip>#include<math.h>using namespace std;struct node{//节点类int i;//节点编号double U,P,Q,delta;//额定电压计算负荷电压相角};struct line{//线路类连接父节点子节点node f_node,s_node;//父节点子节点double R,X,B;//线路参数R X B/2double 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;double q=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){//由前往后推电压double p=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;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=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(int i=0;i<num_n;i++){fin>>u[i];};double *p;//节点有功功率p=new double[num_n];for(int i=0;i<num_n;i++){fin>>p[i];};double *q;//节点无功功率q=new double[num_n];for(int i=0;i<num_n;i++){fin>>q[i];};for(int i=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(int i=0;i<num_l;i++){fin>>r[i];};double *x;//线路电抗x=new double[num_l];for(int i=0;i<num_l;i++){fin>>x[i];};double *b;//线路电纳b=new double[num_l];for(int i=0;i<num_l;i++){fin>>b[i];};for(int i=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(int tc=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;double to_P(0),to_dP(0);for(int i=1;i<num_n;i++){to_P+=nod[i].P;};for(int tc=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(int i=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(int i=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文件3 2113 110 1000 0.17 200 1.7 158.5 1.2220.5 20.20.000282 00 11 2输出为：第1次迭代支路信息:支路1-2:始端功率:1.02165+j0.693296末端功率:1.00434+j0.658675功率损耗:0.0173106+j0.0346211 电压损耗0.274063支路2-3:始端功率:0.5034+j0.3068末端功率:0.5+j0.3功率损耗:0.0034+j0.0068电压损耗0.109支路2-4:始端功率:0.200938+j0.151875 末端功率:0.2+j0.15功率损耗:0.0009375+j0.001875 电压损耗0.0739643节点信息:节点1电压:10.5相角:0节点2电压:10.2259相角:-0.864932节点3电压:10.1169相角:-1.25281节点4电压:10.152相角:-1.072全网信息:总电源有功:1.02165总负荷有功:1总有功损耗:0.0216481网损率:0.0107081最低电压:10.1169最低电压点:3第2次迭代支路信息:支路1-2:始端功率:1.02078+j0.69156末端功率:1.00423+j0.658463功率损耗:0.0165484+j0.0330969电压损耗0.273567支路2-3:始端功率:0.503322+j0.306644末端功率:0.5+j0.3功率损耗:0.00332186+j0.00664371 电压损耗0.108957支路2-4:始端功率:0.20091+j0.151819末端功率:0.2+j0.15功率损耗:0.000909642+j0.00181928 电压损耗0.0739403节点信息:节点1电压:10.5相角:0节点2电压:10.2264相角:-0.86489节点3电压:10.1175相角:-1.25273节点4电压:10.1525相角:-1.07194全网信息:总电源有功:1.02078总负荷有功:1.70434总有功损耗:0.0207799网损率:0.00762533最低电压:10.1175最低电压点:3第3次迭代支路信息:支路1-2:始端功率:1.02078+j0.691556末端功率:1.00423+j0.658462功率损耗:0.0165468+j0.0330936电压损耗0.273566支路2-3:始端功率:0.503322+j0.306643末端功率:0.5+j0.3功率损耗:0.0033215+j0.00664301 电压损耗0.108957支路2-4:始端功率:0.20091+j0.151819末端功率:0.2+j0.15功率损耗:0.000909549+j0.0018191 电压损耗0.0739402节点信息:节点1电压:10.5相角:0节点2电压:10.2264相角:-0.86489节点3电压:10.1175相角:-1.25273节点4电压:10.1525相角:-1.07194全网信息:总电源有功:1.02078总负荷有功:1.70423总有功损耗:0.0207778网损率:0.00762487最低电压:10.1175最低电压点:3例3-2的data：4 310.5 10 10 100 0.3 0.5 0.20 0.2 0.3 0.151.2 1 1.52.4 2 30 0 00 11 21 3输出：第1次迭代支路信息:支路1-2:始端功率:1.02165+j0.693296末端功率:1.00434+j0.658675功率损耗:0.0173106+j0.0346211电压损耗0.275227支路2-3:始端功率:0.5034+j0.3068末端功率:0.5+j0.3功率损耗:0.0034+j0.0068电压损耗0.109244支路2-4:始端功率:0.200938+j0.151875末端功率:0.2+j0.15功率损耗:0.0009375+j0.001875电压损耗0.0740389节点信息:节点1电压:10.5相角:0节点2电压:10.2248相角:-0.864932节点3电压:10.1155相角:-1.2529节点4电压:10.1507相角:-1.07205全网信息:总电源有功:1.02165总负荷有功:1总有功损耗:0.0216481网损率:0.0107081最低电压:10.1155最低电压点:3第2次迭代支路信息:支路1-2:始端功率:1.02078+j0.69157末端功率:1.00423+j0.658465功率损耗:0.0165523+j0.0331045电压损耗0.274734支路2-3:始端功率:0.503323+j0.306646末端功率:0.5+j0.3功率损耗:0.00332278+j0.00664556 电压损耗0.109201支路2-4:始端功率:0.20091+j0.15182末端功率:0.2+j0.15功率损耗:0.000909864+j0.00181973电压损耗0.0740151节点信息:节点1电压:10.5相角:0节点2电压:10.2253相角:-0.86489节点3电压:10.1161相角:-1.25282节点4电压:10.1513相角:-1.07199全网信息:总电源有功:1.02078总负荷有功:1.70434总有功损耗:0.0207849网损率:0.00762714最低电压:10.1161最低电压点:3第3次迭代支路信息:支路1-2:始端功率:1.02078+j0.691566末端功率:1.00423+j0.658464功率损耗:0.0165506+j0.0331013电压损耗0.274733支路2-3:始端功率:0.503322+j0.306645末端功率:0.5+j0.3功率损耗:0.00332243+j0.00664486 电压损耗0.109201支路2-4:始端功率:0.20091+j0.15182末端功率:0.2+j0.15功率损耗:0.000909771+j0.00181954 电压损耗0.074015节点信息:节点1电压:10.5相角:0节点2电压:10.2253相角:-0.86489节点3电压:10.1161相角:-1.25282节点4电压:10.1513相角:-1.07199全网信息:总电源有功:1.02078总负荷有功:1.70423总有功损耗:0.0207828网损率:0.00762669最低电压:10.1161最低电压点:3。