九节点系统潮流计算编程牛N_R法

#include

#include

float divRe(float b1,float b2,float b3,float b4) {

float a1r;

a1r=(b1*b3+b2*b4)/(b3*b3+b4*b4); return(a1r); }

float divIm(float b1,float b2,float b3,float b4) {

float a1i;

a1i=(b2*b3-b1*b4)/(b3*b3+b4*b4); return(a1i); }

float mulRe(float b1,float b2,float b3,float b4) {

float a2r;

a2r=b1*b3-b2*b4; return(a2r); }

float mulIm(float b1,float b2,float b3,float b4) {

float a2i;

a2i=b2*b3+b1*b4; return(a2i);

}

float Max(float a[],int n) {int i; float max; max= fabs(a[0]); for(i=1;imax)max=fabs(a[i]);} return(max); }

void main() {

int i,j,k,h,km; int T=16;

float eps,sumpi1,sumpi2,sumqi1,sumqi2,max,sumir,sumii,I1r,I1i,t,xx,xxx; float pi0[8],qi0[8],detpi[8],detqi[8],Iir0[8],Iii0[8],J0[16][16], detsi[16],detui[16], delta_p[9][9],delta_q[9][9], a[16][32],ni[16][16],H[8][8],N[8][8], J[8][8],L[8][8],ei1[9],fi1[9],sp[9][9],sq[9][9]; static float ybr[9][9]={ {3.3074,-1.3652,0,0,0,-1.9422,0,0,0}, {-1.3652,2.5528,-1.1876,0,0,0,0,0,0}, {0,-1.1876,2.8047,-1.6171,0,0,0,0,0}, {0,0,-1.6171,2.7722,-1.1551,0,0,0,0}, {0,0,0,-1.1551,2.437,-1.282,0,0,0}, {-1.9422,0,0,0,-1.282,3.2242,0,0,0}, {0}, {0}, {0}};

static float ybi[9][9]={ {-39.3089,11.6041,0,0,0,10.5107,0,0,17.3611}, {11.6041,-17.3382,5.9751,0,0,0,0,0,0}, {0,5.9751,-35.4456,13.6980,0,0,16,0,0}, {0,0,13.6980,-23.3033,9.7843,0,0,0,0}, {0,0,0,9.782,-32.1538,5.5882,0,17.0648,0}, {10.5107,0,0,0,5.5882,-15.841,0,0,0}, {0,0,16,0,0,0,-16,0,0}, {0,0,0,0,17.0648,0,0,-17.0648,0}, {17.3611,0,0,0,0,0,0,0,-17.3611}}; static float yd[9][9]={ {0,0.088,0,0,0,0.079,0,0,0}, {0.088,0,0.153,0,0,0,0,0,0}, {0,0.153,0,0.0745,0,0,0,0,0},

{0,0,0.0745,0,0.1045,0,0,0,0}, {0,0,0,0.1045,0,0.179,0,0,0}, {0.079,0,0,0,0.179,0,0,0,0}, {0}, {0}, {0}} ;

float ei0[9]={1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.04}; float fi0[9]={0.0};

float pi[9]={0,-1.25,0,-1.0,0,-0.9,1.63,0.85,0}; float qi[9]={0,-0.5,0,-0.35,0,-0.3,0,0,0}; h=0; km=15;

eps=0.00001; do{ h+=1;

printf(\"\\nNow The %dth times\\n\

for(i=0;i<8;i++)

{printf(\"ei0[%d]=%f\\ printf(\"fi0[%d]=%f\\ }

for(i=0;i<8;i++){

printf(\"pi[%d]=%f\\ printf(\"qi[%d]=%f\\ } sumpi2=0; sumqi2=0;

for(i=0;i<8;i++) {for(j=0;j<9;j++) { sumpi1=(ei0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])+fi0[i]*(ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j])); sumpi2+=sumpi1; }

pi0[i]=sumpi2;

printf(\"pi0[%d]=%f\\ sumpi2=0; }

for(i=0;i<8;i++) {for(j=0;j<9;j++)

{ sumqi1=(fi0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])-ei0[i]*(ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j])); sumqi2+=sumqi1; }

qi0[i]=sumqi2;

printf(\"qi0[%d]=%f\\ sumqi2=0; }

for(i=0;i<8;i++) {

detpi[i]=pi[i]-pi0[i]; detqi[i]=qi[i]-qi0[i]; if(i==6||i==7) { qi0[i]=ei0[i]*ei0[i]+fi0[i]*fi0[i]; detqi[i]=1.051-qi0[i]; }

printf(\"detpi[%d]=%f\\ printf(\"detqi[%d]=%f\\ }

//*************************************************************************************************

//节点的注入电流表达式 for(i=0;i<8;i++) { Iii0[i]=0; Iir0[i]=0; }

for(i=0;i<8;i++)

{for(j=0;j<9;j++) { Iir0[i]+=ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j]; Iii0[i]+=ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j]; }

}

//********************************************************************************************* //求解NHJL矩阵 for(i=0;i<8;i++) {for(j=0;j<8;j++) {if(i==j)

{if(i==6||i==7) {

H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i]; N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i]; J[i][j]=2*fi0[i]; L[i][j]=2*ei0[i]; } else { H[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]+Iii0[i]; N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]+Iir0[i]; J[i][j]=-ybr[i][j]*ei0[j]-ybi[i][i]*fi0[j]+Iir0[i]; L[i][j]=-ybi[i][j]*ei0[j]+ybr[i][j]*fi0[j]-Iii0[i]; } }

else { if(i==6||i==7) { H[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j]; N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]; J[i][j]=0; L[i][j]=0; }

else { H[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j]; N[i][j]=ybr[i][j]*ei0[j]+ybi[i][j]*fi0[j]; J[i][j]=-ybi[i][j]*fi0[j]-ybr[i][j]*ei0[j]; L[i][j]=ybr[i][j]*fi0[j]-ybi[i][j]*ei0[j]; } }

} }

