电力系统潮流程序开发

JU.Add(JT[i]); JJ.Add(i);

//JU[Ucount]=JT[i];

//JJ[Ucount]=i;//非零元素的列号 } }

DT[nl]=DT[nl]/JT[nl];

JI[nl]=count;//非零元素的个数

6.3 计算雅克比矩阵

其中雅可比矩阵的各元素分别为：

??Pi??PiHij?,Nij??ej?fjJij???Qi??Qi,Lij??ej?fj

非对角元素：Hij??Gijfi?BijeiNij??Gijei?BijfiJij?Gijei?Bijfi??NijLij??Gijfi?Bijei?Hij对角元素：Hii????Gijfj?Bijej??Biiei?Giifij?1nn

Nii????Gijej?Bijfj??Giiei?Biifi

j?1Jii????Gijej?Bijfj??Giiei?Biifij?1nLii???Gijfj?Bijej??Biiei?Giifij?1n雅可比矩阵的特点：

1) 雅可比矩阵为一非奇异方阵。传统的，当节点电压以极坐标表示时，该矩阵为2(n-1)-m阶方阵（m为PV节点数）；当节点电压以直角坐标表示时，该矩阵为2(n-1)阶方阵。现在，为了便于编程，一般为经过处理的2n阶。

2) 矩阵元素与节点电压有关，故每次迭代时都要重新计算。

3) 与导纳矩阵具有相似的结构，当Yij=0，Hij、Nij、Jij、Lij均为0，因此也是高度稀疏的矩阵。

4) 具有结构对称性，但数值不对称。

注意：当在计算过程中发生PV节点的无功功率越限时，PV节点要转化为PQ节点。

//形成雅克比矩阵，并进行消去、规格化运算，采取形成一行消去一行的算法

for(int i=0;iGetVoltageP(); double Uid=Nodei->GetVoltageD(); //PQ节点无功不平衡量的雅克比 if(Nodei->GetNodeType()==0) { int nJii,nLii;//对角线元素的位置 double Jii=0.0;//对角线元素初值 double Lii=0.0; //下三角部分 for(int k=0;kGetVoltageP(); double Ujd=Nodek->GetVoltageD(); JT[nk]=YUArray[j].b*Uiv-YUArray[j].g*Uid;//Jij Jii=Jii+YUArray[j].g*Ujd+YUArray[j].b*Ujv;//Jii JT[nk+1]=YUArray[j].g*Uiv+YUArray[j].b*Uid;//Lij Lii=Lii-(YUArray[j].g*Ujv-YUArray[j].b*Ujd);//Lii sprintf_s(sw,sizeof(sw),\,YUArray[j].g,YUArray[j].b,Uiv,Uid,i,k,JT[nk],i,k,JT[nk+1]); fwrite(sw,strlen(sw),1,pf); break; }

} nk+=2; count=count+YIArray[k]; } //JT的对角线元素 nJii=nk++; JT[nJii]=Jii+2.0*YDArray[i].b*Uiv; nLii=nk++; JT[nLii]=Lii+2.0*YDArray[i].b*Uid; //上三角部分 int nl=i;//当前已形成的雅克比元素对应的列号（节点号） for(int j=count;jGetVoltageP(); double Ujd=Nodej->GetVoltageD(); if(Nodej->GetNodeType()!=2) { JT[nk]=YUArray[j].b*Uiv-YUArray[j].g*Uid;//Jij JT[nk+1]=YUArray[j].g*Uiv+YUArray[j].b*Uid;//Lij sprintf_s(sw,sizeof(sw),\,YUArray[j].g,YUArray[j].b,Uiv,Uid,i,nj,JT[nk],i,nj,JT[nk+1]); fwrite(sw,strlen(sw),1,pf); } JT[nJii]=JT[nJii]+YUArray[j].g*Ujd+YUArray[j].b*Ujv;//Jii JT[nLii]=JT[nLii]-(YUArray[j].g*Ujv-YUArray[j].b*Ujd);//Lii } sprintf_s(sw,sizeof(sw),\,nJii,JT[nJii],nLii,JT[nLii],nJii,DT[nJii]); fwrite(sw,strlen(sw),1,pf); Gaosi(nJii); }

//PV节点电压误差的雅克比 else if(Nodei->GetNodeType()==1) { JT[2*i]=-2.0*Uiv; JT[2*i+1]=-2.0*Uid; sprintf_s(sw,sizeof(sw),\,i,i,JT[2*i],i,i,JT[2*i+1],2*i,DT[2

*i]); fwrite(sw,strlen(sw),1,pf); Gaosi(2*i); } nk=0;//JT的下标 count=0;//YU的下标 //有功不平衡量对应的雅克比矩阵行 //临时数组清零 for(int k=0;k<2*Npq+2*Npv;k++) { JT[k]=0.0;//临时数组清零 } int nHii,nNii;//对角线元素的位置 double Hii=0.0;//对角线元素初值 double Nii=0.0; //下三角部分 for(int k=0;k

PNode *Nodek=OpNodeArray[k]; for(int j=count;jGetVoltageP(); double Ujd=Nodek->GetVoltageD(); JT[nk]=-(YUArray[j].g*Uiv+YUArray[j].b*Uid);//Hij Hii=Hii-(YUArray[j].g*Ujv-YUArray[j].b*Ujd);//Hii JT[nk+1]=YUArray[j].b*Uiv-YUArray[j].g*Uid;//Nij Nii=Nii-(YUArray[j].g*Ujd+YUArray[j].b*Ujv);//Nii sprintf_s(sw,sizeof(sw),\,YUArray[j].g,YUArray[j].b,Uiv,Uid,i,k,JT[nk],i,k,JT[nk+1]); fwrite(sw,strlen(sw),1,pf); break; } } nk+=2; count=count+YIArray[k]; } //JT的对角线元素 nHii=nk++;