📄 pf33.cpp
字号:
}
if(Flag==1)break;
}
}
/* cout<<endl<<"Load node new name: number, new node name"<<endl;
for(i=0;i<Num_Load;i++)cout<<setw(5)<<i<<setw(5)<<Load_NodeName[i] \
<<endl;*/
//负荷排序：按照新节点名（号）由小到大的顺序排序，并找出新负荷序号
//对应的旧负荷序号
for(i=0;i<Num_Load;i++)Load_No_NewtoOld[i]=i;
for(i=0;i<Num_Load-1;i++)
{
np=i;
for(j=i+1;j<Num_Load;j++)
{
if(Load_NodeName[j]<Load_NodeName[np])
{
np=j;
}
}
temp=Load_NodeName[np];
Load_NodeName[np]=Load_NodeName[i];
Load_NodeName[i]=temp;
temp=Load_No_NewtoOld[np];
Load_No_NewtoOld[np]=Load_No_NewtoOld[i];
Load_No_NewtoOld[i]=temp;
}
/* cout<<endl<<"Load sequencing: new Load number, new node name, old Load number"<<endl;
for(i=0;i<Num_Load;i++)cout<<setw(5)<<i<<setw(5)<<Load_NodeName[i] \
<<setw(5)<<Load_No_NewtoOld[i]<<endl;*/
//节点类型标志赋初值1（PQ节点）
for(i=0;i<Num_Node;i++)Node_Flag[i]=1; //节点类型标志赋初值1（PQ节点）
//从发电机节点数据中归纳出平衡节点、PV节点、PQ节点的新节点名（号）和对应的旧
//发电机序号，并对平衡节点和PV节点修正其节点类型标志
Nswing=0;
for(i=0;i<Num_Gen;i++)
{
j=Gen_No_NewtoOld[i]; //发电机节点旧顺序号
if(GGen[j].Flag==0)
{
Gen_SWNode[Nswing][0]=Gen_NodeName[i]; //发电机节点名称
Gen_SWNode[Nswing][1]=j;
Node_Flag[Gen_NodeName[i]]=0;
Nswing++;
}
else if(GGen[j].Flag==1)
{
Gen_PQNode[Num_GPQ][0]=Gen_NodeName[i]; //发电机节点名称
Gen_PQNode[Num_GPQ][1]=j;
(Num_GPQ)++;
}
else if(GGen[j].Flag==2)
{
Gen_PVNode[Num_GPV][0]=Gen_NodeName[i]; //发电机节点名称
Gen_PVNode[Num_GPV][1]=j;
Node_Flag[Gen_NodeName[i]]=2;
(Num_GPV)++;
if(Num_GPV>PVMAX)
{
cout<<"Number of PV generators is greater than PVMAX!"<<endl;
exit(7);
}
}
}
}
void Determine(int l) //导纳矩阵上三角非零元素数目限值检查子程
{
if(l>NODEMAX*NODEFACTOR)
{
cout<<"Number of none-zero elements of up_triangle"\
" is gteater than NODEMAX*NODEFACTOR!"<<endl;
exit(8);
}
}
double Y_Diag[NODEMAX][2]; //节点导纳阵的对角元:0-实部;1-虚部
double Y_UpTri[NODEMAX*NODEFACTOR][2]; //节点导纳阵上三角的非零元:0-实部;1-虚部
int Foot_Y_UpTri[NODEMAX*NODEFACTOR]; //上三角按行压缩存储的非零元的列足码
int Num_Y_UpTri[NODEMAX]; //上三角各行非零元素的个数
int No_First_Y_UpTri[NODEMAX]; //上三角各行第一个非零元素在Y_UpTri中的
//顺序号
int Foot_Y_DownTri[NODEMAX*NODEFACTOR]; //下三角按行压缩存储的非零元的列足码
int Num_Y_DownTri[NODEMAX]; //下三角各行非零元素的个数
int No_First_Y_DownTri[NODEMAX]; //下三角各行第一个非零元素在按行压缩存储
//序列中的顺序号
int No_Y_DownTri_RowtoCol[NODEMAX*NODEFACTOR]; //下三角某行非零元所对应的
//按列压缩存储序列中的序号
double Y_BigRtoX[LINEMAX][2]; //特殊处理线路并联支路导纳:0-实部；1-虚部。
int Line_BigRtoX[LINEMAX]; //特殊处理线路对应的原线路号
//形成节点导纳矩阵1子程(不包括线路充电容纳及非标准变比的影响)
void Y_Bus1(int Num_Node,
int Num_Line,
int Num_Load,
int Num_Swing)
{
int i,j,k,k_old,Flag,l;
double R,X,Z,G,B,BK;
for(i=0;i<Num_Node-Num_Swing;i++) //初始化
{
Y_Diag[i][0]=0.0;
Y_Diag[i][1]=0.0;
Num_Y_UpTri[i]=0;
}
l=0;
for(k=0;k<Num_Line;k++)
{
i=Line_NodeName[k][0]; //线路左节点
j=Line_NodeName[k][1]; //线路右节点
if(i>=Num_Node-Num_Swing)
break; //左右节点均为平衡节点，对导纳阵无影响
k_old=Line_No_NewtoOld[k]; //对应的旧线路顺序号
R=LLine[k_old].RXBK[0]; //取线路电阻值
X=LLine[k_old].RXBK[1]; //取线路电抗值
BK=LLine[k_old].RXBK[2]; //取线路容纳的半值或变压器支路的变比值
Z=R*R+X*X;
G=R/Z; //电导
B=-X/Z; //电纳
if(j>=Num_Node-Num_Swing) //左节点为普通节点，右节点为平衡节点
{
if(X>0) //(阻、)感性线路
{
if(R/X>=Resist_to_React) //大R/X线路
Y_Diag[i][1]=Y_Diag[i][1]-1.0/X;
else //小R/X线路
{
Y_Diag[i][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B;
}
}
}
else //左、右节点均为普通节点
{
Flag=0;
if(k>0&&(i==Line_NodeName[k-1][0])&&(j==Line_NodeName[k-1][1]))
Flag=1; //多回线
if(i==j) //接地支路,只需对节点i的自导纳进行处理
{
if(X>0.0) //(阻、)感性线路
{
if(R/X>=Resist_to_React) //大R/X线路
Y_Diag[i][1]=Y_Diag[i][1]-1.0/X;
else //小R/X线路
{
Y_Diag[i][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B;
}
}
}
else //对非接地支路进行处理
{
if(X==0.0) //纯电阻线路
{
Y_Diag[i][1]=Y_Diag[i][1]-1.0/R;
Y_Diag[j][1]=Y_Diag[j][1]-1.0/R;
if(Flag==0) //第一回线
{
Y_UpTri[l][0]=0.0;
Y_UpTri[l][1]=1.0/R;
Foot_Y_UpTri[l]=j;
Num_Y_UpTri[i]++;
l++;
Determine(l);
}
else //多回线
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]+1.0/R;
}
else if(X<0.0) //阻、容性线路
{
Y_Diag[i][1]=Y_Diag[i][1]+1.0/X;
Y_Diag[j][1]=Y_Diag[j][1]+1.0/X;
if(Flag==0) //第一回线
{
Y_UpTri[l][0]=0.0;
Y_UpTri[l][1]=-1.0/X;
Foot_Y_UpTri[l]=j;
Num_Y_UpTri[i]++;
l++;
Determine(l);
}
else //多回线
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-1.0/X;
}
else //阻、感性线路
{
if(R/X>=Resist_to_React) //大R/X线路
{
Y_Diag[i][1]=Y_Diag[i][1]-1.0/X;
Y_Diag[j][1]=Y_Diag[j][1]-1.0/X;
if(Flag==0) //第一回线
{
Y_UpTri[l][0]=0.0;
Y_UpTri[l][1]=1.0/X;
Foot_Y_UpTri[l]=j;
Num_Y_UpTri[i]++;
l++;
Determine(l);
}
else //多回线
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]+1.0/X;
}
else //小R/X线路
{
Y_Diag[i][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B;
Y_Diag[j][0]=Y_Diag[j][0]+G;
Y_Diag[j][1]=Y_Diag[j][1]+B;
if(Flag==0) //第一回线
{
Y_UpTri[l][0]=-G;
Y_UpTri[l][1]=-B;
Foot_Y_UpTri[l]=j;
Num_Y_UpTri[i]++;
l++;
Determine(l);
}
else //多回线
{
Y_UpTri[l-1][0]=Y_UpTri[l-1][0]-G;
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-B;
}
}
}
}
}
}
//将负荷静特性中阻抗成份计入导纳矩阵的对角元中
for(i=0;i<Num_Load;i++)
{
k=Load_No_NewtoOld[i];
if(LLoad[k].Flag==1)
{
j=Load_NodeName[i];
if(j<Num_Node-Num_Swing)
{
Y_Diag[j][0]=Y_Diag[j][0]+LLoad[k].ABC[0];
Y_Diag[j][1]=Y_Diag[j][1]-LLoad[k].ABC[1];
}
}
}
No_First_Y_UpTri[0]=0;
for(i=0;i<Num_Node-Num_Swing;i++)
No_First_Y_UpTri[i+1]=No_First_Y_UpTri[i]+Num_Y_UpTri[i];
//稀疏导纳矩阵上三角及对角元的结果输出
/* cout<<endl<<"对角元素的新节点名(号)、旧节点名(号)及元素值"<<endl;
for(i=0;i<Num_Node-Num_Swing;i++)
cout<<setw(5)<<i<<setw(5)<<Node_Name_NewtoOld[i]<<setw(10) \
<<Y_Diag[i][0]<<setw(10)<<Y_Diag[i][1]<<endl;
cout<<endl<<"上三角 "<<l<<" 个非对角元素的顺序号、列足码及元素值"<<endl;
for(k=0;k<l;k++)cout<<setw(5)<<k<<setw(5)<<Foot_Y_UpTri[k] \
<<setw(10)<<Y_UpTri[k][0]<<setw(10)<<Y_UpTri[k][1]<<endl;
cout<<endl<<"导纳矩阵上三角每行非对角元素的个数：行号及个数"<<endl;
for(i=0;i<Num_Node-Num_Swing;i++)cout<<setw(5)<<i<<setw(5) \
<<Num_Y_UpTri[i]<<endl;
cout<<endl<<"上三角每行第一个非对角元素的顺序号：行号及顺序号"<<endl;
for(i=0;i<Num_Node-Num_Swing+1;i++)cout<<setw(5)<<i<<setw(5) \
<<No_First_Y_UpTri[i]<<endl;
*/
//导纳矩阵下三角按行压缩存储时各行非零元的个数、每行第一个非零元的序号、按行压缩
//存储时非零元的列足码、某一非零元所对应的下三角阵按列压缩存储时相同非零元的序号，
//上述这些值的求取目的在于快速地处理下面迭代过程中修正方程Jacobi矩阵的按行压缩形式
//下的取值、消去和回代运算。
int Row_Down[NODEMAX]; //下三角某列非零元序号的下限工作单元
int Row_Up[NODEMAX]; //下三角某列非零元序号的上限工作单元
int li;
for(i=0;i<Num_Node-Num_Swing;i++) //导纳阵下三角各行非零元个数数组清零
Num_Y_DownTri[i]=0;
for(j=0;j<Num_Node-Num_Swing;j++) //该循环统计下三角各行非零元个数
{
for(k=No_First_Y_UpTri[j];k<No_First_Y_UpTri[j+1];k++)//针对下三角第j列非零元进行处理
{
i=Foot_Y_UpTri[k]; //行足码
Num_Y_DownTri[i]++; //下三角第i行非零元个数增1
}
Row_Down[j]=No_First_Y_UpTri[j];
Row_Up[j]=No_First_Y_UpTri[j+1];
}
No_First_Y_DownTri[0]=0;
for(i=0;i<Num_Node-Num_Swing;i++) //下三角各行第一个非零元的序号
No_First_Y_DownTri[i+1]=No_First_Y_DownTri[i]+Num_Y_DownTri[i];
for(i=1;i<Num_Node-Num_Swing;i++) //该循环确定下三角各行非零元的列足码
{
li=No_First_Y_DownTri[i]; //下三角第i行第一个非零元的序号
for(j=0;j<i;j++) //该循环搜寻下三角第0~i-1列中行号为i的非零元
{
for(k=Row_Down[j];k<Row_Up[j];k++)
{
if(Foot_Y_UpTri[k]==i)
{
Foot_Y_DownTri[li]=j; //记录i行第li个非零元的列足码
No_Y_DownTri_RowtoCol[li]=k;//记录该元素在下三角按列压缩存储序列中的序号
li++; //序号计数器增1，备下次使用
Row_Down[j]++;
}
break;
}
}
}
//稀疏导纳矩阵下三角的结果输出
/* cout<<endl<<"下三角 "<<l<<" 个非对角元素的顺序号、列足码及对应的按列压缩存储时的顺序号"<<endl;
for(k=0;k<l;k++)
cout<<setw(5)<<k<<setw(5)<<Foot_Y_DownTri[k]<<setw(5)\
<<No_Y_DownTri_RowtoCol[k]<<endl;
cout<<endl<<"导纳矩阵下三角每行非对角元素的个数：行号及个数"<<endl;
for(i=0;i<Num_Node-Num_Swing;i++)cout<<setw(5)<<i<<setw(5) \
<<Num_Y_DownTri[i]<<endl;
cout<<endl<<"下三角每行第一个非对角元素的顺序号：行号及顺序号"<<endl;
for(i=0;i<Num_Node-Num_Swing+1;i++)cout<<setw(5)<<i<<setw(5) \
<<No_First_Y_DownTri[i]<<endl; */
}
double Y_Ground[NODEMAX][2]; //部分接地导纳(不纳入最终的节点导纳阵)
//形成节点导纳矩阵2子程(包括线路充电容纳及非标准变比的影响)
void Y_Bus2(int Num_Node,
int Num_Line,
int Num_Load,
int Num_Swing,
int &Num_Line_BigRtoX)
{
int i,j,k,k_old,Flag,l;
int ii,jj; //变压器支路非标准变比侧节点、标准变比侧节点
double R,X,Z,G,B,BK;
for(i=0;i<Num_Node-Num_Swing;i++) //初始化
{
Y_Ground[i][0]=0.0;
Y_Ground[i][1]=0.0;
}
l=0;
for(k=0;k<Num_Line;k++)
{
i=Line_NodeName[k][0]; //线路左节点
j=Line_NodeName[k][1]; //线路右节点
if(i>=Num_Node-Num_Swing)
break; //左右节点均为平衡节点，对导纳阵无影响
k_old=Line_No_NewtoOld[k]; //对应的旧线路顺序号
R=LLine[k_old].RXBK[0]; //取线路电阻值
X=LLine[k_old].RXBK[1]; //取线路电抗值
BK=LLine[k_old].RXBK[2]; //取线路容纳的半值或变压器支路的变比值
Z=R*R+X*X;
G=R/Z; //电导
B=-X/Z; //电纳(数值)
if(j>=Num_Node-Num_Swing) //左节点为普通节点，右节点为平衡节点
{
if(Line_Flag[k]==0) //普通（非变压器）支路
{
Y_Diag[i][1]=Y_Diag[i][1]+BK; //累加线路充电容纳
if(X==0.0) //纯电阻线路
Y_Ground[i][0]=Y_Ground[i][0]+1.0/R;
else if(X<0.0) //阻、容性线路
{
Y_Diag[i][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B;
}
else if(R/X>=Resist_to_React) //大R/X阻、感性线路
{
Y_Ground[i][0]=Y_Ground[i][0]+G;
Y_Ground[i][1]=Y_Ground[i][1]+B+1.0/X;
}
}
else //变压器支路
{
if(Line_Flag[k]==1) //非标准变比在节点i侧
{
if(X==0.0) //纯电阻(实际不存在)
Y_Ground[i][0]=Y_Ground[i][0]+1.0/BK/BK/R;
else if(X<0.0) //三绕组变压器可能发生
{
Y_Diag[i][0]=Y_Diag[i][0]+G/BK/BK;
Y_Diag[i][1]=Y_Diag[i][1]+B/BK/BK;
}
else
{
if(R/X>=Resist_to_React) //大R/X线路
{
Y_Diag[i][1]=Y_Diag[i][1]-(1.0/BK/BK-1.0)/X;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -