pq.cpp

PQ分解法潮流计算程序

         b=p->c;  
		 jy[count]=max(ii,jj);
         uyr[count]=-r;  
		 uyi[count]=-x;
         dyr[ii]+=r; 
		 dyi[ii]+=x+b;
         dyr[jj]+=r;  
		 dyi[jj]+=x+b;
         count++;
      }
      end=tran+num_tran;
      for(p=tran+1;p<=end;p++)
      {
	    ii=p->i;  
		jj=p->j;
        if(min(ii,jj)!=i) continue;
         r=p->a;  
		 x=p->b;  
		 b=r*r+x*x;
	     r/=b;   
		 x/=-b;  
		 b=p->c;
         dyr[ii]+=r;  
		 dyi[ii]+=x;
         jy[count]=max(ii,jj);
	      r/=b; 
		  x/=b;
         uyr[count]=-r;  
		 uyi[count]=-x;
         dyr[jj]+=r/b;   
		 dyi[jj]+=x/b;
         count++;
      }
   }
}

//     PQ 潮流    
void fdlf(int num_node,int num_line,int num_tran,int num_gene,int num_load,
     struct data *line,struct data *tran,struct data *gene,struct data *load,
     double **property,int *iy,int *jy,double *dyr,double *dyi,double *uyr,
     double *uyi,int iter_max,double error_max,int *nto,double *vm,
     double *va,double **pq_node,double **pq_load,int *conv)
{
   int nb,ct,i,j,k,nd,ne,iter,no_error,no_err,*type,*pol,*iuf,*juf,*ius,*jus;
   double error,err,*mis,*duf,*uuf,*dus,*uus;
   nd=num_node+1;  
   ne=nd+1;  
   nb=num_line+num_tran;
   type=(int *)calloc(nd,sizeof(int));
   if(type==NULL)  {  printf("\nno memory in type");   exit(0);  }
   pol=(int *)calloc(nd,sizeof(int));
   if(pol==NULL)  {  printf("\nno memory in pol");   exit(0);  }
   mis=(double *)calloc(nd,sizeof(double));
   if(mis==NULL)  {  printf("\nno memory in mis");   exit(0);  }
   duf=(double *)calloc(nd,sizeof(double));
   if(duf==NULL)   {  printf("\nno memory in duf");   exit(0);  }
   dus=(double *)calloc(nd,sizeof(double));
   if(dus==NULL)   {  printf("\nno memory in dus");   exit(0);  }
   iuf=(int *)calloc(ne,sizeof(int));
   if(iuf==NULL)   {  printf("\nno memory in iuf");   exit(0);  }
   ius=(int *)calloc(ne,sizeof(int));
   if(ius==NULL)   {  printf("\nno memory in ius");   exit(0);  }

   for(i=1;i<=num_gene;i++)
   {
      if(gene[i].j==0)  
	  { 
		  err=gene[i].a/sita; 
		  break;  
	  }
   }
   for(i=1;i<nd;i++)  
   {  
	   type[i]=1; 
	   va[i]=err; 
	   vm[i]=1.0e0; 
   }
   for(j=1;j<=num_gene;j++)
   {
      i=gene[j].i; 
	  k=gene[j].j; 
	  type[i]=k;
      if(k<=0)  vm[i]=gene[j].c;
   }
   for(j=1;j<=num_load;j++)
   {
      i=load[j].i;  
	  k=load[j].j;  
	  pol[i]=j;
      if(k>1)  type[i]=k;
   }
   ne=2*nb;  
   nb/=2;
loop:  ct=0;   printf("\nne=%d",ne);
   juf=(int *)calloc(ne,sizeof(int));
   if(juf==NULL)   {  printf("\nno memory in juf");   exit(0);  }
   uuf=(double *)calloc(ne,sizeof(double));
   if(uuf==NULL)   {  printf("\nno memory in uuf");   exit(0);  }
   fact(0,nd,num_line,num_tran,type,pol,property,line,tran,load,iy,jy,dyi,
 	    uyi,iuf,juf,duf,uuf,ne,&ct);
   if(ct==1) 
   {
	   ne+=nb;  
	   free(juf);  
	   free(uuf);  
	   zero(duf,1,nd);  
	   goto loop; 
   }
   jus=(int *)calloc(ne,sizeof(int));
   if(jus==NULL)   
   {  
	   printf("\nno memory in jus");  
  
	   exit(0); 
   }
   uus=(double *)calloc(ne,sizeof(double));
   if(uus==NULL)   
   {  
	   printf("\nno memory in uus");  
	   exit(0);  
   }
   fact(1,nd,num_line,num_tran,type,pol,property,line,tran,load,iy,jy,dyi,
	    uyi,ius,jus,dus,uus,ne,&ct);
//   开始迭代    
   iter=0;
   for(;;)
   {
      error=0.0e0;
      power_node(0,nd,pq_node,vm,va,iy,jy,dyr,dyi,uyr,uyi);
      mismach(0,nd,num_gene,num_load,property,pq_node,pq_load,gene,load,
              vm,mis,&err,&no_err);
      if(err>error)   { error=err; no_error=no_err; }
      if(err>error_max)  revise(nd,va,mis,iuf,juf,duf,uuf);
      power_node(1,nd,pq_node,vm,va,iy,jy,dyr,dyi,uyr,uyi);
      mismach(1,nd,num_gene,num_load,property,pq_node,pq_load,gene,load,
              vm,mis,&err,&no_err);
      if(err>error)   { error=err; no_error=no_err; }
      if(err>error_max)  revise(nd,vm,mis,ius,jus,dus,uus);
      printf("\niteration:  %5d    max error: %15.7e   at node: %5d",
            iter,error,nto[no_error]);
      iter++;
      if(error<=error_max)  {  *conv=1;  break;  }
      if(iter>iter_max)
      {
         *conv=-1;
	 printf("\n fast decoupled load flow is no convergence");
         break;
      }
   }
//   printf("\nthe top of the heap is: %ld bytes\n",(unsigned long) coreleft());
   free(iuf);  free(juf);  free(duf);  free(uuf);   free(ius);  free(jus);
   free(dus);  free(uus);  free(mis);  free(type);  free(pol);
}

