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];
}