📄 pf30.cpp
字号:
int ix; //因子表上三角元素地址（序号）计数
double Y_Work[NODEMAX]; //工作数组
double Temp1,Temp2; //临时工作单元
int kgpv;
for(i=0;i<Num_Node-Num_Swing;i++)
{
if(IterFlag==1&&Node_Flag[i]==2) //无功迭代对应的PV节点
{
Num_Fact_UpTri[i][IterFlag]=0;
Fact_Diag[i][IterFlag]=0.0;
}
else
{
for(k=i+1;k<Num_Node-Num_Swing;k++)
Y_Work[k]=0.0;
Y_Work[i]=Y_Diag[i][1];
for(k=No_First_Y_UpTri[i];k<No_First_Y_UpTri[i+1];k++)
{
j=Foot_Y_UpTri[k];
Y_Work[j]=Y_UpTri[k][1];
}
if(IterFlag==1)
{
for(kgpv=0;kgpv<Num_GPV;kgpv++)
{
j=Gen_PVNode[kgpv][0]; //PV节点号
Y_Work[j]=0.0;
}
}
ix=0;
for(im=0;im<i;im++)
{
for(k=0;k<Num_Fact_UpTri[im][IterFlag];k++)
{
if(Foot_Fact_UpTri[ix][IterFlag]!=i)
ix++;
else
{
Temp1=Fact_UpTri[ix][IterFlag]\
/Fact_Diag[im][IterFlag];
for(k=k;k<Num_Fact_UpTri[im][IterFlag];k++)
{
j=Foot_Fact_UpTri[ix][IterFlag];
Temp2=Temp1*Fact_UpTri[ix][IterFlag];
Y_Work[j]=Y_Work[j]-Temp2;
ix++;
}
}
}
}
Fact_Diag[i][IterFlag]=1.0/Y_Work[i];
Temp1=Fact_Diag[i][IterFlag];
k=0;
for(j=i+1;j<Num_Node-Num_Swing;j++)
{
if(fabs(Y_Work[j])!=0.0)
{
Fact_UpTri[ix][IterFlag]=Y_Work[j]*Temp1;
Foot_Fact_UpTri[ix][IterFlag]=j;
k++;
ix++;
}
}
Num_Fact_UpTri[i][IterFlag]=k;
}
}
//因子表结果输出
/* cout<<endl<<"因子表对角元素：序号（与节点号同）和元素值"<<endl;
for(i=0;i<Num_Node-Num_Swing;i++)
cout<<setw(5)<<i<<setw(15)<<Fact_Diag[i][IterFlag]<<endl;
cout<<endl<<"因子表上三角非零元素：序号、列足码和元素值"<<endl;
for(i=0;i<ix;i++)
cout<<setw(5)<<i<<setw(5)<<Foot_Fact_UpTri[i][IterFlag]\
<<setw(15)<<Fact_UpTri[i][IterFlag]<<endl;
*/
}
//方程AX=t求解
void Equ_Calculation(int Num_Node,int Num_Swing,
double Power_Dis_Correct[NODEMAX],int IterFlag)
{
int i,j,k,ix; //参见Factor_Table子程说明
double Temp1,Temp2; //临时工作单元
ix=0;
for(i=0;i<Num_Node-Num_Swing;i++) //前代运算开始
{
Temp1=Power_Dis_Correct[i]; //将右端量送入临时工作单元
for(k=0;k<Num_Fact_UpTri[i][IterFlag];k++)
{
j=Foot_Fact_UpTri[ix][IterFlag];
Temp2=Temp1*Fact_UpTri[ix][IterFlag];
Power_Dis_Correct[j]=Power_Dis_Correct[j]-Temp2;
ix++;
}
Power_Dis_Correct[i]=Temp1*Fact_Diag[i][IterFlag];
}
for(i=Num_Node-Num_Swing-1;i>=0;i--) //回代运算开始
{
Temp1=Power_Dis_Correct[i];
for(k=0;k<Num_Fact_UpTri[i][IterFlag];k++)
{
ix--;
j=Foot_Fact_UpTri[ix][IterFlag];
Temp2=Power_Dis_Correct[j]*Fact_UpTri[ix][IterFlag];
Temp1=Temp1-Temp2;
}
Power_Dis_Correct[i]=Temp1; //最终得未知数的解
}
}
double Voltage[NODEMAX][2]; //节点电压:0-相位角;
//1-幅值。
double Current_Const[SWINGMAX*NODEFACTOR][2]; //常数电流:0-实部;
//1-虚部。
int Node_Name_Curr_Const[SWINGMAX*NODEFACTOR]; //常数电流的节点名(号)
double Power_Const[NODEMAX][2]; //各节点注入功率不变部
//分:0-实部;1-虚部。
void Voltage_Initial(int Num_Node,int Num_Swing);
//初始化子程
void Initial(int Num_Node,int Num_Line,int Num_Load,int Num_Swing,
int Num_GPV,int Num_GPQ,int &Num_Cur_Const,
double DVolt[NODEMAX][2],ofstream &of_PQDisOut,
int VolIni_Flag)
{
int i,j,jl,jr,k,kk;
int Flag;
int kg_old,kl_old; //发电机、负荷旧顺序号工作单元
int kl; //负荷计数临时工作单元
int kgpv; //发电机PV节点计数临时工作单元
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 "\
"> 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("vol30.ini",ios::in,filebuf::openprot);
if(if_VIn.fail())
{
cerr<<"Error opening the voltage initial diskette data-file:"\
<<"Vol30.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_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 Ix,Iy; //电流（实、虚部）
double PP,QQ; //有功、无功
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);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -