牛顿拉夫逊迭代法极坐标潮流计算C语言程序

fprintf(fp2,\ \ B[i][k]=-B[i][k]; }

fprintf(fp2,\ ------------------------------------------------------------------------------\\n\ }

/*-------------------------------牛顿——拉夫逊-------------------------------*/ void Calculate_Unbalanced_Para() {

for(i=1;i<=n;i++) {

if(jd[i].ty==1) /*计算PQ节点不平衡量*/ { }

if(jd[i].ty==2) /*计算PV节点不平衡量*/ {

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;

t=jd[i].num; cc[t]=dd[t]=0; for(j=1;j<=n;j++) { }

jd[i].dp=jd[i].p-jd[i].U*cc[t]; jd[i].dq=jd[i].q-jd[i].U*dd[t];

for(a=1;a<=n;a++) { } 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)); Q=jd[a].U*(G[t][j]*sin(jd[i].zkj-jd[a].zkj)-B[t][j]*cos(jd[i].zkj-jd[a].zkj)); cc[t]+=P; dd[t]+=Q;

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));

Q=jd[a].U*(G[t][j]*sin(jd[i].zkj-jd[a].zkj)-B[t][j]*cos(jd[i].zkj-jd[a].zkj));

cc[t]+=P; }

jd[i].dp=jd[i].p-jd[i].U*cc[t];

dd[t]+=Q;

jd[i].q=jd[i].U*dd[t];

for(i=1;i<=pq;i++) /*形成不平衡量矩阵D[M]*/ {

}

D[2*i-1]=jd[i].dp;

D[2*i]=jd[i].dq;

for(i=pq+1;i<=n-1;i++) {

}

fprintf(fp2,\ 不平衡量为:\ /*输出不平衡量*/ for(i=1;i<=pq;i++) { }

for(i=pq+1;i<=n-1;i++) { }

t=jd[i].num;

fprintf(fp2,\ dp[%d]=%f\

t=jd[i].num;

fprintf(fp2,\ dp[%d]=%f\D[pq+i]=jd[i].dp;

fprintf(fp2,\ dq[%d]=%f\

}

void Form_Jacobi_Matric() /*形成雅克比矩阵*/ {

for(i=1;i<=pq;i++) /*形成pq节点子阵*/

for(j=1;j

int j1=jd[j].num; double Ui=jd[i].U; double Uj=jd[j].U; double zi=jd[i].zkj; double zj=jd[j].zkj;

if(j<=pq) /*求j<=pq时的H、N、J、L */

J */ L */

{

if(i!=j) /*求i!=j时的H、N、J、L*/

{

ykb[2*i-1][2*j-1]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */

ykb[2*i-1][2*j]=Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* N */ ykb[2*i][2*j-1]=-ykb[2*i-1][2*j]; /* ykb[2*i][2*j]=ykb[2*i-1][2*j-1]; /*

} {

for(k=1;k<=n;k++) {

int k1=jd[k].num; H=J=0;

J=jd[k].U*(G[i1][k1]*cos(jd[i].zkj-jd[k].zkj)+B[i1][k1]*sin(jd[i].zkj-jd[k].zkj));

else /*求i=j时的H、N、J、L*/

aa[i]=0;bb[i]=0;

H=jd[k].U*(G[i1][k1]*sin(jd[i].zkj-jd[k].zkj)-B[i1][k1]*cos(jd[i].zkj-jd[k].zkj));

aa[i]=aa[i]+H; bb[i]=bb[i]+J;

}

ykb[2*i-1][2*j-1]=-Ui*(aa[i]-Ui*(G[i1][i1]*sin(jd[i].zkj-jd[i].zkj)-B[i1][i1]*cos(jd[i].zkj-jd[i].zkj))); /*H*/

ykb[2*i][2*j-1]=Ui*(bb[i]-Ui*(G[i1][i1]*cos(jd[i].zkj-jd[i].zkj)+B[i1][i1]*sin(jd[i].zkj-jd[i].zkj))); /*J*/ /*N*/ /*L*/

}

}

{ ykb[2*i-1][pq+j]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ }

}

ykb[2*i-1][2*j]=ykb[2*i][2*j-1]+2*Ui*Ui*G[i1][i1];

ykb[2*i][2*j]=-ykb[2*i-1][2*j-1]-2*Ui*Ui*B[i1][i1];

else /*求j>pq时的H、J */

ykb[2*i][pq+j]=-Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* J */

for(i=pq+1;i<=n-1;i++) /*形成pv节点子阵*/ for(j=1;j

{

int i1=jd[i].num;

}

int j1=jd[j].num; double Ui=jd[i].U; double Uj=jd[j].U; double zi=jd[i].zkj; double zj=jd[j].zkj;

if(j<=pq) /*求j<=pq时的H、N */

{ }

else /*求j>pq时的H*/

{

if(i!=j) /*求i!=j时的H*/

ykb[pq+i][pq+j]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ else /*求i=j时的H*/ {

aa[i]=0; {

int k1=jd[k].num;

H=0;

H=jd[k].U*(G[i1][k1]*sin(jd[i].zkj-jd[k].zkj)-B[i1][k1]*cos(jd[i].zkj-jd[k].zkj)); }

ykb[pq+i][2*j-1]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */

ykb[pq+i][2*j]=Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* N */

for(k=1;k<=n;k++)

aa[i]=aa[i]+H;

}

/*------------------------------输出雅克比矩阵--------------------------------*/ for(i=1;i<=(2*pq+pv);i++) { }

for(j=1;j<=2*pq+pv;j++) { }

fprintf(fp2,\ %f\ }

}

ykb[pq+i][pq+j]=-Ui*(aa[i]-Ui*(G[i1][i1]*sin(jd[i].zkj-jd[i].zkj)-B[i1][i1]*cos(jd[i].zkj-jd[i].zkj))); /*H*/

fprintf(fp2,\ 雅克比矩阵为: \

fprintf(fp2,\

void Solve_Equations() /* 求解修正方程组 (LU分解法)*/ {

double l[Nl][Nl]={0}; //定义L矩阵

double u[Nl][Nl]={0}; //定义U矩阵 double y[Nl]={0}; //定义数组Y double x[Nl]={0}; //定义数组X double a[Nl][Nl]={0}; //定义系数矩阵 double b[Nl]={0}; //定义右端项 double sum=0; int i,j,k,s; int n; n=2*pq+pv; for(i=0; i

for(j=0; j

a[i][j]=ykb[i+1][j+1]; } }

for(i=0; i

b[i]=D[i+1]; }

for(i=0; i

for(j=0; j

if(i==j) l[i][j] = 1; } }

for(i=0;i

{ u[0][i]=(float)(a[0][i]);} /*第一步：对矩阵U的首行进行计算*/ for(k=0;k

sum+=l[i][s]*u[s][k]; }

l[i][k]=(float)((a[i][k]-sum)/u[k][k]); }

for(j=k+1;j

sum+=l[k+1][s]*u[s][j]; }

u[k+1][j]=(float)((a[k+1][j]-sum));