电力系统分析潮流计算程序.docx

1、电力系统分析潮流计算程序#include#include#include#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 BranchType int Nbr;/支路号 int Nl;/首节点 int Nr;/末节点 double R;/支路电阻 dou

2、ble X;/支路电抗 double Bn;/对地电抗 double Kt;/支路变比 ;/*全局变量声明*int n;/节点数int nPQ;/PQ节点数int nPV;/PV节点数int nbr;/支路数int ng;/发电机台数int Mark=0;/标记支路参数是否已经转换double *G;/导纳矩阵G部分double *B;/导纳矩阵B部分double *dS;/功率不平衡量double *mid1,*mid2;/求功率不平衡量时的中间变量double *Us;/电压初值double error=1;/误差值double iteration=0.000001;/误差精度double

3、 *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();

4、void Print1(double *,int); void Print2(double *,int,int); 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);/*void LoadData() int i,j; int

5、tN,tType; double te,tf,tPd,tQd,tPs,tQs,tBc;/用于重新排列节点信息的临时变量 FILE *fp;/文件指针 char filename50=; printf(请输入数据文件名：); 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(st

6、ruct 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); /计算PQ节点和PV节点的个数 for(i=0;in;i+) if(Nodei.Type=1) nPQ+; else if(Nodei.Type=2) nPV+; printf(PQ节点个数：%dn,nPQ); prin

7、tf(PV节点个数：%dn,nPV); /重新排列节点参数(冒泡法) 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;No

8、dei+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; /读取支路参数

9、 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+)

10、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;/存储Z2=R2+X2 G=(double *)malloc(sizeof(double *)*n);/为G申请空间 B=(double *)mall

11、oc(sizeof(double *)*n);/为B申请空间 for(i=0;in;i+) Gi=(double *)malloc(sizeof(double)*n); Bi=(double *)malloc(sizeof(double)*n); /初始化G、B 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-=Bra

12、nchi.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; else/变压器支路 GBranchi.Nl-1Branchi.Nr-1-=Branchi.R/Z2/Branchi.Kt; BBranchi.Nl-1Branchi.Nr-1+=Branchi.X/Z2/Branchi.Kt; GBranchi.Nr-1Branchi.Nl-1=G

13、Branchi.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|

14、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; else/变压器支路 BBranchi.Nl-1Branchi.Nl-1-=(1-Branchi.Kt)/Branchi.Kt/Branchi.Kt/Bra

15、nchi.X; BBranchi.Nr-1Branchi.Nr-1-=(Branchi.Kt-1)/Branchi.Kt/Branchi.X; /将并联电容加入到对角元素中 for(i=0;in;i+) Bii=Bii+Nodei.Bc;/*/*求deltaP,deltaQ*void DeltaS()/计算功率不平衡量 int i,j; /为中间变量申请空间 mid1=(double *)malloc(sizeof(double)*n); mid2=(double *)malloc(sizeof(double)*n); /为功率不平衡量申请空间 dS=(double *)malloc(size

16、of(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); else dS2*i+1=Usi*

17、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);

18、i+) for(j=0;j2*(n-1);j+) Jacobij=0; for(j=0;jn-1;j+) /求H,N 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; else Jacob2*i2*i=Bii*Nodei.e-Gii*Nodei.f-mid2i; Jacob2*i2*i+1=-Gii*Nodei.e-Bii*Nodei.f-mid1i; /求J,L for(i=0;inPQ;i+) if(i!=j) Jacob2*i+12*

19、j=Gij*Nodei.e+Bij*Nodei.f; Jacob2*i+12*j+1=Bij*Nodei.e-Gij*Nodei.f; else Jacob2*i+12*i=Gii*Nodei.e+Bii*Nodei.f-mid1i; Jacob2*i+12*i+1=Bii*Nodei.e-Gii*Nodei.f+mid2i; /求R,S 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

20、;/中间变量 /为雅克比矩阵的逆申请空间 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; else invJacij=1; for(i=0;i2*(n-1);i+) for(j=0;j2*(n-1);j+) if(i!=j) temp=Jacobji/Jacobii; fo

21、r(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

22、; 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; /计算PV节点的Q