计算机潮流计算程序代码
c语言编写的牛顿拉夫逊法解潮流程序
c语言编写的牛顿拉夫逊法解潮流程序闲来无事，最近把牛拉法用c语言重写一遍，和matlab相比，c语言编写潮流程序最大的难点在于矩阵求逆，我使用的求逆方法是初等行变换法，程序段如下：#include<stdio.h>#define N 3void main(){int i,j,k;float t;float Jacob[N][N]={{1,2,2},{1,3,4},{2,3,4}};//欲进行求逆的矩阵float inv_J[N][N];//逆矩阵存储于此//初始化inv_J[N][N]for(i=0;i<N;i++)for(j=0;j<N;j++){if(i!=j)inv_J[i][j]=0;elseinv_J[i][j]=1;}//将原矩阵化简为对角阵for(i=0;i<N;i++){for(j=0;j<N;j++){if(i!=j){t=Jacob[j][i]/Jacob[i][i];for(k=0;k<N;k++){Jacob[j][k]-=Jacob[i][k]*t;inv_J[j][k]-=inv_J[i][k]*t;}}}}//原矩阵各对角元素化为1，画出逆矩阵for(i=0;i<N;i++)if(Jacob[i][i]!=1){t=Jacob[i][i];for(j=0;j<N;j++)inv_J[i][j]=inv_J[i][j]/t;}//输出逆矩阵for(i=0;i<N;i++){for(j=0;j<N;j++)printf("%9.4f",inv_J[i][j]);printf("\n");}}整个程序为：//牛拉法解潮流程序#include<stdio.h>#include<math.h>#define N 4 //节点数#define n_PQ 2 //PQ节点数#define n_PV 1 //PV节点数#define n_br 5 //串联支路数void 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}; //无功初值float G[N][N],B[N][N]; //各几点电导电纳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 e[N],f[N],dfe[2*(N-1)]; //e,f存储电压的x轴分量和y轴分量,dfe存储电压修正值float mid1[N],mid2[N],dS[2*(N-1)]; //mid1、mid2存储计算雅克比行列式对角线元素的中间值,dS存储PQU的不平衡量float Jacob[2*(N-1)][2*(N-1)],inv_J[2*(N-1)][2*(N-1)]; //雅克比行列式float dPQU=1.0; //PQU不平衡量最大值int kk=0; //迭代次数int i,j,k;float t;float Pij[n_br]; //存储线路i->j的有功float Qij[n_br]; //存储线路i->j的无功float Pji[n_br]; //存储线路j->i的有功float Qji[n_br]; //存储线路j->i的无功float dPij[n_br]; //存储线路i->j的有功损耗float dQij[n_br]; //存储线路i->j的无功损耗float AA,BB,CC,DD; //存储线路潮流计算时的中间值//形成导纳矩阵--------------------------------------------------------------------------------------------------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);//分离e,ffor(i=0;i<N;i++){e[i]=Us[2*i];f[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]*e[j]-B[i][j]*f[j];mid2[i]=mid2[i]+G[i][j]*f[j]+B[i][j]*e[j];}dS[2*i]=Ps[i]-(e[i]*mid1[i]+f[i]*mid2[i]);if(i<n_PQ)dS[2*i+1]=Qs[i]-(f[i]*mid1[i]-e[i]*mid2[i]);elsedS[2*i+1]=Us[2*i]*Us[2*i]-(e[i]*e[i]+f[i]*f[i]);}dPQU=0;for(i=0;i<2*(N-1);i++){if(dS[i]<0&&dPQU<-dS[i])dPQU=-dS[i];else if(dS[i]>0&&dPQU<dS[i])dPQU=dS[i];}if(dPQU>0.00001){kk++;//形成雅克比行列式-------------------------------------------------------------------------------------------------for(i=0;i<2*(N-1);i++)for(j=0;j<2*(N-1);j++)Jacob[i][j]=0;for(j=0;j<N-1;j++){//求H,Nfor(i=0;i<N-1;i++){if(i!=j){Jacob[2*i][2*j]=B[i][j]*e[i]-G[i][j]*f[i];Jacob[2*i][2*j+1]=-G[i][j]*e[i]-B[i][j]*f[i];}else{Jacob[2*i][2*i]=B[i][i]*e[i]-G[i][i]*f[i]-mid2[i];Jacob[2*i][2*i+1]=-G[i][j]*e[i]-B[i][j]*f[i]-mid1[i];}}//求J,Lfor(i=0;i<n_PQ;i++){if(i!=j){Jacob[2*i+1][2*j]=G[i][j]*e[i]+B[i][j]*f[i];Jacob[2*i+1][2*j+1]=B[i][j]*e[i]-G[i][j]*f[i];}else{Jacob[2*i+1][2*i]=G[i][j]*e[i]+B[i][j]*f[i]-mid1[i];Jacob[2*i+1][2*i+1]=B[i][j]*e[i]-G[i][j]*f[i]+mid2[i];}}//求R,Sfor(i=n_PQ;i<N-1;i++){if(i==j){Jacob[2*i+1][2*i]=-2*f[i];Jacob[2*i+1][2*i+1]=-2*e[i];}}}//雅克比行列式求逆-----------------------------------------------------------------------------------for(i=0;i<2*(N-1);i++)for(j=0;j<2*(N-1);j++){if(i!=j)inv_J[i][j]=0;elseinv_J[i][j]=1;}for(i=0;i<2*(N-1);i++){for(j=0;j<2*(N-1);j++){if(i!=j){t=Jacob[j][i]/Jacob[i][i];for(k=0;k<2*(N-1);k++){Jacob[j][k]-=Jacob[i][k]*t;inv_J[j][k]-=inv_J[i][k]*t;}}}}for(i=0;i<2*(N-1);i++)if(Jacob[i][i]!=1){t=Jacob[i][i];for(j=0;j<2*(N-1);j++)inv_J[i][j]=inv_J[i][j]/t;}//求电压修正值--------------------------------------------------------------------------------------------for(i=0;i<2*(N-1);i++){dfe[i]=0;for(j=0;j<2*(N-1);j++)dfe[i]-=inv_J[i][j]*dS[j];}for(i=0;i<N-1;i++){e[i]+=dfe[2*i+1];f[i]+=dfe[2*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]*e[j]-B[N-1][j]*f[j];mid2[N-1]=mid2[N-1]+G[N-1][j]*f[j]+B[N-1][j]*e[j];}Ps[N-1]=e[N-1]*mid1[N-1]+f[N-1]*mid2[N-1];Qs[N-1]=f[N-1]*mid1[N-1]-e[N-1]*mid2[N-1];for(i=n_PQ;i<N-1;i++)Qs[i]=f[i]*mid1[i]-e[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("\ne=");for(i=0;i<N;i++)printf("%9.4f",e[i]);printf("\nf=");for(i=0;i<N;i++)printf("%9.4f",f[i]);printf("\n");//求线路上的潮流//计算S[i][j]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);AA=-f[ydata[i].nl-1]*ydata[i].Bl+(e[ydata[i].nl-1]-e[ydata[i].nr-1])*ydata[i].R/Z2+(f[ydata[i].nl-1]-f[ydata[i].nr-1])*ydata[i].X/Z2;BB=-e[ydata[i].nl-1]*ydata[i].Bl-(f[ydata[i].nl-1]-f[ydata[i].nr-1])*ydata[i].R/Z2+(e[ydata[i].nl-1]-e[ydata[i].nr-1])*ydata[i].X/Z2;Pij[i]=e[ydata[i].nl-1]*AA-f[ydata[i].nl-1]*BB;Qij[i]=e[ydata[i].nl-1]*BB+f[ydata[i].nl-1]*AA;printf("S[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nl,ydata[i].nr,Pij[i],Qij[i]);dPij[i]=Pij[i]+Pji[i];dQij[i]=Qij[i]+Qji[i];}}printf("\n");//计算S[j][i]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);CC=-f[ydata[i].nr-1]*ydata[i].Br+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].R/Z2+(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].X/Z2;DD=-e[ydata[i].nr-1]*ydata[i].Br-(f[ydata[i].nr-1]-f[ydata[i].nl-1])*ydata[i].R/Z2+(e[ydata[i].nr-1]-e[ydata[i].nl-1])*ydata[i].X/Z2;Pji[i]=e[ydata[i].nr-1]*CC-f[ydata[i].nr-1]*DD;Qji[i]=e[ydata[i].nr-1]*DD+f[ydata[i].nr-1]*CC;printf("S[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nr,ydata[i].nl,Pji[i],Qji[i]);}}printf("\n");//计算dS[i][j]for(i=0;i<n_br;i++){if(ydata[i].nl!=ydata[i].nr){dPij[i]=Pij[i]+Pji[i];dQij[i]=Qij[i]+Qji[i];printf("dS[%d][%d]=%9.4f+j%9.4f\n",ydata[i].nl,ydata[i].nr,dPij[i],dQij[i]);}}printf("\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");}。
潮流计算源程序
潮流计算源程序/***********************************************************文件名称：潮流计算程序.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====================================================================== ========********* 潮流计算结束*********。
电力系统分析潮流计算程序
#include<stdio.h>#include<stdlib.h>#include<math.h>#define PI 3.14159//节点参数结构体structNodeType{int N;//节点号int Type;//节点类型double e;//电压幅值double f;//电压相角double Pd;//负荷有功double Qd;//负荷无功double Ps;//出力有功double Qs;//出力无功double Bc;//并联电容的电抗值};//支路参数结构体structBranchType{intNbr;//支路号intNl;//首节点int Nr;//末节点double R;//支路电阻double X;//支路电抗double Bn;//对地电抗double Kt;//支路变比};//******************************************全局变量声明***************************************int n;//节点数intnPQ;//PQ节点数intnPV;//PV节点数intnbr;//支路数intng;//发电机台数int Mark=0;//标记支路参数是否已经转换double **G;//导纳矩阵G部分double **B;//导纳矩阵B部分double *dS;//功率不平衡量double *mid1,*mid2;//求功率不平衡量时的中间变量double *Us;//电压初值double error=1;//误差值double iteration=0.000001;//误差精度double **Jacob;//雅克比矩阵double **invJac;//雅克比矩阵的逆double *dfe;//节点电压修正值structNodeType *Node;//读入时的节点参数结构体structBranchType *Branch;//读入时的支路参数结构体//*********************************************主程序******************************************void main(){voidLoadData();void FormY();//形成导纳矩阵void DeltaS();//求功率不平衡量void FormJacob();//形成雅克比矩阵void InvJac();//求雅克比矩阵的逆void UpdateU();//修正电压值voidCalculatePQ();void Print1(double *,int);void Print2(double **,int,int);intkk;//迭代次数LoadData();FormY();printf("iteration=%lf\n",iteration);kk=0;DeltaS();while(error>iteration&&kk<50){FormJacob();UpdateU();DeltaS();kk++;}printf("迭代次数为%4d\n",kk);CalculatePQ();printf("error=%e\n",error);}//***************************************************************************** ****************voidLoadData(){inti,j;inttN,tType;double te,tf,tPd,tQd,tPs,tQs,tBc;//用于重新排列节点信息的临时变量FILE *fp;//文件指针char filename[50]={""};printf("请输入数据文件名：");scanf("%s",filename);if((fp=fopen(filename,"r"))==NULL){printf("cannot open the file:data.txt\n");return;}fscanf(fp,"%d",&n);printf("节点个数为：%d\n",n);//为节点参数申请空间Node=(structNodeType *)malloc(sizeof(structNodeType)*n);//读取节点参数printf("调整前的节点参数为：\n");for(i=0;i<n;i++)fscanf(fp,"%d%d%lf%lf%lf%lf%lf%lf%lf",&Node[i].N,&Node[i].Type,&Node[i].e,&Node[i].f,&Node[ i].Pd,&Node[i].Qd,&Node[i].Ps,&Node[i].Qs,&Node[i].Bc);//计算PQ节点和PV节点的个数for(i=0;i<n;i++){if(Node[i].Type==1)nPQ++;else if(Node[i].Type==2)nPV++;}printf("PQ节点个数：%d\n",nPQ);printf("PV节点个数：%d\n",nPV);//重新排列节点参数(冒泡法)for(j=0;j<n-1;j++)for(i=0;i<n-j-1;i++){if(Node[i].Type>Node[i+1].Type){tN=Node[i].N;Node[i].N=Node[i+1].N;Node[i+1].N=tN;tType=Node[i].Type;Node[i].Type=Node[i+1].Type;Node[i+1].Type=tType;te=Node[i].e;Node[i].e=Node[i+1].e;Node[i+1].e=te;tf=Node[i].f;Node[i].f=Node[i+1].f;Node[i+1].f=tf;tPd=Node[i].Pd;Node[i].Pd=Node[i+1].Pd;Node[i+1].Pd=tPd;tQd=Node[i].Qd;Node[i].Qd=Node[i+1].Qd;Node[i+1].Qd=tQd;tPs=Node[i].Ps;Node[i].Ps=Node[i+1].Ps,Node[i+1].Ps=tPs;tQs=Node[i].Qs;Node[i].Qs=Node[i+1].Qs;Node[i+1].Qs=tQs;tBc=Node[i].Bc;Node[i].Bc=Node[i+1].Bc;Node[i+1].Bc=tBc;}}//为电压初值申请空间Us=(double *)malloc(sizeof(double)*(n-1));for(i=0;i<n-1;i++)Us[i]=Node[i].e;//读取支路参数fscanf(fp,"%d",&nbr);printf("支路个数为：%d\n",nbr);//为支路参数申请空间Branch=(structBranchType *)malloc(sizeof(structBranchType)*nbr);//读入的支路参数结构体//读入支路参数for(i=0;i<nbr;i++)fscanf(fp,"%d%d%d%lf%lf%lf%lf",&Branch[i].Nbr,&Branch[i].Nl,&Branch[i].Nr,&Branch[i].R,&Bran ch[i].X,&Branch[i].Bn,&Branch[i].Kt);//支路节点号参数调整for(i=0;i<nbr;i++){Mark=0;for(j=0;j<n;j++){if(Branch[i].Nl==Node[j].N&&Mark==0){Branch[i].Nl=j+1;Mark=1;}}}for(i=0;i<nbr;i++){Mark=0;for(j=0;j<n;j++){if(Branch[i].Nr==Node[j].N&&Mark==0){Branch[i].Nr=j+1;Mark=1;}}}fclose(fp);}//***************************************************************************** ******************************************//*********************************************形成导纳矩阵**************************************************************voidFormY(){inti,j;double Z2;//存储Z^2=R^2+X^2G=(double **)malloc(sizeof(double *)*n);//为G申请空间B=(double **)malloc(sizeof(double *)*n);//为B申请空间for(i=0;i<n;i++){G[i]=(double *)malloc(sizeof(double)*n);B[i]=(double *)malloc(sizeof(double)*n);}//初始化G、Bfor(i=0;i<n;i++)for(j=0;j<n;j++){G[i][j]=0;B[i][j]=0;}//计算非对角元素for(i=0;i<nbr;i++){Z2=Branch[i].R*Branch[i].R+Branch[i].X*Branch[i].X;if(Branch[i].Kt==0)//非变压器支路{G[Branch[i].Nl-1][Branch[i].Nr-1]-=Branch[i].R/Z2;B[Branch[i].Nl-1][Branch[i].Nr-1]+=Branch[i].X/Z2;G[Branch[i].Nr-1][Branch[i].Nl-1]=G[Branch[i].Nl-1][Branch[i].Nr-1];B[Branch[i].Nr-1][Branch[i].Nl-1]=B[Branch[i].Nl-1][Branch[i].Nr-1];}else//变压器支路{G[Branch[i].Nl-1][Branch[i].Nr-1]-=Branch[i].R/Z2/Branch[i].Kt;B[Branch[i].Nl-1][Branch[i].Nr-1]+=Branch[i].X/Z2/Branch[i].Kt;G[Branch[i].Nr-1][Branch[i].Nl-1]=G[Branch[i].Nl-1][Branch[i].Nr-1];B[Branch[i].Nr-1][Branch[i].Nl-1]=B[Branch[i].Nl-1][Branch[i].Nr-1];}}//计算对角元素for(i=0;i<n;i++)for(j=0;j<nbr;j++){Z2=Branch[j].R*Branch[j].R+Branch[j].X*Branch[j].X;if(Branch[j].Kt==0&&(Branch[j].Nl-1==i||Branch[j].Nr-1==i))//非变压器支路{G[i][i]=G[i][i]+Branch[j].R/Z2;B[i][i]=B[i][i]-Branch[j].X/Z2;}else if(Branch[j].Kt!=0&&(Branch[j].Nl-1==i||Branch[j].Nr-1==i))//变压器支路{G[i][i]=G[i][i]+Branch[j].R/Z2/Branch[j].Kt;B[i][i]=B[i][i]-Branch[j].X/Z2/Branch[j].Kt;}}//将对地电纳加入到对角元素中for(i=0;i<nbr;i++){if(Branch[i].Kt==0)//非变压器支路{B[Branch[i].Nl-1][Branch[i].Nl-1]+=Branch[i].Bn;B[Branch[i].Nr-1][Branch[i].Nr-1]+=Branch[i].Bn;}else//变压器支路{B[Branch[i].Nl-1][Branch[i].Nl-1]-=(1-Branch[i].Kt)/Branch[i].Kt/Branch[i].Kt/Branch[i].X;B[Branch[i].Nr-1][Branch[i].Nr-1]-=(Branch[i].Kt-1)/Branch[i].Kt/Branch[i].X;}}//将并联电容加入到对角元素中for(i=0;i<n;i++)B[i][i]=B[i][i]+Node[i].Bc;}//***************************************************************************** ********************//*****************************************求deltaP,deltaQ*****************************************void DeltaS()//计算功率不平衡量{inti,j;//为中间变量申请空间mid1=(double *)malloc(sizeof(double)*n);mid2=(double *)malloc(sizeof(double)*n);//为功率不平衡量申请空间dS=(double *)malloc(sizeof(double)*2*(n-1));//求功率不平衡量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]*Node[j].e-B[i][j]*Node[j].f;mid2[i]=mid2[i]+G[i][j]*Node[j].f+B[i][j]*Node[j].e;}dS[2*i]=Node[i].Ps-Node[i].Pd-(Node[i].e*mid1[i]+Node[i].f*mid2[i]);if(i<nPQ)dS[2*i+1]=Node[i].Qs-Node[i].Qd-(Node[i].f*mid1[i]-Node[i].e*mid2[i]);elsedS[2*i+1]=Us[i]*Us[i]-(Node[i].e*Node[i].e+Node[i].f*Node[i].f);}error=0;for(i=0;i<2*(n-1);i++){if(dS[i]<0&&error<-dS[i])error=-dS[i];else if(dS[i]>0&&error<dS[i])error=dS[i];}}//***************************************************************************** ********************//*********************************************雅克比矩阵******************************************voidFormJacob(){inti,j;//为雅克比行列式申请空间Jacob=(double **)malloc(sizeof(double *)*2*(n-1));for(i=0;i<2*(n-1);i++)Jacob[i]=(double *)malloc(sizeof(double)*2*(n-1));//初始化雅克比行列式for(i=0;i<2*(n-1);i++)for(j=0;j<2*(n-1);j++)Jacob[i][j]=0;for(j=0;j<n-1;j++){//求H,Nfor(i=0;i<n-1;i++){if(i!=j){Jacob[2*i][2*j]=B[i][j]*Node[i].e-G[i][j]*Node[i].f;Jacob[2*i][2*j+1]=-G[i][j]*Node[i].e-B[i][j]*Node[i].f;}else{Jacob[2*i][2*i]=B[i][i]*Node[i].e-G[i][i]*Node[i].f-mid2[i]; Jacob[2*i][2*i+1]=-G[i][i]*Node[i].e-B[i][i]*Node[i].f-mid1[i];}}//求J,Lfor(i=0;i<nPQ;i++){if(i!=j){Jacob[2*i+1][2*j]=G[i][j]*Node[i].e+B[i][j]*Node[i].f;Jacob[2*i+1][2*j+1]=B[i][j]*Node[i].e-G[i][j]*Node[i].f;}else{Jacob[2*i+1][2*i]=G[i][i]*Node[i].e+B[i][i]*Node[i].f-mid1[i]; Jacob[2*i+1][2*i+1]=B[i][i]*Node[i].e-G[i][i]*Node[i].f+mid2[i];}}//求R,Sfor(i=nPQ;i<n-1;i++){if(i==j){Jacob[2*i+1][2*i]=-2*Node[i].f;Jacob[2*i+1][2*i+1]=-2*Node[i].e;}}}}//***************************************************************************** ********************//********************************************雅克比矩阵求逆***************************************voidInvJac(){inti,j,k;double temp;//中间变量//为雅克比矩阵的逆申请空间invJac=(double **)malloc(sizeof(double *)*2*(n-1));for(i=0;i<2*(n-1);i++)invJac[i]=(double *)malloc(sizeof(double)*2*(n-1));//求逆for(i=0;i<2*(n-1);i++)for(j=0;j<2*(n-1);j++){if(i!=j)invJac[i][j]=0;elseinvJac[i][j]=1;}for(i=0;i<2*(n-1);i++){for(j=0;j<2*(n-1);j++){if(i!=j){temp=Jacob[j][i]/Jacob[i][i];for(k=0;k<2*(n-1);k++){Jacob[j][k]-=Jacob[i][k]*temp;invJac[j][k]-=invJac[i][k]*temp;}}}}for(i=0;i<2*(n-1);i++)if(Jacob[i][i]!=1){temp=Jacob[i][i];for(j=0;j<2*(n-1);j++)invJac[i][j]=invJac[i][j]/temp;}}//***************************************************************************** ********************//*********************************************电压修正********************************************voidUpdateU(){void InvJac();//求雅克比矩阵的逆inti,j;dfe=(double *)malloc(sizeof(double)2*(n-1));InvJac();for(i=0;i<2*(n-1);i++){dfe[i]=0;for(j=0;j<2*(n-1);j++)dfe[i]-=invJac[i][j]*dS[j];}for(i=0;i<n-1;i++){Node[i].e+=dfe[2*i+1];Node[i].f+=dfe[2*i];}}voidCalculatePQ(){inti,j;inttN,tType;double te,tf,tPd,tQd,tPs,tQs,tBc;//用于重新排列节点信息的临时变量//计算平衡节点功率mid1[n-1]=0;mid2[n-1]=0;for(j=0;j<n;j++){mid1[n-1]=mid1[n-1]+G[n-1][j]*Node[j].e-B[n-1][j]*Node[j].f;mid2[n-1]=mid2[n-1]+G[n-1][j]*Node[j].f+B[n-1][j]*Node[j].e;}Node[n-1].Ps=Node[n-1].e*mid1[n-1];Node[n-1].Qs=-Node[n-1].e*mid2[n-1];//计算PV节点的Qfor(i=nPQ;i<n-1;i++)Node[i].Qs=Node[i].f*mid1[i]-Node[i].e*mid2[i];//将节点参数排列为按节点号从小到大排列for(j=0;j<n-1;j++)for(i=0;i<n-j-1;i++){if(Node[i].N>Node[i+1].N){tN=Node[i].N;Node[i].N=Node[i+1].N;Node[i+1].N=tN;tType=Node[i].Type;Node[i].Type=Node[i+1].Type;Node[i+1].Type=tType;te=Node[i].e;Node[i].e=Node[i+1].e;Node[i+1].e=te;tf=Node[i].f;Node[i].f=Node[i+1].f;Node[i+1].f=tf;tPd=Node[i].Pd;Node[i].Pd=Node[i+1].Pd;Node[i+1].Pd=tPd;tQd=Node[i].Qd;Node[i].Qd=Node[i+1].Qd;Node[i+1].Qd=tQd;tPs=Node[i].Ps;Node[i].Ps=Node[i+1].Ps,Node[i+1].Ps=tPs;tQs=Node[i].Qs;Node[i].Qs=Node[i+1].Qs;Node[i+1].Qs=tQs;tBc=Node[i].Bc;Node[i].Bc=Node[i+1].Bc;Node[i+1].Bc=tBc;}}//转换为极坐标形式for(i=0;i<n;i++){te=sqrt(Node[i].e*Node[i].e+Node[i].f*Node[i].f);tf=atan(Node[i].f/Node[i].e)*180/PI;Node[i].e=te;Node[i].f=tf;}//输出结果printf("最终结果为：\n");printf(" N Tp Amp DltaPdQd Ps Qs Bc\n");for(i=0;i<n;i++)printf("%4d%4d%9.4lf%9.4lf%9.4lf%9.4lf%9.4lf%9.4lf%9.4lf\n",Node[i].N,Node[i].Type,Node[i].e, Node[i].f,Node[i].Pd,Node[i].Qd,Node[i].Ps,Node[i].Qs,Node[i].Bc);}//***************************************************************************** ********************void Print1(double *Mat,int n){int i;for(i=0;i<n;i++)printf("%9.4lf\n",Mat[i]);printf("\n");}void Print2(double **Mat,intm,int n){inti,j;for(i=0;i<m;i++){for(j=0;j<n;j++)printf("%9.4lf\t",Mat[i][j]);printf("\n");}printf("\n");}。
电力系统分析潮流计算C语言编程-pq分解法2
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");
潮流计算典型程序
function [bus_res,S_res] = PowerFlow_NR_2 % 牛顿-拉夫逊法解潮流方程的主程序bus = [ 1 1 0 0 0 3;2 1 0 -0.4 -0.3 1;3 1 0 -0.5 -0.3 1;4 1 0 -0.147 -0.27 1;5 1 0 -0.344 -0.17 1;6 1 0 -0.8 -0.36 1;7 1 0 -0.3 -0.12 1;8 1 0 -0.6 -0.3 1;9 1 0 -0.18 -0.4 1;10 1 0 1 0 2 ];line=[ 1 20.2450.3990025;230.2330.379001;340.1960.319001;450.1960.319001;560.2210.359 001;590.4530.738001;460.3920.638001;630.5760.938001;670.3430.559001;780.2940.479001;980.2820.459001;710.882 1.436000.04;910.3800.619000.04];[nb,mb]=size(bus);[nl,ml]=size(line); % 计算bus和line矩阵的行数和列数[bus,line,nPQ,nPV,nodenum] = Num_(bus,line); % 对节点重新排序的子程序Y = y_(bus,line); % 计算节点导纳矩阵的子程序myf = fopen('Result.m','w');fprintf(myf,'-- by tongchao ----\n\n');fclose(myf); % 在当前目录下生成“Result.m”文件，写入节点导纳矩阵format longEPS = 1.0e-5; % 设定误差精度for t = 1:100 % 开始迭代计算，设定最大迭代次数为100，以便不收敛情况下及时跳出[dP,dQ] = dPQ_(Y,bus,nPQ,nPV); % 计算功率偏差dP和dQ的子程序J = Jac_(bus,Y,nPQ); % 计算雅克比矩阵的子程序UD = zeros(nPQ,nPQ);for i = 1:nPQUD(i,i) = bus(i,2); % 生成电压对角矩阵enddAngU = J \ [dP;dQ];dAng = dAngU(1:nb-1,1); % 计算相角修正量dU = UD*(dAngU(nb:nb+nPQ-1,1)); % 计算电压修正量bus(1:nPQ,2) = bus(1:nPQ,2) - dU; % 修正电压bus(1:nb-1,3) = bus(1:nb-1,3) - dAng; % 修正相角if (max(abs(dU))<EPS)&(max(abs(dAng))<EPS)breakend % 判断是否满足精度误差，如满足则跳出，否则返回继续迭代endbus = PQ_(bus,Y,nPQ,nPV); % 计算每个节点的有功和无功注入的子程序[bus,line] = ReNum_(bus,line,nodenum); % 对节点恢复编号的子程序YtYm = YtYm_(line); % 计算线路的等效Yt和Ym的子程序，以计算线路潮流bus_res = bus_res_(bus); % 计算节点数据结果的子程序S_res = S_res_(bus,line,YtYm); % 计算线路潮流的子程序myf = fopen('Result.m','a');fprintf(myf,'----牛顿－拉夫逊法潮流计算结果------\n\n 节点计算结果：\n节点节点电压节点相角（角度）节点注入功率\n');for i = 1:nbfprintf(myf,'%f ',bus_res(i,1));fprintf(myf,'%f ',bus_res(i,2));fprintf(myf,'%f ',bus_res(i,3));fprintf(myf,'%f + j %f\n',real(bus_res(i,4)),imag(bus_res(i,4)));endfprintf(myf,'\n线路计算结果：\n节点I 节点J 线路功率S(I,J) 线路功率S(J,I) 线路损耗dS(I,J)\n');for i = 1:nlfprintf(myf,'%f ',S_res(i,1));fprintf(myf,'%f ',S_res(i,2));fprintf(myf,'%f + j %f ',real(S_res(i,3)),imag(S_res(i,3)));fprintf(myf,'%f + j %f ',real(S_res(i,4)),imag(S_res(i,4)));fprintf(myf,'%f + j %f\n',real(S_res(i,5)),imag(S_res(i,5)));endfclose(myf); % 迭代结束后继续在“Result.m”写入节点计算结果和线路计算结果程序结束function S_res = S_res_(bus,line,YtYm)[nl,ml]=size(line);S_res = zeros(nl,5); % S_res矩阵储存着线路潮流计算结果S_res(:,1:2) = line(:,1:2); % 前两列为节点编号for k=1:nlI = S_res(k,1);J = S_res(k,2);if (J~=0)&(I~=0)S_res(k,3)=bus(I,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,4)))-bus(I,2)*bus(J,2)*(cos(bus(I,3))+j*sin(bus(I,3)))*(conj(cos(bus(J,3))+j*sin(bus(J,3))))*conj(YtYm(k,3 ));S_res(k,4)=bus(J,2)^2*(conj(YtYm(k,3))+conj(YtYm(k,5)))-bus(I,2)*bus(J,2)*(conj(cos(bus(I,3))+j*sin(bus(I,3))))*(cos(bus(J,3))+j*sin(bus(J,3)))*conj(YtYm(k,3 ));S_res(k,5)=S_res(k,3) + S_res(k,4); % 利用公式计算非接地支路的潮流elseif J==0S_res(k,3)=bus(I,2)^2*conj(YtYm(k,4));S_res(k,5)=S_res(k,3)+S_res(k,4);elseS_res(k,4)=bus(J,2)^2*conj(YtYm(k,5));S_res(k,5)=S_res(k,3)+S_res(k,4); % 利用公式计算接地支路的潮流endendfunction bus_res = bus_res_(bus);[nb,mb]=size(bus);bus_res = zeros(nb,4); % bus_res矩阵储存着节点计算结果bus_res(:,1:2) = bus(:,1:2);bus_res(:,3) = bus(:,3) *180 / pi; % 相角采用角度制bus_res(:,4) = bus(:,4) + (sqrt(-1))*bus(:,5); % 注入功率function YtYm = YtYm_(line)[nl,ml]=size(line);YtYm = zeros(nl,5); % 对YtYm矩阵赋初值0YtYm(:,1:2) = line(:,1:2); % 矩阵前两列为线路两段节点编号，后三列分别为线路等效Yt，i侧的等效Ym，j侧的等效Ymj = sqrt(-1);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4);if Zt~=0Yt=1/Zt;endYm=line(k,5)+j*line(k,6);K=line(k,7);if (K==0)&(J~=0) % 普通线路: K=0YtYm(k,3) = Yt;YtYm(k,4) = Ym;YtYm(k,5) = Ym;endif (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0YtYm(k,4) = Ym;endif K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧YtYm(k,3) = Yt/K;YtYm(k,4) = Ym+Yt*(K-1)/K;YtYm(k,5) = Yt*(1-K)/K/K;endif K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i 侧YtYm(k,3) = -Yt*K;YtYm(k,4) = Ym+Yt*(1+K);YtYm(k,5) = Yt*(K^2+K);endendfunction [bus,line] = ReNum_(bus,line,nodenum)[nb,mb]=size(bus);[nl,ml]=size(line);bus_temp = zeros(nb,mb); % bus_temp矩阵用于临时存放bus矩阵的数据k = 1;for i = 1 :nbfor j = 1 : nbif bus(j,1) == kbus_temp(k,:) = bus(j,:);k = k + 1;endendend % 利用bus矩阵的首列编号重新对bus矩阵排序并存入bus_temp矩阵中bus = bus_temp; % 重新赋值回bus，完成bus矩阵的重新编号for i=1:nlfor j=1:2for k=1:nbif line(i,j)==nodenum(k,1)line(i,j)=nodenum(k,2);breakendendendend % 恢复line的编号function bus = PQ_(bus,Y,nPQ,nPV)n = nPQ+nPV+1; % n为总节点数for i = nPQ+1:n-1for j = 1:nbus(i,5)=bus(i,5)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endend % 利用公式计算PV节点的无功注入for j =1:nbus(n,4)=bus(n,4)+bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));bus(n,5)=bus(n,5)+bus(n,2)*bus(j,2)*(real(Y(n,j))*sin(bus(n,3)-bus(j,3))-imag(Y(n,j))*cos(bus(n,3)-bus(j,3)));end % 计算平衡节点的无功及有功注入function J = Jac_(bus,Y,nPQ)[nb,mb]=size(bus);H = zeros(nb-1,nb-1);N = zeros(nb-1,nPQ);K = zeros(nPQ,nb-1);L = zeros(nPQ,nPQ); % 将雅克比矩阵分块，即：J = [H N; K L]，并初始化Qi = zeros(nb-1,1);Pi = zeros(nb-1,1);for i = 1:nb-1for j = 1:nbPi(i,1)=Pi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));Qi(i,1)=Qi(i,1)+bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endend % 初始化并计算Qi和Pifor i = 1:nb-1for j = 1:nb-1if i~=jH(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));elseH(i,j)=Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);end % 分别计算H矩阵的对角及非对角元素if j < nPQ+1if i~=jN(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));elseN(i,j)=-Pi(i,1)-real(Y(i,j))*((bus(i,2))^2);endend % 分别计算N矩阵的对角及非对角元素if i < nPQ+1if i~=jK(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));elseK(i,j)=-Pi(i,1)+real(Y(i,j))*((bus(i,2))^2);end % 分别计算K矩阵的对角及非对角元素if j < nPQ+1if i~=jL(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));elseL(i,j)=-Qi(i,1)+imag(Y(i,j))*((bus(i,2))^2);endend % 分别计算L矩阵的对角及非对角元素endendendJ = [H N; K L]; % 生成雅克比矩阵function [dP,dQ] =dPQ_(Y,bus,nPQ,nPV) % nPQ、nPV为相应节点个数n = nPQ + nPV +1; % 总节点个数dP = bus(1:n-1,4);dQ = bus(1:nPQ,5); % 对dP和dQ赋初值PV节点不需计算dQ 平衡节点不参与计算for i = 1:n-1for j = 1:ndP(i,1) = dP(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(bus(i,3)-bus(j,3))+imag(Y(i,j))*sin(bus(i,3)-bus(j,3)));if i<nPQ+1dQ(i,1) = dQ(i,1)-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(bus(i,3)-bus(j,3))-imag(Y(i,j))*cos(bus(i,3)-bus(j,3)));endendend % 利用循环计算求取dP和dQfunction Y = y_(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb); % 对导纳矩阵赋初值0for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+j*line(k,4); % 读入线路参数if Zt~=0Yt=1/Zt; % 接地支路不计算YtendYm=line(k,5)+j*line(k,6);K=line(k,7);if (K==0)&(J~=0) % 普通线路: K=0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt+Ym;Y(I,J)=Y(I,J)-Yt;Y(J,I)=Y(I,J);endif (K==0)&(J==0) % 对地支路: K=0,J=0,R=X=0Y(I,I)=Y(I,I)+Ym;endif K>0 % 变压器线路: Zt和Ym为折算到i侧的值,K在j 侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+Yt/K/K;Y(I,J)=Y(I,J)-Yt/K;Y(J,I)=Y(I,J);endif K<0 % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=Y(J,J)+K*K*Yt;Y(I,J)=Y(I,J)+K*Yt;Y(J,I)=Y(I,J);endendfunction [bus,line,nPQ,nPV,nodenum] = Num_(bus,line)[nb,mb]=size(bus);[nl,ml]=size(line);nSW = 0; % nSW为平衡节点个数nPV = 0; % nPV为PV节点个数nPQ = 0; % nPQ为PQ节点个数for i = 1:nb, % nb为总节点数type= bus(i,6);if type == 3,nSW = nSW + 1;SW(nSW,:)=bus(i,:);% 计算并储存平衡节点elseif type == 2,nPV = nPV +1;PV(nPV,:)=bus(i,:);% 计算并储存PV节点elsenPQ = nPQ + 1;PQ(nPQ,:)=bus(i,:);% 计算并储存PQ节点endendbus=[PQ;PV;SW]; % 对bus矩阵按PQ、PV、平衡节点的顺序重新排序nodenum=[[1:nb]' bus(:,1)];% 生成新旧节点对照表for i=1:nlfor j=1:2for k=1:nbif line(i,j)==nodenum(k,2)line(i,j)=nodenum(k,1);breakendendendend % 按排序以后的节点顺序对line矩阵重新编号。
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语言程序
＃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语言程序，两行，大家可以看看，仔细研究，然后在这个基础上修改。
谢谢#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<<",迭代结束。
潮流计算程序
%潮流计算fprintf('开始潮流计算\n');fprintf('请输入待求网络的相应参数\n');%参数输入部分n=input('网络中的节点数:n=');L=input('网络中的支路数:L=');ss=input('平衡节点ss=');pr=input('误差精度:pr=');X1=input('支路参数:X1=');X2=input('节点参数:X2=');X=input('节点号和对地参数:X=');fprintf('参数输入部分结束\n\n');Y=zeros(n);%置迭代次数mm=1;%创建节点导纳矩阵for i=1:Lif X1(i,6)==0 %不含变压器的支路p=X1(i,1);q=X1(i,2);Y(p,q)=Y(p,q)-1/X1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/X1(i,3)+0.5*X1(i,4); Y(q,q)=Y(q,q)+1/X1(i,3)+0.5*X1(i,4); else %含有变压器的支路p=X1(i,1);q=X1(i,2);Y(p,q)=Y(p,q)-1/(X1(i,3)*X1(i,5)); Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/X1(i,3);Y(q,q)=Y(q,q)+1/(X1(i,5)^2*X1(i,3)); endendY;OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化%创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理if i~=ss&&X2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(X2(i,3))*(real(Y(i,j))*real(X2(j,3)) -imag(Y(i,j))*imag(X2(j,3)))+imag(X2(i,3))*(real(Y(i,j))*imag(X2(j,3) )+imag(Y(i,j))*real(X2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(X2(i,3))*(real(Y(i,j))*real(X2(j,3))-ima g(Y(i,j))*imag(X2(j,3)))-real(X2(i,3))*(real(Y(i,j))*imag(X2(j,3))+im ag(Y(i,j))*real(X2(j,3)));endendendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0if i~=ss&&X2(i,6)==0h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(X2(i,3))*(real(Y(i,j))*real(X2(j,3)) -imag(Y(i,j))*imag(X2(j,3)))+imag(X2(i,3))*(real(Y(i,j))*imag(X2(j,3) )+imag(Y(i,j))*real(X2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(X2(i,3))*(real(Y(i,j))*real(X2(j,3))-ima g(Y(i,j))*imag(X2(j,3)))-real(X2(i,3))*(real(Y(i,j))*imag(X2(j,3))+im ag(Y(i,j))*real(X2(j,3)));endendendOrgS;%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif X2(i,6)==0t=t+1;PVU(t,1)=X2(i,3);endendPVU;%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理if i~=ss&&X2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(X2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(X2(i,2))-OrgS(2*h,1);endendt=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0if i~=ss&&X2(i,6)==0h=h+1;t=t+1;DetaS(2*h-1,1)=real(X2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(X2(i,3))^2-imag(X 2(i,3))^2;endendDetaS;%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=ssh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(X2(i,3));endendI;%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if X2(i,6)==2h=h+1;for j=1:nif j~=ssk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) )+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) );Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if X2(i,6)==0h=h+1;for j=1:nif j~=ssk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) )+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(X2(i,3));Jacbi(2*h,2*k)=2*real(X2(i,3));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) );Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendJacbi;%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU;%修正节点电压j=0;for i=1:n %对PQ节点处理if X2(i,6)==2j=j+1;X2(i,3)=X2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n %对PV节点的处理if X2(i,6)==0j=j+1;X2(i,3)=X2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendX2;%开始循环********************************************************************* *while abs(max(DetaU))>prOrgS=zeros(2*n-2,1); %!!!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i~=ss&&X2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(X2(i,3))*(real(Y(i,j))*real(X2(j,3)) -imag(Y(i,j))*imag(X2(j,3)))+imag(X2(i,3))*(real(Y(i,j))*imag(X2(j,3) )+imag(Y(i,j))*real(X2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(X2(i,3))*(real(Y(i,j))*real(X2(j,3))-ima g(Y(i,j))*imag(X2(j,3)))-real(X2(i,3))*(real(Y(i,j))*imag(X2(j,3))+im ag(Y(i,j))*real(X2(j,3)));endendendfor i=1:nif i~=ss&&X2(i,6)==0h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(X2(i,3))*(real(Y(i,j))*real(X2(j,3)) -imag(Y(i,j))*imag(X2(j,3)))+imag(X2(i,3))*(real(Y(i,j))*imag(X2(j,3) )+imag(Y(i,j))*real(X2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(X2(i,3))*(real(Y(i,j))*real(X2(j,3))-ima g(Y(i,j))*imag(X2(j,3)))-real(X2(i,3))*(real(Y(i,j))*imag(X2(j,3))+im ag(Y(i,j))*real(X2(j,3)));endendendOrgS;%创建DetaSh=0;for i=1:nif i~=ss&&X2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(X2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(X2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=ss&&X2(i,6)==0h=h+1;t=t+1;DetaS(2*h-1,1)=real(X2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(X2(i,3))^2-imag(X 2(i,3))^2;endendDetaS;%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=ssh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(X2(i,3));endendI;%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif X2(i,6)==2h=h+1;for j=1:nif j~=ssk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) )+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1)); elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) );Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif X2(i,6)==0h=h+1;for j=1:nif j~=ssk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) )+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(X2(i,3));Jacbi(2*h,2*k)=2*real(X2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(X2(i,3))+real(Y(i,j))*imag(X2(i ,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(X2(i,3))+imag(Y(i,j))*imag(X2(i,3) );Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:nif X2(i,6)==2j=j+1;X2(i,3)=X2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:nif X2(i,6)==0j=j+1;X2(i,3)=X2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendX2mm=mm+1; %迭代次数加1endmm。
潮流计算程序
潮流程序说明书1.程序原始数据输入格式如下：B1是支路参数矩阵，第一列和第二列为节点编号，按照从小到大编写；第三列为支路串联阻抗参数；第四列为对地导纳参数；第五列为变压器变比；第六列为是否含有变压器，“1”为含有变压器，“2”为不含变压器。
B2是节点参数矩阵，第一列为节点注入发电功率参数；第二列为负荷功率参数；第三列为节点电压参数；其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点。
X是节点号和对地参数矩阵，第一列为节点号；第二列为对地情况，“0”表示没有对地。
2.基本原理： 采用牛顿-—拉夫逊潮流算法（直角坐标系）主要变量定义：B1—支路参数矩阵；B2—节点参数矩阵；X—节点号和对地参数；S—节点功率数值，序号为奇数的储存有功功率，偶数的储存无功功率；DS—序号为奇数的储存有功功率不平衡量，对PQ节点，序号为偶数的储存无功功率不平衡量，对PV节点，序号为偶数的储存电压不平衡量；I—节点电流数值；DetaU—节点电压修正值；Number—迭代次数；dS—线路上损耗的功率；3.程序流程框图如下：启动输入原始数据形成节点导纳矩阵Y设置节点电压初始值B2（i,3）设置迭代次数Number=1根据给定初始值计算功率数值PV 节点计算有功功率功率不平衡量DS 和电压不平衡量DetaUPQ 节点计算有功功率和无功功率不平衡量DSJacobi 矩阵是否全部形成，i>n计算Jacobi 矩阵的对角元素与非对角元素增大节点号i=i+1置节点号i=1解修正方程求出功率增量和各节点电压增量求出电压增量最大值max(DetaU)迭代是否收敛，最大值是否小于ac求出各节点功率S、各条支路首末端功率Si 、Sj,节点电压U 、节点电流I停止增大迭代次数Number=Number+1计算各节点电压值是否是否4.计算结果如下（完整计算结果见文本文档）其中：3节点运行结果如下：请输入节点数:n=3请输入支路数:n1=3请输入平衡节点号:bl=3请输入误差精度:ac=0.000001请输入支路参数：B1=[1 2 0.15i 0 1 1;2 3 0.93i 0 1 0;2 3 0.5i 0 1 0]请输入节点参数：B2=[0.9 0 1 0 0 3;0 0 1 0 0 2;0 0 0.9001 1.05 0 1]节点号和对地参数:X=[1 0;2 0;3 0]Y =0 + 6.6667i0 - 6.6667i0 + 3.0753i0 + 6.6667i0 - 9.7419i0 - 3.0753i0 0 + 3.0753iU =0.96941.00000.9001各节点的功率S为（按节点号排列）：0 - 0.1976i0 + 0.5110i0 - 0.2765i各条支路的首端功率Si为（同B1顺序一样）：0 + 0.2038i0 + 0.1074i0 + 0.1998i各条支路的末端功率Sj为（同B1顺序一样）：0 - 0.1976i0 - 0.0967i0 - 0.1798i各条支路的功率损耗dS为（同B1顺序一样）：0 + 0.0062i0 + 0.0200i0 + 0.0200i9节点运行结果如下：请输入节点数:n=9请输入支路数:n1=9请输入平衡节点号:bl=1请输入误差精度:ac=0.001请输入支路参数：B1=[1 4 0+0.0576i 0 1 1;2 7 0+0.0625i 0 1 1;3 9 0+0.0586i 0 1 1;4 5 0.01+0.085i 0.176i 1 0;4 6 0.017+0.092i 0.158i 1 0;5 7 0.032+0.161i 0.306i 1 0;6 9 0.039+0.17i 0.358i 1 0;7 0.0085+0.072i 0.149i 1 0;8 9 0.0119+0.1008i 0.209i 1 0]请输入节点参数：B2=[0 0 1.04 1.05 0 1;1.63+0i 0 1.025 1.05 0 3;0.85+0i 0 1.025 1.05 0 3;0 0 1 0 0 2;0 -1.25-0.50i 1 0 0 2;0 -0.90-0.30i 1 0 0 2;0 0 1 0 0 2;0 -1-0.35i 1 0 0 2;0 0 1 0 0 2]节点号和对地参数:X=[1 0;2 0;3 0;4 0;5 0;6 0;7 0;8 0;9 0]Y =Columns 1 through 50 -17.3611i 0 0 0 +17.3611i0 0 -16.0000i 0 0 00 0 0 -17.0648i 0 00 +17.3611i 0 0 3.3074 -39.3089i-1.3652 +11.6041i0 0 0 -1.3652 +11.6041i2.5528 -17.3382i0 0 0 -1.9422 +10.5107i0 0 +16.0000i 0 0-1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i0 0 0-1.9422 +10.5107i0 00 -1.1876 + 5.9751i0 0 -1.2820 + 5.5882i3.2242 -15.8409i-1.6171 +13.6980i0 2.8047 -35.4456i-1.1551 + 9.7843i2.7722 -23.3032i0 -1.6171 +13.6980i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i -0.0168 -0.1328iU =1.04003.2481 - 0.9555i3.2076 - 0.8593i3.2428 - 0.5596i3.2235 - 0.7453i1.8842 - 0.3587i2.3853 - 0.6321i2.4410 - 0.8344i各节点的功率S为（按节点号排列）：1.0e+002 *0.1010 - 0.3977i-0.0361 + 0.4978i-0.0571 + 0.0252i0.0687 + 1.7360i-0.0214 + 0.1542i-0.0808 - 0.4329i0.0531 - 0.5394i-0.1579 - 0.1256i0.2074 + 0.4176i-0.00180.0071-0.00710.00410.0044各条支路的首端功率Si为（同B1顺序一样）： -1.0105e+001 +1.2945e+002i3.6128 -36.2000i5.7087 - 2.3123i7.0539 - 2.2715i9.9160 +46.4206i4.8627 +14.6652i-1.8285 -15.5465i5.7060 - 4.7477i-10.1578 -16.9569i各条支路的末端功率Sj为（同B1顺序一样）： 10.1046 -39.7723i-3.6128 +49.7840i-5.7087 + 2.5192i-7.0064 + 0.7595i-6.2529 -27.7428i-4.0131 -12.9970i4.2137 +23.3616i-5.6348 + 4.4012i10.8152 +20.7070i各条支路的功率损耗dS为（同B1顺序一样）： -0.0000 +89.6771i0 +13.5840i0 + 0.2069i3.6631 +18.6778i0.8496 + 1.6681i2.3852 + 7.8152i0.0712 - 0.3465i0.6574 + 3.7501i39节点运行结果如下：请输入节点数:n=39请输入支路数:n1=46请输入平衡节点号:bl=31请输入误差精度:ac=0.1请输入支路参数：B1=[1 2 0.00350+0.04110i 0.6987i 1 0;1 39 0.00100+0.02500i 0.7500i 1 0;2 3 0.00130+0.01510i 0.2572i 1 0;2 25 0.00700+0.00860i 0.1460i 1 0;30 2 0.00000+0.01810i 0.00 102.5 1;3 4 0.00130+0.02130i 0.2214i 1 0;3 18 0.00110+0.01330i 0.2138i 1 0;4 5 0.00080+0.01280i 0.1342i 1 0;4 14 0.00080+0.01290i 0.1382i 1 0;5 6 0.00020+0.00260i -0.02170 1 0;5 8 0.00080+0.01120i 0.1476i 1 0;6 7 0.00060+0.00920i 0.1130i 1 0;6 11 0.00070+0.00820i 0.1389i 1 0;31 6 0.00000+0.02500i 0.00 107.000 1;7 8 0.00040+0.00460i 0.0780i 1 0;8 9 0.00230+0.03630i 0.3804i 1 0;9 39 0.00100+0.02500i 1.2000i 1 0;10 11 0.00040+0.00430i 0.0729i 1 0;10 13 0.00040+0.00430i 0.0729i 1 0;32 10 0.00000+0.02000i 0.00 107.000 1;11 12 0.00160+0.04350i 0.00 100.600 1;13 12 0.00160+0.04350i 0.00 100.600 1;13 14 0.00090+0.01010i 0.1723i 1 0;14 15 0.00180+0.02170i 0.3660i 1 0;15 16 0.00090+0.00940i 0.3420i 1 0;16 17 0.00070+0.00890i 0.1342i 1 0;16 19 0.00160+0.01950i 0.3040i 1 0;16 21 0.00080+0.01350i -0.12740 1 0;16 24 0.00030+0.00590i 0.0680i 1 0;17 18 0.00070+0.00820i 0.1319i 1 0;17 27 0.00130+0.01730i 0.3216i 1 0;20 19 0.00070+0.01380i 0.00 106.000 1;33 19 0.00070+0.01420i 0.00 107.000 1;34 20 0.00090+0.01800i 0.00 100.900 1;21 22 0.00080+0.01400i 0.2565i 1 0;22 23 0.00060+0.00960i 0.1846i 1 0;35 22 0.00000+0.01430i 0.00 102.500 1;23 24 0.00220+0.03500i 0.3610i 1 0;36 23 0.00050+0.02720i 0.00 100.000 1;25 26 0.00320+0.03230i 0.5130i 1 0;37 25 0.00060+0.02320i 0.00 102.500 1;26 27 0.00140+0.01470i 0.2396i 1 0;26 28 0.00430+0.04740i 0.7802i 1 0;26 29 0.00570+0.06250i 1.0290i 1 0;28 29 0.00140+0.01510i 0.2490i 1 0;38 29 0.00080+0.01560 0.00 102.500 1;] 请输入节点参数：B2=[0 0 1 0 0 2;0 0 1 0 0 2;1 0 0 2;0 -3.22-0.024i1 0 0 2;0 -5.00-1.84i0 0 1 0 0 2;0 0 1 0 0 2;1 0 0 2;0 -2.338-0.84i1 0 0 2;0 -5.22-1.76i0 0 1 0 0 2;0 0 1 0 0 2;0 0 1 0 0 2;1 0 0 2;0 -0.085-0.88i0 0 1 0 0 2;0 0 1 0 0 2;1 0 0 2;0 -3.20-1.53i1 0 0 2;0 -3.29-0.323i0 0 1 0 0 2;0 -1.58-0.30i1 0 0 2;0 0 1 0 0 2;1 0 0 2;0 -6.80-1.03i1 0 0 2;0 -2.74-1.15i0 0 1 0 0 2;1 0 0 2;0 -2.475-0.846i0 -3.086+0.922i 1 0 0 2;1 0 0 2;0 -2.24-0.472i1 0 0 2;0 -1.39-0.17i1 0 0 2;0 -2.81-0.755i1 0 0 2;0 -2.06-0.276i0 -2.835-1.269i1 0 0 2;2.50+1.3621i 0 1.04750 0 0 3;0 0 1.050 1.05 0 1;6.50+1.7590i 0 0.98310 1.05 0 3;6.32+1.0335i 0 0.99720 1.05 0 3;5.08+1.6440i 0 1.01230 1.05 0 3;6.50+2.0484i 0 1.04930 1.05 0 3;5.60+0.9688i 0 1.06350 1.05 0 3;5.40-0.0444i 0 1.02780 1.05 0 3;8.30+0.1939i 0 1.02650 1.05 0 3;1.03 1.05 0 3]10.0+0.68460i -11.04-2.50i节点号和对地参数:X=[1 0;2 0;3 0;4 0;5 0;6 0;7 0;8 0;9 0;10 0;11 0;12 0;13 0;14 0;15 0;16 0;17 0;18 0;19 0;20 0;21 0;22 0;23 0;24 0;25 0;26 0;27 0;28 0;29 0;30 0;31 0;32 0;33 0;34 0;35 0;36 0;37 0;38 0;39 0]Y =1.0e+002 *Columns 1 through 50.0365 - 0.6337i -0.0206 + 0.2416i 0 0 00 0-0.0206 + 0.2416i-0.0566 + 0.6574i0.6465 - 1.5929i-0.0285 + 0.4677i0.1469 - 1.8684i0 -0.0566 + 0.6574i-0.04860 0 -0.0285 + 0.4677i0.1251 - 2.0157i+ 0.7782i0 0 0 -0.0486 + 0.7782i0.4061 - 5.4887i0 0 0 0-0.2941 + 3.8235i0 0 0 0 00 0 0 0-0.0635 + 0.8883i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 -0.0479 + 0.7722i0 0 0 0 00 0 0 0 00 0 -0.0618 + 0.7468i 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 -0.5693 + 0.6994i 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 + 0.0054i 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.0160 + 0.3994i 0 0 0 0Columns 6 through 100 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.2941 + 3.8235i 0 -0.0635 + 0.8883i 0 00.4679 - 6.1153i -0.0706 + 1.0824i 0 0 00 0-0.0706 + 1.0824i-0.1876 + 2.1576i0.2582 - 3.2390i-0.0174 + 0.2744i0.2685 - 3.3173i0 -0.1876 + 2.1576i0 0 -0.0174 + 0.2744i 0.0334 - 0.6658i0 0 0 00.4290 - 4.6106i-0.1034 + 1.2107i 0 0 0-0.2145 + 2.3056i0 0 0 0 00 0 0 0-0.2145 + 2.3056i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 + 0.0037i 0 0 0 00 0 0 0 0 + 0.0047i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 -0.0160 + 0.3994iColumns 11 through 150 0 0 0 00 0 0 0 00 0 0 -0.0479 + 0.7722i0 0 0 0 0-0.1034 + 1.2107i 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.2145 + 2.3056i 0 -0.2145 + 2.3056i 0 00.3263 - 3.7448i -0.0001 + 0.0023i 0 0 00 0-0.0001 + 0.0023i-0.0001 + 0.0023i0.0000 - 0.0000i-0.0875 + 0.9823i0 -0.0001 + 0.0023i0.3105 - 3.5163i-0.03800.1734 - 2.2088i0 0 -0.0875 + 0.9823i+ 0.4577i0 0 0 -0.0380 + 0.4577i0.1389 - 1.5083i0 0 0 0-0.1009 + 1.0542i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0Columns 16 through 200 0 0 0 00 0 0 0 00 0 -0.0618 + 0.7468i 0 00 0 0 0 00 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.1009 + 1.0542i 0 0 0 00.3596 - 5.1047i0 -0.0418 + 0.5094i-0.0878 + 1.1167i-0.1034 + 1.2107i0 00.2344 - 2.8992i-0.0878 + 1.1167i0 -0.1034 + 1.2107i 0.1651 - 1.9557i 0 00 0 0.0418 - 0.5080i-0.0003-0.0418 + 0.5094i+ 0.0068i0 0 0 -0.0003 + 0.0068i0.0367 - 0.7228i-0.0437 + 0.7381i 0 0 0 00 0 0 0 00 0 0 0 0-0.0860 + 1.6905i 0 0 0 00 0 0 0 00 0 0 0 00 -0.0432 + 0.5748i 0 0 00 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 -0.0003 + 0.0066i0 0 0 0-0.0003 + 0.0055i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0Columns 21 through 250 0 0 0 00 0 0 0-0.5693 + 0.6994i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.0437 + 0.7381i 0 0 -0.0860 + 1.6905i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00.0838 - 1.4488i -0.0407 + 0.7120i 0 0 0-0.0649 + 1.0376i0 00.1055 - 1.7474i-0.0407 + 0.7120i-0.0179 + 0.2846i0.0827 - 1.3195i0 -0.0649 + 1.0376i0 0 -0.0179 + 0.2846i 0.1038 - 1.9730i0 0 0 00.5997 - 1.0027i0 0 0 0-0.0304 + 0.3066i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 + 0.0068i 0 0 00 0 -0.0001 + 0.0037i 0 00 0 0 0-0.0001 + 0.0042i0 0 0 0 00 0 0 0 0Columns 26 through 300 0 0 0 00 0 0 0 0 + 0.0054i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 00 0 0 0 00 0 0 0 00 -0.0432 + 0.5748i 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0-0.0304 + 0.3066i 0 0 0 0-0.0190 + 0.2092i-0.0145 + 0.1587i0.1280 - 1.3359i-0.0642 + 0.6742i-0.0642 + 0.6742i 0.1074 - 1.2461i 0 0 0-0.0609 + 0.6566i0 0.0799 - 0.8607i-0.0190 + 0.2092i0.0754 - 0.8089i0 -0.0609 + 0.6566i-0.0145 + 0.1587i0 0 0 0 0 - 0.5525i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 -0.00590 0 0 0 0Columns 31 through 350 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 + 0.0037i 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 + 0.0047i 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 -0.0003 + 0.0066i 00 0 0 -0.0003 + 0.0055i0 0 0 0 00 0 0 0 0 + 0.0068i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 00 - 0.4000i 0 0 0 00 0 - 0.5000i 0 0 00 0 0.0346 - 0.7025i 0 00 0 0 0.0277 - 0.5542i0 0 0 0 0 - 0.6993i0 0 0 0 00 0 0 0 00 0 0 0 00 0 0 0 0Columns 36 through 390 0 0 -0.0160 + 0.3994i 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0160 + 0.3994i 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0-0.0001 + 0.0037i0 0 0 00 00 -0.0001 + 0.0042i0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0059 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 00.0068 - 0.3675i0 00 0.0111 - 0.4307i0 0 0.6098 0 0 0 0 0.0319 - 0.7890i U =1.0e+003 *0.0915 - 0.1339i0.0712 - 0.1171i0.0618 - 0.0724i0.0680 - 0.0540i0.0666 - 0.0494i0.0745 - 0.0550i0.0792 - 0.0574i0.1659 - 0.0776i0.0535 - 0.0308i0.0555 - 0.0296i3.8741 +4.1652i0.0514 - 0.0316i0.0546 - 0.0614i0.0497 - 0.1055i0.0496 - 0.1237i0.0547 - 0.1336i0.0605 - 0.1277i0.0504 - 0.1256i-0.0002 + 0.0015i0.0466 - 0.1265i0.0461 - 0.1286i0.0455 - 0.1290i0.0484 - 0.1255i0.0744 - 0.1436i0.0516 - 0.1896i0.0520 - 0.1647i0.0334 - 0.2574i0.0283 - 0.2765i0.0016 - 0.0779i0.00110.0014 - 0.0726i-0.0072 - 0.0790i-0.0062 - 0.0964i0.0011 - 0.0800i-0.0023 - 0.0823i-0.0031 - 0.0825i-4.3537 + 0.6080i0.2219 - 0.0873i各节点的功率S为（按节点号排列）： 1.0e+009 *-0.0000 + 0.0000i-0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0001i 0.0002 + 0.0001i 0.0000 - 0.0001i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i -0.0000 + 0.0003i 0.0000 + 0.0000i -0.0000 + 0.0003i 0.0000 + 0.0004i 0.0000 + 0.0005i -0.0000 + 0.0004i 0.0000 + 0.0002i0.0000 + 0.0003i1.1785 + 0.0007i 0.0005 + 0.0008i 0.00000.0000-0.00000.0000-0.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000各条支路的首端功率Si为（同B1顺序一样）： 3.5400e+005 +2.4736e+005i-3.5581e+005 -2.3628e+005i1.0071e+005 +2.6013e+005i2.3682e+005 -1.6375e+005i3.7264e+003 +1.4471e+006i-8.0159e+004 +2.7957e+005i1.4664e+005 -5.0894e+004i-1.1893e+005 +8.1179e+004i-6.8733e+003 +9.6236e+004i-8.1580e+004 +1.3584e+005i-3.9332e+004 -8.2485e+004i-7.5059e+004 +2.1608e+005i-1.9406e+001 +2.7536e+005i-2.5689e+004 -1.0282e+005i-1.0682e+005 -2.1610e+005i-1.2535e+005 -4.1708e+005i-2.9483e+004 -1.4331e+004i2.6881e+004 +1.7017e+004i1.7943e+003 +1.8913e+005i2.7401e+007 +7.4282e+008i2.7399e+007 +7.4282e+008i1.3064e+005 -1.2139e+005i1.1489e+005 -1.2327e+005i7.7164e+004 -2.1257e+005i-2.8359e+004 -1.6599e+005i-1.7425e+003 -1.7302e+004i3.5700e+004 -1.7516e+004i3.8835e+004 -3.2323e+004i-1.2909e+005 +6.6903e+004i1.0067e+005 -2.4294e+005i6.7162e+004 +1.3247e+006i6.6365e+004 +1.2810e+006i-7.0330e+000 +2.1637e+002i9.7641e+003 -2.0069e+004i1.0245e+004 -4.5881e+003i2.4215e+003 +1.2980e+006i-1.4470e+004 +6.4486e+003i1.4051e+004 +6.8360e+005i1.9061e+005 -1.7744e+005i3.1770e+004 +1.1214e+006i-6.2319e+004 +3.1983e+005i1.2274e+005 -2.7770e+005i1.1898e+005 -2.7519e+005i9.8516e+004 -3.3183e+005i4.8839e+006 -7.0589e+005i各条支路的末端功率Sj为（同B1顺序一样）： -3.3804e+005 -8.4052e+004i3.5993e+005 +3.0184e+005i-9.6779e+004 -2.2023e+005i-2.1492e+005 +1.8682e+005i-3.7264e+003 -5.6637e+003i8.6098e+004 -1.8535e+005i-1.4524e+005 +6.3681e+004i1.2077e+005 -5.2863e+004i8.2090e+004 -1.2718e+005i4.0208e+004 +9.3498e+004i8.3567e+003 +9.6235e+004i8.0400e+004 -1.5426e+005i19.4058 -26.1530i2.6210e+004 +1.0809e+005i1.2059e+005 +4.2529e+005i1.3051e+005 +4.9198e+005i2.9595e+004 +1.5259e+004i-2.6774e+004 -1.6139e+004i-1.7943e+003 -1.0552e+003i-7.9694e+004 -1.7952e+004i-7.7356e+004 -1.2573e+004i-1.2279e+005 +2.0853e+005i-1.0739e+005 +2.0990e+005i-7.3846e+004 +2.4186e+005i2.9462e+004 +1.7743e+005i1.7620e+003 +1.2057e+004i-3.7913e+004 +1.8782e+004i-3.8792e+004 +3.1940e+004i1.2980e+005 -6.1207e+004i-9.6457e+004 +2.9084e+005i3.9717e+001 +1.3853e+002i-3.5171e+003 -6.0815e+003i2.4385e+001 +1.3068e+002i-9.7460e+003 +1.5660e+004i-1.0242e+004 +1.1971e+003i-2.4215e+003 -7.0055e+003i1.4506e+004 -1.2520e+004i-1.5552e+003 -3.8100e+003i-1.8260e+005 +2.4172e+005i-2.8930e+003 -4.7791e+003i6.6279e+004 -2.8646e+005i-1.1337e+005 +3.3958e+005i-1.0726e+005 +3.4408e+005i-9.6140e+004 +3.3945e+005i2.8547e+005 +7.0589e+005i各条支路的功率损耗dS为（同B1顺序一样）： 1.5956e+004 +1.6331e+005i4.1140e+003 +6.5558e+004i3.9344e+003 +3.9904e+004i2.1901e+004 +2.3079e+004i-7.8671e-011 +1.4414e+006i1.4002e+003 +1.2787e+004i 1.8394e+003 +2.8317e+004i 8.3252e+002 +1.2332e+004i 5.0953e+002 +8.6579e+003i 8.7681e+002 +1.1013e+004i 6.5898e+002 +9.2305e+003i 5.3409e+003 +6.1812e+004i -1.3866e-011 +2.7534e+005i 5.2034e+002 +5.2760e+003i 1.3774e+004 +2.0919e+005i 5.1660e+003 +7.4900e+004i 1.1258e+002 +9.2738e+002i 1.0692e+002 +8.7808e+002i 7.5033e-012 +1.8808e+005i 2.7321e+007 +7.4280e+008i 2.7322e+007 +7.4281e+008i 7.8440e+003 +8.7133e+004i 7.4944e+003 +8.6625e+004i3.3182e+003 +2.9296e+004i 1.1030e+003 +1.1435e+004i 1.9504e+001 -5.2453e+003i -2.2133e+003 +1.2653e+003i4.2504e+001 -3.8293e+002i 7.1612e+002 +5.6968e+003i 4.2119e+003 +4.7899e+004i6.7202e+004 +1.3248e+006i 6.2848e+004 +1.2749e+006i 1.7352e+001 +3.4704e+002i 1.8043e+001 -4.4090e+003i3.6386e+000 -3.3909e+003i4.4565e-011 +1.2910e+006i 3.5978e+001 -6.0715e+003i 1.2496e+004 +6.7979e+005i 8.0137e+003 +6.4281e+004i2.8877e+004 +1.1166e+006i3.9594e+003 +3.3374e+004i 9.3637e+003 +6.1880e+004i1.1719e+004 +6.8892e+004i2.3756e+003 +7.6167e+003i 5.1694e+006。
潮流计算的代码。。
运行环境：虚拟机Windows XP Professional SP2编译器：Microsoft Visual Studio 6.0（周登文老师简易版）程序主文件如下：黑色粗体表示文件名。
/**************************************************************/// stdafx.cpp : source file that includes just the standard includes// flow.pch will be the pre-compiled header// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"// TODO: reference any additional headers you need in STDAFX.H// and not in this file/**************************************************************/ // 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.Value(0)=5;ob1.Value(1)=4;ob1.Run();printf("x1=%f\n",ob1.Value(0));printf("x2=%f\n",ob1.Value(1));}void GetData()//Read the data{FILE *fp;int i;if((fp=fopen("C:\\Documents and Settings\\Administrator\\桌面\\程序模板\\data\\data.txt","r"))==NULL){printf("无法打开文件：'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].Voltage,&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,bus1,bus2;float r,x,d,g,b,g1,b1,g2,b2,g3,b3;FILE *fp;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++){if(gLine[i].k==0){bus1=gLine[i].No_I-1;bus2=gLine[i].No_J-1;r=gLine[i].R;x=gLine[i].X;d=r*r+x*x;g=r/d;b=-x/d;gY_G[bus1][bus1]=gY_G[bus1][bus1]+g;gY_G[bus2][bus2]=gY_G[bus2][bus2]+g;gY_G[bus1][bus2]=gY_G[bus1][bus2]-g;gY_G[bus2][bus1]=gY_G[bus2][bus1]-g;gY_B[bus1][bus1]=gY_B[bus1][bus1]+b+gLine[i].B; gY_B[bus2][bus2]=gY_B[bus2][bus2]+b+gLine[i].B; gY_B[bus1][bus2]=gY_B[bus1][bus2]-b;gY_B[bus2][bus1]=gY_B[bus2][bus1]-b;}else{bus1=gLine[i].No_I-1;bus2=gLine[i].No_J-1;r=gLine[i].R;x=gLine[i].X;d=r*r+x*x;g=r/d;b=-x/d;g1=g/gLine[i].k;b1=b/gLine[i].k;g2=g*(1-gLine[i].k)/(gLine[i].k*gLine[i].k);b2=b*(1-gLine[i].k)/(gLine[i].k*gLine[i].k);g3=g*(gLine[i].k-1)/gLine[i].k;b3=b*(gLine[i].k-1)/gLine[i].k;gY_G[bus1][bus1]=gY_G[bus1][bus1]+g1+g2;gY_G[bus2][bus2]=gY_G[bus2][bus2]+g1+g3;gY_G[bus1][bus2]=gY_G[bus1][bus2]-g1;gY_G[bus2][bus1]=gY_G[bus2][bus1]-g1;gY_B[bus1][bus1]=gY_B[bus1][bus1]+b1+b2;gY_B[bus2][bus2]=gY_B[bus2][bus2]+b1+b3;gY_B[bus1][bus2]=gY_B[bus1][bus2]-b1;gY_B[bus2][bus1]=gY_B[bus2][bus1]-b1;}}// output the Y matrixif((fp=fopen("C:\\Documents and Settings\\Administrator\\桌面\\程序模板\\data\\ymatrix.txt","w"))==NULL){printf("无法打开文件：'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 和else的嵌套，就留给大家了，很简单的哦~{gf[i]=gBus[i].Voltage*sin(gBus[i].Phase);ge[i]=gBus[i].Voltage*cos(gBus[i].Phase);}}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].Voltage*gBus[i+1].Voltage-(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]+g Y_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("C:\\Documents and Settings\\Administrator\\桌面\\程序模板\\data\\unbalance.txt","w"))==NULL){printf("无法打开文件：'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("C:\\Documents and Settings\\Administrator\\桌面\\程序模板\\data\\jaccobi.txt","w"))==NULL){printf("无法打开文件：'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.Value(i)=gDelta_PQ[i];ob1.Run();for(i=0;i<=Bus_Num-2;i++){gDelta_f[i]=ob1.Value(2*i);gDelta_e[i]=ob1.Value(2*i+1);gDelta_fe[2*i]=gDelta_f[i];gDelta_fe[2*i+1]=gDelta_e[i];}if((fp=fopen("C:\\Documents and Settings\\Administrator\\桌面\\程序模板\\data\\revised.txt","w"))==NULL){printf("无法打开文件：'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 GetNewValue(){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("C:\\Documents and Settings\\Administrator\\桌面\\程序模板\\data\\newvalue.txt","w"))==NULL){printf("无法打开文件：'newvalue.txt' \n");return ;}fprintf(fp,"---New Value---\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 maxValue;test();GetData();GetYMatrix();SetInitial();for(Count_Num=0;Count_Num<=100;Count_Num++){GetUnbalance();GetJaccobi();GetRevised();//GetNewValue();maxValue=fabs(gDelta_fe[0]);for(i=1;i<=2*(Bus_Num-1)-1;i++){if(maxValue<fabs(gDelta_fe[i])){maxValue=fabs(gDelta_fe[i]);}}if(maxValue<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;}/**************************************************************/ 各头文件如下：/**************************************************************/// stdafx.h : include file for standard system include files,// or project specific include files that are used frequently, but// are changed infrequently//#if !defined(AFX_STDAFX_H__B4FD5F15_56CF_4BD1_9653_6EC14D29A916__INCLUDED_)#define AFX_STDAFX_H__B4FD5F15_56CF_4BD1_9653_6EC14D29A916__INCLUDED_#if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers#include <stdio.h>// TODO: reference additional headers your program requires here//{{AFX_INSERT_LOCATION}}// Microsoft Visual C++ will insert additional declarations immediately before the previous line.#endif // !defined(AFX_STDAFX_H__B4FD5F15_56CF_4BD1_9653_6EC14D29A916__INCLUDED_) /**************************************************************/// NEquation.h: interface for the NEquation class.////by bhack#if !defined(AFX_NEQUATION_H__7D6DCBC7_B370_4F2E_A720_4FB8BA7B55D9__INCLUDED_)#define AFX_NEQUATION_H__7D6DCBC7_B370_4F2E_A720_4FB8BA7B55D9__INCLUDED_#if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000#define CHECKERRORclass NEquation{public:NEquation();virtual ~NEquation();public:void SetSize(int size);float & Data(int lhs, int rhs);double & Value(int lhs);int Run();protected:int m_nNumber;float * m_nDataBuffer;double * m_nValue;void Factorial();void Forward();void Backward();virtual float ZEROD();virtual double ZEROV();};NEquation::NEquation(){m_nDataBuffer = NULL;m_nValue = NULL;m_nNumber = 0;}NEquation::~NEquation(){if(m_nDataBuffer != NULL) delete []m_nDataBuffer;if(m_nValue != NULL) delete []m_nValue;}float NEquation::ZEROD(){return 0.0;}double NEquation::ZEROV(){return 0.0;}void NEquation::SetSize(int size){if(size < 1)return;if(m_nDataBuffer != NULL) delete []m_nDataBuffer;if(m_nValue != NULL) delete []m_nValue;m_nDataBuffer = new float[size * size];m_nValue = new double[size];for(int i=0; i< size*size; i++)m_nDataBuffer[i] = ZEROD();for(i=0; i< size; i++)m_nValue[i] = ZEROV();m_nNumber = size;}float & NEquation::Data(int lhs, int rhs){#if defined(CHECKERROR)if((lhs<m_nNumber)&&(lhs>=0)&&(rhs<m_nNumber)&&(rhs>=0)) #endifreturn m_nDataBuffer[lhs * m_nNumber + rhs];#if defined(CHECKERROR)else// AfxMessageBox("Error");cout << "Error" << endl;return m_nDataBuffer[0];#endif}double & NEquation::Value(int lhs){#if defined(CHECKERROR)if((lhs<m_nNumber)&&(lhs>=0))#endifreturn m_nValue[lhs];#if defined(CHECKERROR)else// AfxMessageBox("Error");cout << "Error" << endl;return m_nValue[0];#endif}int NEquation::Run(){Factorial();Forward();Backward();return 1;}float reverse( float ff){return (float)1.0/ff;};void NEquation::Factorial(){int i, j, k;for(i=0; i<m_nNumber; i++){// 规格化format line;Data(i, i) = reverse(Data(i, i));for(j = i+1; j< m_nNumber; j++) Data(i, j) = Data(i, i) * Data(i, j); // 消去for(j = i+1; j<m_nNumber; j++){for(k=i+1; k<m_nNumber; k++){Data(j,k) = Data(j, k) - Data(j,i) * Data(i, k);}}}}void NEquation::Forward(){int i, j;for(i=0; i< m_nNumber; i++){Value(i) = Data(i, i) * Value(i);for(j = i+1; j < m_nNumber; j++)Value(j) = Value(j) - Data(j, i) * Value(i);}}void NEquation::Backward(){int i, j;for(i = 1; i < m_nNumber; i++){for(j = 1; j < i+1; j++)Value(m_nNumber - i -1) -= Data(m_nNumber -i-1, m_nNumber -j) * Value(m_nNumber -j); }}int fangcheng(int num, float *A, float *B){NEquation ob;ob.SetSize(num);for(int i=0;i<num; i++){for(int j=0;j<num; j++)ob.Data(i, j) = A[i*num+j];ob.Value(i) = B[i];}ob.Run();for(i=0;i<num; i++){B[i] = ob.Value(i);}return 0;}#endif// !defined(AFX_NEQUATION_H__7D6DCBC7_B370_4F2E_A720_4FB8BA7B55D9__INCLUDED_) /**************************************************************///config.h 所有头文件名都是可以随意命名的，只要在cpp文件中调用的时候别搞错了。
潮流计算程序
潮流计算程序clear;clc;%读入原始数据cl1=textread('C:\Users\adm\Desktop\IEEE 5\LoadFlow_In.txt','','headerlines',1);lo_i=textread('C:\Users\adm\Desktop\IEEE 5\Load_Control_In.txt');%考虑某两个节点之间有并联支路的情况cfjd=0;cl=cl1;for m=2:cl(1,3)-cfjdfor n=m+1:cl(1,3)+1-cfjdif cl(m,2)==cl(n,2)if cl(m,3)==cl(n,3)cfjd=cfjd+1;cl(m,4)=cl(m,4)/2;%简化两条数值相同的支路cl(m,5)=cl(m,5)/2;cl(n,:)=[];endelseif cl(m,2)==cl(n,3)if cl(m,3)==cl(n,2)cl(m,4)=cl(m,4)/2;%简化两条数值相同的支路cl(m,5)=cl(m,5)/2;cl(n,:)=[];endendendendcl(1,3)=cl(1,3)-cfjd;%形成节点导纳矩阵Y，以及B'和B''Y=zeros(cl(1,1),cl(1,1));B1=zeros(cl(1,1),cl(1,1));B2=zeros(cl(1,1),cl(1,1));for m=2:cl(1,3)+1if cl(m,6)>=0 %线路%用于计算功率误差的完整的节点导纳矩阵Y(cl(m,2),cl(m,3))=-1/(cl(m,4)+cl(m,5)*1j);Y(cl(m,3),cl(m,2))=-1/(cl(m,4)+cl(m,5)*1j);Y(cl(m,2),cl(m,2))=Y(cl(m,2),cl(m,2))-Y(cl(m,2),cl(m,3))+cl(m,6)*1j;Y(cl(m,3),cl(m,3))=Y(cl(m,3),cl(m,3))-Y(cl(m,2),cl(m,3))+cl(m,6)*1j;%用于迭代的B',即Y的虚部B1(cl(m,2),cl(m,3))=imag(Y(cl(m,2),cl(m,3)));B1(cl(m,3),cl(m,2))=imag(Y(cl(m,3),cl(m,2)));B1(cl(m,2),cl(m,2))=B1(cl(m,2),cl(m,2))-B1(cl(m,2),cl(m,3));%imag(Y(cl(m,2),cl(m,2)))-cl(m,6);B1(cl(m,3),cl(m,3))=B1(cl(m,3),cl(m,3))-B1(cl(m,2),cl(m,3));%imag(Y(cl(m,3),cl(m,3)))-cl(m,6);%用于迭代的B''，忽略电阻B2(cl(m,2),cl(m,3))=1/cl(m,5);B2(cl(m,3),cl(m,2))=1/cl(m,5);B2(cl(m,2),cl(m,2))=B2(cl(m,2),cl(m,2))-B2(cl(m,2),cl(m,3))+cl(m,6);B2(cl(m,3),cl(m,3))=B2(cl(m,3),cl(m,3))-B2(cl(m,2),cl(m,3))+cl(m,6);else %变压器支路%用于计算功率误差的完整的节点导纳矩阵Y(cl(m,2),cl(m,3))=1/((cl(m,4)+cl(m,5)*1j)*cl(m,6));Y(cl(m,3),cl(m,2))=1/((cl(m,4)+cl(m,5)*1j)*cl(m,6));Y(cl(m,2),cl(m,2))=Y(cl(m,2),cl(m,2))-Y(cl(m,2),cl(m,3))+(1+cl(m,6))/(cl(m,6)*(cl(m,4)+cl(m,5)*1j));Y(cl(m,3),cl(m,3))=Y(cl(m,3),cl(m,3))-Y(cl(m,2),cl(m,3))+(1+cl(m,6))/(cl(m,6)^2*(cl(m,4)+cl(m,5)*1j));%用于迭代的B'B1(cl(m,2),cl(m,3))=imag(1/((cl(m,4)+cl(m,5)*1j)*cl(m,6)));%cl (m,5)/(cl(m,4)^2+cl(m,5)^2);B1(cl(m,3),cl(m,2))=imag(1/((cl(m,4)+cl(m,5)*1j)*cl(m,6)));%cl (m,5)/(cl(m,4)^2+cl(m,5)^2);B1(cl(m,2),cl(m,2))=B1(cl(m,2),cl(m,2))-B1(cl(m,2),cl(m,3));B1(cl(m,3),cl(m,3))=B1(cl(m,3),cl(m,3))-B1(cl(m,2),cl(m,3));%用于迭代的B''，忽略电阻B2(cl(m,2),cl(m,3))=-1/(cl(m,5)*cl(m,6));B2(cl(m,3),cl(m,2))=-1/(cl(m,5)*cl(m,6));B2(cl(m,2),cl(m,2))=B2(cl(m,2),cl(m,2))-B2(cl(m,2),cl(m,3))+(1+cl(m,6))/(cl(m,6)*cl(m,5));%%%%%%%%%%%这里做了变化，看看对不对B2(cl(m,3),cl(m,3))=B2(cl(m,3),cl(m,3))-B2(cl(m,2),cl(m,3))+(1+cl(m,6))/(cl(m,6)^2*cl(m,5));endend%考虑并联补偿装置，只影响自导纳Y和B2的对角元素forn=cl(1,3)+cl(1,4)+cl(1,5)+cl(1,6)+2:cl(1,3)+cl(1,4)+cl(1,5)+cl(1,6) +cl(1,7)+1Y(cl(n,1),cl(n,1))=Y(cl(n,1),cl(n,1))+cl(n,2)*1j;B2(cl(n,1),cl(n,1))=B2(cl(n,1),cl(n,1))+cl(n,2);end%将B'中的平衡节点去掉适应迭代B1(cl(1,2),:)=[];B1(:,cl(1,2))=[];%将B''中平衡节点和PV节点去掉%先判断一下PV节点中有没有平衡节点t=0;form=cl(1,3)+cl(1,4)+cl(1,5)+2:cl(1,3)+cl(1,4)+cl(1,5)+cl(1,6)+1 if cl(m,1)~=cl(1,2)t=t+1;endend%将PV节点和平衡节点都放在数组pvjd里边if t==cl(1,6)pvjd=zeros(cl(1,6)+1,1);for n=1:cl(1,6)pvjd(n,1)=cl(cl(1,3)+cl(1,4)+cl(1,5)+1+n,1);endpvjd(cl(1,6)+1,1)=cl(1,2);elsepvjd=zeros(cl(1,6),1);for n=1:cl(1,6)pvjd(n,1)=cl(cl(1,3)+cl(1,4)+cl(1,5)+1+n,1);endend%将pvjd里面的数据按照从小到大的顺序进行排序for c=1:length(pvjd)-1min=pvjd(c,1);for d=c+1:length(pvjd)if min>pvjd(d,1)min=pvjd(d,1);pvjd(d,1)=pvjd(c,1);pvjd(c,1)=min;endendend%将B''中的PV节点和平衡节点一起去掉s=0;%作计数用for n=1:length(pvjd)B2(pvjd(n,1)-s,:)=[];B2(:,pvjd(n,1)-s)=[];s=s+1;end%应该把电压幅值和电压相角都生成两个，方便对有功和无功误差进行计算并修正%为各节点送初始电压值，PQ节点送电压平均值1，PV节点送要维持的电压值%为方便去掉平衡节点，刚开始将平衡节点送入，之后再去掉dyfz=ones(cl(1,1),1);for c=cl(1,3)+2:cl(1,3)+cl(1,4)+cl(1,5)+1dyfz(cl(c,1),1)=cl(c,4);enddyxj=zeros(cl(1,1),1);%送入各节点有功功率和无功功率p=zeros(cl(1,1),1);q=zeros(cl(1,1),1);for d=cl(1,3)+2:cl(1,3)+cl(1,4)+cl(1,5)+1p(cl(d,1),1)=p(cl(d,1),1)+cl(d,2);q(cl(d,1),1)=q(cl(d,1),1)+cl(d,3);end%判断有功功率和无功功率的误差是否满足精度要求pmax=1;qmax=1;ddcs=0;lf_m=zeros(1,3);while(max(pmax,qmax)>lo_i(2)&&ddcs<50)%计算有功功率误差p1=zeros(cl(1,1),1);pv1=zeros(cl(1,1),1);for g=1:cl(1,1)for h=1:cl(1,1)p1(g,1)=p1(g,1)-dyfz(g,1)*dyfz(h,1)*(real(Y(g,h))*cos(dyxj(g,1)-dyxj(h,1))+imag(Y(g,h))*sin(dyxj(g,1)-dyxj(h,1)));endp1(g,1)=p(g,1)+p1(g,1);%计算有功功率误差除于电压幅值pv1(g,1)=p1(g,1)/dyfz(g,1);end%将每次运行的最大功率误差和相应的节点号求出来，并保存%%%%%为防止是平衡节点，将平衡节点置0p2=p1;p2(cl(1,2),1)=0;for r=1:cl(1,1)if p2(r,1)<0p2(r,1)=-p2(r,1);endend[pmax,pjdh]=max(p2);%将迭代次数、最大有功误差以及对应的节点编号存放在矩阵lf_m中lf_m1=zeros(1,3);lf_m1(1,1)=ddcs;lf_m1(1,2)=pjdh;lf_m1(1,3)=pmax;lf_m=[lf_m(1:ddcs+1,:);lf_m1];%将pv1的平衡节点去掉，以方便求取电压相角的修正量pv1(cl(1,2),:)=[];xj1=B1\pv1;%将电压幅值的的平衡节点去掉dyfz1=dyfz;dyfz1(cl(1,2),:)=[];xj2=zeros(cl(1,1)-1,1);%xj1除以相应的电压幅值，得到相角修正量for h=1:cl(1,1)-1xj2(h,1)=xj1(h,1)/dyfz1(h,1);end%将电压相角修正量还原成总节点xj=xj2;%判断平衡节点是不是第一个节点if cl(1,2)==1xj=[0;xj(1:cl(1,1)-1,:)];elsexj=[xj(1:cl(1,2)-1,:);0;xj(cl(1,2):cl(1,1)-1,:)];end%修正电压相角for m=1:cl(1,1)dyxj(m,1)=dyxj(m,1)-xj(m,1);end%计算无功误差q1=zeros(cl(1,1),1);qv1=zeros(cl(1,1),1);for g=1:cl(1,1)for h=1:cl(1,1)q1(g,1)=q1(g,1)-dyfz(g,1)*dyfz(h,1)*(real(Y(g,h))*sin(dyxj(g,1)-dyxj(h,1))-imag(Y(g,h))*cos(dyxj(g,1)-dyxj(h,1)));endq1(g,1)=q(g,1)+q1(g,1);%计算无功功率除以电压值qv1(g,1)=q1(g,1)/dyfz(g,1);end%将每次运行的最大功率误差和相应的节点号求出来，并保存q2=q1;for m=1:length(pvjd)q2(pvjd(m,1),1)=0;endfor r=1:cl(1,1)-length(pvjd)if q2(r,1)<0q2(r,1)=-q2(r,1);endend[qmax,qjdh]=max(q2);%将qv1的PV节点和平衡节点去掉，以方便计算电压幅值的修正量s=0;for n=1:length(pvjd)qv1(pvjd(n,1)-s,:)=[];s=s+1;endfz1=B2\qv1;%将电压幅值修正量还原成总节点数fz=fz1;%判断平衡和PV节点是不是第一个节点for n=1:length(pvjd)if pvjd(n,1)==1fz=[0;fz(1:length(fz),:)];elseif pvjd(n,1)>length(fz)fz=[fz(1:length(fz),:);0];elsefz=[fz(1:pvjd(n,1)-1,:);0;fz(pvjd(n,1):length(fz),:)];endend%修正电压幅值for m=1:cl(1,1)dyfz(m,1)=dyfz(m,1)-fz(m,1);end%判断迭代次数是否达到规定的最高次数ddcs=ddcs+1;endlf_m(1,:)=[];%将电压用向量表示dyxl=zeros(cl1(1,1),1);for m=1:cl1(1,1)dyxl(m,1)=dyfz(m,1)*cos(dyxj(m,1))+dyfz(m,1)*sin(dyxj(m,1)) *1i;end%求各条支路的电流zldl=zeros(cl1(1,3),1);zlgl1=zeros(cl1(1,3),1);zlgl2=zeros(cl1(1,3),1);for n=1:cl1(1,3)%支路为线路if cl1(n+1,6)>=0zldl(n,1)=(dyxl(cl1(n+1,2),1)-dyxl(cl1(n+1,3),1))/(cl1(n+1,4)+cl1(n+1,5)*1i);%求输电线路两端的功率zlgl1(n,1)=dyxl(cl1(n+1,2),1)*conj(zldl(n,1));zlgl2(n,1)=-dyxl(cl1(n+1,3),1)*conj(zldl(n,1));%减去线路的充电功率zlgl1(n,1)=zlgl1(n,1)-cl1(n+1,6)*dyfz(cl1(n+1,2),1)^2*1i;zlgl2(n,1)=zlgl2(n,1)-cl1(n+1,6)*dyfz(cl1(n+1,3),1)^2*1i;%支路为变压器elsezldl(n,1)=(dyxl(cl1(n+1,2),1)+dyxl(cl1(n+1,3),1)/cl1(n+1,6))/( cl1(n+1,4)+cl1(n+1,5)*1i);%求输电线路两端的功率zlgl1(n,1)=dyxl(cl1(n+1,2),1)*conj(zldl(n,1));zlgl2(n,1)=dyxl(cl1(n+1,3),1)/cl1(n+1,6)*conj(zldl(n,1));endend%求各条支路的功率损耗zlsh=zeros(cl1(1,3),1);for n=1:cl1(1,3)zlsh(n,1)=zlgl1(n,1)+zlgl2(n,1);end%求取各节点的有功、无功出力（发电机节点）和有功、无功负荷（负荷节点）%PQ节点有功和无功功率已知，PV节点需要求无功功率，平衡节点需要求有功和无功功率%计算平衡节点功率jdgl2=zeros(cl1(1,1),1);for n=1:cl1(1,1)for m=1:cl1(1,1)jdgl2(n,1)=jdgl2(n,1)+dyxl(n,1)*conj(Y(m,n))*conj(dyxl(m,1));endend%判断这些节点是单纯的PQ节点、单纯的PV节点（平衡节点也包括在内）、既是PV又是PQ节点jdgl=zeros(cl1(1,1),4);for s=cl1(1,3)+cl1(1,4)+2:cl1(1,3)+cl1(1,4)+cl1(1,5)+1jdgl(cl1(s,1),3)=-cl1(s,2);jdgl(cl1(s,1),4)=-cl1(s,3);endfor t=1:cl(1,1)jdgl(t,1)=real(jdgl2(t,1))+jdgl(t,3);jdgl(t,2)=imag(jdgl2(t,1))+jdgl(t,4);end%计算网络中的总有功、无功出力和总有功、无功负荷以及总的网络损耗zgl=zeros(8,1);for c=1:4for k=1:cl1(1,1)zgl(c,1)=zgl(c,1)+jdgl(k,c);endendfor n=1:cl1(1,3)zgl(5,1)=zgl(5,1)+zlsh(n,1);endzgl(6,1)=imag(zgl(5,1));zgl(5,1)=real(zgl(5,1));zgl(7,1)=zgl(1,1)/sqrt(zgl(1,1)^2+zgl(2,1)^2);zgl(8,1)=zgl(3,1)/sqrt(zgl(3,1)^2+zgl(4,1)^2);%将要输出的节点电压幅值、相角以及出力、负荷都集中在一个矩阵zjdgl中zjdgl=zeros(cl1(1,1),7);for s=1:cl1(1,1)zjdgl(s,1)=s;zjdgl(s,2)=dyfz(s,1);zjdgl(s,3)=dyxj(s,1)*180/3.1415926;zjdgl(s,4)=jdgl(s,1);zjdgl(s,5)=jdgl(s,2);zjdgl(s,6)=jdgl(s,3);zjdgl(s,7)=jdgl(s,4);end%将要输出的支路数据i、j侧的数据以及支路损耗zzlsj=zeros(cl1(1,3),9);for m=1:cl1(1,3)for n=1:3zzlsj(m,n)=cl1(m+1,n);endzzlsj(m,4)=real(zlgl1(m,1));zzlsj(m,5)=imag(zlgl1(m,1));zzlsj(m,6)=real(zlgl2(m,1));zzlsj(m,7)=imag(zlgl2(m,1));zzlsj(m,8)=real(zlsh(m,1));zzlsj(m,9)=imag(zlsh(m,1));%将程序运行结果输出fid=fopen('C:\Users\weepoer\Desktop\IEEE5\loadflow_out2.txt','wt');fprintf(fid,'潮流方案名： RST%d节点系统\n',cl1(1,1));fprintf(fid,'节点数： %d\n',cl1(1,1));fprintf(fid,'平衡节点编号： %d\n',cl1(1,2));fprintf(fid,'支路数： %d\n',cl1(1,3));fprintf(fid,'发电节点数： %d\n',cl1(1,4));fprintf(fid,'负荷节点数： %d\n',cl1(1,5));fprintf(fid,'PV节点数： %d\n',cl1(1,6));fprintf(fid,'并联补偿装置数： %d\n\n',cl1(1,7));fprintf(fid,'总有功发电： %6f\n',zgl(1,1));fprintf(fid,'总无功发电： %6f\n',zgl(2,1));fprintf(fid,'Costg： %6f\n',zgl(7,1));fprintf(fid,'总有功负荷： %6f\n',zgl(3,1));fprintf(fid,'总无功负荷： %6f\n',zgl(4,1));fprintf(fid,'Costl： %6f\n',zgl(8,1));fprintf(fid,'总有功网损： %6f\n',zgl(5,1));fprintf(fid,'总无功网损： %6f\n\n',zgl(6,1));fprintf(fid,'节点名称电压幅值电压相角有功出力无功出力有功损耗无功损耗\n');for m=1:cl(1,1)fprintf(fid,'%d ',zjdgl(m,1));for n=2:6fprintf(fid,'%6f ',zjdgl(m,n));endfprintf(fid,'%6f\n',zjdgl(m,7));fprintf(fid,' \n');fprintf(fid,'支路编号 I侧节点名称 J侧节点名称 I侧注入有功 I侧注入无功 J侧注入有功 J侧注入无功支路有功损耗支路无功损耗\n');for m=1:cl(1,3)for n=1:3fprintf(fid,'%d ',zzlsj(m,n));endfprintf(fid,'%6f ',zzlsj(m,4));fprintf(fid,'%6f ',zzlsj(m,5));fprintf(fid,'%6f ',zzlsj(m,6));fprintf(fid,'%6f ',zzlsj(m,7));fprintf(fid,'%6f ',zzlsj(m,8));fprintf(fid,'%6f\n',zzlsj(m,9));endfclose(fid);%输出迭代中间过程fid=fopen('C:\Users\weepoer\Desktop\IEEE5\loadflow_media2.txt','wt');fprintf(fid,'迭代次数节点名称最大误差\n');for m=1:ddcs-1for n=1:2fprintf(fid,'%u ',lf_m(m,n));endfprintf(fid,'%6f\n',lf_m(m,3));endfclose(fid);。
潮流计算程序及计算结果
附表1：计算机计算潮流程序：%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号clear;n=13;%input('请输入节点数：n=');nl=13;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[];%input('请输入由支路参数形成的矩阵：B1=');B2=[];%input('请输入各节点参数形成的矩阵：B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); %－－－－－－－修改部分－－－－－－－－－－－－ym=0;SB=100;UB=220;%ym=input('您输入的参数是标么值？（若不是则输入一个不为零的数值）'); if ym~=0%SB=input('请输入功率基准值:SB=');%UB=input('请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1矩阵B1=');disp(B1)for i=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2矩阵B2=');disp(B2)end% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于低压侧p=B1(i,1);q=B1(i,2);elsep=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵disp('导纳矩阵Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部for i=1:n %给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV节点电压给定模值endfor i=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%=========================================================== ========P=real(S);Q=imag(S);ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1:nif i~=isb %非平衡节点C(i)=0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求P',Q'V2=e(i)^2+f(i)^2; %电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素=========if B2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算=================%================= 求取Jacobi矩阵===================for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df=dQ/deX3=X2; % X2=dp/df X3=dQ/deX4=-X1; % X1=dP/de X4=dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素===========DP=P(i)-P1; % PV节点有功误差DV=V(i)^2-V2; % PV节点电压误差for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/dfX5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;endendendendend%========= 以上为求雅可比矩阵的各个元素===================== for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%======================================================== ====k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行前的k4行消去for k2=k1:N1 % k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endelsefor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endendend%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵=====for k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); %修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%=============== 计算各输出量=========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);VV=V*UB;disp(VV);disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);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('各节点的功率S为(节点号从小到大排列)：');disp(S);disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');SS=S*SB;disp(SS);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Si(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))));Siz(i)=Si(p,q);elseSi(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))));Siz(i)=Si(p,q);enddisp(Si(p,q));SSi(p,q)=Si(p,q)*SB;ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);%disp(SSi(p,q));disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(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))));Sjy(i)=Sj(q,p);elseSj(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))));Sjy(i)=Sj(q,p);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p)*SB;ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);%disp(SSj(q,p));disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);%disp(DDS(i));disp('-----------------------------------------------------');endfigure(1);subplot(2,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(2,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;subplot(2,2,3);bar(real(S));ylabel('节点注入有功');grid on;subplot(2,2,4);bar(Siz);ylabel('支路首端无功');grid on;1.冬季最大运行方式潮流计算结果:计算机运行的B1,B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.88+0.545*i 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.88+0.545*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.2+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：节点号电压值相角值支路标号首端功率末端功率支路功率损耗1 242.0000 0 1-2 148.9391+3.327637i -143-27.45476i 5.93913-24.1271i2 231.0000 -3.0341 1-11 89.1956+37.193i -88.19803-61.4687i 0.997519-24.2757i3 227.4154 -4.4383 11-3 88.19803+61.4687i -88-54.5i 0.19803+6.9687i4 226.3304 -4.9004 1-12 77.8364+36.4186i -77.2404-55.7062i 0.596039-19.2876i5 229.8318 -3.9547 12-4 77.2404+55.7062i -77-47.72i 0.24036+7.9862i6 217.7226 -8.2720 1-5 63.5438+12.2204i -61.8773-32.0321i 1.66656-19.8116i7 234.9326 -3.1047 5-6 77.2597+56.3501i -77-47.72i 0.25974+8.6301i8 226.0991 -6.2429 5-7 15.3825-24.3181i 15.5354+0.274108i 0.152961-24.044i9 231.0000 0.5851 7-8 88.20085+61.55009i -88-54.5i 0.20085+7.0501i10 231.0000 3.4950 1-7 40.806-6.05737i -40.0647-22.0555i 0.741264-28.1128i11 236.1931 -1.3348 7-13 -63.6716-39.7687i 64.0965+17.5832i 0.424934-22.1855i12 237.9346 -0.8892 13-9 -64.0965-17.5832i 64.2+21.0187i 0.10348+3.4355i13 238.1868 -2.2028 1-10 79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果图形：2.冬季最小运行方式潮流计算结果：计算机运行的B1B2矩阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.616+0.3817*i 1 0 0 20 0.539+0.3817*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.1+0.06197*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]电压调整措施：变电所1、4变压器变比：+2.5% 水电厂变压器变比：-2.5%5 234.2241 -2.0051 (12--4) 54.0234+42.2706i -53.9-38.17i 0.12342+4.1006i6 226.2159 -4.8577 (1--5) 34.1669+4.22856 -33.6427-28.276i 0.524165-24.0474i7 235.1128 -0.4737 (5--6) 54.0179+37.3169i -53.9-33.4i 0.11789+3.9169i8 223.9941 -2.4603 (5--7) -20.3752-9.04094i 20.5371-15.4558i 0.161947-24.4968i9 231.0000 3.4429 (1--7) 11.0013+1.45186i -10.8258-31.4513i 0.175442-29.9995i10 231.0000 3.9321 (7--8) 53.9768+36.0957i -53.9-33.4i 0.076797+2.6957i11 238.4187 -0.9561 (7--13) -63.6881+10.8115i 64.0874-32.7763i 0.39932-21.9648i12 239.1742 -0.6185 (13--9) -64.0874+32.7763i 64.2-29.0386i 0.11258+3.7377i13 234.9468 0.6942 (1--10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机运行结果的图形：3.夏季最大运行方式计算机计算结果：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.776+0.481*i 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.776+0.481*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：10 231.0000 3.4950 (7,8) 77.7585+53.6634i -77.6-48.1i0.1585+5.5634i11 237.0938 -1.1858 (7,13) -113.4984+9.476406i 114.6926-28.38885i 1.19422-18.9124i12 239.1742 -0.6185 (13,9) -114.6926+28.38885 115-18.18369i0.307384+10.2052i13 233.3912 3.2303 (1,10) -79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果如图：4.夏季最小运行方式：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.543+0.336*i 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.543+0.336*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]10 231.0000 3.9321 (7,8) 54.3735+36.2509i -54.3-33.67i 0.073528+2.5809i11 239.0099 -0.8518 (7,13) -113.4428+24.39707i 114.6754-43.79301i 1.23255-19.3959i12 239.8726 -0.5689 (13,9) -114.6754+43.79301i 115-33.01613i 0.324605+10.7769i13 235.9565 4.4510 (1,10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机计算结果如图：5.夏季故障运行状态:调压及无功补偿措施如下:变电所3的变压器变比为-2.5%,无功补偿容量为20Mvar。
FORTRAN语言潮流计算程序
$DEBUG$LARGEc for 2002 studensimplicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1 /BNDA/ZB(80,5),ZN(50,7),PV A(30,2) /GBDA/GBII(50,4),1 GBIJ(80,3) /FACT/JF(2,100),DJ(100),UJ(2,2000)OPEN(1,FILE='lf02.res')CALL INPUT(SB,U0,REFA)CALL YMCALL DEVR(U0)DO 18 I=1,NPVK=PV A(I,1)18 V(1,K)=PV A(I,2)VS=V(1,NN)V(1,NN)=VS*DCOS(REFA/180*3.1415926)V(2,NN)=VS*DSIN(REFA/180*3.1415926)KT=027 CALL PQ(1,KT)IF(KT.GT.20) STOPIF(DABS(AXM).LE.ES) GOTO 58CALL JACO(1)CALL SBW(1,1)KT=KT+1GOTO 2758 CONTINUEWRITE(1,*)'ITERATING NUMBER: ',KTWRITE(1,7)'MAXMUM ERROR',F(1,KKK),F(2,KKK), ' AT NODE',KKK 7 FORMAT(1X,A,2F18.11,A,I5)CALL PQ(2,1)close(1)STOPEND*----------------------------------------------------------------SUBROUTINE INPUT(SB,U0,REFA)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1 /BNDA/ZB(80,5),ZN(50,7),PV A(30,2) /NNOO/NEW(150),NOLD(50,2)COMMON /GBDA/GBII(50,4),GBIJ(80,3)OPEN(7,FILE='lf02.DA T',STA TUS='OLD')REWIND(7)READ(7,*)NN,NB,NPV,NSWRITE(1,*)'NN,NB,NPV,NS:'WRITE(1,1)NN,NB,NPV,NS WRITE(1,*)'ES,U0,SB,REFA:'C READ(7,*)ES,U0,SB,REFAU0=1.0SB=1.0REFA=0.0READ(7,*)ESWRITE(1,2)ES,U0,SB,REFA DO 8 I=1,NBREAD(7,*)(ZB(I,J),J=1,5)8 CONTINUEDO 18 I=1,NN18 READ(7,*)(ZN(I,J),J=1,5)DO 38 I=1,NPV38 READ(7,*)(PV A(I,J),J=1,2)DO 48 I=1,3*NN48 NEW(I)=0DO 58 I=1,NBK=DABS(ZB(I,1))L=DABS(ZB(I,2))NEW(K)=-158 NEW(L)=-1DO 68 I=1,NPVK=PV A(I,1)68 NEW(K)=-2NEW(NS)=-3KL=1DO 338 I=1,3*NNIF(NEW(I).NE.-1) GOTO 338NEW(I)=KLKL=KL+1338 CONTINUEDO 348 I=1,3*NNIF(NEW(I).NE.-2) GOTO 348NEW(I)=KLKL=KL+1348 CONTINUENEW(NS)=KLDO 358 I=1,3*NNIF(NEW(I).EQ.0) GOTO 358K=NEW(I)NOLD(K,1)=I358 CONTINUEDO 78 I=1,NBK=ABS(ZB(I,1))L=ABS(ZB(I,2))ZB(I,1)=ZB(I,1)/K*NEW(K)78 ZB(I,2)=ZB(I,2)/L*NEW(L)DO 98 I=1,NNK=DABS(ZN(I,1))98 ZN(I,1)=NEW(K)DO 128 I=1,NPVK=PV A(I,1)128 PV A(I,1)=NEW(K)DO 148 I=1,NBK=DABS(ZB(I,1))L=DABS(ZB(I,2))IF(K.LE.L) GOTO 148T=ZB(I,1)ZB(I,1)=ZB(I,2)ZB(I,2)=T148 CONTINUECALL REGDA(ZB,80,5,2,NB)CALL REGDA(ZN,50,7,1,NN)CALL REGDA(PV A,30,2,1,NPV)C RETURNWRITE(1,*)' (ZB(I,J),J=1,5):'DO 28 I=1,NBWRITE(1,2)(ZB(I,J),J=1,5)28 CONTINUEWRITE(1,*)' (ZN(I,J),J=1,5):'DO 168 I=1,NN168 WRITE(1,2)(ZN(I,J),J=1,5)WRITE(1,*)' (PV A(I,J),J=1,2):'DO 188 I=1,NPV188 WRITE(1,2)(PV A(I,J),J=1,2)1 FORMAT(10I5)2 FORMAT(6F12.6)END*---------------------------------------------------------------- SUBROUTINE REGDA(ARY,KR,KC,K2,KNB) implicit real*8 (a-h,p-z)real*8 ARY(KR,KC)DO 38 K=1,KNB-1I=dabs(ARY(K,1))I1=dabs(ARY(K,K2))II=III1=I1DO 48 K1=K+1,KNBJ=dabs(ARY(K1,1))J1=dabs(ARY(K1,K2))IF(I.GT.J.OR.(I.EQ.J.AND.I1.GT.J1)) THENI=JI1=J1KK=K1ENDIF48 CONTINUEIF(I.EQ.II.AND.I1.EQ.II1) GOTO 38DO 58 K1=1,KCAI=ARY(KK,K1)ARY(KK,K1)=ARY(K,K1)58 ARY(K,K1)=AI38 CONTINUEEND*----------------------------------------------------------------SUBROUTINE DEVR(U0)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1 /BNDA/ZB(80,5),ZN(50,7),PV A(30,2)DO 118 I=1,NNJ=ZN(I,1)V(1,J)=U0V(2,J)=0.D0P(1,J)=ZN(I,2)-ZN(I,4)P(2,J)=ZN(I,3)-ZN(I,5)C WRITE(1,'(F5.0,6F12.6)') (ZN(I,K),K=1,5),P(1,J),P(2,J)118 CONTINUEEND*----------------------------------------------------------SUBROUTINE YMimplicit real*8 (a-h,p-z)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKKCOMMON /GBDA/GBII(50,4),GBIJ(80,3)1 /BNDA/ZB(80,5),ZN(50,7),PV A(30,2)DO 8 I=1,NNDO 8 J=1,48 GBII(I,J)=0.D0DO 58 I=1,NBNC=0N=1IOLD=0JOLD=0DO 18 K=1,NBI=ZB(K,1)J=ZB(K,2)AK=ZB(K,5)TI=1.0TJ=1.0IF(I.LT.0) TI=AKIF(J.LT.0) TJ=AKIF(I.LT.0.OR.J.LT.0) AK=0I=dabs(ZB(K,1))J=dabs(ZB(K,2))IF(I.EQ.J) THENNC=NC+1GBII(I,2)=GBII(I,2)+ZB(K,5)GOTO 18ELSEGIJ=ZB(K,3)BIJ=ZB(K,4)RR=GIJ*GIJ+BIJ*BIJGIJ=GIJ/RRBIJ=-BIJ/RRGBII(I,1)=GBII(I,1)+GIJ/TI/TIGBII(I,2)=GBII(I,2)+BIJ/TI/TI+AKGBII(J,1)=GBII(J,1)+GIJ/TJ/TJGBII(J,2)=GBII(J,2)+BIJ/TJ/TJ+AKGBII(I,3)=GBII(I,3)+1IF(I.EQ.IOLD.AND.J.EQ.JOLD) THENN=N-1GBII(I,3)=GBII(I,3)-1ENDIFGBIJ(N,1)=GBIJ(N,1)-GIJ/TI/TJGBIJ(N,2)=GBIJ(N,2)-BIJ/TI/TJGBIJ(N,3)=JENDIFIF(I.NE.J) N=N+1IOLD=IJOLD=J18 CONTINUEDO 38 I=1,NN38 GBII(I+1,4)=GBII(I,3)+GBII(I,4)C RETURN15 FORMAT(3X,6F12.6)WRITE(1,*)' (GBII(I,J),J=1,4):'DO 108 I=1,NNWRITE(1,15)(GBII(I,J),J=1,4)108 CONTINUEKI=GBII(NN+1,4)-1WRITE(1,*)' (GBIJ(I,J),J=1,3):'DO 118 I=1,KI118 WRITE(1,15)(GBIJ(I,J),J=1,3)END*----------------------------------------------------------------SUBROUTINE PQ(III,KT)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1 /GBDA/GBII(50,4),GBIJ(80,3) /NNOO/NEW(150),NOLD(50,2)1 /BNDA/ZB(80,5),ZN(50,7),PV A(30,2)IF(III.EQ.2) GOTO 30DO 18 I=1,NNF(1,I)=0.D018 F(2,I)=0.D0DO 28 I=1,NNei=v(1,i)fi=v(2,i)VI=EI*EI+FI*FIF(1,I)=F(1,I)+VI*GBII(I,1)F(2,I)=F(2,I)-VI*GBII(I,2)IF(I.EQ.NN) GOTO 28M=GBII(I,4)N=GBII(I+1,4)-1DO 38 L=M,NJ=GBIJ(L,3)EJ=V(1,J)FJ=V(2,J)GIJ=GBIJ(L,1)BIJ=GBIJ(L,2)F(1,I)=F(1,I)+EI*(GIJ*EJ-BIJ*FJ)+FI*(GIJ*FJ+BIJ*EJ)F(2,I)=F(2,I)+FI*(GIJ*EJ-BIJ*FJ)-EI*(GIJ*FJ+BIJ*EJ)F(1,J)=F(1,J)+EJ*(GIJ*EI-BIJ*FI)+FJ*(GIJ*FI+BIJ*EI)F(2,J)=F(2,J)+FJ*(GIJ*EI-BIJ*FI)-EJ*(GIJ*FI+BIJ*EI)38 CONTINUE28 CONTINUEL=1AXM=0.D0KKK=0DO 48 I=1,NN-1F(1,I)=P(1,I)-F(1,I)F(2,I)=P(2,I)-F(2,I)IPV=PV A(L,1)IF(I.EQ.IPV) THENP(2,I)=P(2,I)-F(2,I)EI=V(1,I)FI=V(2,I)F(2,I)=PV A(L,2)*PV A(L,2)-EI*EI-FI*FIL=L+1ENDIFIF(DABS(F(1,I)).GT.DABS(AXM)) THENAXM=F(1,I)KKK=IENDIFIF(DABS(F(2,I)).GT.DABS(AXM)) THENAXM=F(2,I)KKK=IENDIF48 CONTINUEP(1,NN)=F(1,NN)P(2,NN)=F(2,NN)F(1,NN)=0.D0F(2,NN)=0.D0WRITE(1,5)'NO.=',KT,': ERR=',F(1,KKK)1 ,F(2,KKK),'AT BUS:',NOLD(KKK,1)return5 FORMA T(1X,A,I2,A,2F18.11,5X,A,I4)WRITE(1,*)'F(I)---DP(I)'goto 5030 CONTINUEDO 58 I=1,NNAI=V(1,I)BI=V(2,I)F(1,I)=DSQRT(AI*AI+BI*BI)BI=BI/AIF(2,I)=dATAN(BI)58 CONTINUEWRITE(1,*)'P(I)=Pi+j Qi:'WRITE(1,15)(P(1,NEW(I)),P(2,NEW(I)),I=1,NN)WRITE(1,*)'Vi=Ei+jFi'WRITE(1,17)(F(1,NEW(KL)),F(2,NEW(KL)),KL=1,NN)DO 68 I=1,NNF(2,I)=F(2,I)*180./3.141592668 CONTINUEWRITE(1,*)'V(I)=Vi,angle:'50 WRITE(1,16)(F(1,NEW(KL)),F(2,NEW(KL)),KL=1,NN)IF(III.EQ.1) WRITE(1,*)'V(I)=Ei+j Fi:'IF(III.EQ.1) WRITE(1,15)(V(1,NEW(I)),V(2,NEW(I)),I=1,NN)15 FORMAT(1X,8F9.4)16 FORMAT(1X,8F9.4)17 FORMAT(1X,4(F9.4,F9.6))END*----------------------------------------------------------------SUBROUTINE JACO(KK)implicit real*8 (a-h,p-z)real*8 H(100),E(100)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F COMMON /GBDA/GBII(50,4),GBIJ(80,3)1 /BNDA/ ZB(80,5),ZN(50,7),PV A(30,2)1 /FACT/ JF(2,100),DJ(100),UJ(2,2000)N2=2*NNII1=3II2=3JF(1,1)=2JF(2,1)=2DO 408 K=1,NN-1RXI=0.D0RXJ=0.D0DO 18 I=1,N2H(I)=0.D018 E(I)=0.D0L=0DO 28 I=1,NPVKPV=PV A(I,1)IF(KPV.NE.K) GOTO 28L=1GOTO 3028 CONTINUE30 R=0.D0X=0.D0DO 38 I=1,K-1M=GBII(I,4)N=GBII(I+1,4)-1DO 48 J=M,NKJ=GBIJ(J,3)IF(KJ.NE.K) GOTO 48CALL JAA(H,E,I,J,K,N2,L,KK)RXI=RXI+GBIJ(J,1)*V(1,I)-GBIJ(J,2)*V(2,I)RXJ=RXJ+GBIJ(J,1)*V(2,I)+GBIJ(J,2)*V(1,I)GOTO 3848 CONTINUE38 CONTINUEM=GBII(K,4)N=GBII(K+1,4)-1DO 58 J=M,NI=GBIJ(J,3)CALL JAA(H,E,I,J,K,N2,L,KK)RXI=RXI+GBIJ(J,1)*V(1,I)-GBIJ(J,2)*V(2,I)RXJ=RXJ+GBIJ(J,1)*V(2,I)+GBIJ(J,2)*V(1,I) 58 CONTINUERXI=RXI+GBII(K,1)*V(1,K)-GBII(K,2)*V(2,K)RXJ=RXJ+GBII(K,1)*V(2,K)+GBII(K,2)*V(1,K) GII=GBII(K,1)BII=GBII(K,2)EI=V(1,K)FI=V(2,K)J=2*KH(J-1)=BII*EI-GII*FI-RXJH(J)=-GII*EI-BII*FI-RXIIF(L.NE.1.OR.KK.NE.1) THENE(J-1)=GII*EI+BII*FI-RXIE(J)=BII*EI-GII*FI+RXJELSEE(J-1)=-2.*FIE(J)=-2.*EIENDIFR=F(1,K)X=F(2,K)IF(KK.EQ.2) THENR=0.D0X=0.D0ENDIFH(N2-1)=RE(N2-1)=XI1=2*K-2CALL JBB(H,N2,II1,II2,I1,KK)I1=2*K-1CALL JBB(E,N2,II1,II2,I1,KK)408 CONTINUERETURN5 FORMAT(3X,8I9)6 FORMAT(3X,8F9.4)WRITE(1,*)'JF,UJ'WRITE(1,5)(JF(1,KL),KL=1,10)WRITE(1,6)(UJ(1,KL),KL=1,50)WRITE(1,*) ' Factor table pointer is',ii1END*----------------------------------------------------------------SUBROUTINE JAA(H,E,I,J,K,N2,L,KK)implicit real*8 (a-h,p-z)real*8 H(N2),E(N2)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F COMMON /GBDA/GBII(50,4),GBIJ(80,3)G=GBIJ(J,1)B=GBIJ(J,2)MM=2*IH(MM-1)=B*V(1,K)-G*V(2,K)H(MM)=-G*V(1,K)-B*V(2,K)IF (L.EQ.1.AND.KK.EQ.1) GOTO 10E(MM-1)=-H(MM)E(MM)=H(MM-1)10 CONTINUERETURNEND*----------------------------------------------------------------SUBROUTINE JBB(H,N2,II1,II2,K,KK)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50),H(N2)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F COMMON /FACT/ JF(2,100),DJ(100),UJ(2,2000)K1=0IF(II1.EQ.3) GOTO 40DO 38 I=1,KIF (dabs(H(I)).LT.1.0E-15) GOTO 38IF(KK.EQ.2) THENUJ(2,II2)=H(I)UJ(2,II2+1)=III2=II2+2K1=K1+1IF(II2.GE.2000-10) WRITE(1,*) ' Factor table pointer is',ii2IF(II2.GE.2000) WRITE(1,*) ' Add UJ(2,2000) 'IF(II2.GE.2000) stopENDIFM=JF(1,I)N=JF(1,I+1)-1DO 18 J=M,NLL=UJ(1,2*J)H(LL)=H(LL)-H(I)*UJ(1,2*J-1)18 CONTINUE38 CONTINUE40 K=K+1S=H(K)IF(KK.EQ.2) THENJF(2,K+1)=JF(2,K)+K1DJ(K)=SWRITE(1,*) k,sENDIFK1=0DO 48 I=K+1,N2-1IF (dabs(H(I)).LT.1.0E-15) GOTO 48UJ(1,II1)=H(I)/SUJ(1,II1+1)=IK1=K1+1II1=II1+248 CONTINUEIF(II1.GE.2000-10) WRITE(1,*) ' Factor table pointer is',ii1IF(II1.GE.2000) WRITE(1,*) ' Add UJ(1,2000) 'IF(II1.GE.2000) stopJF(1,K+1)=JF(1,K)+K1RETURNEND*----------------------------------------------------------------SUBROUTINE SBW(KK,KK1)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F COMMON /FACT/ JF(2,100),DJ(100),UJ(2,2000)1 /BNDA/ZB(80,5),ZN(50,7),PV A(30,2) /BH/ H(100)IF(KK.EQ.2) THENDO 68 I=KK1,2*NN-2K=JF(2,I)DO 78 J=K,JF(2,I+1)-1LL=UJ(2,2*J)78 H(I)=H(I)-H(LL)*UJ(2,2*J-1)H(I)=H(I)/DJ(I)68 CONTINUEELSEDO 18 I=1,2*NN18 H(I)=0.D0ENDIFK1=2*NN-2DO 38 I=K1,1,-1K=JF(1,I)K2=JF(1,I+1)-1LL=UJ(1,2*K2)IF(LL.NE.2*NN-1.or.k.gt.k2) THENK2=K2+1GOTO 30ELSEH(I)=UJ(1,2*K2-1)ENDIF30 DO 28 J=K,K2-1L=UJ(1,2*J)28 H(I)=H(I)-UJ(1,2*J-1)*H(L)38 CONTINUEIF(KK.EQ.2) GOTO 60DO 48 I=1,NN-1K=2*IV(1,I)=V(1,I)-H(K)V(2,I)=V(2,I)-H(K-1)48 CONTINUE60 CONTINUEEND*----------------------------------------------------------------RESNN,NB,NPV,NS:5 7 1 1ES,U0,SB,REFA:0.000200 1.000000 1.000000 0.000000(ZB(I,J),J=1,5):1.0000002.000000 0.060000 0.180000 0.0000001.000000 3.000000 0.060000 0.180000 0.0000001.000000 4.000000 0.040000 0.120000 0.0000001.000000 5.000000 0.020000 0.060000 0.0000002.0000003.000000 0.010000 0.030000 0.0000002.000000 5.000000 0.080000 0.240000 0.0000003.0000004.000000 0.080000 0.240000 0.000000(ZN(I,J),J=1,5):1.000000 0.200000 0.200000 0.000000 0.0000002.000000 0.000000 0.000000 0.450000 0.1500003.000000 0.000000 0.000000 0.400000 0.0500004.000000 0.000000 0.000000 0.600000 0.1000005.000000 0.000000 0.000000 0.000000 0.000000(PV A(I,J),J=1,2):5.000000 1.060000(GBII(I,J),J=1,4):10.833333 -32.500000 4.000000 1.00000012.916667 -38.750000 2.000000 5.00000012.916667 -38.750000 1.000000 7.0000003.750000 -11.250000 0.000000 8.0000006.250000 -18.750000 0.000000 8.000000(GBIJ(I,J),J=1,3):-1.666667 5.000000 2.000000-1.666667 5.000000 3.000000-2.500000 7.500000 4.000000-5.000000 15.000000 5.000000-10.000000 30.000000 3.000000-1.250000 3.750000 5.000000-1.250000 3.750000 4.000000NO.= 0: ERR= 0.50000000000 1.10000000000 AT BUS: 2 NO.= 1: ERR= -0.0770******* -0.022******** A T BUS: 2 NO.= 2: ERR= -0.00000341934 -0.00084234440 A T BUS: 5 NO.= 3: ERR= -0.00000002401 -0.00000014864 A T BUS: 5 ITERATING NUMBER: 3MAXMUM ERROR -0.00000002401 -0.00000014864 AT NODE 4 P(I)=Pi+j Qi:1.2982 0.2445 0.2000 0.2000 -0.4500 -0.1500 -0.4000 -0.0500-0.6000 -0.1000Vi=Ei+jFi1.0600 0.000000 1.0365-0.046069 1.0088-0.083907 1.0073-0.0896091.0016-0.104414V(I)=Vi,angle:1.0600 0.0000 1.0365 -2.6396 1.0088 -4.8075 1.0073 -5.13421.0016 -5.9825DA TE5,7,1,10.0002,1.00,1,0.01,2,0.02,0.06,0.001,3,0.08,0.24,0.002,3,0.06,0.18,0.002,4,0.06,0.18,0.002,5,0.04,0.12,0.003,4,0.01,0.03,0.004,5,0.08,0.24,0.001,0.0,0.0,0.0,0.02,0.2,0.2,0.0,0.03,0.0,0.0,0.45,0.154,0.0,0.0,0.4,0.055,0.0,0.0,0.6,0.11,1.06610.833333 -32.500000 1.00000012.916667 -38.750000 5.00000012.916667 -38.750000 7.0000003.750000 -11.250000 8.0000006.250000 -18.750000 8.0000000.000000 0.000000 8.0000007-1.666667 5.000000 2.000000-1.666667 5.000000 3.000000-2.500000 7.500000 4.000000-5.000000 15.000000 5.000000-10.000000 30.000000 3.000000-1.250000 3.750000 5.000000-1.250000 3.750000 4.000000NN,NB,NPV,NS:5 7 1 1ES,U0,SB,REFA:.000200 1.000000 1.000000 .000000(GBII(I,J),J=1,3):10.833333 -32.500000 1.00000012.916667 -38.750000 5.00000012.916667 -38.750000 7.0000003.750000 -11.250000 8.0000006.250000 -18.750000 8.000000.000000 .000000 8.000000(GBIJ(I,J),J=1,3):-1.666667 5.000000 2.000000-1.666667 5.000000 3.000000-2.500000 7.500000 4.000000-5.000000 15.000000 5.000000-10.000000 30.000000 3.000000-1.250000 3.750000 5.000000-1.250000 3.750000 4.000000(ZB(I,J),J=1,5):1.0000002.000000 .060000 .180000 .0000001.000000 3.000000 .060000 .180000 .0000001.000000 4.000000 .040000 .120000 .0000001.000000 5.000000 .020000 .060000 .0000002.0000003.000000 .010000 .030000 .0000002.000000 5.000000 .080000 .240000 .0000003.0000004.000000 .080000 .240000 .000000(ZN(I,J),J=1,5):1.000000 .200000 .200000 .000000 .0000002.000000 .000000 .000000 .450000 .1500003.000000 .000000 .000000 .400000 .0500004.000000 .000000 .000000 .600000 .1000005.000000 .000000 .000000 .000000 .000000(PV A(I,J),J=1,2):5.000000 1.060000NO.= 0: ERR= .50000100000 1.10000000000 AT BUS: 2 NO.= 1: ERR= -.0770******* -.022******** A T BUS: 2 NO.= 2: ERR= -.00000341921 -.00084234361 A T BUS: 5 NO.= 3: ERR= -.00000002401 -.00000014864 A T BUS: 5 ITERATING NUMBER: 3MAXMUM ERROR -.00000002401 -.00000014864 AT NODE 4 P(I)=Pi+j Qi:。