//     形成因子表     
void fact(int kpq,int nd,int num_line,int num_tran,int *type,
     int *pol,double **property,struct data *line,struct data *tran,
     struct data *load,int *iy,int *jy,double *dyi,double *uyi,int *iu,
     int *ju,double *du,double *uu,int ne,int *ct)
{
   int tk,i,j,k,fa,ii,jj,kk;
   double *temp,x,a,b;
   struct data *p,*end;
   temp=(double *)calloc(nd,sizeof(double));
   if(temp==NULL)  {  printf("\nno memory in temp");   exit(0);  }
   for(fa=1,k=1;k<nd;k++)
   {
      tk=type[k];  
	  iu[k]=fa;  
	  zero(temp,k+1,nd);
      if(kpq==0)
      {
         for(end=line+num_line,p=line+1;p<=end;p++)
         {
	        ii=p->i;  
	 	    jj=p->j;
            if(min(ii,jj)!=k||ii==jj)  continue;
            x=1.0e0/p->b; 
			du[ii]-=x; 
			du[jj]-=x;
            kk=max(ii,jj); 
			if(type[kk]!=0)  temp[kk]=x;
         }
         for(end=tran+num_tran,p=tran+1;p<=end;p++)
         {
            ii=p->i; 
			jj=p->j;
            if(min(ii,jj)!=k)  continue;
            x=1.0e0/p->b; 
			du[ii]-=x;  
			du[jj]-=x;
            kk=max(ii,jj); 
			if(type[kk]!=0)  temp[kk]=x;
         }
         if(tk==0)  continue;
      }
      else
      {
         du[k]=dyi[k];
         if(tk<=0)  continue;
         if(tk>1)
		 {
            kk=pol[k]; 
			a=load[kk].c; 
			a*=a;  
			x=load[kk].b/a;
            du[k]=dyi[k]-x*(2.0e0*property[tk][3]+property[tk][4]);
         }
         for(kk=iy[k+1],j=iy[k];j<kk;j++)
         {
            i=jy[j];
            if(type[i]>0)  temp[i]=uyi[j];
         }
      }
      x=du[k];
      for(i=1;i<k;i++)
      {
         for(kk=0,j=iu[i+1],ii=iu[i];ii<j;ii++)
         {
	        jj=ju[ii];
            if(kk==0&&jj>k)  break;
            if(jj==k) 
			{  
				a=uu[ii];
				b=a*du[i]; 
				x-=a*b; 
				kk=1;  
				continue; 
			}
            if(kk==1)  temp[jj]-=b*uu[ii];
		 }
      }
      du[k]=x;
      for(j=k+1;j<nd;j++)
      {
         a=temp[j];  
		 if(a==0.0e0)  continue;
	     uu[fa]=a/x;  
		 ju[fa++]=j;
	    if(fa>=ne) 
		{ 
			*ct=1; 
			goto loop;  
		}
      }
   }
loop:   free(temp);
}

