📄 flow.cpp
字号:
}
}
//-------------2004.10.27，处理小阻抗和R/X比值较大的支路而增加的类
if(bUseParallelCompenTech && nAbnormalBranchNum>0)
{
ParallelCompensBranch * pcb;
Complex Yij;
for(i=0;i<nAbnormalBranchNum;i++)
{
pcb=ParalCompensBranches+i;
ni=OptimisedNode(pcb->NodeI-1);
nj=OptimisedNode(pcb->NodeJ-1);
tii=GetY(ni,ni); //获得导纳矩阵原来元素的值
tij=GetY(ni,nj);
tji=GetY(nj,ni);
tjj=GetY(nj,nj);
Yij=Complex(pcb->Gn,pcb->Bn);
if(pcb->BranchType==LN)
{
tii=tii+Yij;
tjj=tjj+Yij;
tii.Imag+=pcb->B_2;
tjj.Imag+=pcb->B_2;
tij=tij-Yij;
tji=tji-Yij;
}
else //变压器支路
{
Complex y10,y12,y20;
y12=Yij/tb->K;
y20=y12/(tb->K);
y10=Yij;
tii=tii+y20;
tjj=tjj+y10;
tij=tij-y12;
tji=tji-y12;
}
SetY(ni,ni,tii); //设置导纳矩阵的值
SetY(ni,nj,tij);
SetY(nj,ni,tji);
SetY(nj,nj,tjj);
}
}
/*
FILE * AdmitFile=fopen("Admittance.dat","w");
fprintf(AdmitFile,"\n\n节点导纳阵如下：\n");
for(i=0;i<NodeNum;i++)
{
fprintf(AdmitFile,"\n");
for(int j=0;j<NodeNum;j++)
{
if(j>0 && (j%5==0))fprintf(AdmitFile,"\n");
fprintf(AdmitFile,"%lf+j%lf ",GetY(i,j).Real,GetY(i,j).Imag);
}
fprintf(AdmitFile,"\n");
}
fprintf(AdmitFile,"\n\n");
fclose(AdmitFile);
*/
}
//--------------------------------------------------------------------------------
inline Complex Flow::GetY(int Row, int Col)
{
return Y.GetElement(Row,Col);
}
//--------------------------------------------------------------------------------
inline void Flow::SetY(int Row, int Col, Complex Value)
{
Y.SetElement(Row,Col,Value);
}
//--------------------------------------------------------------------------------
void Flow::PrepareACLoadFlow( Method method, bool refresh )
{
if (refresh) {
SwitchPallelCompenTech(false);
NodeSerialOptimize1();
StatisAbnormalBranches();
BuildMatrixY();
if ( method == PQ_BX_METHOD || method == PQ_XB_METHOD ) {
BuildMatrixB();
BuildMatrixX();
}
}
return;
}
//--------------------------------------------------------------------------------
bool Flow::NR1() //牛拉法解潮流,极坐标形式
{
JacNR Jacobian;
int * PQNode=new int[SumOfPQ]; //PQ节点编号表
int i,j;
int nsl=OptimisedNode(SlackNodeOrder-1);
Phase[nsl]=pNodeData[SlackNodeOrder-1].Fa*PI/180; //平衡节点角度
for(i=0;i<NodeNum;i++) //设置节点电压和相位初值
{
int ti;
ti=OriginalNode(i);
if(pNodeData[ti].NodeType==PQ)
U[i]=1.0; //设置PQ电压幅值为1
else
U[i]=pNodeData[ti].DesireVolt; //设置PV节点和平衡节点电压幅值
Phase[i]=Phase[nsl]; //设置节点电压相位
}
int pqcount=0; //用于统计PQ节点个数
for(i=0;i<NodeNum;i++) //获得PQ节点优化后的编号列表
{
int fi;
fi=OriginalNode(i); //获得未优化的节点编号
if(pNodeData[fi].NodeType==PQ)
{
PQNode[pqcount]=i;
pqcount++;
}
}
bPFconverge=false;
int ITER; //迭代次数
double MaxError; //max mismatch
int MaxErrNode; //最大误差点
double * PQError;
int oldOrder;
double Gij,Bij;
double Pij;
double temp;
Complex Yij;
int JacRows;
JacRows=NodeNum-1+SumOfPQ;
PQError=new double[JacRows]; //用于存放功率误差
Jacobian.Initialize(JacRows);
for(ITER=1;ITER<=MaxIterNum;ITER++)
{
Jacobian.SetJacStatus(false); //设置状态
MaxError=0.0;
CalcNodalPowerP();
for(i=0;i<NodeNum-1;i++) //注意，平衡节点号在最后
{
oldOrder=OriginalNode(i);
PQError[i]=pNodeData[oldOrder].P-NodalPowerP[i];
if(fabs(PQError[i])>MaxError)
{
MaxError=fabs(PQError[i]);
MaxErrNode=oldOrder;
}
}
/*
#ifndef NDEBUG
printf("\nThe Error of P\n");
for(i=0;i<NodeNum-1;i++)
printf("%lf ",PQError[i]);
#endif
*/
CalcNodalPowerQ();
for(i=0;i<SumOfPQ;i++)
{
oldOrder=OriginalNode(PQNode[i]);
PQError[i+NodeNum-1]=pNodeData[oldOrder].Q-NodalPowerQ[PQNode[i]];
if(fabs(PQError[i+NodeNum-1]) > MaxError)
{
MaxError=fabs(PQError[i+NodeNum-1]);
MaxErrNode=oldOrder;
}
}
/*
#ifndef NDEBUG
printf("\nThe ERROR of Q\n");
for(i=0;i<SumOfPQ;i++)
printf("%lf ",PQError[i+NodeNum-1]);
#endif
*/
//#ifndef NDEBUG
std::cout<<setw(4)<<ITER<<setw(14)<<MaxError<<std::endl;
// getch();
//#endif
if(MaxError<eps)
{
bPFconverge=true;
break;
}
//形成雅可比矩阵
//形成H元素
for(i=0;i<NodeNum-1;i++)
{
for(j=0;j<NodeNum-1;j++)
{
if(i==j)
{
temp=U[i]*U[i]*GetY(i,i).Imag;
temp+=NodalPowerQ[i];
Jacobian.SetJacElement(i,i,temp);
/* temp=0.0;
int tempj;
for(tempj=0;tempj<NodeNum;tempj++)
{
if(tempj!=i && Adjacency.GetElement(i,tempj))
{
Yij=GetY(i,tempj);
Gij=Yij.Real;
Bij=Yij.Imag;
Pij=Phase[i]-Phase[tempj];
temp+=U[tempj]*(Gij*sin(Pij)-Bij*cos(Pij));
}
}
temp=temp*U[i];
Jacobian.SetJacElement(i,j,temp);
*/
}
// else if(Adjacency.GetElement(i,j))
else if(Adjacency[i][j])
{
Yij=GetY(i,j);
Gij=Yij.Real;
Bij=Yij.Imag;
Pij=Phase[i]-Phase[j];
temp=-U[i]*U[j]*(Gij*sin(Pij)-Bij*cos(Pij));
Jacobian.SetJacElement(i,j,temp);
}
}
}
//形成N元素
for(i=0;i<NodeNum-1;i++)
{
for(int tj=0;tj<SumOfPQ;tj++)
{
j=PQNode[tj];
temp=0.0;
// if(j!=i && Adjacency.GetElement(i,j))
if(j!=i && Adjacency[i][j])
{
Yij=GetY(i,j);
Gij=Yij.Real;
Bij=Yij.Imag;
Pij=Phase[i]-Phase[j];
temp=-U[i]*U[j]*(Gij*cos(Pij)+Bij*sin(Pij));
Jacobian.SetJacElement(i,tj+NodeNum-1,temp);
}
else if(i==j)
{
double Gii=GetY(i,i).Real;
temp=-U[i]*U[i]*Gii-pNodeData[OriginalNode(i)].P;
Jacobian.SetJacElement(i,tj+NodeNum-1,temp);
}
}
}
//形成J元素
int ti;
for(ti=0;ti<SumOfPQ;ti++)
{
i=PQNode[ti];
for(j=0;j<NodeNum-1;j++)
{
// if(j!=i && Adjacency.GetElement(i,j))
if(j!=i && Adjacency[i][j])
{
Yij=GetY(i,j);
Gij=Yij.Real;
Bij=Yij.Imag;
Pij=Phase[i]-Phase[j];
temp=U[i]*U[j]*(Gij*cos(Pij)+Bij*sin(Pij));
Jacobian.SetJacElement(ti+NodeNum-1,j,temp);
}
else if(i==j)
{
double Gii=GetY(i,i).Real;
temp=U[i]*U[i]*Gii-pNodeData[OriginalNode(i)].P;
Jacobian.SetJacElement(ti+NodeNum-1,j,temp);
}
}
}
//形成L元素
for(ti=0;ti<SumOfPQ;ti++)
{
i=PQNode[ti];
int tj;
for(tj=0;tj<SumOfPQ;tj++)
{
j=PQNode[tj];
// if(j!=i && Adjacency.GetElement(i,j))
if(j!=i && Adjacency[i][j])
{
Yij=GetY(i,j);
Gij=Yij.Real;
Bij=Yij.Imag;
Pij=Phase[i]-Phase[j];
temp=-U[i]*U[j]*(Gij*sin(Pij)-Bij*cos(Pij));
Jacobian.SetJacElement(ti+NodeNum-1,tj+NodeNum-1,temp);
}
else if(i==j)
{
double Bii=GetY(i,i).Imag;
temp=U[i]*U[i]*Bii-pNodeData[OriginalNode(i)].Q;
Jacobian.SetJacElement(ti+NodeNum-1,tj+NodeNum-1,temp);
}
}
}
//输出雅可比矩阵
/*
FILE * fmidresults=fopen("MideResults.dat","a");
for(i=0;i<JacRows;i++)
{
for(j=0;j<JacRows;j++)
fprintf(fmidresults,"%12.6f",Jacobian.GetJacElement(i,j));
fprintf(fmidresults,"\n");
}
fprintf(fmidresults,"\n");
fclose(fmidresults);
*/
//解雅可比矩阵
Jacobian.SetJacStatus(true);
if(ITER==1)
{
Jacobian.CountElement();
Jacobian.BuildFactorTable();
}
/*
#ifndef NDEBUG
printf("\nArray of PQError before resolve\n");
for(i=0;i<NodeNum-1+SumOfPQ;i++)
printf("%lf ",PQError[i]);
#endif
*/
Jacobian.ResolveEquation(PQError);
/*
#ifndef NDEBUG
printf("\nArray of PQError after resolve\n");
for(i=0;i<NodeNum-1+SumOfPQ;i++)
printf("%lf ",PQError[i]);
#endif
*/
//修正电压变量
for(i=0;i<NodeNum-1;i++)
{
Phase[i]=Phase[i]-PQError[i];
}
for(i=0;i<SumOfPQ;i++)
{
int ti=PQNode[i];
U[ti]=U[ti]-PQError[i+NodeNum-1]*U[ti];
}
}
delete [] PQNode;
delete [] PQError;
return true;
}
//--------------------------------------------------------------------------------
//---------------2008.10.28-------------------------------------------------------
//将雅可比中的元素按dP,dQ的顺序排放
bool Flow::NR2() //牛拉法解潮流,极坐标形式
{
#ifndef NDEBUG
std::cout << "\n迭代次数 最大偏差\n";
#endif
JacNR Jacobian;
int i,j;
int nsl=OptimisedNode(SlackNodeOrder-1);
Phase[nsl]=pNodeData[SlackNodeOrder-1].Fa*PI/180; //平衡节点角度
for(i=0;i<NodeNum;i++) //设置节点电压和相位初值
{
int ti;
ti=OriginalNode(i);
if(pNodeData[ti].NodeType==PQ)
U[i]=1.0; //设置PQ电压幅值为1
else
U[i]=pNodeData[ti].DesireVolt; //设置PV节点和平衡节点电压幅值
Phase[i]=Phase[nsl]; //设置节点电压相位
}
bPFconverge=false;
int ITER; //迭代次数
double MaxError; //max mismatch
int MaxErrNode; //最大误差点
double * PQError;
int oldOrder;
double Gij,Bij,Gii,Bii;
double Pij;
Complex Yij;
double Hij,Nij,Jij,Lij;
int JacRows;
JacRows=(NodeNum-1)*2;
PQError=new double[JacRows]; //用于存放功率误差
Jacobian.Initialize(JacRows);
for(ITER=1;ITER<=MaxIterNum;ITER++)
{
Jacobian.SetJacStatus(false); //设置状态
MaxError=0.0;
CalcNodalPowerP();
CalcNodalPowerQ();
for(i=0;i<NodeNum-1;i++) //注意，平衡节点号在最后
{
oldOrder=OriginalNode(i);
PQError[2*i]=pNodeData[oldOrder].P-NodalPowerP[i];
if(fabs(PQError[2*i])>MaxError)
{
MaxError=fabs(PQError[2*i]);
MaxErrNode=oldOrder;
}
if(pNodeData[oldOrder].NodeType==PQ)
{
PQError[2*i+1]=pNodeData[oldOrder].Q-NodalPowerQ[i];
if(fabs(PQError[2*i+1])>MaxError)
{
MaxError=fabs(PQError[2*i+1]);
MaxErrNode=oldOrder;
}
}
else
PQError[2*i+1]=0;
}
/*
#ifndef NDEBUG
printf("\nThe Error of P\n");
for(i=0;i<NodeNum-1;i++)
printf("%lf ",PQError[i]);
#endif
*/
/*
#ifndef NDEBUG
printf("\nThe ERROR of Q\n");
for(i=0;i<SumOfPQ;i++)
printf("%lf ",PQError[i+NodeNum-1]);
#endif
*/
#ifndef NDEBUG
std::cout<<setw(4)<<ITER<<setw(14)<<MaxError<<std::endl;
#endif
if(MaxError<eps)
{
bPFconverge=true;
break;
}
//形成雅可比矩阵
for(i=0;i<NodeNum-1;i++)
{
for(j=0;j<NodeNum-1;j++)
{
if(i==j)
{
Yij=GetY(i,i);
Gii=Yij.Real;
Bii=Yij.Imag;
oldOrder=OriginalNode(i);
Hij=U[i]*U[i]*Bii+NodalPowerQ[i];
Jacobian.SetJacElement(2*i,2*i,Hij);
if(pNodeData[oldOrder].NodeType==PQ)
{
Jij=U[i]*U[i]*Gii-NodalPowerP[i]; // pNodeData[oldOrder].P;
Nij=-U[i]*U[i]*Gii-NodalPowerP[i]; // pNodeData[oldOrder].P;
Lij=U[i]*U[i]*Bii-NodalPowerQ[i]; // pNodeData[oldOrder].Q;
Jacobian.SetJacElement(2*i+1,2*i,Jij);
Jacobian.SetJacElement(2*i,2*i+1,Nij);
Jacobian.SetJacElement(2*i+1,2*i+1,Lij);
}
else
{
Jacobian.SetJacElement(2*i+1,2*i+1,1.0);
}
}
// else if(Adjacency.GetElement(i,j))
else if(Adjacency[i][j])
{
int oldOrderI,oldOrderJ;
oldOrderI=OriginalNode(i);
oldOrderJ=OriginalNode(j);
Yij=GetY(i,j);
Gij=Yij.Real;
Bij=Yij.Imag;
Pij=Phase[i]-Phase[j];
Hij=-U[i]*U[j]*(Gij*sin(Pij)-Bij*cos(Pij));
Jacobian.SetJacElement(2*i,2*j,Hij);
if(pNodeData[oldOrderI].NodeType==PQ && pNodeData[oldOrderJ].NodeType==PQ)
{
Jij=U[i]*U[j]*(Gij*cos(Pij)+Bij*sin(Pij));
Nij=-Jij;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -