牛顿拉夫讯

baoruzhong · 发表于 2011-5-19 14:41:45

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有账号？立即加入

×

谁能帮忙看看这个牛顿拉夫讯怎么改进的，能用于配电网潮流计算？
牛顿拉夫讯
double PowerFlow(int iNode, int iBalc, ROAD *pRoad, NodePQ *pNodePQ, bool vOutputVoltage = false)
{
int M=2*(iNode-1);
int flag1[N]; // PV节点标志

double g[N][N]; // 电导
double b[N][N]; // 电纳
double dvv[N]; // 电压微分
double angle[N]; // 相角
double p[N],dp[N]; // 有功及其微分
double q[N],dq[N]; // 无功及其微分
double e[N]; // 节点电压实部
double f[N]; // 节点电压虚部
double *jcbb = new double[M];
double *jcb = new double  [M][M];

int i,j,i1,j1,n,m,t,k1,balc,node/*,flag,npv*/;
double rr,xx,bb,k,com1,com2,sump=0.0,sumq=0.0;
double lossp=0.0, lossq=0.0;
node = iNode; // 节点数
balc = iBalc - 1; // 平衡节点
for(i=0;i<node;i++){
  flag1[i]=0;
  dp[i]=0.0;
  dq[i]=0.0;
  dvv[i]=0.0;
}

for(i=0;i<node;i++) {
  for(j=0;j<node;j++){
g[i][j]=0.0;
b[i][j]=0.0;
  }
}
for(n = 0; n < node-1; n++){
  i = pRoad[n].head;
  j = pRoad[n].tail;
  rr = pRoad[n].r;
  xx = pRoad[n].x;
//  printf("%d\t%d\t%.3f\t%.3f\n", pRoad[n].head, pRoad[n].tail, pRoad[n].r, pRoad[n].x);
  bb = 0;
  k = 1;

  g[i-1][j-1]+=-(rr/(rr*rr+xx*xx)/k);
  b[i-1][j-1]+=-(-xx/(rr*rr+xx*xx)/k);
  g[j-1][i-1]+=-(rr/(rr*rr+xx*xx)/k);
  b[j-1][i-1]+=-(-xx/(rr*rr+xx*xx)/k);
  g[i-1][i-1]+=rr/((rr*rr+xx*xx)*k*k);
  g[j-1][j-1]+=rr/(rr*rr+xx*xx);
  b[i-1][i-1]+=bb+(-xx/((rr*rr+xx*xx)*k*k));
  b[j-1][j-1]+=bb+(-xx/(rr*rr+xx*xx));
}
for(i=0;i<node;i++){
  angle[i] = 0;
  v[i] = 1;
  p[i] = pNodePQ[i].p;
  q[i] = pNodePQ[i].q;
//  printf("P:%.3f\t Q:%.3f\n", pNodePQ[i].p, pNodePQ[i].q);
  e[i]=v[i]*cos(angle[i]*PI/180);
  f[i]=v[i]*sin(angle[i]*PI/180);
}
// flag1[1]=1;

// printf("共有n=%d条支路\n",n);
// printf("节点导纳矩阵如下：\n");

// for(i=0;i<node;i++){
//  for(j=0;j<node;j++)
// printf("%lf+j%lf ", g[i][j],b[i][j]);
//  printf("\n\n");}
//
// printf("初始值如下：\n");
// printf("节点号电压    相角    有功    无功       e       f\n");
// for(i=0;i<node;i++)
//  printf("%d %lf  %lf  %lf  %lf  %lf  %lf\n", (i+1), v[i], angle[i], p[i], q[i], e[i], f[i]);

for(t=0;t<50;t++){
//  printf("\n第%d次迭代开始：\n", t);
  for(i=0;i<node;i++){
if(i==balc) /*如果是平衡节点*/
continue;
else{
com1=0.0, com2=0.0;
for(j=0;j<node;j++){
   com1+=g[i][j]*e[j]-b[i][j]*f[j];
   com2+=g[i][j]*f[j]+b[i][j]*e[j];
}
if(flag1[i]!=1){                      /*如果是PQ节点*/
   dp[i]=p[i]-e[i]*com1-f[i]*com2;
   dq[i]=q[i]-f[i]*com1+e[i]*com2;
}
else{
   dp[i]=p[i]-e[i]*com1-f[i]*com2;
   dvv[i]=v[i]*v[i]-(e[i]*e[i]+f[i]*f[i]);
}
}
  }

  for(i=0,j=0;i<node;i++){
if(i==balc) continue;
else if(flag1[i]!=1){
jcbb[j++]=dp[i];
jcbb[j++]=dq[i];
// printf("\ndp[%d]=%lf\ndq[%d]=%lf",(i+1),dp[i],(i+1),dq[i]);
}
else{
jcbb[j++]=dp[i];
jcbb[j++]=dvv[i];
// printf("\ndp[%d]=%lf\ndvv[%d]=%lf",(i+1),dp[i],(i+1),dvv[i]);
}
  }
//  printf("\n");

  for(i=0;i<node;i++) {
if(fabs(dp[i])>=E||fabs(dq[i])>=E||fabs(dvv[i])>=E)
break;
  }
  if(i>=node)
break;


  for(i=0,i1=0;i<node;i++){
if(i==balc)
continue;
i1++;
n=2*i1;
for(j=0,j1=0;j<node;j++){
if(j==balc) continue;
j1++;
m=j1*2;
if(i!=j){
   jcb[n-2][m-2]=-(g[i][j]*e[i]+b[i][j]*f[i]);
   jcb[n-2][m-1]=b[i][j]*e[i]-g[i][j]*f[i];
   if(flag1[i]!=1)
   {
   jcb[n-1][m-2]=jcb[n-2][m-1];
   jcb[n-1][m-1]=-jcb[n-2][m-2];
   }
   else
   jcb[n-1][m-2]=jcb[n-1][m-1]=0.0;
}
else{
   com1=0.0, com2=0.0;
   for(k1=0;k1<node;k1++)
   {
   com1+=g[i][k1]*e[k1]-b[i][k1]*f[k1];
   com2+=g[i][k1]*f[k1]+b[i][k1]*e[k1];
   }
   jcb[n-2][m-2]=-com1-g[i][i]*e[i]-b[i][i]*f[i];
   jcb[n-2][m-1]=-com2+b[i][i]*e[i]-g[i][i]*f[i];
   if(flag1[i]!=1)
   {
   jcb[n-1][m-2]=com2+b[i][i]*e[i]-g[i][i]*f[i];
   jcb[n-1][m-1]=-com1+g[i][i]*e[i]+b[i][i]*f[i];
   }
   else
   {
   jcb[n-1][m-2]=-2*e[i];
   jcb[n-1][m-1]=-2*f[i];
   }
}
}
  }

//  printf("\n第%d次输出-J矩阵如下：\n", t);
  for(i=0;i<2*(node-1);i++){
for(j=0;j<2*(node-1);j++){
jcb[i][j]*=-1.0;
// printf("%lf ", jcb[i*M+j]);
}
// printf("\n\n");
  }

  /*算线形方程组*/
  if (agaus(&jcb[0],&jcbb[0],2*(node-1))!=0) {
for (i=0,j=0;i<2*(node-1);i++) {
if(j==balc)
   j++;

if(i%2==0){
   e[j]+=jcbb[i];
//    printf("de[%d]=%lf e[%d]=%lf\n", (j+1),jcbb[i],(j+1),e[j]);
}
else{
   f[j]+=jcbb[i];
//    printf("df[%d]=%lf f[%d]=%lf\n",(j+1),jcbb[i],(j+1),f[j]);
   j++;
}
}
  }
}

// printf("共经过了%d次迭代\n", t);

/*求平衡节点P，Q*/
i=balc;
com1=0.0, com2=0.0;
for(j=0;j<node;j++){
  com1+=g[i][j]*e[j]-b[i][j]*f[j];
  com2+=g[i][j]*f[j]+b[i][j]*e[j];
}
p[i]=e[i]*com1+f[i]*com2;
q[i]=f[i]*com1-e[i]*com2;

/*求其他各节点P，Q*/
for(i=0;i<node;i++){
  if(i==balc) continue;
  else
  {
com1=0.0, com2=0.0;
for(j=0;j<node;j++)
{
com1+=g[i][j]*e[j]-b[i][j]*f[j];
com2+=g[i][j]*f[j]+b[i][j]*e[j];
}
p[i]=e[i]*com1+f[i]*com2;
q[i]=f[i]*com1-e[i]*com2;
  }
}

// printf("最终结果为：\n");
// printf("节点号电压    相角    有功    无功       e       f\n");
for(i=0;i<node;i++)
{
  v[i]=sqrt(e[i]*e[i]+f[i]*f[i]);
  angle[i]=atan(f[i]/e[i])*180/PI;
  sump+=p[i];
  sumq+=q[i];
//  printf("%d %lf  %lf  %lf  %lf  %lf  %lf\n", (i+1), v[i], angle[i], p[i], q[i], e[i], f[i]);
}

账号		自动登录	找回密码
密码			立即加入

牛顿拉夫讯

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

相关帖子