📄 nr_powerflow.cpp
字号:
int nCount2 = bpNR->nJacobiCount;
double *pJacobi = bpNR->pJacobi;
double *pPQ_Uef = bpNR->pPQ_Uef;
double aii,bii;//引入节点注入电流Ii=aii+j(bii)
int i=0,j=0,k=0,l=0;//引入变量k,l的目的是形成【nCount2×nCount2】阶矩阵
for(i=0,k=0;i<nCount1;i++)
{
if(pType[i]==0) continue;//平衡节点则进入下一次循环
else if(pType[i]==1)//PQ节点
{
for(j=0,l=0;j<nCount1;j++)
{
if(i!=j)
{
if(pType[j]==0) continue;//平衡节点则进入下一次循环
pJacobi[(2*k)*nCount2+l*2+0]=-pB[i*nCount1+j]*pEn[i]+pG[i*nCount1+j]*pFn[i];//Hij
pJacobi[(2*k)*nCount2+l*2+1]=pG[i*nCount1+j]*pEn[i]+pB[i*nCount1+j]*pFn[i];//Nij
pJacobi[((2*k)+1)*nCount2+l*2+0]=-pJacobi[(2*k)*nCount2+l*2+1];//Jij=-Nij
pJacobi[((2*k)+1)*nCount2+l*2+1]=pJacobi[(2*k)*nCount2+l*2+0];//Lij=Hij
l++;//形成两列
}
else//i=j时
{//if(pType[j]==0) continue;//平衡节点则进入下一次循环
aii=0;bii=0;
for(int m=0;m<nCount1;m++)
{
aii+=(pG[i*nCount1+m]*pEn[m]-pB[i*nCount1+m]*pFn[m]);
bii+=(pG[i*nCount1+m]*pFn[m]+pB[i*nCount1+m]*pEn[m]);
}
pJacobi[(2*k)*nCount2+l*2+0]=bii-pB[i*nCount1+i]*pEn[i]+pG[i*nCount1+i]*pFn[i];
pJacobi[(2*k)*nCount2+l*2+1]=aii+pG[i*nCount1+i]*pEn[i]+pB[i*nCount1+i]*pFn[i];
pJacobi[((2*k)+1)*nCount2+l*2+0]=aii-pG[i*nCount1+i]*pEn[i]-pB[i*nCount1+i]*pFn[i];
pJacobi[((2*k)+1)*nCount2+l*2+1]=-bii-pB[i*nCount1+i]*pEn[i]+pG[i*nCount1+i]*pFn[i];
l++;//形成两列
}
}
k++;//形成两行
}//if(pType[i]==1)
else if(pType[i]==2)//PV节点
{
for(j=0,l=0;j<nCount1;j++)
{
if(i!=j)
{
if(pType[j]==0) continue;//平衡节点则进入下一次循环
pJacobi[(2*k)*nCount2+l*2+0]=-pB[i*nCount1+j]*pEn[i]+pG[i*nCount1+j]*pFn[i];//Hij
pJacobi[(2*k)*nCount2+l*2+1]=pG[i*nCount1+j]*pEn[i]+pB[i*nCount1+j]*pFn[i];//Nij
pJacobi[((2*k)+1)*nCount2+l*2+0]=0;//Rij=0
pJacobi[((2*k)+1)*nCount2+l*2+1]=0;//Sij=0
l++;//形成两列
}
else//i=j时
{//if(pType[j]==0) continue;//平衡节点则进入下一次循环
aii=0;bii=0;
for(int m=0;m<nCount1;m++)
{
aii+=(pG[i*nCount1+m]*pEn[m]-pB[i*nCount1+m]*pFn[m]);
bii+=(pG[i*nCount1+m]*pFn[m]+pB[i*nCount1+m]*pEn[m]);
}
pJacobi[(2*k)*nCount2+l*2+0]=bii-pB[i*nCount1+i]*pEn[i]+pG[i*nCount1+i]*pFn[i];//Hii
pJacobi[(2*k)*nCount2+l*2+1]=aii+pG[i*nCount1+i]*pEn[i]+pB[i*nCount1+i]*pFn[i];//Nii
pJacobi[((2*k)+1)*nCount2+l*2+0]=2*pFn[i];//Rii=2*pFn[i]
pJacobi[((2*k)+1)*nCount2+l*2+1]=2*pEn[i];//Sii=2*pEn[i]
l++;//形成两列
}
}
k++;//形成两行
}// else if(pType[i]==2)
}//for(i=0,k=0;i<nCount1;i++)
//**********For Test**********
CStdioFile fp;
CString szText;
/*fp.Open("c:\\test.txt", CFile::modeCreate | CFile::modeWrite);
fp.WriteString("************************************************\n");
for(i=0; i<nCount2; i++)
{
for(j=0; j<nCount2; j++)
{
if( i == j )
{
szText.Format("%6.2lf\n", pJacobi[i*nCount2+j]);
fp.WriteString(szText);
}
}
}
fp.Close();*/
}
/*
编写者：侯学勇
编写日期:2003-10-22
功能：全选主元,高斯消去法迭代求解电压修正量
输入参数：
pJacobi：指向雅克比矩阵（一维数组）的首地址；
pResult：解空间首地址；
nCount：求解维数
返回参数：
整型（1:成功，0:失败）
备注:
根据雅克比矩阵和迭代后的功率及电压平方误差量解出各节点的电压虚部和实部的修正量
输入参数：雅克比矩阵首地址，修正量首地址，修正量
总数（除平衡节点的节点数，即PQ＋PV节点数）
返回参数整型，1表示求解成功，0表示失败
*/
int CNR_PowerFlowApp::GaussJordan(NodeParam *bpNode, NRParam *bpNR)
{
int nCount = (bpNode->PQCount + bpNode->PVCount)*2;
double *pJacobi = bpNR->pJacobi;
double *pResult = bpNR->pPQ_Uef;
if(nCount<=0)
{
return 0;
}
int *js,i,j,k,is,u,v;
double d,t;
js=new int[nCount];
for (k=0; k<=nCount-1; k++)
{ d=0.0;
for (i=k; i<=nCount-1; i++)
for (j=k; j<=nCount-1; j++)
{ t=fabs(pJacobi[i*nCount+j]);
if (t>d)
{
d=t;
js[k]=j;
is=i;
}
}
if (d+1.0==1.0)
{
delete []js;
return(0);
}
if (is!=k)
{ for (j=k; j<=nCount-1; j++)
{ u=k*nCount+j; v=is*nCount+j;
t=pJacobi[u]; pJacobi[u]=pJacobi[v]; pJacobi[v]=t;
}
t=pResult[k]; pResult[k]=pResult[is]; pResult[is]=t;
}
if (js[k]!=k)
for (i=0; i<=nCount-1; i++)
{ u=i*nCount+k; v=i*nCount+js[k];
t=pJacobi[u]; pJacobi[u]=pJacobi[v]; pJacobi[v]=t;
}
t=pJacobi[k*nCount+k];
for (j=k+1; j<=nCount-1; j++)
{ u=k*nCount+j;
if (pJacobi[u]!=0.0) pJacobi[u]=pJacobi[u]/t;
}
pResult[k]=pResult[k]/t;
for (j=k+1; j<=nCount-1; j++)
{ u=k*nCount+j;
if (pJacobi[u]!=0.0)
{ for (i=0; i<=nCount-1; i++)
{ v=i*nCount+k;
if ((i!=k)&&(pJacobi[v]!=0.0))
{
is=i*nCount+j;
pJacobi[is]=pJacobi[is]-pJacobi[v]*pJacobi[u];
}
}
}
}
for (i=0; i<=nCount-1; i++)
{ u=i*nCount+k;
if ((i!=k)&&(pJacobi[u]!=0.0))
pResult[i]=pResult[i]-pJacobi[u]*pResult[k];
}
}
for (k=nCount-1; k>=0; k--)
if (k!=js[k])
{
t=pResult[k];
pResult[k]=pResult[js[k]];
pResult[js[k]]=t;
}
delete []js;
return(1);
}
/*
编写者：侯学勇
编写日期:2003-10-22
功能：检测是否收敛
输入参数：
nJacobiCount：修正量个数及雅克比矩阵阶数；
pType：节点类型数组；
返回参数：
整型（收敛返回0，不收敛返回1）
备注:
电压修正量检测，修正量最大值小于给定阈值返回0，否则返回1
*/
int CNR_PowerFlowApp::DeltaCheck(NRParam *bpNR)
{
int nJacobiCount = bpNR->nJacobiCount;
double *pPQ_Uef = bpNR->pPQ_Uef;
double duEps = bpNR->duEps;
int i = 0, k = 0 ;
double maxdelta=0;//初始化最大不平衡量
//求电压修正量最大值
for( i=0;i<nJacobiCount;i++)
{
maxdelta = maxdelta < fabs(pPQ_Uef[i]) ? fabs(pPQ_Uef[i]) : maxdelta;
}
// TRACE("本次迭代最大误差：%.5f\n",maxdelta);
//如迭代已收敛后，先看看PV节点无功是否越限，如果越限则是假收敛，需进行PV转PQ处理进入下次迭代
if( maxdelta < duEps )
{
return 0;
}
else
{
return 1;
}
}
/*
编写者：侯学勇
编写日期:2003-10-22
功能：牛拉法叠代过程
输入参数：无
返回参数：
整型（0表示已收敛，1表示未收敛或PV节点无功越界,进入下次迭代）
备注:
牛拉法叠代过程（结束条件为迭代达到最大次数，或已收敛）
*/
int CNR_PowerFlowApp::didai(NodeParam *bpNode,BranchParam *bpBranch,NRParam *bpNR)
{
int flag = 1;//迭代收敛标志
while ( ( flag == 1 )&&( RemainTimes > 0 ) )
{
TRACE("修正耗时：");
DWORD time1=::GetTickCount();
xiuzhen(bpNode, bpNR);/*电压修正及有功无功电压平方不平衡量形成*/
DWORD time2=::GetTickCount();
TRACE("%d\n",time2-time1);
TRACE("雅克比矩阵形成耗时：");
time1=::GetTickCount();
int nJocbiNumber = bpNR->nJacobiCount;
JACOBI(bpNode, bpNR);/*雅克比矩阵形成*/
double MAX = fabs(bpNR->pJacobi[0]);
double MIN = fabs(bpNR->pJacobi[0]);
for(int i=1;i<nJocbiNumber;i++)
{
if(fabs(bpNR->pJacobi[i*nJocbiNumber+i])>MAX)
{
MAX=fabs(bpNR->pJacobi[i*nJocbiNumber+i]);
//TRACE("pfJocbi[%d][%d]:%4.3lf\t\t",i,i,pfJocbi[i*nJocbiNumber+i]);
}
if(fabs(bpNR->pJacobi[i*nJocbiNumber+i])<MIN)
{
MIN=fabs(bpNR->pJacobi[i*nJocbiNumber+i]);
//TRACE("pfJocbi[%d][%d]:%4.3lf\t\t",i,i,pfJocbi[i*nJocbiNumber+i]);
}
}
time2=::GetTickCount();
TRACE("%d\n",time2-time1);
TRACE("求解耗时：");
time1=::GetTickCount();
GaussJordan(bpNode, bpNR);/*全选主元高斯消去法求解电压修正量*/
time2=::GetTickCount();
TRACE("%d\n",time2-time1);
flag=DeltaCheck(bpNR);/*检验修正量，判断是否收敛*/
RemainTimes--;
}
if( flag == 0 )
{
return 0;
}
else
{
return 1;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -