电力系统潮流计算 注:这是一个基于N-R法的潮流计算通用程序,仅提供了子程序,需要做些处理才能成为一个可运行的计算程序!此程序非我原创,仅与大家共享!!! /******************************************************************* * 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是 * * Newton_Raphson法. * * 程序中所用的变量说明如下: * * N:网络节点总数. M:网络的PQ节点数. * * L:网络的支路总数. N0:雅可比矩阵的行数. * * N1:N0+1 K:打印开关.K=1,则打印;否则,不打印.* * K1:子程序PLSC中判断输入电压的形式.K1=1,则为极座标形式.否则* * 为直角坐标形式. * * D:有功及无功功率误差的最大值. * * G(I,J):Ybus的电导元素(实部). * * B(I,J):Ybus的电纳元素(虚部). * * G1(I) :第I支路的串联电导. B1(I):第I支路的串联电纳. * * C1(I) :第I支路的pie型对称接地电纳. * * C(I,J):第I节点J支路不对称接地电纳. * * CO(I) :第I节点的接地电纳. * * S1(I) :第I节点的起始节点号. E1(I):第I节点的终止节点号. * * P(I) :第I节点的注入有功功率. Q(I):第I节点的注入无功功率.* * P0(I) :第I节点有功功率误差. Q0(I):第I节点无功功率误差. * * V0(I) :第I节点(PV节点)的电压误差(平方误差). * * V(I) :第I节点的电压误差幅值. * * E(I) :第I节点的电压的实部. F(I):第I节点的电压的虚部. * * JM(I,J):Jacoby矩阵的第I行J列元素. * * A(I,J):修正方程的增广矩阵,三角化矩阵的第I行J列元素,运算结 * * 束后A矩阵的最后一列存放修正的解. * * P1(I) :第I支路由S1(I)节点注入的有功功率. * * Q1(I) :第I支路由S1(I)节点注入的无功功率. * * P2(I) :第I支路由E1(I)节点注入的有功功率. * * Q2(I) :第I支路由E1(I)节点注入的无功功率. * * P3(I) :第I支路的有功功率损耗. * * Q3(I) :第I支路的无功功率损耗. * * ANGLE(I):第I节点电压的角度. * *******************************************************************/ #include #include #define f1(i) (i-1) /* 把习惯的一阶矩阵的下标转化为C语言数组下标*/ #define f2(i,j,n) ((i-1)*(n)+j-1) /* 把习惯的二阶矩阵的下标转化为C语言数组下标*/ /**************************************************** * 本子程序根据所给的支路导纳及有关信息,形成结点 * * 导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B * ****************************************************/ void ybus(int n,int l,int m,float *g,float *b,float *g1,float *b1,float *c1,\ float *c,float *co,int k,int *s1,int *e1) { extern FILE *file4; FILE *fp; int i,j,io,i0; int pos1,pos2; int st,en; if(file4==NULL) { fp=stdout; } else { fp=file4; /* 输出到文件 */ } /* 初始化矩阵G,B */ for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { pos2=f2(i,j,n); g[pos2]=0;b[pos2]=0; } } /* 计算支路导纳 */ for(i=1;i<=l;i++) { /* 计算对角元 */ pos1=f1(i); st=s1[pos1];en=e1[pos1]; pos2=f2(st,st,n); g[pos2]+=g1[pos1]; b[pos2]+=b1[pos1]+c1[pos1]; pos2=f2(en,en,n); g[pos2]+=g1[pos1]; b[pos2]+=b1[pos1]+c1[pos1]; /* 计算非对角元 */ pos2=f2(st,en,n); g[pos2]-=g1[pos1]; b[pos2]-=b1[pos1]; g[f2(en,st,n)]=g[f2(st,en,n)]; b[f2(en,st,n)]=b[f2(st,en,n)]; } /* 计算接地支路导纳 */ for(i=1;i<=n;i++) { /* 对称部分 */ b[f2(i,i,n)]+=co[f1(i)]; /* 非对称部分 */ for(j=1;j<=l;j++) { b[f2(i,i,n)]+=c[f2(i,j,l)]; } } if(k!=1) { return; /* 如果K不为 1,则返回;否则,打印导纳矩阵 */ } fprintf(fp,"\n BUS ADMITTANCE MATRIX Y(BUS):"); fprintf(fp,"\n ******************* ARRAY G ********************"); for(io=1;io<=n;io+=5) { i0=(io+4)>n?n:(io+4); fprintf(fp,"\n"); for(j=io;j<=i0;j++) { fprintf(fp,"%13d",j); } for(i=1;i<=n;i++) { fprintf(fp,"\n%2d",i); for(j=io;j<=i0;j++) { fprintf(fp,"%13.6f",g[f2(i,j,n)]); } } fprintf(fp,"\n"); } fprintf(fp,"\n ******************* ARRAY B ********************"); for(io=1;io<=n;io+=5) { i0=(io+4)>n?n:(io+4); fprintf(fp,"\n"); for(j=io;j<=i0;j++) { fprintf(fp,"%13d",j); } for(i=1;i<=n;i++) { fprintf(fp,"\n%2d",i); for(j=io;j<=i0;j++) { fprintf(fp,"%13.6f",b[f2(i,j,n)]); } } fprintf(fp,"\n"); } fprintf(fp,"\n************************************************"); } /******************************************* * 本子程序根据所给的功率及电压等数据 * * 求出功率及电压误差量,并返回最大有功功率 * * 以用于与给定误差比较.如打印参数K=1,则输 * * 出P0,Q0(对PQ结点),V0(对PV结点). * *******************************************/ void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,int m,\ int n,float *e,float *f,int k,float *g,float *b,float *dd) { extern FILE *file4; FILE *fp; int i,j,l; int pos1,pos2; float a1,a2,d1,d; if(file4==NULL) { fp=stdout; /* 输出到屏幕 */ } else { fp=file4; /* 输出到文件 */ } l=n-1; if(k==1) { fprintf(fp,"\n CHANGE OF P0,V**2,P0(I),Q0(I),V0(I) "); fprintf(fp,"\n I P0(I) Q0(I)"); } for(i=1;i<=l;i++) { a1=0;a2=0; pos1=f1(i); for(j=1;j<=n;j++) { /* a1,a2对应课本p185式(4-86)中括号内的式子 */ pos2=f2(i,j,n); a1+=g[pos2]*e[f1(j)]-b[pos2]*f[f1(j)]; a2+=g[pos2]*f[f1(j)]+b[pos2]*e[f1(j)]; } /* 计算式(4-86)(4-87)中的deltaPi */ p0[pos1]=p[pos1]-e[pos1]*a1-f[pos1]*a2; if(i <= m) { /* 计算PQ结点中的deltaQi */ q0[pos1]=q[pos1]-f[pos1]*a1+e[pos1]*a2; } else { /* 计算PV结点中的deltaVi平方 */ v0[pos1]=v[pos1]*v[pos1]-e[pos1]*e[pos1]-f[pos1]*f[pos1]; } /* 输出结果 */ if(k==1) { if(i fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]); } else if(i==m) { fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],q0[pos1]); fprintf(fp,"\n I P0(I) V0(I)"); } else { fprintf(fp,"\n %2d %15.6e %15.6e",i,p0[pos1],v0[pos1]); } } } /* 找到deltaP和deltaQ中的最小者,作为收敛指标, 存在dd中 */ d=0; for(i=1;i<=l;i++) { pos1=f1(i); d1=p0[pos1]>0 ? p0[pos1] : -p0[pos1]; if(d d=d1; } if(i<=m) { d1=q0[pos1]>0?q0[pos1]:-q0[pos1]; if(d d=d1; } } } (*dd)=d; } /*************************************************** * 本子程序根据节点导纳及电压求Jacoby矩阵,用于求* * 电压修正量,如打印参数K=1,则输出Jacoby矩阵. * * 对应于课本P186式(4-89)(4-90) * ***************************************************/ void jmcc(int m,int n,int n0,float *e,float *f,float *g,float *b,float *jm,int k) { extern FILE *file4; FILE *fp; int i,j,i1,io,i0,ns; int pos1,pos2; if(file4==NULL) { fp=stdout; } else { fp=file4; } /* 初始化矩阵jm */ for(i=1;i<=n0;i++) { for(j=1;j<=n0;j++) { jm[f2(i,j,n0)]=0; } } ns=n-1; /* 去掉一个平衡结点 */ /* 计算式(4-89)(4-90) */ for(i=1;i<=ns;i++) { /* 计算式(4-90) */ for(i1=1;i1<=n;i1++) { /* pos1是式(4-90)中的j */ pos1=f1(i1); /* pos2是式(4-90)中的ij */ pos2=f2(i,i1,n); if(i<=m) /* i是PQ结点 */ { /* 计算式(4-90)中的Jii等式右侧第一部分 */ jm[f2(2*i-1,2*i-1,n0)]+=g[pos2]*f[pos1]+b[pos2]*e[pos1]; /* 计算式(4-90)中的Lii等式右侧第一部分 */ jm[f2(2*i-1,2*i,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1]; } /* 计算式(4-90)中的Hii等式右侧第一部分 */ jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]+b[pos2]*f[pos1]; /* 计算式(4-90)中的Nii等式右侧第一部分 */ jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]-b[pos2]*e[pos1]; } /* pos2是式(4-90)中的ii */ pos2=f2(i,i,n); /* pos1是式(4-90)中的i */ pos1=f1(i); if(i<=m) /* i是PQ结点 */ { /* 计算式(4-90)中的Jii */ jm[f2(2*i-1,2*i-1,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1]; /* 计算式(4-90)中的Lii */ jm[f2(2*i-1,2*i,n0)]+=g[pos2]*e[pos1]+b[pos2]*f[pos1]; } /* 计算式(4-90)中的Hii */ jm[f2(2*i,2*i-1,n0)]+=-g[pos2]*e[pos1]-b[pos2]*f[pos1]; /* 计算式(4-90)中的Jii */ jm[f2(2*i,2*i,n0)]+=-g[pos2]*f[pos1]+b[pos2]*e[pos1]; if(i>m) /* PV结点 */ { /* 计算式(4-90)中的Rii */ jm[f2(2*i-1,2*i-1,n0)]=-2*e[pos1]; /* 计算式(4-90)中的Sii */ jm[f2(2*i-1,2*i,n0)]=-2*f[pos1]; } /* 计算式(4-89) */ for(j=1;j<=ns;j++) { if(j!=i) { /* pos1是式(4-89)中的i */ pos1=f1(i); /* pos2是式(4-89)中的ij */ pos2=f2(i,j,n); /* 计算式(4-89)中的Nij */ jm[f2(2*i,2*j,n0)]=b[pos2]*e[pos1]-g[pos2]*f[pos1]; /* 计算式(4-89)中的Hij */ jm[f2(2*i,2*j-1,n0)]=-g[pos2]*e[pos1]-b[pos2]*f[pos1]; if(i<=m) /* i是PQ结点 */ { /* 计算式(4-89)中的Lij (=-Hij) */ jm[f2(2*i-1,2*j,n0)]=-jm[f2(2*i,2*j-1,n0)]; /* 计算式(4-89)中的Jij (=Nij) */ jm[f2(2*i-1,2*j-1,n0)]=jm[f2(2*i,2*j,n0)]; } else /* i是PV结点 */ { /* 计算式(4-89)中的Rij (=0) */ jm[f2(2*i-1,2*j-1,n0)]=0; /* 计算式(4-89)中的Sij (=0) */ jm[f2(2*i-1,2*j,n0)]=0; } } } } if(k!=1) { return; } /* 输出Jacoby矩阵 */ fprintf(fp,"\n J MATRIX(C)"); for(io=1;io<=n0;io+=5) { i1=(io+4)>n0?n0:(io+4); fprintf(fp,"\n"); for(j=io;j<=i1;j++) { fprintf(fp,"%10d",j); } for(i=1;i<=n0;i++) { fprintf(fp,"\n%2d",i); for(j=io;j<=i1;j++) { fprintf(fp,"%12.6f",jm[f2(i,j,n0)]); } } } fprintf(fp,"\n"); } /********************************************** * 本子程序用选列主元素的高斯消元法求解组 * * 性方程组求各结点电压修正量,如打印参数K=1,则* * 输出增广矩阵变换中的上三角及电压修正量.如果* * 无唯一解,则给出信息,并停止程序运行. * **********************************************/ void sevc ( float a[], int n0, int k, int n1) { extern FILE *file4; FILE *fp; int i,j,l,n2,n3,n4,i0,io,j1,i1; float t0,t,c; if(file4==NULL) fp=stdout; else fp=file4; for(i=1;i<=n0;i++) { l=i; for(j=i;j<=n0;j++) { if( fabs(a[f2(j,i,n1)]) > fabs(a[f2(l,i,n1)]) ) { l=j; /* 找到这行中的最大元 */ } } if(l!=i) { /* 行列交换 */ for (j=i;j<=n1;j++) { t=a[f2(i,j,n1)]; a[f2(i,j,n1)]=a[f2(l,j,n1)]; a[f2(l,j,n1)]=t; } } if (fabs(a[f2(i,i,n1)]-0)<1e-10) { /* 对角元近似于0, 无解 */ printf("\nNo Solution\n"); exit (1); } t0=a[f2(i,i,n1)]; for(j=i;j<=n1;j++) { /* 除对角元 */ a[f2(i,j,n1)]/=t0; } if(i==n0) { /* 最后一行,不用消元 */ continue; } /* 消元 */ j1=i+1; for(i1=j1;i1<=n0;i1++) { c=a[f2(i1,i,n1)]; for(j=i;j<=n1;j++) { a[f2(i1,j,n1)] -= a[f2(i,j,n1)] *c; } } } if(k==1) { /* 输出上三角矩阵 */ fprintf(fp,"\nTrianglar Angmentex Matrix "); for(io=1;io<=n1;io+=5) { i0=(io+4)>n1?n1:(io+4); fprintf(fp,"\n"); fprintf(fp," "); for(i=io;i<=i0;i++) { fprintf(fp,"%12d",i); } for(i=1;i<=n0;i++) { fprintf(fp,"\n"); fprintf(fp,"%2d",i); for(j=io;j<=i0;j++) { fprintf(fp,"%15.6f", a[f2(i,j,n1)]); } } } } /* 回代求方程解 */ n2=n1-2; for(i=1;i<=n2;i++) { n3=n1-i; for(i1=n3;i1<=n0;i1++) { n4=n0-i; a[f2(n4,n1,n1)] -= a[f2(i1,n1,n1)]*a[f2(n4,i1,n1)]; } } if(k!=1) { return; } /* 输出电压修正值 */ fprintf(fp,"\nVoltage correction E(i), F(i) :"); for(io=1;io<=n0;io+=4) { i1=(io+1)/2; i0=((io+3)/2)>(n0/2)?(n0/2):((io+3)/2); fprintf(fp,"\n"); for(j=i1;j<=i0;j++) { fprintf(fp,"%16d%16d",j,j); } i1 = 2*i0; fprintf(fp,"\n"); for(i=io;i<=i1;i++) { fprintf(fp,"%15.6f", a[f2(i,n1,n1)]); } } } #define Pi 3.1415927/180 void plsc(int n,int l,int m,float g[],float b[],float e[],float f[],\ int e1[],int s1[],float g1[],float b1[],float c1[],float c[],\ float co[],float p1[],float q1[],float p2[],float q2[],float p3[],\ float q3[],float p[],float q[],float v[],float angle[],int k1) { extern FILE *file4; FILE *fp; float t1,t2,st,en,cm,x,y,z,x1,x2,y1,y2; int i,i1,j,m1,ns,pos1,pos2,km; ns=n-1; if(file4==NULL) { fp=stdout; } else { fp=file4; } fprintf(fp,"\nTHE RESULT ARE:"); if(k1==1) { for(i=0;i angle[i]*=Pi; e[i]=v[i]*cos(angle[i]); f[i]=v[i]*sin(angle[i]); } } t1=0.0;t2=0.0; for(i=1;i<=n;i++) { pos1=f1(i);pos2=f2(n,i,n); t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1]; t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1]; } pos1=f1(n); p[pos1]=t1*e[pos1]; q[pos1]=-t2*e[pos1]; m1=m+1; for(i1=m1;i1<=ns;i1++) { t1=0;t2=0; for(i=1;i<=n;i++) { pos1=f1(i);pos2=f2(i1,i,n); t1+=g[pos2]*e[pos1]-b[pos2]*f[pos1]; t2+=g[pos2]*f[pos1]+b[pos2]*e[pos1]; } pos1=f1(i1); q[pos1]=f[pos1]*t1-e[pos1]*t2; } for(i=0;i cm=co[i]; if(cm!=0) { q[i]-=(e[i]*e[i]+f[i]*f[i])*cm; } } fprintf(fp,"\nBUS DATA"); fprintf(fp,"\nBUS VOLTAGE ANGLE(DEGS.) BUS P BUS Q"); for(i=0;i v[i]=sqrt(e[i]*e[i]+f[i]*f[i]); x=e[i]; y=f[i]; z=y/x; angle[i]=atan(z); angle[i]/=Pi; fprintf(fp,"\n%3d%13.5e%15.5f%15.5e%15.5e",i+1,v[i],angle[i],p[i],q[i]); } fprintf(fp,"\n LINE FLOW "); for(i=1;i<=l;i++) { pos1=f1(i); st=s1[pos1]; en=e1[pos1]; x1=e[f1(st)]*e[f1(st)]+f[f1(st)]*f[f1(st)]; x2=e[f1(en)]*e[f1(en)]+f[f1(en)]*f[f1(en)]; y1=e[f1(st)]*e[f1(en)]+f[f1(st)]*f[f1(en)]; y2=f[f1(st)]*e[f1(en)]-e[f1(st)]*f[f1(en)]; p1[pos1]=(x1-y1)*g1[pos1]-y2*b1[pos1]; q1[pos1]=-x1*(c1[pos1]+b1[pos1])+y1*b1[pos1]-y2*g1[pos1]; p2[pos1]=(x2-y1)*g1[pos1]+y2*b1[pos1]; q2[pos1]=-x2*(c1[pos1]+b1[pos1])+y1*b1[pos1]+y2*g1[pos1]; for(j=1;j<=n;j++) { cm=c[f2(j,i,l)]; if(cm!=0.0) { km=1; if(en==j) { km=2; } if(km==1) { q1[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm; } else { q2[pos1]-=(e[f1(j)]*e[f1(j)]+f[f1(j)]*f[f1(j)])*cm; } } } p3[pos1]=p1[pos1]+p2[pos1] ; q3[pos1]=q1[pos1]+q2[pos1] ; fprintf(fp,"\n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d%11d%13.6e%13.6e",\ i,s1[pos1],e1[pos1],p1[pos1],q1[pos1],p3[pos1],q3[pos1],\ e1[pos1],s1[pos1],p2[pos1],q2[pos1]); } } |
|小黑屋|最新主题|手机版|微赢网络技术论坛 ( 苏ICP备08020429号 )
GMT+8, 2024-9-30 03:33 , Processed in 0.235917 second(s), 12 queries , Gzip On, MemCache On.
Powered by Discuz! X3.5
© 2001-2023 Discuz! Team.