1、.电力系统潮流计算注：这是一个基于 NR 法的潮流计算通用程序，仅提供了子程序，需要做些处理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享！/* 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是 * Newton_Raphson 法.* 程序中所用的变量说明如下:* N:网络节点总数. M:网络的 PQ 节点数.* L:网络的支路总数. N0:雅可比矩阵的行数.* N1:N0+1 K:打印开关.K=1,则打印;否则,不打印.* K1:子程序 PLSC 中判断输入电压的形式.K1=1,则为极座标形式.否则* 为直角坐标形式.* D:有功及无功功率误差的最大值.* G(I,
2、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) :第
3、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)节点
4、注入的无功功率.* 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
5、 *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;elsefp=file4; /* 输出到文件 */* 初始化矩阵 G,B */for(i=1;in?n:(io+4);fprintf(fp,“n“);for(j=io;jn?n:(io+4);fprintf(fp,“n“);for(j=io;j0 ?
6、 p0pos1 : -p0pos1;if(d0?q0pos1:-q0pos1;if(dm) /* PV 结点 */* 计算式(4-90)中的 Rii */jmf2(2*i-1,2*i-1,n0)=-2*epos1;/* 计算式(4-90)中的 Sii */jmf2(2*i-1,2*i,n0)=-2*fpos1;/* 计算式(4-89) */for(j=1;jn0?n0:(io+4);fprintf(fp,“n“);for(j=io;j fabs(af2(l,i,n1) )l=j; /* 找到这行中的最大元 */if(l!=i) /* 行列交换 */for (j=i;jn1?n1:(io+4);
7、fprintf(fp,“n“);fprintf(fp,“ “);for(i=io;i(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“, af2(i,n1,n1);#define Pi 3.1415927/180void plsc(int n,int l,int m,float g,float b,float e,float f,int e1,int s1
8、,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;elsefp=file4;fprintf(fp,“nTHE
9、 RESULT ARE:“);if(k1=1)for(i=0;in;i+)anglei*=Pi;ei=vi*cos(anglei);fi=vi*sin(anglei);t1=0.0;t2=0.0;for(i=1;i=n;i+)pos1=f1(i);pos2=f2(n,i,n);t1+=gpos2*epos1-bpos2*fpos1;t2+=gpos2*fpos1+bpos2*epos1;pos1=f1(n);ppos1=t1*epos1;qpos1=-t2*epos1;m1=m+1;for(i1=m1;i1=ns;i1+)t1=0;t2=0;for(i=1;i=n;i+)pos1=f1(i);
10、pos2=f2(i1,i,n);t1+=gpos2*epos1-bpos2*fpos1;t2+=gpos2*fpos1+bpos2*epos1;pos1=f1(i1);qpos1=fpos1*t1-epos1*t2;.for(i=0;in; i+)cm=coi;if(cm!=0)qi-=(ei*ei+fi*fi)*cm;fprintf(fp,“nBUS DATA“);fprintf(fp,“nBUS VOLTAGE ANGLE(DEGS.) BUS P BUS Q“);for(i=0;in;i+)vi=sqrt(ei*ei+fi*fi);x=ei;y=fi;z=y/x;anglei=atan(
11、z);anglei/=Pi;fprintf(fp,“n%3d%13.5e%15.5f%15.5e%15.5e“,i+1,vi,anglei,pi,qi);fprintf(fp,“n LINE FLOW “);for(i=1;i=l;i+)pos1=f1(i);st=s1pos1;en=e1pos1;x1=ef1(st)*ef1(st)+ff1(st)*ff1(st);x2=ef1(en)*ef1(en)+ff1(en)*ff1(en);y1=ef1(st)*ef1(en)+ff1(st)*ff1(en);y2=ff1(st)*ef1(en)-ef1(st)*ff1(en);p1pos1=(x1
12、-y1)*g1pos1-y2*b1pos1;q1pos1=-x1*(c1pos1+b1pos1)+y1*b1pos1-y2*g1pos1;p2pos1=(x2-y1)*g1pos1+y2*b1pos1;q2pos1=-x2*(c1pos1+b1pos1)+y1*b1pos1+y2*g1pos1;for(j=1;j=n;j+)cm=cf2(j,i,l);if(cm!=0.0)km=1;if(en=j).km=2;if(km=1)q1pos1-=(ef1(j)*ef1(j)+ff1(j)*ff1(j)*cm;elseq2pos1-=(ef1(j)*ef1(j)+ff1(j)*ff1(j)*cm;p3pos1=p1pos1+p2pos1 ;q3pos1=q1pos1+q2pos1 ;fprintf(fp,“n%2d%8d%11d%13.6e%13.6e%13.6e%13.6e%17d%11d%13.6e%13.6e“,i,s1pos1,e1pos1,p1pos1,q1pos1,p3pos1,q3pos1,e1pos1,s1pos1,p2pos1,q2pos1);