//形成jacobian矩阵 for(i=0;i<16;i++) for(j=0;j<16;j++){ if(i%2==0&&j%2==0) J0[i][j]=H[i/2][j/2]; else if(i%2==0&&j%2!=0) J0[i][j]=N[i/2][(j-1)/2]; else if(i%2!=0&&j%2==0) J0[i][j]=J[(i-1)/2][j/2]; else J0[i][j]=L[(i-1)/2][(j-1)/2]; }

//for(i=0;i<16;i++) //for(j=0;j<16;j++)

// printf(\"J0[%d][%d]=%f\\

//************************************************************************************** //求detui

//*************************************************************************************

for(i=0;i<16;i++) {if(i%2==0) detsi[i]=detpi[i/2]; else detsi[i]=detqi[(i-1)/2]; }//将detp和detq用一个数组表示 for(i=0;ifor(j=0;j<(2*T);j++) { xx=a[i][j]*t; a[k][j]=a[k][j]-xx; } } } } for(i=0;ifor(i=0;ifor(j=0;j{xxx=xxx+ni[i][j]*detsi[j]; }

detui[i]=xxx;}

//检测detui满足要求与否 max=Max(detui,16);

printf(\"max=%f\\n\

//****************************************************************** //算第n次迭代后的u for(i=0;i//********************************************* for(i=0;i<8;i++)//下一次迭代赋初值 {ei0[i]=ei1[i]; fi0[i]=fi1[i]; }

for(i=0;i<8;i++)

{printf(\"ei0[%d]=%f\\ printf(\"fi0[%d]=%f\\ }

for(i=0;i<8;i++)

{pi[i]=detpi[i]+pi0[i]; qi[i]=detqi[i]+qi0[i]; }

}while(max>eps&&h//********************************************** for(i=0;i<9;i++)//平衡节点功率计算 {

I1r=mulRe(ybr[8][i],-ybi[8][i],ei0[i],-fi0[i]); I1i=mulIm(ybr[8][i],-ybi[8][i],ei0[i],-fi0[i]); sumir+=I1r; sumii+=I1i; }

pi[8]=mulRe(ei0[8],fi0[8],sumir,sumii); qi[8]=mulIm(ei0[8],fi0[8],sumir,sumii); printf(\"S9=%f+j%f\\n\

sumpi1=0; sumpi2=0; sumqi1=0; sumqi2=0;

for(i=0;i<9;i++)

{for(j=0;j<9;j++)

if(i!=j&&ybi[i][j]!=0) {

sumpi1=mulRe(ei0[i],fi0[i],0.0,yd[i][j]); sumqi1=mulIm(ei0[i],fi0[i],0.0,yd[i][j]);

sumpi2=mulRe(ei0[i]-ei0[j],fi0[i]-fi0[j],-ybr[i][j],-ybi[i][j]); sumqi2=mulIm(ei0[i]-ei0[j],fi0[i]-fi0[j],-ybr[i][j],-ybi[i][j]);

sumpi1+=sumpi2; sumqi1+=sumqi2;

sp[i][j]=mulRe(ei0[i],fi0[i],sumpi1,-sumqi1); sq[i][j]=mulIm(ei0[i],fi0[i],sumpi1,-sumqi1);

printf(\"S%d=%f+j%f\\n\ } }

for(i=0;i<9;i++) { for(j=0;j<9;j++) if(j !=i&&ybi[i][j]!=0) { delta_p[i][j]=sp[i][j]+sp[j][i]; delta_q[i][j]=sq[i][j]+sq[j][i]; if( (delta_p[i][j] !=0) || (delta_q[i][j] !=0)) { printf(\"dS%d=%f+j%f\\n\ } } }

ei1[8]=ei0[8]; fi1[8]=fi0[8]; for(i=0;i<9;i++)

{printf(\"u%d=%f<%f\\n\}