//      计算节点功率      
void power_node(int kpq,int nd,double **pq_node,double *vm,double *va,\
     int *iy,int *jy,double *dyr,double *dyi,double *uyr,double *uyi)
{
   int i,j,k,m,ip;
   double a,b,vi,vij,cij,si,co,aa;
   for(i=1;i<nd;i++) 
	   pq_node[i][kpq]=0.0e0;
   for(i=1;i<nd;i++)
   {
      vi=vm[i];
      if(kpq==0) 
		  a=dyr[i];
      else 
		  a=-dyi[i];
      pq_node[i][kpq]+=vi*vi*a;
      for(m=iy[i+1],k=iy[i];k<m;k++)
      {
         if(kpq==0)
		 {  
			 a=uyr[k];  
			 b=uyi[k];
		 }
         else  
		 { 
			 a=-uyi[k]; 
			 b=uyr[k];
		 }
         j=jy[k];  
		 vij=vi*vm[j]; 
		 cij=va[i]-va[j];
         aa=cij*cij; 
		 si=cij-cij*aa/6.0e0; 
		 co=1.0e0-0.5e0*aa+aa*aa/24.0e0;
         a=a*vij*co;  
		 b=b*vij*si;
         pq_node[i][kpq]+=a+b;   
		 pq_node[j][kpq]+=a-b;
      }
   }
}

//       计算节点功率误差      
void mismach(int kpq,int nd,int num_gene,int num_load,double **property,
     double **pq_node,double **pq_load,struct data *gene,struct data *load,
     double *vm,double *mis,double *err,int *no_err)
{
   int i,j;
   double v,a,b,c,pow;
   struct data *p,*end;
   for(i=1;i<nd;i++) 
	   mis[i]=-pq_node[i][kpq];
   for(end=load+num_load,p=load+1;p<=end;p++)
   {
      i=p->i;  j=p->j;
      if(kpq==0)  pow=p->a;
      else  pow=p->b;
      if(j>1)
      {
         v=vm[i]/p->c;
         if(kpq==0) { a=property[j][0];  b=property[j][1];  c=property[j][2]; }
         else       { a=property[j][3];  b=property[j][4];  c=property[j][5]; }
         pow*=c+v*(b+v*a);
      }
      pq_load[i][kpq]=pow;  mis[i]-=pow;
   }
   for(end=gene+num_gene,p=gene+1;p<=end;p++)
   {
      i=p->i;  j=p->j;
      if(kpq==0)
      {
         if(j==0)  mis[i]=0.0e0;
         else      mis[i]+=p->a;
      }
      else
      {
         if(j<=0)  mis[i]=0.0e0;
         else      mis[i]+=p->b;
      }
   }
   for(*err=0.0e0,i=1;i<nd;i++)
   {
      pow=mis[i];  mis[i]=pow/vm[i];  pow=fabs(pow);
      if(pow>*err)  { *err=pow;  *no_err=i;  }
   }
}

//      修正电压幅值和角度      
void revise(int nd,double *v,double *mis,int *iu,int *ju,double *du,
     double *uu)
{
   int i,j,k,m;
   for(i=1;i<nd;i++)
   {
      for(m=iu[i+1],j=iu[i];j<m;j++)  
	  { 
		  k=ju[j];  mis[k]-=uu[j]*mis[i]; 
	  }
      mis[i]/=du[i];
   }
   for(i=nd-2;i>0;i--)
   {
      for(m=iu[i+1],j=iu[i];j<m;j++)
      {
	 k=ju[j];  mis[i]-=uu[j]*mis[k];
      }
   }
   for(i=1;i<nd;i++)  v[i]-=mis[i];
}


