c语言编写的牛顿拉夫逊法解潮流程序

c语言编写的牛顿拉夫逊法解潮流程序闲来无事，最近把牛拉法用c语言重写一遍，和matlab相比，c语言编写潮流程序最大的难点在于矩阵求逆，我使用的求逆方法是初等行变换法，程序段如下：#include<stdio.h>#defineN3voidmain(){inti,j,k;floatt;floatJacob[N][N]={{122},{1,3,4},{2,3,4}};/澈进行求逆的矩阵floatinv_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");}}整个程序为：//牛拉法解潮流程序〃节点数//PQ〃节点数//PQ节点数//PV节点数〃串联支路数#definen_PV1#definen_br5voidmain(){〃矩阵显示函数〃电压初值〃矩阵显示函数〃电压初值〃有功初值〃无功初值//各几点电导//阻floatUs[2*N]={1.0,0,1.0,0,1.05,0,1.05,0};floatPs[N]={0,-0.5,0.2};floatQs[N]={0,-0.3};floatG[N][N],B[N][N];电纳struct抗参数{intnl;//左节点intnr;//右节点floatR;//串联电阻值floatX;//串联电抗值floatBl;〃左节点并联电导floatBr;//右节点并联电纳}ydata[n_br]={{120,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}};

floatZ2;//Z^2=RA2+XA2各串联阻抗值的平方floate[N],f[N],dfe[2*(N-1)]; //e,f存储电压的x轴分量和y轴分量,dfe存储电压修正值floatmid1[N],mid2[N],dS[2*(N-1)]; //mid1、mid2存储计算雅克比行列式对角线元素的中间值dS存储PQU的不平衡量floatJacob[2*(N-1)][2*(N-1)],inv_J[2*(N-1)][2*(N-1)]; //雅克比行列式float dPQU=1.0;//PQU不平衡量最大值int kk=0;〃迭代次数inti,j,k;floatt;floatPij[n_br]; //存储线路i->j的有功floatQij[n_br]; 〃存储线路i->j的无功floatPji[n_br]; 〃存储线路j->i的有功floatQji[n_br]; 〃存储线路j->i的无功floatdPij[n_br]; 〃存储线路i->j的有功损耗floatdQij[n_br]; 〃存储线路i->j的无功损耗floatAA,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];elseif(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项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])*yda