牛拉法计算C程序

#include#include#include#define PI 3.14159struct NodeTypeint N;int Type;double e;double f;double Pd;double Qd;double Ps;double Qs;double Bc;struct BranchTypeint Nbr;int Nl;int Nr;double R;double X;double Bn;double Kt;int n;int nPQ;int nPV;int nPH;int nbr;int ng;int Mark=0;double *G;double *B;double *dS;double *mid1,*mid2;double *Us;double error=1;double iteration=0.000001;double *Jacob;double *invJac;double *dfe;struct NodeType *Node;struct BranchType *Branch;void main()void LoadData();void FormY();void DeltaS();void FormJacob();void InvJac();void UpdateU();void CalculatePQ();int kk;LoadData();FormY();printf( iteration=%lfn,iteration);kk=0;DeltaS();while(erroriteration&kk50)FormJacob();UpdateU();DeltaS();kk+;printf(%4dn,kk);CalculatePQ();printf( error=%en,error);printf(n);printf(n);printf( :nNTp1PQ2PV3nAmpDltaPdQdnPsQsBcn);void LoadData()int i,j;int tN,tType;double te,tf,tPd,tQd,tPs,tQs,tBc;FILE *fp;char filename50= ;printf(TXT);scanf(%s,filename);if(fp=fopen(filename,r)=NULL)printf(cannot open the file:data.txtn);return;fscanf(fp,%d,&n);printf(%dn,n);Node=(struct NodeType *)malloc(sizeof(struct NodeType)*n);printf(n);for(i=0;in;i+)fscanf(fp,%d%d%lf%lf%lf%lf%lf%lf%lf,&Nodei.N,&Nodei.Type,&Nodei.e,&Nodei.f,&Nodei.Pd,&Nodei.Qd,&Nodei.Ps,&Nodei.Qs,&Nodei.Bc);for(i=0;in;i+)if(Nodei.Type=1)nPQ+;else if(Nodei.Type=2)nPV+;else if(Nodei.Type=3) nPH+;printf(PQ%dn,nPQ);printf(PV%dn,nPV);printf(%dn,nPH);for(j=0;jn-1;j+)for(i=0;iNodei+1.Type)tN=Nodei.N;Nodei.N=Nodei+1.N;Nodei+1.N=tN;tType=Nodei.Type;Nodei.Type=Nodei+1.Type;Nodei+1.Type=tType;te=Nodei.e;Nodei.e=Nodei+1.e;Nodei+1.e=te;tf=Nodei.f;Nodei.f=Nodei+1.f;Nodei+1.f=tf;tPd=Nodei.Pd;Nodei.Pd=Nodei+1.Pd;Nodei+1.Pd=tPd;tQd=Nodei.Qd;Nodei.Qd=Nodei+1.Qd;Nodei+1.Qd=tQd;tPs=Nodei.Ps;Nodei.Ps=Nodei+1.Ps,Nodei+1.Ps=tPs;tQs=Nodei.Qs;Nodei.Qs=Nodei+1.Qs;Nodei+1.Qs=tQs;tBc=Nodei.Bc;Nodei.Bc=Nodei+1.Bc;Nodei+1.Bc=tBc;Us=(double *)malloc(sizeof(double)*(n-1);for(i=0;in-1;i+)Usi=Nodei.e;fscanf(fp,%d,&nbr);printf(%dn,nbr);Branch=(struct BranchType *)malloc(sizeof(struct BranchType)*nbr);for(i=0;inbr;i+)fscanf(fp,%d%d%d%lf%lf%lf%lf,&Branchi.Nbr,&Branchi.Nl,&Branchi.Nr,&Branchi.R,&Branchi.X,&Branchi.Bn,&Branchi.Kt);for(i=0;inbr;i+)Mark=0;for(j=0;jn;j+)if(Branchi.Nl=Nodej.N&Mark=0)Branchi.Nl=j+1;Mark=1;for(i=0;inbr;i+)Mark=0;for(j=0;jn;j+)if(Branchi.Nr=Nodej.N&Mark=0)Branchi.Nr=j+1;Mark=1;fclose(fp);void FormY()int i,j;double Z2;G=(double *)malloc(sizeof(double *)*n);B=(double *)malloc(sizeof(double *)*n);for(i=0;in;i+)Gi=(double *)malloc(sizeof(double)*n);Bi=(double *)malloc(sizeof(double)*n);for(i=0;in;i+)for(j=0;jn;j+)Gij=0;Bij=0;for(i=0;inbr;i+)Z2=Branchi.R*Branchi.R+Branchi.X*Branchi.X;if(Branchi.Kt=0)GBranchi.Nl-1Branchi.Nr-1-=Branchi.R/Z2;BBranchi.Nl-1Branchi.Nr-1+=Branchi.X/Z2;GBranchi.Nr-1Branchi.Nl-1=GBranchi.Nl-1Branchi.Nr-1;BBranchi.Nr-1Branchi.Nl-1=BBranchi.Nl-1Branchi.Nr-1;elseGBranchi.Nl-1Branchi.Nr-1-=Branchi.R/Z2/Branchi.Kt;BBranchi.Nl-1Branchi.Nr-1+=Branchi.X/Z2/Branchi.Kt;GBranchi.Nr-1Branchi.Nl-1=GBranchi.Nl-1Branchi.Nr-1;BBranchi.Nr-1Branchi.Nl-1=BBranchi.Nl-1Branchi.Nr-1;for(i=0;in;i+)for(j=0;jnbr;j+)Z2=Branchj.R*Branchj.R+Branchj.X*Branchj.X;if(Branchj.Kt=0&(Branchj.Nl-1=i|Branchj.Nr-1=i)Gii=Gii+Branchj.R/Z2;Bii=Bii-Branchj.X/Z2;else if(Branchj.Kt!=0&(Branchj.Nl-1=i|Branchj.Nr-1=i)Gii=Gii+Branchj.R/Z2/Branchj.Kt;Bii=Bii-Branchj.X/Z2/Branchj.Kt;for(i=0;inbr;i+)if(Branchi.Kt=0)BBranchi.Nl-1Branchi.Nl-1+=Branchi.Bn;BBranchi.Nr-1Branchi.Nr-1+=Branchi.Bn;elseBBranchi.Nl-1Branchi.Nl-1-=(1-Branchi.Kt)/Branchi.Kt/Branchi.Kt/Branchi.X;BBranchi.Nr-1Branchi.Nr-1-=(Branchi.Kt-1)/Branchi.Kt/Branchi.X;for(i=0;in;i+)Bii=Bii+Nodei.Bc;void DeltaS()int i,j;mid1=(double *)malloc(sizeof(double)*n);mid2=(double *)malloc(sizeof(double)*n);dS=(double *)malloc(sizeof(double)*2*(n-1);for(i=0;in-1;i+)mid1i=0;mid2i=0;for(j=0;jn;j+)mid1i=mid1i+Gij*Nodej.e-Bij*Nodej.f;mid2i=mid2i+Gij*Nodej.f+Bij*Nodej.e;dS2*i=Nodei.Ps-Nodei.Pd-(Nodei.e*mid1i+Nodei.f*mid2i);if(inPQ)dS2*i+1=Nodei.Qs-Nodei.Qd-(Nodei.f*mid1i-Nodei.e*mid2i);elsedS2*i+1=Usi*Usi-(Nodei.e*Nodei.e+Nodei.f*Nodei.f);error=0;for(i=0;i2*(n-1);i+)if(dSi0&error0&errordSi)error=dSi;void FormJacob()int i,j;Jacob=(double *)malloc(sizeof(double *)*2*(n-1);for(i=0;i2*(n-1);i+)Jacobi=(double *)malloc(sizeof(double)*2*(n-1);for(i=0;i2*(n-1);i+)for(j=0;j2*(n-1);j+)Jacobij=0;for(j=0;jn-1;j+)for(i=0;in-1;i+)if(i!=j)Jacob2*i2*j=Bij*Nodei.e-Gij*Nodei.f;Jacob2*i2*j+1=-Gij*Nodei.e-Bij*Nodei.f;elseJacob2*i2*i=Bii*Nodei.e-Gii*Nodei.f-mid2i;Jacob2*i2*i+1=-Gii*Nodei.e-Bii*Nodei.f-mid1i;for(i=0;inPQ;i+)if(i!=j)Jacob2*i+12*j=Gij*Nodei.e+Bij*Nodei.f;Jacob2*i+12*j+1=Bij*Nodei.e-Gij*Nodei.f;elseJacob2*i+12*i=Gii*Nodei.e+Bii*Nodei.f-mid1i;Jacob2*i+12*i+1=Bii*Nodei.e-Gii*Nodei.f+mid2i;for(i=nPQ;in-1;i+)if(i=j)Jacob2*i+12*i=-2*Nodei.f;Jacob2*i+12*i+1=-2*Nodei.e;void InvJac()int i,j,k;double temp;invJac=(double *)malloc(sizeof(double *)*2*(n-1);for(i=0;i2*(n-1);i+)invJaci=(double *)malloc(sizeof(double)*2*(n-1);for(i=0;i2*(n-1);i+)for(j=0;j2*(n-1);j+)if(i!=j)invJacij=0;elseinvJacij=1;for(i=0;i2*(n-1);i+)for(j=0;j2*(n-1);j+)if(i!=j)temp=Jacobji/Jacobii;for(k=0;k2*(n-1);k+)Jacobjk-=Jacobik*temp;invJacjk-=invJacik*temp;for(i=0;i2*(n-1);i+)if(Jacobii!=1)temp=Jacobii;for(j=0;j2*(n-1);j+)invJacij=invJacij/temp;void UpdateU()void InvJac();int i,j;dfe=(double *)malloc(sizeof(double)2*(n-1);InvJac();for(i=0;i2*(n-1);i+)dfei=0;for(j=0;j2*(n-1);j+)dfei-=invJacij*dSj;for(i=0;in-1;i+)Nodei.e+=dfe2*i+1;Nodei.f+=dfe2*i;void CalculatePQ()int i,j;int tN,tType;double te,tf,tPd,tQd,tPs,tQs,tBc;mid1n-1=0;mid2n-1=0;for(j=0;jn;j+)mid1n-1=mid1n-1+Gn-1j*Nodej.e-Bn-1j*Nodej.f;mid2n-1=mid2n-1+Gn-1j*Nodej.f+Bn-1j*Nodej.e;Noden-1.Ps=Noden-1.e*mid1n-1;Noden-1.Qs=-Noden-1.e*mid2n-1;for(i=nPQ;in-1;i+)Nodei.Qs=Nodei.f*mid1i-Nodei.e*mid2i;for(j=0;jn-1;j+)for(i=0;iNodei+1.N)tN=Nodei.N;Nodei.N=Nodei+1.N;Nodei+1.N=tN;tType=Nodei.Type;Nodei.Type=Nodei+1.Type;