yakebi.txt

最优潮流计算中形成雅克比矩阵（用c原理）

void form_ykb()                      /* 形成雅可比矩阵 */
 { float ei,ej,fi,fj,a=0,b=0;
   int i1,j1,k1;
   float dP[M],dQ[M],dU2[M];

   
   for(i=1;i<=2*(pq+pv);i++)           /*初始化矩阵*/
      { dP[i]=0;
	dQ[i]=0;
	dU2[i]=0;
	D[i]=0;
	for(j=1;j<=2*(pq+pv)+1;j++)
	   ykb[i][j]=0;
      }


    for(i=1;i<=pq;i++)                   /*PQ节点功率不平衡量*/
      {  ei=jiedian[i].e;
	 fi=jiedian[i].f;
	 for(j=1;j<=n;j++)
	  { i1=jiedian[i].num;
	    j1=jiedian[j].num;
	    ej=jiedian[j].e;
	    fj=jiedian[j].f;
	    dP[i]+=ei*(G[i1][j1]*ej-B[i1][j1]*fj)
		  +fi*(G[i1][j1]*fj+B[i1][j1]*ej);
	    dQ[i]+=fi*(G[i1][j1]*ej-B[i1][j1]*fj)
		  -ei*(G[i1][j1]*fj+B[i1][j1]*ej);
	    }
	 dP[i]=jiedian[i].p-dP[i];
	 dQ[i]=jiedian[i].q-dQ[i];
	}


   for(i=pq+1;i<=pq+pv;i++)                /*PV节点有功功率不平衡量和电压不平衡量*/
	{  ei=jiedian[i].e;
	   fi=jiedian[i].f;
	   for(j=1;j<=n;j++)
	    { i1=jiedian[i].num;
	      j1=jiedian[j].num;
	      ej=jiedian[j].e;
	      fj=jiedian[j].f;
	      dP[i]+=ei*(G[i1][j1]*ej-B[i1][j1]*fj)
		    +fi*(G[i1][j1]*fj+B[i1][j1]*ej);
	     }
	   dP[i]=jiedian[i].p-dP[i];
	   dU2[i]=jiedian[i].v*jiedian[i].v-ei*ei-fi*fi;
	}


   for(i=1;i<=pq;i++)                                   /*形成不平衡量矩阵D[M]*/
      {  D[2*i-1]=dP[i];
	 D[2*i]=dQ[i]; }
   for(i=pq+1;i<=pq+pv;i++)
      {  D[2*i-1]=dP[i];
	 D[2*i]=dU2[i];  }

                              /*形成pq节点子阵*/
   for(i=1;i<=pq;i++)                                 
      for(j=1;j<n;j++)
	 {  i1=jiedian[i].num;
	    j1=jiedian[j].num;
	    ei=jiedian[i].e;
	    ej=jiedian[j].e;
	    fi=jiedian[i].f;
	    fj=jiedian[j].f;
	 
		if(i!=j)                                     /*求i!=j时的H、N、J、L*/
	    {  ykb[2*i-1][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi;      /* H */
	       ykb[2*i-1][2*j]=G[i1][j1]*ei+B[i1][j1]*fi;         /* N */
	       ykb[2*i][2*j-1]=-G[i1][j1]*ei-B[i1][j1]*fi;        /* J */
	       ykb[2*i][2*j]=-B[i1][j1]*ei+G[i1][j1]*fi;          /* L */
	     }
	  
		
		else                                         /*求i=j时的H、N、J、L*/
	    { a=0;b=0;
	      for(k=1;k<=n;k++)
		  if(k!=i)
		    { k1=jiedian[k].num;
		      a=a+G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;
		      b=b+G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;
		     }
	      ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].f+b;                         /*H*/
	      ykb[2*i][2*j]=-2*B[i1][i1]*jiedian[i].e-b;                              /*L*/
	      ykb[2*i-1][2*i]=2*G[i1][i1]*jiedian[i].e+a;                           /*N*/
	      ykb[2*i][2*i-1]=-2*B[i1][i1]*jiedian[i].f+a;                             /*J*/
	    }
	  }


   for(i=pq+1;i<=pq+pv;i++)                                  /*形成pv节点子阵*/
      for(j=1;j<n;j++)
	 {  i1=jiedian[i].num;
	    j1=jiedian[j].num;
	    ei=jiedian[i].e;
	    ej=jiedian[j].e;
	    fi=jiedian[i].f;
	    fj=jiedian[j].f;
	  if(i!=j)                                           /*求i!=j时的H、N*/                                          
	    { ykb[2*i-1][2*j-1]=-B[i1][j1]*ei+G[i1][j1]*fi;     /*H*/
	      ykb[2*i-1][2*j]=G[i1][j1]*ei+B[i1][j1]*fi;        /*N*/
	     }
	  
	  
	  else                                               /*求i=j时的H、N、R、S*/
	    { a=0;b=0;
	       for(k=1;k<=n;k++)
		  if(k!=i)
		    { k1=jiedian[k].num;
		      a+=G[i1][k1]*jiedian[k].e-B[i1][k1]*jiedian[k].f;
		      b+=G[i1][k1]*jiedian[k].f+B[i1][k1]*jiedian[k].e;
		     }
	       ykb[2*i-1][2*j-1]=2*G[i1][i1]*jiedian[i].f+b;    /*H*/
	       ykb[2*i-1][2*i]=2*G[i1][i1]*jiedian[i].e+a;      /*N*/
	       ykb[2*i][2*j-1]=2*fi;                            /*R*/
	       ykb[2*i][2*j]=2*ei;                              /*S*/
	      }
	 }

   for(i=1;i<=2*(pq+pv);i++)                          /*把不平衡量矩阵编入修正方程*/
      ykb[i][2*(pq+pv)+1]=D[i];


  }