//       求解线性方程组      
void eliminate(int m,double *x,double *s,int n,int *fa,int *iu,int *ju,
     double *uu)
{
   int i,j,k;
   double aa,bb;
   iu[m]=*fa;
   for(i=1;i<m;i++)
   {
      aa=x[i];  
	  if(aa==0.0e0)  continue;
      for(k=iu[i+1],j=iu[i];j<k;j++) 
		  x[ju[j]]-=aa*uu[j];
          s[m]-=aa*s[i];
   }
   aa=x[m];
   for(i=m+1;i<n;i++)
   {
      bb=x[i];  
	  if(bb==0.0e0)  continue;
      ju[*fa]=i;  
	  uu[*fa]=bb/aa; 
	  (*fa)++;
   }
   s[m]/=aa;
}

/*finding (i,j) element in addmittance matrix------------------------------*/
void czm(int i,int j,int *iy,int *jy,int *k)
{
   int m,ih,jh,ib;
   *k=0;  ih=i;  jh=j;
   if(ih>jh) {  m=ih;  ih=jh;  jh=m;  }
   for(ib=iy[ih+1],m=iy[ih];m<ib;m++)
   {
      if(jy[m]==jh)   {   *k=m;  break;  }
   }
}

/*array *a set up zero from start to end-----------------------------------*/
void zero(double *a,int start,int end)
{
   int i;
   for(i=start;i<end;i++)   a[i]=0.0e0;
}

//     计算支路功率      
void cbr(int num_line,int num_tran,struct data *line,struct data *tran,
	 int *nto,double *vm,double *va)
{
   int i,j,i2,j2;
   double r,x,t,b,d,cd,sd,ri,rj,xi,xj,vi,vj,vij,pij,qij,pji,qji,dpb,dqb,ph,qh;
   struct data *p,*end;
   fprintf(fou,"\nno.i  no.j     pij        qij        pji        qji        dph        dqh\n");
   ph=0.0e0;  qh=0.0e0;
   for(end=line+num_line,p=line+1;p<=end;p++)
   {
      i=p->i;  j=p->j;
      r=p->a;  x=p->b;  b=r*r+x*x;
      if(i==j)
      {
	 vi=vm[i];  b=vi*vi/b;
	 pij=r*b;  qij=x*b;  pji=0.0e0;  qji=0.0e0;
	 dpb=pij;  ph+=dpb;  dqb=qij;    qh+=dqb;
      }
      else
      {
	 r/=b;  x/=-b;  b=p->c;  d=va[i]-va[j];
	 vi=vm[i];  vj=vm[j];  vij=vi*vj;  vi*=vi;  vj*=vj;
	 cd=vij*cos(d);  sd=vij*sin(d);
	 pij=vi*r-r*cd-x*sd;  pji=vj*r-r*cd+x*sd;  dpb=pij+pji;  ph+=dpb;
	 qij=-vi*(b+x)+x*cd-r*sd;  qji=-vj*(b+x)+x*cd+r*sd;  dqb=qij+qji;  qh+=dqb;
      }
      fprintf(fou,"%4d%6d%10.4f%11.4f%11.4f%11.4f%11.4f%11.4f\n",
	      nto[i],nto[j],pij,qij,pji,qji,dpb,dqb);
   }
   for(end=tran+num_tran,p=tran+1;p<=end;p++)
   {
      i=p->i;  j=p->j;  r=p->a;  x=p->b;  t=p->c;
      b=t*(r*r+x*x);  r/=b;  x/=-b;
      b=t-1.0e0;  ri=r*b; xi=x*b;  rj=-ri/t;  xj=-xi/t;
      vi=vm[i];  vj=vm[j];  vij=vi*vj;  vi*=vi;  vj*=vj;
      d=va[i]-va[j];  cd=vij*cos(d);  sd=vij*sin(d);
      pij=vi*(ri+r)-r*cd-x*sd;  pji=vj*(rj+r)-r*cd+x*sd;  dpb=pij+pji;  ph+=dpb;
      qij=-vi*(xi+x)+x*cd-r*sd;  qji=-vj*(xj+x)+x*cd+r*sd;  dqb=qij+qji;  qh+=dqb;
      fprintf(fou,"%4d%6d%10.4f%11.4f%11.4f%11.4f%11.4f%11.4f\n",
	      nto[i],nto[j],pij,qij,pji,qji,dpb,dqb);
   }
   fprintf(fou,"\n\n 系统总损耗:");
   fprintf(fou,"\n    有功损耗= %f",ph);
   fprintf(fou,"\n    无功损耗= %f",qh);
}