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");}。
