电力系统潮流计算C程序.doc

电力系统潮流计算注：这是一个基于NR法的潮流计算通用程序，仅提供了子程序，需要做些处理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享!！/* 这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是 * 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 math.h#include stdio.h#define f1（i) （i1）/ 把习惯的一阶矩阵的下标转化为C语言数组下标/define f2（i，j,n） （(i1)（n)+j1)/ 把习惯的二阶矩阵的下标转化为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); gpos2=0；bpos2=0; / 计算支路导纳 */for(i=1;i=l;i+） / 计算对角元 */ pos1=f1（i）; st=s1pos1；en=e1pos1； pos2=f2(st，st,n）； gpos2+=g1pos1; bpos2+=b1pos1+c1pos1; pos2=f2（en,en，n); gpos2+=g1pos1; bpos2+=b1pos1+c1pos1; /* 计算非对角元 */ pos2=f2（st,en，n); gpos2=g1pos1; bpos2-=b1pos1； gf2（en,st,n)=gf2(st，en,n)； bf2(en,st,n)=bf2（st,en，n)；/ 计算接地支路导纳 */for(i=1;in？n：(io+4）； fprintf（fp,”n”）； for（j=io;j=i0;j+） fprintf(fp，13d，j）； for（i=1;i=n；i+） fprintf(fp,n2d”,i); for(j=io；jn？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”,bf2(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=n1;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式（486）中括号内的式子 */ pos2=f2(i，j，n）； a1+=gpos2ef1(j）bpos2ff1(j）; a2+=gpos2*ff1(j)+bpos2*ef1(j)； /* 计算式(4-86)（487）中的deltaPi*/ p0pos1=ppos1epos1*a1-fpos1*a2； if（i = m) / 计算PQ结点中的deltaQi/ q0pos1=qpos1fpos1*a1+epos1a2; else / 计算PV结点中的deltaVi平方*/ v0pos1=vpos1*vpos1-epos1*epos1fpos1fpos1; / 输出结果 / if（k=1） if(im) fprintf（fp,n %2d 15.6e %15.6e”，i，p0pos1,q0pos1）； else if（i=m） fprintf(fp,n %2d 15.6e 15。6e,i，p0pos1,q0pos1)； fprintf(fp,n I P0(I) V0(I); else fprintf（fp,n %2d 15.6e %15。6e，i，p0pos1,v0pos1）; /* 找到deltaP和deltaQ中的最小者，作为收敛指标， 存在dd中 /d=0;for（i=1；i0 ? p0pos1 : p0pos1； if(dd1) d=d1； if（i=m) d1=q0pos10?q0pos1:-q0pos1； if(dd1) 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+) jmf2(i,j,n0）=0； ns=n-1; / 去掉一个平衡结点 /* 计算式(489）（4-90） */for(i=1；i=ns;i+) / 计算式(4-90） / for（i1=1;i1=n；i1+） / pos1是式（490)中的j / pos1=f1（i1）; /* pos2是式（490)中的ij */ pos2=f2（i,i1,n）； if（im) /* PV结点 */ / 计算式(4-90)中的Rii */ jmf2(2*i1,2i-1,n0）=2epos1； / 计算式(490）中的Sii / jmf2（2i1,2*i，n0)=-2*fpos1; /* 计算式（4-89） */ for（j=1;j=ns；j+) if（j!=i） /* pos1是式（489）中的i / pos1=f1（i)； /* pos2是式（489)中的ij / pos2=f2（i,j，n); /* 计算式（4-89)中的Nij */ jmf2（2i,2*j，n0）=bpos2epos1gpos2*fpos1; /* 计算式(489）中的Hij / jmf2(2*i，2*j1，n0)=gpos2epos1-bpos2*fpos1; if（i=m) / i是PQ结点 */ / 计算式(489）中的Lij (=-Hij) / jmf2(2*i1，2*j，n0）=-jmf2(2i，2j1，n0）； /* 计算式（4-89）中的Jij （=Nij) */ jmf2（2i1,2*j1,n0)=jmf2（2i,2*j，n0)； else / i是PV结点 */ / 计算式(489）中的Rij (=0） */ jmf2（2i-1，2*j-1,n0）=0; /* 计算式（489）中的Sij (=0) */ jmf2（2*i1,2*j，n0)=0; if（k！=1） return；/ 输出Jacoby矩阵 */fprintf（fp,”n J MATRIX(C)）;for（io=1；ion0？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,jmf2(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（af2(j,i，n1)） fabs(af2(l,i,n1）) ） l=j; /* 找到这行中的最大元 */ if（l！=i） / 行列交换 */ for （j=i;j=n1;j+） t=af2(i,j,n1)； af2(i,j,n1)=af2（l，j,n1）； af2(l，j,n1)=t； if (fabs（af2(i,i，n1)0）1e-10） / 对角元近似于0, 无解 */ printf（”nNo Solutionn）； exit （1)； t0=af2(i，i，n1）; for(j=i；j=n1；j+） / 除对角元 / af2(i,j，n1）/=t0; if(i=n0） /* 最后一行，不用消元 / continue； / 消元 */ j1=i+1; for(i1=j1；i1=n0；i1+） c=af2(i1，i,n1）； for（j=i；j=n1;j+） af2（i1，j，n1） -= af2（i,j,n1） c； if(k=1) /* 输出上三角矩阵 / fprintf(fp，nTrianglar Angmentex Matrix )； for(io=1；ion1？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, af2（i,j，n1）)； / 回代求方程解 */n2=n1-2;for（i=1;i=n2;i+） n3=n1i； for（i1=n3;i1=n0；i1+) n4=n0i； af2(n4,n1,n1) = af2（i1,n1，n1)*af2（n4，i1,n1); if(k!=1) return;/ 输出电压修正值 */fprintf(fp,nVoltage correction E（i）, F（i) ：）；for（io=1;io（n0/2）?（n0/2）:(（io+3）/2)； fprintf(fp,”n”）; for（j=i1;j=i0;j+) fprintf（fp,”%16d%16d,j,j）; i1 = 2i0; 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，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=n1；if（file4=NULL) fp=stdout；else fp=file4;fprintf（fp，nTHE RESULT ARE：)；if（k1=1） for(i=0；in;i+） anglei*=Pi; ei=vi*cos(anglei）; fi=visin(anglei)； t1=0.0；t2=0.0；for（i=1；i=n；i+） pos1=f1(i）;pos2=f2(n，i，n）; t1+=gpos2epos1bpos2*fpos1; t2+=gpos2*fpos1+bpos2epos1;pos1=f1（n）；ppos1=t1epos1;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)；pos2=f2(i1，i，n)； t1+=gpos2epos1bpos2fpos1; t2+=gpos2fpos1+bpos2*epos1; pos1=f1(i1); qpos1=fpos1t1epos1*t2；for（i=0；in； i+） cm=coi； if(cm！=0) qi-=(eiei+fifi）*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（z)； anglei/=Pi； fprintf(fp，n%3d13.5e15.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-y1)g1pos1-y2*b1pos1; q1pos1=-x1*(c1pos1+b1pos1)+y1*b1pos1y2*g1pos1; p2pos1=(x2y1)g1pos1+y2b1pos1； q2pos1=-x2（c1pos1+b1pos1)+y1*b1pos1+y2g1pos1; 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； else q2pos1-=(ef1（j）ef1(j)+ff1(j）*ff1(j）