搜索
/*利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。/*可计算最大节点数为100，可计算PQ,PV,平衡节点*//*可计算非标准变比和平行支路*/#includestdio.h#includemath.h#includestdlib.h#defineM100/*最大矩阵阶数*/#defineNl100/*迭代次数*/inti,j,k,a,b,c;/*循环控制变量*/intt,l;doubleP,Q,H,J;/*中间变量*/intn,/*节点数*/m,/*支路数*/pq,/*PQ节点数*/pv;/*PV节点数*/doubleeps;/*迭代精度*/doubleaa[M],bb[M],cc[M],dd[M],max,rr,tt;/*中间变量*/doublemo,c1,d1,c2,d2;/*复数运算函数的返回值*/doubleG[M][M],B[M][M],Y[M][M];/*节点导纳矩阵中的实部、虚部及其模方值*/doubleykb[M][M],D[M],d[M],dU[M];/*雅克比矩阵、不平衡量矩阵*/structjd/*节点结构体*/{intnum,ty;/*num为节点号，ty为节点类型*/doublep,q,S,U,zkj,dp,dq,du,dj;/*节点有功、无功功率，功率模值，电压模值,阻抗角牛顿--拉夫逊中功率不平衡量、电压修正量*/}jd[M];structzl/*支路结构体*/{intnumb;/*numb为支路号*/intp1,p2;/*支路的两个节点*/doublekx;/*非标准变比*/doubler,x;/*支路的电阻与电抗*/}zl[M];FILE*fp1,*fp2;voiddata()/*读取数据*/{inth,number;fp1=fopen(input.txt,r);fscanf(fp1,%d,%d,%d,%d,%lf\n,&n,&m,&pq,&pv,&eps);/*输入节点数,支路数,PQ节点数,PV节点数和迭代精度*/for(i=1;i=n;i++)/*输入节点编号、类型、输入功率和电压初值*/{fscanf(fp1,%d,%d,&number,&h);if(h==1)/*类型h=1是PQ节点*/{fscanf(fp1,%lf,%lf,%lf,%lf\n,&jd[i].p,&jd[i].q,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;}if(h==2)/*类型h=2是pv节点*/{fscanf(fp1,,%lf,%lf,%lf\n,&jd[i].p,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;jd[i].q=-1.567;}if(h==3)/*类型h=3是平衡节点*/{fscanf(fp1,,%lf,%lf\n,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;}}for(i=1;i=m;i++)/*输入支路阻抗*/fscanf(fp1,%d,%lf,%d,%d,%lf,%lf\n,&zl[i].numb,&zl[i].kx,&zl[i].p1,&zl[i].p2,&zl[i].r,&zl[i].x);fclose(fp1);if((fp2=fopen(output.txt,w))==NULL){printf(cannotopenfile!\n);exit(0);}fprintf(fp2,电力系统潮流计算\n);fprintf(fp2,**********原始数据*********\n);fprintf(fp2,================================================================================\n);fprintf(fp2,节点数:%d支路数:%dPQ节点数:%dPV节点数:%d精度:%f\n,n,m,pq,pv,eps);fprintf(fp2,------------------------------------------------------------------------------\n);for(i=1;i=pq;i++)fprintf(fp2,PQ节点:节点%dP[%d]=%fQ[%d]=%f\n,jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].q);for(i=pq+1;i=pq+pv;i++)fprintf(fp2,PV节点:节点%dP[%d]=%fU[%d]=%f初值Q[%d]=%f\n,jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].U,jd[i].num,jd[i].q);fprintf(fp2,平衡节点:节点%de[%d]=%ff[%d]=%f\n,jd[n].num,jd[n].num,jd[n].U,jd[n].num,jd[n].zkj);fprintf(fp2,-------------------------------------------------------------------------------\n);for(i=1;i=m;i++)fprintf(fp2,支路%d:相关节点:%d,%d非标准变比:kx=%fR=%fX=%f\n,i,zl[i].p1,zl[i].p2,zl[i].kx,zl[i].r,zl[i].x);fprintf(fp2,==============================================================================\n);}/*------------------------------------复数运算函数--------------------------------------*/doublemozhi(doublea0,doubleb0)/*复数求模值函数*/{mo=sqrt(a0*a0+b0*b0);returnmo;}doubleji(doublea1,doubleb1,doublea2,doubleb2)/*复数求积函数a1为电压模值，a2为阻抗角，a3为导纳实部，a4为导纳虚部*/{a1=a1*cos(b1);b1=a1*tan(b1);c1=a1*a2-b1*b2;d1=a1*b2+a2*b1;returnc1;returnd1;}doubleshang(doublea3,doubleb3,doublea4,doubleb4)/*复数除法求商函数*/{c2=(a3*a4+b3*b4)/(a4*a4+b4*b4);d2=(a4*b3-a3*b4)/(a4*a4+b4*b4);returnc2;returnd2;}/*--------------------------------计算节点导纳矩阵----------------------------------*/voidForm_Y(){for(i=1;i=n;i++)/*节点导纳矩阵元素附初值*/for(j=1;j=n;j++)G[i][j]=B[i][j]=0;for(i=1;i=n;i++)/*节点导纳矩阵的主对角线上的元素,非对地导纳加入相应的主对角线元素中(考虑非标准变比）*/for(j=1;j=m;j++)if(zl[j].p1==i){if(zl[j].kx==1){mozhi(zl[j].r,zl[j].x);if(mo==0)continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2;B[i][i]+=d2;}else{mozhi(zl[j].r,zl[j].x);if(mo==0)continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2/zl[j].kx+c2*(1-zl[j].kx)/(zl[j].kx*zl[j].kx);B[i][i]+=d2/zl[j].kx+d2*(1-zl[j].kx)/(zl[j].kx*zl[j].kx);}}elseif(zl[j].p2==i){if(zl[j].kx==1){mozhi(zl[j].r,zl[j].x);if(mo==0)continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2;B[i][i]+=d2;}else{mozhi(zl[j].r,zl[j].x);if(mo==0)continue;shang(1,0,zl[j].r,zl[j].x);G[i][i]+=c2/zl[j].kx+c2*(zl[j].kx-1)/zl[j].kx;B[i][i]+=d2/zl[j].kx+d2*(zl[j].kx-1)/zl[j].kx;}}for(k=1;k=m;k++)/*节点导纳矩阵非主对角线上（考虑非标准变比）的元素*/if(zl[k].kx==1){i=zl[k].p1;j=zl[k].p2;mozhi(zl[k].r,zl[k].x);if(mo==0)continue;shang(1,0,zl[k].r,zl[k].x);G[i][j]-=c2;B[i][j]-=d2;G[j][i]=G[i][j];B[j][i]=B[i][j];}else{i=zl[k].p1;j=zl[k].p2;mozhi(zl[k].r,zl[k].x);if(mo==0)continue;shang(1,0,zl[k].r,zl[k].x);G[i][j]-=c2/zl[k].kx;B[i][j]-=d2/zl[k].kx;G[j][i]=G[i][j];B[j][i]=B[i][j];}/*--------------------------输出节点导纳矩阵------------------------------*/fprintf(fp2,\n\n*********潮流计算过程*********\n);fprintf(fp2,\n==============================================================================\n);fprintf(fp2,\n节点导纳矩阵为:);for(i=1;i=n;i++){fprintf(fp2,\n);for(k=1;k=n;k++){fprintf(fp2,%f,G[i][k]);if(B[i][k]=0){fprintf(fp2,+j);fprintf(fp2,%f,B[i][k]);}else{B[i][k]=-B[i][k];fprintf(fp2,-j);fprintf(fp2,%f,B[i][k]);B[i][k]=-B[i][k];}}}fprintf(fp2,\n------------------------------------------------------------------------------\n);}/*-------------------------------牛顿——拉夫逊-------------------------------*/voidCalculate_Unbalanced_Para(){for(i=1;i=n;i++){if(jd[i].ty==1)/*计算PQ节点不平衡量*/{t=jd[i].num;cc[t]=dd[t]=0;for(j=1;j=n;j++){for(a=1;a=n;a++){if(jd[a].num==j)break;}P=Q=0;P=jd[a].U*(G[t][j]*cos(jd[i].zkj-jd[a].zkj)+B[t][j]*sin(jd[i].zkj-jd[a].zkj))