📄 pf33.cpp
字号:
Y_Ground[i][0]=Y_Ground[i][0]+G/BK/BK;
Y_Ground[i][1]=Y_Ground[i][1]+(B+1.0/X)/BK/BK;
}
else //小R/X线路
{
Y_Diag[i][0]=Y_Diag[i][0]+(1.0/BK/BK-1.0)*G;
Y_Diag[i][1]=Y_Diag[i][1]+(1.0/BK/BK-1.0)*B;
}
}
}
else //非标准变比在节点j侧
{
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 //左、右节点均为普通节点
{
Flag=0;
if(k>0&&(i==Line_NodeName[k-1][0])&&(j==Line_NodeName[k-1][1]))
Flag=1; //多回线
if(Line_Flag[k]==0) //普通（非变压器）支路
{
if(i==j) //接地支路
{
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 //对非接地支路进行处理
{
Y_Diag[i][1]=Y_Diag[i][1]+BK; //累加线路充电容纳
Y_Diag[j][1]=Y_Diag[j][1]+BK;
if(X==0.0) //纯电阻线路
{
Line_BigRtoX[Num_Line_BigRtoX]=k; //记录线路号
Y_BigRtoX[Num_Line_BigRtoX][0]=1.0/R; //线路支路导纳
Y_BigRtoX[Num_Line_BigRtoX][1]=1.0/R;
(Num_Line_BigRtoX)++;
}
else if(X<0.0) //阻、容性线路
{
Line_BigRtoX[Num_Line_BigRtoX]=k; //记录线路号
Y_BigRtoX[Num_Line_BigRtoX][0]=G; //线路支路导纳
Y_BigRtoX[Num_Line_BigRtoX][1]=B-1.0/X;
(Num_Line_BigRtoX)++;
}
else if(R/X>=Resist_to_React) //大R/X阻、感性线路
{
Line_BigRtoX[Num_Line_BigRtoX]=k; //记录线路号
Y_BigRtoX[Num_Line_BigRtoX][0]=G; //线路支路导纳
Y_BigRtoX[Num_Line_BigRtoX][1]=B+1.0/X;
(Num_Line_BigRtoX)++;
}
if(Flag==0)l++;
}
}
else //变压器支路(此种情形不存在直接接地的支路)
{
if(Line_Flag[k]==1) //非标准变比在节点i侧
{
ii=i;
jj=j;
}
else //非标准变比在节点j侧
{
ii=j;
jj=i;
}
if(X==0.0) //纯电阻线路
{
Y_Diag[ii][1]=Y_Diag[ii][1]-(1.0/BK/BK-1.0)/R;
Y_Ground[ii][0]=Y_Ground[ii][0]+(1.0/BK-1.0)/BK/R;
Y_Ground[ii][1]=Y_Ground[ii][1]+(1.0/BK-1.0)/BK/R;
Y_Ground[jj][0]=Y_Ground[jj][0]+(1.0-1.0/BK)/R;
Y_Ground[jj][1]=Y_Ground[jj][1]+(1.0-1.0/BK)/R;
Line_BigRtoX[Num_Line_BigRtoX]=k; //记录线路号
Y_BigRtoX[Num_Line_BigRtoX][0]=1.0/R/BK; //线路支路导纳
Y_BigRtoX[Num_Line_BigRtoX][1]=1.0/R/BK;
(Num_Line_BigRtoX)++;
if(Flag==0) //第一回线
{
Y_UpTri[l][1]=Y_UpTri[l][1]+(1.0/BK-1.0)/R;
l++;
}
else //多回线
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]+(1.0/BK-1.0)/R;
}
else if(X<0.0) //阻、容性线路
{
Y_Diag[ii][1]=Y_Diag[ii][1]+(1.0/BK/BK-1.0)/X;
Y_Ground[ii][0]=Y_Ground[ii][0]+(1.0/BK-1.0)/BK*G;
Y_Ground[ii][1]=Y_Ground[ii][1]+(1.0/BK-1.0)/BK*(B-1.0/X);
Y_Ground[jj][0]=Y_Ground[jj][0]+(1.0-1.0/BK)*G;
Y_Ground[jj][1]=Y_Ground[jj][1]+(1.0-1.0/BK)*(B-1.0/X);
Line_BigRtoX[Num_Line_BigRtoX]=k; //记录线路号
Y_BigRtoX[Num_Line_BigRtoX][0]=G/BK; //线路支路导纳
Y_BigRtoX[Num_Line_BigRtoX][1]=(B-1.0/X)/BK;
(Num_Line_BigRtoX)++;
if(Flag==0) //第一回线
{
Y_UpTri[l][1]=Y_UpTri[l][1]-(1.0/BK-1.0)/X;
l++;
}
else //多回线
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-(1.0/BK-1.0)/X;
}
else //阻、感性线路
{
if(R/X>=Resist_to_React) //大R/X线路
{
Y_Diag[ii][1]=Y_Diag[ii][1]-(1.0/BK/BK-1.0)/X;
Y_Ground[ii][0]=Y_Ground[ii][0]+(1.0/BK-1.0)/BK*G;
Y_Ground[ii][1]=Y_Ground[ii][1]+(1.0/BK-1.0)/BK*(B+1.0/X);
Y_Ground[jj][0]=Y_Ground[jj][0]+(1.0-1.0/BK)*G;
Y_Ground[jj][1]=Y_Ground[jj][1]+(1.0-1.0/BK)*(B+1.0/X);
Line_BigRtoX[Num_Line_BigRtoX]=k; //记录线路号
Y_BigRtoX[Num_Line_BigRtoX][0]=G/BK; //线路支路导纳
Y_BigRtoX[Num_Line_BigRtoX][1]=(B+1.0/X)/BK;
(Num_Line_BigRtoX)++;
if(Flag==0) //第一回线
{
Y_UpTri[l][1]=Y_UpTri[l][1]+(1.0/BK-1.0)/X;
l++;
}
else //多回线
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]+(1.0/BK-1.0)/X;
}
else //小R/X线路
{
Y_Diag[ii][0]=Y_Diag[ii][0]+(1.0/BK/BK-1.0)*G;
Y_Diag[ii][1]=Y_Diag[ii][1]+(1.0/BK/BK-1.0)*B;
if(Flag==0) //第一回线
{
Y_UpTri[l][0]=Y_UpTri[l][0]-(1.0/BK-1.0)*G;
Y_UpTri[l][1]=Y_UpTri[l][1]-(1.0/BK-1.0)*B;
l++;
}
else //多回线
{
Y_UpTri[l-1][0]=Y_UpTri[l-1][0]-(1.0/BK-1.0)*G;
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-(1.0/BK-1.0)*B;
}
}
}
}
}
}
//稀疏导纳矩阵上三角及对角元的结果输出
/* 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<<"接地导纳的新节点名(号)、旧节点名(号)及元素值"<<endl;
for(i=0;i<Num_Node-Num_Swing;i++)
cout<<setw(5)<<i<<setw(5)<<Node_Name_NewtoOld[i]<<setw(10) \
<<Y_Ground[i][0]<<setw(10)<<Y_Ground[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;
*/
}
//复数A和B相加得C:极坐标形式（0-相位、1-幅值）
void Comp_Plus(double C[2],double A[2],double B[2])
{
double an1,an2,r,x;
an1=A[0];
an2=B[0];
r=A[1]*cos(an1)+B[1]*cos(an2);
x=A[1]*sin(an1)+B[1]*sin(an2);
C[0]=atan2(x,r);
C[1]=sqrt(r*r+x*x);
}
//复数A和B相减得C:极坐标形式（0-相位、1-幅值）
void Comp_Minus(double C[2],double A[2],double B[2])
{
double an1,an2,r,x;
an1=A[0];
an2=B[0];
r=A[1]*cos(an1)-B[1]*cos(an2);
x=A[1]*sin(an1)-B[1]*sin(an2);
C[0]=atan2(x,r);
C[1]=sqrt(r*r+x*x);
}
//复数A和B相乘得C:直角坐标形式（0-实部、1-需部）
void Comp_Mul(double C[2],double A[2],double B[2])
{
C[0]=A[0]*B[0]-A[1]*B[1];
C[1]=A[0]*B[1]+A[1]*B[0];
}
//复数A和B相除得C:直角坐标形式（0-实部、1-需部）
void Comp_Div(double C[2],double A[2],double B[2])
{
double t[2];
double tt;
tt=B[0]*B[0]+B[1]*B[1];
if(tt==0.0)
{
cout<<"Divided by zero!"<<endl;
exit(0);
}
else
{
t[0]=B[0];
t[1]=-B[1];
Comp_Mul(C,A,t);
C[0]=C[0]/tt;
C[1]=C[1]/tt;
}
}
double Fact_Diag[NODEMAX][2]; //因子表对角元素
double Fact_UpTri[NODEMAX*NODEFACTOR][2]; //因子表上三角非零元素
int Foot_Fact_UpTri[NODEMAX*NODEFACTOR][2]; //因子表上三角非零元素列足码
int Num_Fact_UpTri[NODEMAX][2]; //因子表上三角各行非零非对角
//元素的个数
//形成节点导纳矩阵因子表
void Factor_Table(int Num_Node,int Num_Swing,int Num_GPV,int IterFlag)
{
int i; //因子表正在形成的所在行号
int im; //因子表正在消去的所在行号
int j; //列足码暂存单元
int k; //临时记数单元
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节点计数临时工作单元
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -