九节点系统潮流计算编程牛N-R法

如图所示系统，试计算潮流分布，相关数据见《版潮流计算用户手册》P121。

#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;i

{if(fabs(a[i])>max)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]={

{,,0,0,0,,0,0,0},

{,,,0,0,0,0,0,0},

{0,,,,0,0,0,0,0},

{0,0,,,,0,0,0,0},

{0,0,0,,,,0,0,0},

{,0,0,0,,,0,0,0},

{0},

{0},

{0}};

static float ybi[9][9]={

{,,0,0,0,,0,0,},

{,,,0,0,0,0,0,0},

{0,,,,0,0,16,0,0},

{0,0,,,,0,0,0,0},

{0,0,0,,,,0,,0},

{,0,0,0,,,0,0,0},

{0,0,16,0,0,0,-16,0,0},

{0,0,0,0,,0,0,,0},

{,0,0,0,0,0,0,0,}};

static float yd[9][9]={

{0,,0,0,0,,0,0,0},

{,0,,0,0,0,0,0,0},

{0,,0,,0,0,0,0,0},

{0,0,,0,,0,0,0,0},

{0,0,0,,0,,0,0,0},

{,0,0,0,,0,0,0,0},

{0},

{0},

{0}} ;

float ei0[9]={,,,,,,,,};

float fi0[9]={};

float pi[9]={0,,0,,0,,,,0};

float qi[9]={0,,0,,0,,0,0,0};

h=0;

km=15;

eps=;

do{

h+=1;

printf("\nNow The %dth times\n",h);

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

{printf("ei0[%d]=%f\t",i,ei0[i]);

printf("fi0[%d]=%f\t",i,fi0[i]);

}

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

printf("pi[%d]=%f\t",i,pi[i]);

printf("qi[%d]=%f\t",i,qi[i]);

}

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\t",i,pi0[i]);

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\t",i,qi0[i]);

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]=[i];

}

printf("detpi[%d]=%f\t",i,detpi[i]);

printf("detqi[%d]=%f\t",i,detqi[i]);

}

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

//节点的注入电流表达式

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

{ Iii0[i]=0;

Iir0[i]=0;

}

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

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

{