📄 flow.cpp
字号:
Lij=Hij;
Jacobian.SetJacElement(2*i+1,2*j,Jij);
Jacobian.SetJacElement(2*i,2*j+1,Nij);
Jacobian.SetJacElement(2*i+1,2*j+1,Lij);
}
else if(pNodeData[oldOrderI].NodeType==PQ && pNodeData[oldOrderJ].NodeType==PV)
{
Jij=U[i]*U[j]*(Gij*cos(Pij)+Bij*sin(Pij));
Jacobian.SetJacElement(2*i+1,2*j,Jij);
}
else if (pNodeData[oldOrderI].NodeType==PV && pNodeData[oldOrderJ].NodeType==PQ)
{
Nij=-U[i]*U[j]*(Gij*cos(Pij)+Bij*sin(Pij));
Jacobian.SetJacElement(2*i,2*j+1,Nij);
}
}
}
}
//解雅可比矩阵
Jacobian.SetJacStatus(true);
if(ITER==1)
{
Jacobian.CountElement();
Jacobian.BuildFactorTable();
}
/* 2008.3.27 Disable unnecessary IO operations.
//输出雅可比矩阵
fmidresults<<"………………第"<<ITER<<"次迭代的雅可比矩阵………………"<<std::endl;
for(i=0;i<JacRows;i++)
{
for(j=0;j<JacRows;j++)
{
double temp=Jacobian.GetJacElement(i,j);
if(fabs(temp)<SMALLEST)
continue;
fmidresults<<"("<<i+1<<","<<j+1<<") "<<temp<<" ";
}
fmidresults<<std::endl;
}
fmidresults<<std::endl;
*/
/*
#ifndef NDEBUG
fprintf(fmidresults,"\nArray of PQError before resolve\n");
for(i=0;i<NodeNum-1;i++)
fprintf(fmidresults,"%lf %lf\n",PQError[2*i],PQError[2*i+1]);
fprintf(fmidresults,"\n");
#endif
*/
Jacobian.ResolveEquation(PQError);
/*
#ifndef NDEBUG
fprintf(fmidresults,"\nArray of PQError after resolve\n");
for(i=0;i<NodeNum-1;i++)
fprintf(fmidresults,"%lf %lf\n",PQError[2*i],PQError[2*i+1]);
fprintf(fmidresults,"\n");
fclose(fmidresults);
#endif
*/
//修正电压变量
for(i=0;i<NodeNum-1;i++)
{
oldOrder=OriginalNode(i);
Phase[i]=Phase[i]-PQError[2*i];
if(pNodeData[oldOrder].NodeType==PQ)
U[i]=U[i]-PQError[2*i+1]*U[i];
}
}
delete [] PQError;
if(ITER>=MaxIterNum) return false;
return true;
}
///-------------------------------------------------------------------------------
///20008.5.29
///单独形成雅克比矩阵并保存
bool Flow::BackupJacobian()
{
if(bBackupedJacobian) return false;
return bBackupedJacobian;
}
//----------------------------计算并列补偿支路的等效注入
void Flow::CalcParalCompensBranchesPQ()
{
if(bUseParallelCompenTech && nAbnormalBranchNum>0)
{
ParallelCompensBranch * pcb;
int i;
int ni,nj;
double PhaseIJ;
for(i=0;i<nAbnormalBranchNum;i++)
{
pcb=ParalCompensBranches+i;
ni=OptimisedNode(pcb->NodeI-1);
nj=OptimisedNode(pcb->NodeJ-1);
PhaseIJ=Phase[ni]-Phase[nj];
pcb->Pi=pcb->Gc*U[ni]*U[ni]-U[ni]*U[nj]*(pcb->Gc*cos(PhaseIJ)+pcb->Bc*sin(PhaseIJ));
pcb->Qi=-pcb->Bc*U[ni]*U[ni]+U[ni]*U[nj]*(pcb->Bc*cos(PhaseIJ)-pcb->Gc*sin(PhaseIJ));
pcb->Pj=pcb->Gc*U[nj]*U[nj]-U[ni]*U[nj]*(pcb->Gc*cos(PhaseIJ)-pcb->Bc*sin(PhaseIJ));
pcb->Qj=-pcb->Bc*U[nj]*U[nj]+U[ni]*U[nj]*(pcb->Bc*cos(PhaseIJ)+pcb->Gc*sin(PhaseIJ));
}
}
}
//--------------------------------------------------------------------------------
double Flow::GetPError(int i) //获得功率误差，参数为优化后的节点编号
{
double sum=0;
double VI=U[i];
double PI=pNodeData[OriginalNode(i)].P;
double AngleIJ,G,B;
Complex Yij;
double VJ;
for(int j=0;j<NodeNum;j++)
{
// if(j==i || Adjacency.GetElement(i,j))
if(Adjacency[i][j])
{
AngleIJ=Phase[i]-Phase[j];
Yij=GetY(i,j);
G=Yij.Real;
B=Yij.Imag;
VJ=U[j];
sum+=VJ*(G*cos(AngleIJ)+B*sin(AngleIJ));
}
}
sum*=VI;
return (PI-sum);
}
///-------------------------------------------------------------
void Flow::CalcNodalPowerP()
{
int i,j;
double A,B;
Complex Yij;
double AngleIJ;
for(i=0;i<NodeNum;i++)
{
NodalPowerP[i]=0;
for(j=0;j<NodeNum;j++)
{
// if(j==i || Adjacency.GetElement(i,j))
if(Adjacency[i][j])
{
AngleIJ=Phase[i]-Phase[j];
Yij=GetY(i,j);
A=Yij.Real*cos(AngleIJ);
B=Yij.Imag*sin(AngleIJ);
NodalPowerP[i]+=U[j]*(A+B);
}
}
NodalPowerP[i]*=U[i];
}
//-------------2008.10.27，处理小阻抗和R/X比值较大的支路而增加
if(bUseParallelCompenTech && nAbnormalBranchNum>0)
{
ParallelCompensBranch * pcb;
int ni,nj;
for(i=0;i<nAbnormalBranchNum;i++)
{
pcb=ParalCompensBranches+i;
ni=OptimisedNode(pcb->NodeI-1);
nj=OptimisedNode(pcb->NodeJ-1);
NodalPowerP[ni]+=pcb->Pi;
NodalPowerP[nj]+=pcb->Pj;
}
}
/* //-------------------[3/27/2009 ZQfalcon更改] 处理BPA格式的并联导纳负荷-------------------------------------
if(NodeNum==3141)
{
for(i=0;i<NodeNum;i++)
{
j=OptimisedNode(i);
NodalPowerP[j]+=U[j]*U[j]*ShuntLoad[i].ShuntP; //累加并联有功负荷
}
}
*/ //-------------------[3/27/2009 ZQfalcon更改] 处理BPA格式的并联导纳负荷-------------------------------------
BOOST_FOREACH( DCBus const& bus, dc_buses ) {
double coef = bus.GetType() == DCBus::RECTIFIER ? +1.0 : -1.0;
double Pd = bus.GetPower();
int index = OptimisedNode(bus.GetID()-1);
NodalPowerP[index] -= coef * Pd;
}
}
double Flow::GetQError(int i) //获得功率误差，参数为优化后的节点编号
{
double sum=0;
double VI=U[i];
double QI=0;
Complex Yij;
double AngleIJ,G,B;
double VJ;
QI=pNodeData[OriginalNode(i)].Q;
for(int j=0;j<NodeNum;j++)
{
// if(j==i || Adjacency.GetElement(i,j))
if(Adjacency[i][j])
{
AngleIJ=Phase[i]-Phase[j];
Yij=GetY(i,j);
G=Yij.Real;
B=Yij.Imag;
VJ=U[j];
sum+=VJ*(G*sin(AngleIJ)-B*cos(AngleIJ));
}
}
sum*=VI;
return (QI-sum);
}
//--------------------------------------------------------------------------------
void Flow::CalcNodalPowerQ() //获得求最后发电机的无功出力，参数为优化后的节点编号
{
int i,j;
double A,B;
Complex Yij;
double AngleIJ;
for(i=0;i<NodeNum;i++)
{
NodalPowerQ[i]=0;
for(j=0;j<NodeNum;j++)
{
// if(j==i || Adjacency.GetElement(i,j))
if(Adjacency[i][j])
{
AngleIJ=Phase[i]-Phase[j];
Yij=GetY(i,j);
A=Yij.Real*sin(AngleIJ);
B=Yij.Imag*cos(AngleIJ);
NodalPowerQ[i]+=U[j]*(A-B);
}
}
NodalPowerQ[i]*=U[i];
}
//-------------2008.10.27，处理小阻抗和R/X比值较大的支路而增加
if(bUseParallelCompenTech && nAbnormalBranchNum>0)
{
ParallelCompensBranch * pcb;
int ni,nj;
for(i=0;i<nAbnormalBranchNum;i++)
{
pcb=ParalCompensBranches+i;
ni=OptimisedNode(pcb->NodeI-1);
nj=OptimisedNode(pcb->NodeJ-1);
NodalPowerQ[ni]+=pcb->Qi;
NodalPowerQ[nj]+=pcb->Qj;
}
}
/* //-------------------[3/27/2009 ZQfalcon更改] 处理BPA格式的并联导纳负荷-------------------------------------
if(NodeNum==3141)
{
for(i=0;i<NodeNum;i++)
{
j=OptimisedNode(i);
NodalPowerQ[j]-=U[j]*U[j]*ShuntLoad[i].ShuntQ; //累加并联无功负荷
}
}
*/ //-------------------[3/27/2009 ZQfalcon更改] 处理BPA格式的并联导纳负荷-------------------------------------
BOOST_FOREACH( DCBus const& bus, dc_buses ) {
double coef = bus.GetType() == DCBus::RECTIFIER ? +1.0 : -1.0;
double Qd = bus.GetPower() * std::tan(std::acos(bus.GetFactor()));
int index = OptimisedNode(bus.GetID()-1);
NodalPowerQ[index] -= coef * Qd;
}
}
//--------------------------------------------------------------------------------
double Flow::GetQOutput1(int i) //获得求最后发电机节点的综合无功出力，参数为优化后的节点编号
{
double sum=0;
double VI=U[i];
double QI=0;
double AngleIJ,G,B;
double VJ;
QI=pNodeData[OriginalNode(i)].Lq;
for(int j=0;j<NodeNum;j++)
{
// if(j==i || Adjacency.GetElement(i,j))
if(Adjacency[i][j])
{
AngleIJ=Phase[i]-Phase[j];
G=GetY(i,j).Real;
B=GetY(i,j).Imag;
VJ=U[j];
sum+=VJ*(G*sin(AngleIJ)-B*cos(AngleIJ));
}
}
sum*=VI;
return (-sum);
}
//--------------------------------------------------------------------------------
void Flow::BuildMatrixX()
{
int i;
ACBranch * tb;
int ni,nj; //优化后的节点编号
double tii,tij,tji,tjj,tc;
X.DeleteAllElements(); // 2008.05.21 Clear all non-zero elements.
for(i=0;i<BranchNum;i++)
{
tb=pBranchData+i;
ni=OptimisedNode(tb->NodeI-1);
if ( tb->BranchType!=CAP) nj=OptimisedNode(tb->NodeJ-1);
tii=GetX(ni,ni); //获得导纳矩阵原来元素的值
if(tb->BranchType!=CAP)
{
tij=GetX(ni,nj);
tji=GetX(nj,ni);
tjj=GetX(nj,nj);
}
switch(tb->BranchType) //不考虑对地电容
{
case LN:
tc=-1.0/tb->X;
tii=tii+tc;
tjj=tjj+tc;
tij=tij-tc;
tji=tji-tc;
break;
case TS:
{
double y12;
tc=tb->X;
y12=-1.0/tc;
tii=tii+y12;
tjj=tjj+y12;
tij=tij-y12;
tji=tji-y12;
}
break;
case CAP:
break;
case TSP:
break;
default:
assert(false);
}
SetX(ni,ni,tii); //设置X矩阵的值
if(tb->BranchType!=CAP){
SetX(ni,nj,tij);
SetX(nj,ni,tji);
SetX(nj,nj,tjj);}
}
//-------------2008.10.27，处理小阻抗和R/X比值较大的支路而增加的类
if(bUseParallelCompenTech && nAbnormalBranchNum>0)
{
ParallelCompensBranch * pcb;
double Y2,X;
for(i=0;i<nAbnormalBranchNum;i++)
{
pcb=ParalCompensBranches+i;
ni=OptimisedNode(pcb->NodeI-1);
nj=OptimisedNode(pcb->NodeJ-1);
tii=GetX(ni,ni); //获得导纳矩阵原来元素的值
tij=GetX(ni,nj);
tji=GetX(nj,ni);
tjj=GetX(nj,nj);
Y2=pcb->Gn*pcb->Gn+pcb->Bn*pcb->Bn;
X=-pcb->Bn/Y2;
if(pcb->BranchType==LN)
{
tc=-1.0/X;
tii=tii+tc;
tjj=tjj+tc;
tij=tij-tc;
tji=tji-tc;
}
else //变压器支路
{
double y12;
y12=-1.0/X/tb->K;
tii=tii+y12;
tjj=tjj+y12;
tij=tij-y12;
tji=tji-y12;
}
SetX(ni,ni,tii); //设置导纳矩阵的值
SetX(ni,nj,tij);
SetX(nj,ni,tji);
SetX(nj,nj,tjj);
}
}
/*
#ifndef NDEBUG
printf("\n\nX Matrix");
for(i=0;i<NodeNum;i++)
{
printf("\n");
for(int j=0;j<NodeNum;j++)
{
printf("%lf ",GetX(i,j));
}
}
#endif
*/
}
void Flow::BuildMatrixB()
{
int i;
ACBranch * pb;
int ni,nj; //优化后的节点编号
B.DeleteAllElements(); // 2008.05.21 Clear all non-zero elements.
for(i=0;i<NodeNum;i++)
B.SetElement(i,i,Y.GetElement(i,i).Imag);
for(i=0;i<BranchNum;i++)
{
pb=pBranchData+i;
ni=OptimisedNode(pb->NodeI-1);
nj=OptimisedNode(pb->NodeJ-1);
B.SetElement(ni,nj,Y.GetElement(ni,nj).Imag);
B.SetElement(nj,ni,Y.GetElement(nj,ni).Imag);
}
/* Branch * pb;
int ni,nj; //优化后的节点编号
double R,X;
double Z2;
double tii,tij,tji,tjj,tc;
for(i=0;i<BranchNum;i++)
{
pb=pBranchData+i;
ni=OptimisedNode(pb->NodeI-1);
if ( pb->BranchType!=CAP) nj=OptimisedNode(pb->NodeJ-1);
tii=B.GetElement(ni,ni); //获得导纳矩阵原来元素的值
if(pb->BranchType!=CAP)
{
tij=B.GetElement(ni,nj);
tji=B.GetElement(nj,ni);
tjj=B.GetElement(nj,nj);
R=pb->R;
X=pb->X;
Z2=R*R+X*X;
}
switch(pb->BranchType)
{
case LN:
tc=-X/Z2;
tii=tii+tc+pb->B_2;
tjj=tjj+tc+pb->B_2;
tij=tij-tc;
tji=tji-tc;
break;
case TS:
{
tc=-X/Z2;
double y12=tc/pb->K;
tii=tii+tc;
tjj=tjj+y12/pb->K;
tij=tij-y12;
tji=tji-y12;
}
break;
case CAP: //不考虑对地电容
break;
case TSP:
break;
default:
assert(false);
}
B.SetElement(ni,ni,tii); //设置X矩阵的值
if(pb->BranchType!=CAP){
B.SetElement(ni,nj,tij);
B.SetElement(nj,ni,tji);
B.SetElement(nj,nj,tjj);}
}
//建立导纳矩阵还需要加入各节点导纳
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -