电力系统潮流计算机算法牛顿.doc

电力系统潮流计算机算法牛顿拉夫逊法 /* -设计思想：- 电力系统潮流计算机算法,解非线性方程采用的方法是牛顿拉夫逊法, 修正方程以直角坐标表示，解修正方程采用的是列主元素消去法 -程序功能- 程序可计算带PQ节点、PV节点、平行支路、接地支路、非标准变比变压器支路 的电力系统潮流，通用性良好*/#include #include /*主要变量说明如下*/#define M 50 /*矩阵阶数*/#define N 50 /*迭代次数*/int i,j,k,l,z; /* 循环变量 */int n, /* 节点数 */ m, /* 支路数 */ byqn, /* 变压器数 */ dd, /*对地支路数*/ pq, /* PQ节点数 */ pv; /* PV节点数 */float eps, /* 精度 */ max;char filename10;static float GMM,BMM,DM,ykbMM;struct jiedian /* 节点结构体 */ int num,s; /* num为节点号，s为节点类型*/ float p,q,e,f,v; jiedianM;struct zhilu /* 支路结构体 */ int num; int p1,p2; /*支路的两个节点*/ float r,x; /*支路的电阻与电抗*/ zhiluM; struct byq /*变压器结构体*/ int num; int p1,p2;float r,x,k;float s; byqM; FILE *fp1,*fp2; /*input的输入要求 input.c的输入要求： 节点数n, 支路数m, 对地支路数dd, 变压器数byqn，PQ节点数pq, PV节点数pv, 精度eps 节点编号numb,节点类型h,节点输入功率（实数部分），节点输入功率（虚数部分），节电电压初值（实数部分）e，节电电压初值（虚数部分）f 平衡节点编号，节点类型，节电电压初值（实数部分）e，节电电压初值（虚数部分）f (继续输入节点信息） 支路的编号,支路一端节点的编号,支路另一端节点的编号,支路等效阻抗的实数部分,支路等效阻抗的虚数部分 (继续输入支路信息） 变压器编号num，所连节点p1，p2，阻抗r，x，变比k，容量s */void data1() /* 读取数据 */ int h,numb,queren; printf(程序可计算带PQ节点、PV节点、平行支路、接地支路、非标准变比变压器支路的系统潮流n); printf(输入必须符合程序中input的输入要求（如下），若不符合，可能导致错误n); printf( *输入要求*n); printf(-n); printf(节点数n,支路数m, 对地支路数dd, 变压器数byqn，PQ节点数pq, PV节点数pv, 精度epsn); printf(节点编号,节点类型h, 输入功率p, q， 节电电压初值e，fn); printf(（平衡节点不必输入p,q）n); printf(继续输入节点信息）n); printf(支路编号,支路一端节点的编号p1,另一端节点的编号p2,支路等效阻抗r,xn); printf(n); printf(继续输入支路信息）n); printf(变压器编号，所连节点p1，p2，阻抗r，x，变比k，容量sn); printf(-n); printf(n *本程序典型的 测试 *n); printf(测试带PQ节点系统（我的原题目），输入文件名为t1.c, 输出文件名为j1.c,n); printf( 测试带pq、Pv节点系统，输入文件名为t2.c, 输出文件名为j2.c,n); printf( 测试带对地支路系统， 输入文件名为t3.c, 输出文件名为j3.c,n); printf( 测试带平行支路系统， 输入文件名为t4.c, 输出文件名为j4.c,n); printf( 测试带变压器支路系统，输入文件名为t5.c, 输出文件名为j5.c,nnn); printf( *若要 自行 输入其他系统计算，请确认输入符合输入要求*n); printf( 输入文件名为input.c,输出文件名为output.c,n); printf(请键入 输入文件名n); scanf(%s,&filename); if(fp1=fopen(filename,r)=NULL) printf( can not open the %s file,please check it!n,filename); exit(0); fscanf(fp1,%d,%d,%d,%d,%d,%d,%fn,&n,&m,&dd,&byqn,&pq,&pv,&eps); /*输入节点数n,支路数m,对地支路数dd,变压器数byqn，PQ节点数pq,PV节点数pv,精度eps*/ j=1;k=pq+1; for(i=1;i=n;i+) fscanf(fp1,%d,%d,&numb,&h); /*输入节点编号numb,节点类型h*/if(h=1) /*类型h=1是PQ节点*/ fscanf(fp1,%f,%f,%f,%fn, &jiedianj.p, /*输入节点的输入功率的实数部分*/ &jiedianj.q, /*输入节点的输入功率的虚数部分*/ &jiedianj.e, /*输入节电电压初值实数部分*/ &jiedianj.f); /*输入节电电压初值虚数部分*/ jiedianj.num=numb; /*赋予节点编号numb,节点类型h*/ jiedianj.s=h; j+; if(h=2) /*类型h=2是pv节点*/ fscanf(fp1,%f,%f,%f,%fn, &jiediank.p, /*输入节点的输入功率的实数部分*/ &jiediank.v, /*输入节点的输入功率的虚数部分*/ &jiediank.e, /*输入节电电压初值实数部分*/ &jiediank.f); /*输入节电电压初值虚数部分*/ jiediank.num=numb; /*赋予节点编号numb,节点类型h*/ jiediank.s=h; k+; if(h=3) /*类型h=3是平衡节点*/ fscanf(fp1,%f,%fn,&jiediann.e,&jiediann.f);/*输入节电电压初值jiediann.e,jiediann.f*/ jiediann.num=numb; /*赋予节点编号numb,节点类型h*/ jiediann.s=h; for(i=1;i=m;i+) /*输入支路阻抗，最大支路数为m*/ fscanf(fp1,%d,%d,%d,%f,%fn, &zhilui.num, /*输入支路的编号*/ &zhilui.p1, /*输入支路一端节点的编号*/ &zhilui.p2, /*输入支路另一端节点的编号*/ &zhilui.r, /*输入支路等效阻抗的实数部分*/ &zhilui.x); /*输入支路等效阻抗的虚数部分*/ for(i=1;i=byqn;i+) fscanf(fp1,%d,%d,%d,%f,%f,%f,%fn,&byqi.num,&byqi.p1,&byqi.p2,&byqi.r,&byqi.x,&byqi.k,&byqi.s); /*变压器编号num，所连节点p1，p2，阻抗r，x，变比k，容量s*/ fclose(fp1);printf(请键入 输出文件名n); scanf(%s,&filename); if(fp2=fopen(filename,w)=NULL) /* 判断输出文件是否可以打开 */ printf( can not write the %s file,please check !n,filename); exit(0); fprintf(fp2,n 电自042班 刘振宗 学号200401030210n); fprintf(fp2,n 原始数据 n); /*显示原始数据*/ fprintf(fp2, 节点数:%2d 支路数:%2d 对地支路数:%2d 变压器数:%2d PQ节点数:%2d PV节点数:%2d 精度:%fn, n,m,dd,byqn,pq,pv,eps); for(i=1;i=pq;i+) fprintf(fp2, PQ节点 节点%2d P%d=%8.5f Q%d=%8.5fn, jiediani.num,jiediani.num,jiediani.p,jiediani.num,jiediani.q); for(i=pq+1;i=pq+pv;i+) fprintf(fp2, PV节点 节点%2d P%d=%8.5f V%d=%8.5fn, jiediani.num,jiediani.num,jiediani.p,jiediani.num,jiediani.v); fprintf(fp2, 平衡节点 节点%2d e%d=%8.5f f%d=%8.5fn, jiediann.num,jiediann.num,jiediann.e,jiediann.num,jiediann.f); fprintf(fp2,-n); for(i=1;i=m;i+) fprintf(fp2, 支路%2d 相关节点:%2d,%2d R=8.5f, X=%8.5fn, i,zhilui.p1,zhilui.p2,zhilui.r,zhilui.x); for(i=1;i=byqn;i+) fprintf(fp2,变压器编号%2d 相关节点:%2d,%2d R=%8.5f X=%8.5f 变比K=%8.5fn,byqi.num, byqi.p1,byqi.p2,byqi.r,byqi.x,byqi.k); fprintf(fp2,n);void form_y() /* 形成节点导纳矩阵 */ float S; /*分母有理化中的分母*/ float g,b,g1,b1,g2,b2; int bn; /*计算节点导纳矩阵的主对角线上的导纳*/ for(i=1;i=m;i+) /* m, 支路数 */ for(j=0;j=n;j+) /*n, 节点数 */ if(zhilui.p1=j)|(zhilui.p2=j) /*zhilui.p1为支路一端节点的编号,zhilui.p2为支路另一端节点的编号*/ S=zhilui.r*zhilui.r+zhilui.x*zhilui.x; /*分母有理化*/ if(S=0) continue; Gjj+=zhilui.r/S; /*得到导纳实数部分 */ Bjj+=-zhilui.x/S; /*得到导纳虚数部分 */ for(bn=1;bn=byqn;bn+) if(byqbn.p1=j) S=byqbn.r*byqbn.r+byqbn.x*byqbn.x; g=byqbn.r/S; b=-byqbn.x/S; g2=g/byqbn.k; b2=b/byqbn.k; g1=g*(1-byqbn.k)/(byqbn.k*byqbn.k); b1=b*(1-byqbn.k)/(byqbn.k*byqbn.k); Gjj=Gjj+g2+g1; Bjj=Bjj+b2+b1; if(byqbn.p2=j) S=byqbn.r*byqbn.r+byqbn.x*byqbn.x; g=byqbn.r/S; b=-byqbn.x/S; g2=g/byqbn.k; b2=b/byqbn.k; g1=g*(byqbn.k-1)/byqbn.k; b1=b*(byqbn.k-1)/byqbn.k; Gjj=Gjj+g2+g1; Bjj=Bjj+b2+b1; for(j=1;j=m;j+) /*对地导纳*/ if(zhiluj.p1=0)|(zhiluj.p2=0) S=zhiluj.r*zhiluj.r+zhiluj.x*zhiluj.x; if(S=0) Gzhiluj.p1zhiluj.p2=Bzhiluj.p1zhiluj.p2=Gzhiluj.p2zhiluj.p1=Bzhiluj.p2zhiluj.p1=0;continue; Gzhiluj.p1zhiluj.p2+=zhilui.r/S; Bzhiluj.p1zhiluj.p2+=-zhilui.x/S; Gzhiluj.p2zhiluj.p1=Gzhiluj.p1zhiluj.p2; Bzhiluj.p2zhiluj.p1=Bzhiluj.p1zhiluj.p2; for(k=1;k=m;k+) /*节点导纳矩阵非主对角线上的导纳*/ i=zhiluk.p1;j=zhiluk.p2;S=zhiluk.r*zhiluk.r+zhiluk.x*zhiluk.x;if(S=0) continue;Gij+=-zhiluk.r/S; /*得到导纳实数部分 */Bij+=zhiluk.x/S; /*得到导纳虚数部分 */Gji=Gij; /*利用节点导纳矩阵的对称性*/Bji=Bij; for(bn=1;bn=byqn;bn+) i=byqbn.p1;j=byqbn.p2; S=byqbn.r*byqbn.r+byqbn.x*byqbn.x; Gij+=-(byqbn.r/S)/byqbn.k; Bij+=(byqbn.x/S)/byqbn.k; Gji=Gij; Bji=Bij; /*输出节点导纳矩阵*/fprintf(fp2,nnn 计算结果n); fprintf(fp2,n 节点导纳矩阵为:); for(i=1;i=n;i+) fprintf(fp2,n ); for(j=1;j=n;j+) fprintf(fp2,%8.5f+j%8.5f ,Gij,Bij); fprintf(fp2,n ); void form_ykb() /* 形成雅可比矩阵 */ float ei,ej,fi,fj,a=0,b=0; int i1,j1,k1; float dPM,dQM,dU2M; for(i=1;i=2*(pq+pv);i+) /*初始化矩阵*/ dPi=0;dQi=0;dU2i=0;Di=0;for(j=1;j=2*(pq+pv)+1;j+) ykbij=0; for(i=1;i=pq;i+) /*PQ节点功率不平衡量*/ ei=jiediani.e; fi=jiediani.f; i1=jiediani.num; for(j=1;j=n;j+) j1=jiedianj.num; ej=jiedianj.e; fj=jiedianj.f; dPi+=ei*(Gi1j1*ej-Bi1j1*fj) +fi*(Gi1j1*fj+Bi1j1*ej); dQi+=fi*(Gi1j1*ej-Bi1j1*fj) -ei*(Gi1j1*fj+Bi1j1*ej); dPi=jiediani.p-dPi; dQi=jiediani.q-dQi; for(i=pq+1;i=pq+pv;i+) /*PV节点有功不平衡量和电压不平衡量*/ ei=jiediani.e; fi=jiediani.f; for(j=1;j=n;j+) i1=jiediani.num; j1=jiedianj.num; ej=jiedianj.e; fj=jiedianj.f; dPi+=ei*(Gi1j1*ej-Bi1j1*fj) +fi*(Gi1j1*fj+Bi1j1*ej); dPi=jiediani.p-dPi; dU2i=jiediani.v*jiediani.v-ei*ei-fi*fi; for(i=1;i=pq;i+) /*形成不平衡量矩阵DM*/ D2*i-1=dPi; D2*i=dQi; /*pq节点部分，基数位为dP,偶数位为dQ*/ for(i=pq+1;i=pq+pv;i+) D2*i-1=dPi; D2*i=dU2i; /*pv节点部分，基数位为dP,偶数位为dV2*/ /*形成pq节点子阵*/ for(i=1;i=pq;i+) for(j=1;jn;j+) i1=jiediani.num; j1=jiedianj.num; ei=jiediani.e; ej=jiedianj.e; fi=jiediani.f; fj=jiedianj.f; if(i!=j) /*求i!=j时的H、N、J、L*/ ykb2*i-12*j-1=-Bi1j1*ei+Gi1j1*fi; /* H 位于奇数行奇数列*/ ykb2*i-12*j=Gi1j1*ei+Bi1j1*fi; /* N 位于奇数行偶数列*/ ykb2*i2*j-1=-Gi1j1*ei-Bi1j1*fi; /* J 位于偶数行奇数列*/ ykb2*i2*j=-Bi1j1*ei+Gi1j1*fi; /* L 位于偶数行偶数列*/ else /*求i=j时的H、N、J、L*/ a=0;b=0; for(k=1;k=n;k+) if(k!=i) k1=jiediank.num; a=a+Gi1k1*jiediank.e-Bi1k1*jiediank.f; b=b+Gi1k1*jiediank.f+Bi1k1*jiediank.e; ykb2*i-12*j-1=2*Gi1i1*jiediani.f+b; /*H*/ ykb2*i-12*i=2*Gi1i1*jiediani.e+a; /*N*/ ykb2*i2*i-1=-2*Bi1i1*jiediani.f+a; /*J*/ ykb2*i2*j=-2*Bi1i1*jiediani.e-b; /*L*/ for(i=pq+1;i=pq+pv;i+) /*形成pv节点子阵*/ for(j=1;jn;j+) i1=jiediani.num; j1=jiedianj.num; ei=jiediani.e; ej=jiedianj.e; fi=jiediani.f; fj=jiedianj.f; if(i!=j) /*求i!=j时的H、N*/ ykb2*i-12*j-1=-Bi1j1*ei+Gi1j1*fi; /*H*/ ykb2*i-12*j=Gi1j1*ei+Bi1j1*fi; /*N*/ /*省略Rij=0及Sij=0*/ else /*求i=j时的H、N、R、S*/ a=0;b=0; for(k=1;k=n;k+) if(k!=i) k1=jiediank.num; a+=Gi1k1*jiediank.e-Bi1k1*jiediank.f; b+=Gi1k1*jiediank.f+Bi1k1*jiediank.e; ykb2*i-12*j-1=2*Gi1i1*jiediani.f+b; /*H*/ ykb2*i-12*i=2*Gi1i1*jiediani.e+a; /*N*/ ykb2*i2*j-1=2*fi; /*R*/ ykb2*i2*j=2*ei; /*S*/ for(i=1;i=2*(pq+pv);i+) /*把不平衡量矩阵编入修正方程*/ ykbi2*(pq+pv)+1=Di; /*输出迭代次数和雅可比矩阵*/ void print_ykb(int z) int i,j; fprintf(fp2,输出雅可比矩阵和迭代次数n); fprintf(fp2,n 迭代的次数为 %2dn, z); for(i=1;i=2*(pq+pv);i+) fprintf(fp2,n); for(j=1;jfabs(d) /*列选主元*/ d=ykbik; /*将列主元素存入d*/ l=i; /*用l记录列主元素所在行*/ i+; while(i=2*(pq+pv); if(l!=k) for(j=k;j=2*(pq+pv)+1;j+) /*交换行*/ t=ykblj; ykblj=ykbkj; ykbkj=t; for(j=k+1;j=2*(pq+pv)+1;j+) /*正消*/ ykbkj/=ykbkk; for(i=k+1;i=2*(pq+pv);i+) for(j=k+1;j=2*(pq+pv)+1;j+) ykbij-=ykbik*ykbkj; k+; while(k=1;i-) /*回代*/ t=0; for(j=i+1;j=2*(pq+pv);j+) t+=ykbij*Dj; Di=ykbi2*(pq+pv)+1-t; for(i=1;i=(pq+pv);i+) jiediani.e+=D2*i; jiediani.f+=D2*i-1; max=fabs(D1); for(i=1;imax)max=fabs(Di); /*迭代过程中的电压值*/ void voltage() int i,j; fprintf(fp2,nn输出 df,de:n); for(i=1;i=(pq+pv);i+) fprintf(fp2,n); fprintf(fp2,节点为 %2d df=%8.5f de=%8.5fn,i,D2*i-1,D2*i); fprintf(fp2,n输出迭代过程中的电压值:n); for(i=1;i=(pq+pv);i+) fprintf(fp2,n 节点为 %2d f=%8.5f e=%8.5fn,i,jiediani.f,jiediani.e); fprintf(fp2, n-n); float mul_Re(x1,y1,x2,y2)float x1,x2,y1,y2; float x; x=x1*x2-y1*y2; return(x); float mul_Im(x1,y1,x2,y2)float x1,x2,y1,y2; float y; y=x1*y2+x2*y1; return(y); void data2() /* 潮流计算结果 */ float sp=0,sq=0,mo1,mo2,mo3,mo4,mo5,mo6,ff; /*mo1线路Sij功率计算中使用中间量、mo2线路Sji功率计算中使用中间量*/ /*mo3、mo4、mo5、mo6平行线路功率计算中使用中间量*/ int i1,j1,i11=0,j11=0; static float PMM,QMM; fprintf(fp2,nn 各节点电压为:); for(i=1;i=n;i+) fprintf(fp2,n U%d=%8.6f + j %8.6f, jiediani.num,jiediani.e,jiediani.f); fprintf(fp2,nn 平衡节点功率为:); for(i=1;i=n;i+) i1=jiediann.num;j1=jiediani.num;sp+=mul_Re(Gi1j1,-Bi1j1,jiediani.e,-jiediani.f);sq+=mul_Im(Gi1j1,-Bi1j1,jiediani.e,-jiediani.f); jiediann.p=mul_Re(jiediann.e,jiediann.f,sp,sq); jiediann.q=mul_Im(jiediann.e,jiediann.f,sp,sq); fprintf(fp2,n S%d=%8.6f+j %8.6f,jiediann.num,jiediann.p,jiediann.q); fprintf(fp2,nn 线路功率如下:n); for(k=1;k=(m-dd);k+) i1=zhiluk.p1; j1=zhiluk.p2; for(l=1;l=n;l+) if(jiedianl.num=i1)i=l; if(jiedianl.num=j1)j=l; sp=mul_Re(jiediani.e-jiedianj.e,-jiediani.f+jiedianj.f,-Gi1j1,Bi1j1); sq=mul_Im(jiediani.e-jiedianj.e,-jiediani.f+jiedianj.f,-Gi1j1,Bi1j1); mo1=jiediani1.e*jiediani1.e+jiediani1.f*jiediani1.f; Pi1j1=mul_Re(jiediani.e,jiediani.f,sp,sq)+mo1*Gi10; Qi1j1=mul_Im(jiediani.e,jiediani.f,sp,sq)-mo1*Bj10; sp=mul_Re(jiedianj.e-jiediani.e,-jiedianj.f+jiediani.f,-Gj1i1,Bj1i1); sq=mul_Im(jiedianj.e-jiediani.e,-jiedianj.f+jiediani.f,-Gj1i1,Bj1i1); mo2=jiedianj1.e*jiedianj1.e+jiedianj1.f*jiedianj1.f; Pj1i1=mul_Re(jiedianj.e,jiedianj.f,sp,sq)+mo2*Gj10; Qj1i1=mul_Im(jiedianj.e,jiedianj.f,sp,sq)-mo2*Bj10; if(i1!=i11)|(j1!=j11) fprintf(fp2,n 线路%d-%d的功率:,i1,j1); fprintf(fp2, %8.6f + j %8.6f,Pi1j1,Qi1j1); fprintf(fp2,n