📄 pf33.cpp
字号:
int kgpq; //发电机PQ节点计数临时工作单元
double R,X,Z,Ang; //线路参数及平衡节点电压相角临时单元
double yij[2],V_Temp[2],I_Temp[2]; //临时工作单元
for(i=0;i<Num_Node-Num_Swing;i++)
{
Voltage[i][0]=0.0; //电压相位的初值：0.0(弧度)
Voltage[i][1]=1.0; //电压幅值的初值：1.0
DVolt[i][0]=0.0; //电压相量变化量赋初值0.0
DVolt[i][1]=0.0;
Power_Const[i][0]=0.0; //各节点注入功率赋初值
Power_Const[i][1]=0.0;
}
if(VolIni_Flag==1)Voltage_Initial(Num_Node,Num_Swing);
for(kgpv=0;kgpv<Num_GPV;kgpv++) //发电机PV节点的电压幅值=VG
{
i=Gen_PVNode[kgpv][0];
kg_old=Gen_PVNode[kgpv][1];
Voltage[i][1]=GGen[kg_old].PQV[1];
}
for(i=0;i<Num_Line;i++) //常数电流信息求解
{
jl=Line_NodeName[i][0]; //线路左节点
jr=Line_NodeName[i][1]; //线路右节点
if(jl<Num_Node-Num_Swing&&jr>=Num_Node-Num_Swing) //jl为普通节点
{ //jr为平衡节点
k=Line_No_NewtoOld[i]; //对应的旧线路号
R=LLine[k].RXBK[0];
X=LLine[k].RXBK[1];
Z=R*R+X*X;
if(Line_Flag[i]==0)
{
R=R/Z;
X=-X/Z;
}
else
{
R=R/Z/LLine[k].RXBK[2];
X=-X/Z/LLine[k].RXBK[2];
}
yij[0]=R;
yij[1]=X; //至此求得支路导纳除以非标准变比
Flag=0;
for(k=0;k<Num_Swing;k++)
{
if(Gen_SWNode[k][0]==jr)
{
kk=Gen_SWNode[k][1];
Ang=GGen[kk].PQV[1]*Deg_to_Rad;
V_Temp[0]=GGen[kk].PQV[0]*cos(Ang);
V_Temp[1]=GGen[kk].PQV[0]*sin(Ang);
Flag=1;
}
if(Flag==1)break;
} //至此求得相应平衡节点电压的实虚部
Flag=0;
for(j=0;j<Num_Cur_Const;j++)
{
if(Node_Name_Curr_Const[j]==jl)Flag=1;
if(Flag==1)break;
}
if(Flag==0) //新增常数电流节点
{
Node_Name_Curr_Const[Num_Cur_Const]=jl;
Comp_Mul(Current_Const[Num_Cur_Const],yij,V_Temp);
Num_Cur_Const++;
if(Num_Cur_Const>SWINGMAX*NODEFACTOR)
{
cout<<"Number of constant-current nodes is"\
" gteater than SWINGMAX*NODEFACTOR!"<<endl;
exit(10);
}
}
else //该常数电流节点已经出现过
{
Comp_Mul(I_Temp,yij,V_Temp);
Current_Const[j][0]=Current_Const[j][0]+I_Temp[0];
Current_Const[j][1]=Current_Const[j][1]+I_Temp[1];
}
}
}
//输出常数电流节点数据
/* cout<<endl<<"Constant current:Number, new node name, current value"<<endl;
for(i=0;i<Num_Cur_Const;i++)cout<<setw(5)<<i<<setw(5)<< \
Node_Name_Curr_Const[i]<<setw(10)<<Current_Const[i][0]<<setw(10) \
<<Current_Const[i][1]<<endl;*/
//各节点注入功率不变部分求值处理
kgpv=0;
kgpq=0;
kl=0;
for(i=0;i<Num_Node-Num_Swing;i++)
{
if(kgpv<Num_GPV&&i==Gen_PVNode[kgpv][0]) //发电机PV节点
{
kg_old=Gen_PVNode[kgpv][1]; //发电机旧顺序号
if(kl<Num_Load&&i==Load_NodeName[kl]) //负荷节点
{
kl_old=Load_No_NewtoOld[kl]; //负荷旧顺序号
Power_Const[i][0]=Power_Const[i][0]+GGen[kg_old].PQV[0] \
-LLoad[kl_old].ABC[4]; //有功部分加入PG和C1
kl++;
}
else //非负荷节点
{
Power_Const[i][0]=Power_Const[i][0]+GGen[kg_old].PQV[0];
//有功部分加入PG
}
kgpv++;
}
else if(kgpq<Num_GPQ&&i==Gen_PQNode[kgpq][0]) //发电机PQ节点
{
kg_old=Gen_PQNode[kgpq][1]; //发电机旧顺序号
if(kl<Num_Load&&i==Load_NodeName[kl]) //负荷节点
{
kl_old=Load_No_NewtoOld[kl]; //负荷旧顺序号
Power_Const[i][0]=Power_Const[i][0]+GGen[kg_old].PQV[0] \
-LLoad[kl_old].ABC[4]; //有功部分加入PG和C1
Power_Const[i][1]=Power_Const[i][1]+GGen[kg_old].PQV[1] \
-LLoad[kl_old].ABC[5]; //无功部分加入QG和C2
kl++;
}
else //非负荷节点
{
Power_Const[i][0]=Power_Const[i][0]+GGen[kg_old].PQV[0];
//有功部分加入PG
Power_Const[i][1]=Power_Const[i][1]+GGen[kg_old].PQV[1];
//无功部分加入QG
}
kgpq++;
}
else //既非发电机PV节点，又非发电机PQ节点
{
if(kl<Num_Load&&i==Load_NodeName[kl]) //负荷节点
{
kl_old=Load_No_NewtoOld[kl]; //负荷旧顺序号
Power_Const[i][0]=Power_Const[i][0]-LLoad[kl_old].ABC[4];
//有功部分加入C1
Power_Const[i][1]=Power_Const[i][1]-LLoad[kl_old].ABC[5];
//无功部分加入C2
kl++;
}
}
}
//各节点注入功率不变部分处理结果输出：新节点名（号），有功功率，无功功率
/* cout<<endl<<"Node input constant power:new node name, real and reactive power"<<endl;
for(i=0;i<Num_Node-Num_Swing;i++)cout<<setw(5)<<i<<setw(10)<< \
Power_Const[i][0]<<setw(10)<<Power_Const[i][1]<<endl;*/
//功率失配量输出部分的标题输出
of_PQDisOut.open("pqdis.dat",ios::out,filebuf::openprot);//打开失配量输出磁盘文件
if(of_PQDisOut.fail()) //判断打开是否有错
{
cerr<<"Error opening the diskette data-file:"<<"pqdis.dat"<<endl;
exit(10);
}
of_PQDisOut<<setiosflags( ios::fixed);
of_PQDisOut<<setw(15)<<"Iterating No"<<setw(15)<<"P_Dis_Max"<<setw(6)<<"Node"\
<<setw(15)<<"Q_Dis_Max"<<setw(6)<<"Node"<<endl;
of_PQDisOut<<setw(15)<<"=============="<<setw(15)<<"=============="<<setw(6)<<"====="\
<<setw(15)<<"=============="<<setw(6)<<"====="<<endl;
cout<<setiosflags( ios::fixed);
cout<<setw(15)<<"Iterating No"<<setw(15)<<"P_Dis_Max"<<setw(6)<<"Node"\
<<setw(15)<<"Q_Dis_Max"<<setw(6)<<"Node"<<endl;
cout<<setw(15)<<"=============="<<setw(15)<<"=============="<<setw(6)<<"====="\
<<setw(15)<<"=============="<<setw(6)<<"====="<<endl;
}
void Voltage_Initial(int Num_Node,int Num_Swing)//读取电压初始化值(按内部节点号)子程
{
ifstream if_VIn("vol33.ini",ios::in,filebuf::openprot);
if(if_VIn.fail())
{
cerr<<"Error opening the voltage initial diskette data-file:Vol33.ini"<<endl;
exit(0);
}
for(int i=0;i<Num_Node-Num_Swing;i++)
if_VIn>>Voltage[i][0]>>Voltage[i][1];
if_VIn.close(); //关闭磁盘文件
}
//求节点功率失配量子程(PV节点的dQi=0)
void PowerDis_Comp(int Num_Node,
int Num_Load,
int Num_Swing,
int Num_GPV,
int Num_Line_BigRtoX,
int Num_Cur_Const,
double Power_Dis[NODEMAX][2],
double Pij_Sum[NODEMAX][2],
double DVolt[NODEMAX][2],
ofstream of_PQDisOut,
int Num_Iteration,
double &Power_Dis_Max)
{
int i,j,k,li,kl_old,kl,kgpv,ki,k1;
double V,Ang; //节点i电压幅值和相位
double VV[2]; //节点i电压实部和虚部
double V_Temp[2]; //节点电压临时工作单元（实、虚部）
double Cur_Count[2]; //注入电流统计单元（实、虚部）
double Cur_Temp[2]; //注入电流临时工作单元（实、虚部）
double G,B; //大R/X线路并联处理部分的支路导纳
double Vi,Vj; //电压幅值
double Ix,Iy; //电流（实、虚部）
double PP,QQ; //有功、无功
double CS,SN; //余弦、正弦变量
int ipmax,iqmax; //最大有功、无功失配量对应的节点号
double P_Dis_Max,Q_Dis_Max; //最大有功、无功失配量
kl=0;
kgpv=0;
ki=0;
for(i=0;i<Num_Node-Num_Swing;i++)
{
Power_Dis[i][0]=Power_Const[i][0]; //将节点i注入功率不变部分送入失配量单元
Power_Dis[i][1]=Power_Const[i][1];
Ang=Voltage[i][0]; //上次迭代后节点i电压相位和幅值
V=Voltage[i][1];
VV[0]=V*cos(Ang); //上次迭代后节点i电压实部和虚部
VV[1]=V*sin(Ang);
//下面求Pij和Qij及从节点i发出的所有Pij和Qij之和
Cur_Count[0]=0.0;
Cur_Count[1]=0.0;
for(k=No_First_Y_DownTri[i];k<No_First_Y_DownTri[i+1];k++) //下三角第i行非零元循环
{
j=Foot_Y_DownTri[k]; //下三角第i行当前非零元的列足码
V_Temp[0]=Voltage[j][1]*cos(Voltage[j][0]);
V_Temp[1]=Voltage[j][1]*sin(Voltage[j][0]);
li=No_Y_DownTri_RowtoCol[k]; //对应的按列压缩存储序列中同一非零元的序号
Comp_Mul(Cur_Temp,Y_UpTri[li],V_Temp);
Cur_Count[0]=Cur_Count[0]+Cur_Temp[0];
Cur_Count[1]=Cur_Count[1]+Cur_Temp[1];
}
Comp_Mul(Cur_Temp,Y_Diag[i],VV);
Cur_Count[0]=Cur_Count[0]+Cur_Temp[0];
Cur_Count[1]=Cur_Count[1]+Cur_Temp[1];
for(k=No_First_Y_UpTri[i];k<No_First_Y_UpTri[i+1];k++)
{
j=Foot_Y_UpTri[k];
V_Temp[0]=Voltage[j][1]*cos(Voltage[j][0]);
V_Temp[1]=Voltage[j][1]*sin(Voltage[j][0]);
Comp_Mul(Cur_Temp,Y_UpTri[k],V_Temp);
Cur_Count[0]=Cur_Count[0]+Cur_Temp[0];
Cur_Count[1]=Cur_Count[1]+Cur_Temp[1];
}
Cur_Count[1]=-Cur_Count[1]; //电流取共軛
Comp_Mul(Pij_Sum[i],VV,Cur_Count); //至此，求得从节点i发出的所有Pij和Qij之和
if(kgpv<Num_GPV&&i==Gen_PVNode[kgpv][0]) //发电机PV节点
{
if(kl<Num_Load&&i==Load_NodeName[kl]) //负荷节点
{
kl_old=Load_No_NewtoOld[kl]; //负荷旧顺序号
if(LLoad[kl_old].Flag==1) //计及负荷静特性
Power_Dis[i][0]=Power_Dis[i][0]-LLoad[kl_old].ABC[2]*V;
kl++;
}
if(ki<Num_Cur_Const&&i==Node_Name_Curr_Const[ki])
{ //含常数电流节点
Power_Dis[i][0]=Power_Dis[i][0]+V*(Current_Const[ki][0]*cos(Ang)+\
Current_Const[ki][1]*sin(Ang));
ki++;
}
Power_Dis[i][0]=Power_Dis[i][0]-Pij_Sum[i][0]-Y_Ground[i][0]*V*V;
kgpv++;
}
else //PQ节点(包括发电机、负荷及联络节点）
{
if(kl<Num_Load&&i==Load_NodeName[kl]) //负荷节点
{
kl_old=Load_No_NewtoOld[kl]; //负荷旧顺序号
if(LLoad[kl_old].Flag==1) //计及负荷静特性
{
Power_Dis[i][0]=Power_Dis[i][0]-LLoad[kl_old].ABC[2]*V;
Power_Dis[i][1]=Power_Dis[i][1]-LLoad[kl_old].ABC[3]*V;
}
kl++;
}
if(ki<Num_Cur_Const&&i==Node_Name_Curr_Const[ki])
{ //含常数电流节点
Power_Dis[i][0]=Power_Dis[i][0]+V*(Current_Const[ki][0]*cos(Ang)+\
Current_Const[ki][1]*sin(Ang));
Power_Dis[i][1]=Power_Dis[i][1]+V*(Current_Const[ki][0]*sin(Ang)-\
Current_Const[ki][1]*cos(Ang));
ki++;
}
Power_Dis[i][0]=Power_Dis[i][0]-Pij_Sum[i][0]-Y_Ground[i][0]*V*V;
Power_Dis[i][1]=Power_Dis[i][1]-Pij_Sum[i][1]+Y_Ground[i][1]*V*V;
}
}
for(k=0;k<Num_Line_BigRtoX;k++) //特殊线路产生的注入累加处理
{
k1=Line_BigRtoX[k]; //对应的新线路号
i=Line_NodeName[k1][0]; //新线路左节点号
j=Line_NodeName[k1][1]; //新线路右节点号
G=Y_BigRtoX[k][0];
B=Y_BigRtoX[k][1];
Ang=Voltage[i][0]-Voltage[j][0]; //节点i,j电压相位之差
Vi=Voltage[i][1];
Vj=Voltage[j][1];
Power_Dis[i][0]=Power_Dis[i][0]+Vi*(G*(Vj*cos(Ang)-Vi)+B*Vj*sin(Ang));
Power_Dis[j][0]=Power_Dis[j][0]+Vj*(G*(Vi*cos(Ang)-Vj)-B*Vi*sin(Ang));
if(Node_Flag[i]==1) //节点i为PQ节点
Power_Dis[i][1]=Power_Dis[i][1]+Vi*(G*Vj*sin(Ang)-B*(Vj*cos(Ang)-Vi));
if(Node_Flag[j]==1) //节点j为PQ节点
Power_Dis[j][1]=Power_Dis[j][1]-Vj*(G*Vi*sin(Ang)+B*(Vi*cos(Ang)-Vj));
}
P_Dis_Max=0.0;
Q_Dis_Max=0.0;
for(i=0;i<Num_Node-Num_Swing;i++)
{
if(fabs(Power_Dis[i][0])>P_Dis_Max)
{
P_Dis_Max=fabs(Power_Dis[i][0]);
ipmax=i;
}
if(fabs(Power_Dis[i][1])>Q_Dis_Max)
{
Q_Dis_Max=fabs(Power_Dis[i][1]);
iqmax=i;
}
}
Power_Dis_Max=__max(P_Dis_Max,Q_Dis_Max);
cout<<setw(15)<<Num_Iteration<<setw(15)<<setprecision(8)<<P_Dis_Max\
<<setw(6)<<Node_Name_NewtoOld[ipmax]\
<<setw(15)<<setprecision(8)<<Q_Dis_Max\
<<setw(6)<<Node_Name_NewtoOld[iqmax]<<endl;
of_PQDisOut<<setw(15)<<Num_Iteration<<setw(15)<<setprecision(8)<<P_Dis_Max\
<<setw(6)<<Node_Name_NewtoOld[ipmax]\
<<setw(15)<<setprecision(8)<<Q_Dis_Max\
<<setw(6)<<Node_Name_NewtoOld[iqmax]<<endl;
//节点功率一次偏差修正项
/* for(k=0;k<Num_Cur_Const;k++) //常数电流修正项处理
{
i=Node_Name_Curr_Const[k]; //节点号
if(i<Num_Node-Num_Swing)
{
Ix=Current_Const[k][0];
Iy=Current_Const[k][1];
Ang=Voltage[i][0];
PP=(-Ix*sin(Ang)+Iy*cos(Ang))*DVolt[i][0]+(Ix*cos(Ang)+Iy*sin(Ang))*DVolt[i][1];
QQ=(Ix*cos(Ang)+Iy*sin(Ang))*DVolt[i][0]+(Ix*sin(Ang)-Iy*cos(Ang))*DVolt[i][1];
Power_Dis[i][0]=Power_Dis[i][0]-PP;
Power_Dis[i][1]=Power_Dis[i][1]-QQ;
}
}
for(i=0;i<Num_Node-Num_Swing;i++)//i节点所有出线的Pij、Qij之和、接地导纳修正项处理
{
V=Voltage[i][1];
PP=Pij_Sum[i][1]*DVolt[i][0]-Pij_Sum[i][0]*DVolt[i][1];
QQ=-Pij_Sum[i][0]*DVolt[i][0]-Pij_Sum[i][1]*DVolt[i][1];
// Power_Dis[i][0]=Power_Dis[i][0]-PP/V;
// Power_Dis[i][1]=Power_Dis[i][1]-QQ/V;
Power_Dis[i][0]=Power_Dis[i][0]-PP;
Power_Dis[i][1]=Power_Dis[i][1]-QQ;
G=Y_Ground[i][0];
B=Y_Ground[i][1];
PP=-2.0*G*V*DVolt[i][1];
QQ=2.0*B*V*DVolt[i][1];
Power_Dis[i][0]=Power_Dis[i][0]-PP;
Power_Dis[i][1]=Power_Dis[i][1]-QQ;
}
for(k=0;k<Num_Load;k++) //负荷静特性之修正项处理
{
i=Load_NodeName[k];
k1=Load_No_NewtoOld[k];
if(i<Num_Node-Num_Swing&&LLoad[k1].Flag==1)
{
PP=-LLoad[k1].ABC[2]*DVolt[i][1];
QQ=-LLoad[k1].ABC[3]*DVolt[i][1];
Power_Dis[i][0]=Power_Dis[i][0]-PP;
Power_Dis[i][1]=Power_Dis[i][1]-QQ;
}
}
for(k=0;k<Num_Line_BigRtoX;k++) //特殊线路修正项处理
{
k1=Line_BigRtoX[k]; //对应的新线路号
i=Line_NodeName[k1][0]; //左节点
j=Line_NodeName[k1][1]; //右节点
Vi=Voltage[i][1]; //左节点电压幅值
Vj=Voltage[j][1]; //右节点电压幅值
Ang=Voltage[i][0]-Voltage[j][0]; //左右节点电压相位差
CS=cos(Ang);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -