c語(yǔ)言編寫(xiě)的牛頓拉夫遜法解潮流程序

c語(yǔ)言編寫(xiě)的牛頓拉夫遜法解潮流程序閑來(lái)無(wú)事，最近把牛拉法用c語(yǔ)言重寫(xiě)一遍，和matlab相比，c語(yǔ)言編寫(xiě)潮流程序最大的難點(diǎn)在于矩陣求逆，我使用的求逆方法是初等行變換法，程序段如下：#include<stdio.h>#defineN3voidmain(){inti,j,k;floatt;floatJacob[N][N]={{122},{1,3,4},{2,3,4}};/澈進(jìn)行求逆的矩陣floatinv_J[N][N];//逆矩陣存儲(chǔ)于此〃初始化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;}〃將原矩陣化簡(jiǎn)為對(duì)角陣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;}}}}//原矩陣各對(duì)角元素化為1，畫(huà)出逆矩陣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");}}整個(gè)程序?yàn)椋?/牛拉法解潮流程序〃節(jié)點(diǎn)數(shù)//PQ〃節(jié)點(diǎn)數(shù)//PQ節(jié)點(diǎn)數(shù)//PV節(jié)點(diǎn)數(shù)〃串聯(lián)支路數(shù)#definen_PV1#definen_br5voidmain(){〃矩陣顯示函數(shù)〃電壓初值〃矩陣顯示函數(shù)〃電壓初值〃有功初值〃無(wú)功初值//各幾點(diǎn)電導(dǎo)//阻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抗參數(shù){intnl;//左節(jié)點(diǎn)intnr;//右節(jié)點(diǎn)floatR;//串聯(lián)電阻值floatX;//串聯(lián)電抗值floatBl;〃左節(jié)點(diǎn)并聯(lián)電導(dǎo)floatBr;//右節(jié)點(diǎn)并聯(lián)電納}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各串聯(lián)阻抗值的平方floate[N],f[N],dfe[2*(N-1)]; //e,f存儲(chǔ)電壓的x軸分量和y軸分量,dfe存儲(chǔ)電壓修正值floatmid1[N],mid2[N],dS[2*(N-1)]; //mid1、mid2存儲(chǔ)計(jì)算雅克比行列式對(duì)角線元素的中間值dS存儲(chǔ)PQU的不平衡量floatJacob[2*(N-1)][2*(N-1)],inv_J[2*(N-1)][2*(N-1)]; //雅克比行列式float dPQU=1.0;//PQU不平衡量最大值int kk=0;〃迭代次數(shù)inti,j,k;floatt;floatPij[n_br]; //存儲(chǔ)線路i->j的有功floatQij[n_br]; 〃存儲(chǔ)線路i->j的無(wú)功floatPji[n_br]; 〃存儲(chǔ)線路j->i的有功floatQji[n_br]; 〃存儲(chǔ)線路j->i的無(wú)功floatdPij[n_br]; 〃存儲(chǔ)線路i->j的有功損耗floatdQij[n_br]; 〃存儲(chǔ)線路i->j的無(wú)功損耗floatAA,BB,CC,DD; 〃存儲(chǔ)線路潮流計(jì)算時(shí)的中間值//形成導(dǎo)納矩陣 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);〃串聯(lián)阻抗等效導(dǎo)納值〃非對(duì)角元素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;〃對(duì)角元素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);〃并聯(lián)導(dǎo)納等效導(dǎo)納值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){〃計(jì)算功率不平衡量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項(xiàng)i][j]*dS[j];}for(i=0;i<N-1;i++){e[i]+=dfe[2*i+1];f[i]+=dfe[2*i];}}elsebreak;}〃循環(huán)結(jié)束 //求平衡節(jié)點(diǎn)功率 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];// //顯示輸出結(jié)果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");//求線路上的潮流//計(jì)算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");//計(jì)算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