牛拉法计算C程序

1、#include<stdio.h> #include<stdlib.h> #include<math.h> #define PI 3.14159 struct NodeType int 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

2、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 *i nvJac;double *dfe;struct NodeType *Node;struct Bran chType *Bra nch;void mai n()void LoadData();void FormY();void DeltaS();void Fo

3、rmJacob();void In vJac();void UpdateU();void CalculatePQ();int kk;LoadData();FormY();printf(" 误差？度 iteration=%lfn",iteration);kk=O;DeltaS();while(error>iteratio n&&kk<50)FormJacob();UpdateU();DeltaS();kk+;printf("迭代次数为 %4dn",kk);CalculatePQ();prin tf("误差值 error

4、=%en",error);prin tf("n");prin tf("n");printf("注:nN为节点号，Tp为节点类型(其中1为PQ节点，2为PV节点，3为平衡节点？c) , nAmp为电压大小，Dlta为相？角，Pd为节点输出的有功功率， Qd为节点输出的无功 功率，nPs为发电”涑鲇泄 B剩琎s为发电”涑鑫薰 B剩珺c为并联电容的电抗 值n");void LoadData()int i,j;int tN,tType;double te,tf,tPd,tQd,tPs,tQs,tBc;FILE *fp; char f

5、ilename50=" "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(struc

6、t NodeType)*n)；printf(" 调整前的节点参数为： n")； for(i=0；i<n；i+)fscanf(fp,"%d%d%lf%lf%lf%lf%lf%lf%lf",&Nodei.N,&Nodei.Type,&Nodei.e,&Node i.f,&Nodei.Pd,&Nodei.Qd,&Nodei.Ps,&Nodei.Qs,&Nodei.Bc)； for(i=0；i<n；i+)if(Nodei.Type=1)nPQ+；else if(Nodei.Typ

7、e=2)nPV+；else if(Nodei.Type=3)nPH+；printf("PQ节点个数： %dn",nPQ)；printf("PV节点个数： %dn",nPV)；printf(" 平衡节点个数： %dn",nPH)；for(j=0；j<n-1；j+)for(i=0；i<n-j-1；i+)if(Nodei.Type>Nodei+1.Type)tN=Nodei.N；Nodei.N=Nodei+1.N；Nodei+1.N=tN； tType=Nodei.Type；Nodei.Type=Nodei+1.Type；

8、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

9、=Nodei.Bc;Nodei.Bc=Nodei+1.Bc;Nodei+1.Bc=tBc;Us=(double *)malloc(sizeof(double)*(n-1); for(i=0;i<n-1;i+)Usi=Nodei.e;fscanf(fp,"%d",&nbr);printf(" 支路个数为： %dn",nbr);Branch=(struct BranchType *)malloc(sizeof(struct BranchType)*nbr);for(i=0;i<nbr;i+)fscanf(fp,"%d%d%d%l

10、f%lf%lf%lf",&Branchi.Nbr,&Branchi.Nl,&Branchi.Nr,&Branc hi.R,&Branchi.X,&Branchi.Bn,&Branchi.Kt);for(i=0;i<nbr;i+)Mark=0;for(j=0;j<n;j+)if(Branchi.Nl=Nodej.N&&Mark=0)Branchi.Nl=j+1;Mark=1; for(i=0;i<nbr;i+)Mark=0;for(j=0;j<n;j+)if(Branchi.Nr=Nodej

11、.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;i<n;i+)Gi=(double *)malloc(sizeof(double)*n);Bi=(double *)malloc(sizeof(double)*n);for(i=0;i<n;i+)for(j=0;j<n;j+)Gij=0;Bij=

12、0;for(i=0;i<nbr;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-=Bra

13、nchi.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;i<n;i+)for(j=0;j<nbr;j+)Z2=Branchj.R*Branchj.R+Branchj.X*Branchj.X;if(Branchj.Kt=0&&(Branchj.Nl-1=i|B

14、ranchj.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;i<nbr;i+)if(Branchi.Kt=0)BBranchi.Nl-1Branchi.Nl-1+=Branchi.Bn;BBranchi.Nr-1Branchi.Nr-1+=Branchi.Bn;else

15、BBranchi.Nl-1Branchi.Nl-1-=(1-Branchi.Kt)/Branchi.Kt/Branchi.Kt/Br anchi.X;BBranchi.Nr-1Branchi.Nr-1-=(Branchi.Kt-1)/Branchi.Kt/Branchi.X; for(i=0;i<n;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(siz

16、eof(double)*2*(n-1);for(i=0;i<n-1;i+) mid1i=0;mid2i=0;for(j=0;j<n;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(i<nPQ)dS2*i+1=Nodei.Qs-Nodei.Qd-(Nodei.f*mid1i-Nodei.e*mid2i); elsedS2*i+1=Usi*Usi-(Nodei.e

17、*Nodei.e+Nodei.f*Nodei.f);error=0;for(i=0;i<2*(n-1);i+)if(dSi<0&&error<-dSi)error=-dSi;else if(dSi>0&&error<dSi)error=dSi;void FormJacob()int i,j;Jacob=(double *)malloc(sizeof(double *)*2*(n-1); for(i=0;i<2*(n-1);i+)Jacobi=(double *)malloc(sizeof(double)*2*(n-1);for

18、(i=0;i<2*(n-1);i+)for(j=0;j<2*(n-1);j+)Jacobij=0;for(j=0;j<n-1;j+)for(i=0;i<n-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;i<nPQ;i+)if(i!=j)Jac

19、ob2*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;i<n-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

20、*)malloc(sizeof(double *)*2*(n-1);for(i=0;i<2*(n-1);i+)invJaci=(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)invJacij=0;elseinvJacij=1;for(i=0;i<2*(n-1);i+)for(j=0;j<2*(n-1);j+)if(i!=j)temp=Jacobji/Jacobii;for(k=0;k<2*(n-1);k+)Jacobjk-=Ja

21、cobik*temp; invJacjk-=invJacik*temp;for(i=0;i<2*(n-1);i+)if(Jacobii!=1)temp=Jacobii;for(j=0;j<2*(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;i<2*(n-1);i+)dfei=0;for(j=0;j<2*(n-1);j+)dfei-=invJacij*dSj;for

22、(i=0;i<n-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;j<n;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;i<n-1;i+)Nodei.Qs=Nodei.f*mid1i-Nodei.e*mid2i;for(j=0;j<n-1;j+)for(i=0;i<n-j-1;i+)if(Nodei.N>Nodei+1.N)tN=Nodei.N;Nodei.N=Nodei+1.N;Nodei+1.N=