📄 nflow_sub.cpp
字号:
// netbus_->PLAST_LN2[j]=(1000.0/1.732)*sqrt(netbus_->W_LN2[j]*netbus_->W_LN2[j]+netbus_->R_LN2[j]*netbus_->R_LN2[j])/netbus_->KV_BS[netbus_->I$BS_ND[netparm_->I$ND_LN[i]-1]-1];
// }
// }
// if(netbus_->KV_BS[netbus_->I$BS_ND[netparm_->Z$ND_LN[i]-1]-1]>20.001){
// if(fabs(netbus_->W_LN2[k])>0.01||fabs(netbus_->R_LN2[k])>0.01){
// netbus_->PLAST_LN2[k]=(1000.0/1.732)*sqrt(netbus_->W_LN2[k]*netbus_->W_LN2[k]+netbus_->R_LN2[k]*netbus_->R_LN2[k])/netbus_->KV_BS[netbus_->I$BS_ND[netparm_->Z$ND_LN[i]-1]-1];
// }
// }
// }
//
// for(i=0;i<netmom_->LV$XF;i++) {
// j=netparm_->I$XF2_XF[i]-1;
// k=netparm_->Z$XF2_XF[i]-1;
// netbus_->PLAST_XF2[j]=0.0;
// netbus_->PLAST_XF2[k]=0.0;
// if(netbus_->KV_BS[netbus_->I$BS_ND[netparm_->I$ND_XF[i]-1]-1]>2.001){
// if(fabs(netbus_->W_XF2[j])>0.01||fabs(netbus_->R_XF2[j])>0.01){
// netbus_->PLAST_XF2[j]=(1000.0/1.732)*sqrt(netbus_->W_XF2[j]*netbus_->W_XF2[j]+netbus_->R_XF2[j]*netbus_->R_XF2[j])/netbus_->KV_BS[netbus_->I$BS_ND[netparm_->I$ND_XF[i]-1]-1];
// }
// }
// if(netbus_->KV_BS[netbus_->I$BS_ND[netparm_->Z$ND_XF[i]-1]-1]>2.001){
// if(fabs(netbus_->W_XF2[k])>0.01||fabs(netbus_->R_XF2[k])>0.01){
// netbus_->PLAST_XF2[k]=(1000.0/1.732)*sqrt(netbus_->W_XF2[k]*netbus_->W_XF2[k]+netbus_->R_XF2[k]*netbus_->R_XF2[k])/netbus_->KV_BS[netbus_->I$BS_ND[netparm_->Z$ND_XF[i]-1]-1];
// }
// }
// }
//
return RET_FA;
}
void pf_setdef()
{
SBS=100.0;
IT=0;//ITER_NUM
eps_set=0.0005;//0.000005;//收敛偏差
eps_max=5000.0;//发散计算偏差
FLATSTRT=1;
DYN_FLAG=1;//1;//动态潮流,0:静态潮流
OPT_FLAG=0;
SYS_DIM=2;
NEWT_FLAG=1;
}
int pf_businjections(int NOWISLAND)
{
pf_form_bsinjection(NOWISLAND);//BS注入,有名
pf_form_businjection(NOWISLAND);//BUS注入,标幺
pf_set_initvalue(NOWISLAND);//设置计算值初始值
pf_form_busvolt(NOWISLAND);//根据调压设置的固定电压设置参考节点和PV节点的V的计算初始值
return RET_OK;
}
int pf_form_bsinjection(int NOWISLAND)//BS注入,有名
{
int I;
// printf("begin ps_form_bsinjection()\n");
// printf("NOWISLAND=%d\n",NOWISLAND);
for(I=0;I<lv_bs;I++)
{
if(map_bs[I].AC_island==NOWISLAND)
{
map_bs[I].p=0.0;
map_bs[I].q=0.0;
}
}
// print_debug_info(16);
// print_debug_info(17);
for(I=0;I<lv_map_ACun;I++)
{
// printf("map_ACun[%d]:ibs==%d,p==%f,q==%f\n",I,map_ACun[I].ibs,map_ACun[I].p,map_ACun[I].q);
if(map_ACun[I].ibs>0)
{
// printf("map_ACun[%d].AC_island==%d\n",I,map_bs[map_ACun[I].ibs-1].AC_island);
if(map_bs[map_ACun[I].ibs-1].AC_island==NOWISLAND)
{
map_bs[map_ACun[I].ibs-1].p=map_bs[map_ACun[I].ibs-1].p+map_ACun[I].p;//老程序为+CW_UN[i]??
map_bs[map_ACun[I].ibs-1].q=map_bs[map_ACun[I].ibs-1].q+map_ACun[I].q;//
}
}
}
// printf("\n\n\n\n");
for(I=0;I<lv_map_ACld;I++)
{
// printf("map_ACld[%d]:ibs==%d,p==%f,q==%f\n\n",I,map_ACun[I].ibs,map_ACld[I].p,map_ACld[I].q);
if(map_ACld[I].ibs>-1)
{
// printf("map_ACld[%d].AC_island==%d\n",I,map_bs[map_ACld[I].ibs-1].AC_island);
if(map_bs[map_ACld[I].ibs-1].AC_island==NOWISLAND)
{
map_bs[map_ACld[I].ibs-1].p=map_bs[map_ACld[I].ibs-1].p-map_ACld[I].p;
map_bs[map_ACld[I].ibs-1].q=map_bs[map_ACld[I].ibs-1].q-map_ACld[I].q;
}
}
}
//虚负荷XLD和电容RC
return RET_OK;
}
int pf_form_businjection(int NOWISLAND)//BUS注入,标幺
{
int I;
for(I=0;I<lv_bus;I++)
{
map_bus[I].Pi=0.0;
map_bus[I].Qi=0.0;
}
for(I=0;I<lv_bs;I++)
{
if(map_bs[I].AC_island==NOWISLAND)
{
map_bus[map_bs[I].ibus].Pi=map_bus[map_bs[I].ibus].Pi+map_bs[I].p/SBS;//为什么不是直接==?
map_bus[map_bs[I].ibus].Qi=map_bus[map_bs[I].ibus].Qi+map_bs[I].q/SBS;
}
}
return RET_OK;
}
int pf_set_initvalue(int NOWISLAND)//设置计算值初始值
{
int I;
if(FLATSTRT==1)
{//平值启动,取初始值
for(I=0;I<lv_bus;I++)
{
map_bus[I].v_pu=1.0;
map_bus[I].a_pu=0.0;
}
//ADD BY DHW 07.07.05 BEGIN
for(I=0;I<lv_map_ACun;I++)
{
printf("map_ACun[%d]:ibs==%d,p==%f,q==%f\n",I,map_ACun[I].ibs,map_ACun[I].p,map_ACun[I].q);
if(map_ACun[I].vset>0.5&&map_ACun[I].vset<1.5)
{
map_bus[map_bs[map_ACun[I].ibs-1].ibus].v_pu=map_ACun[I].vset;
printf("v_set==%f\n",map_ACun[I].vset);
}
}
//ADD BY DHW 07.07 END
}
else if(FLATSTRT==0)
{//取上次计算值,作为初始值开始计算,需要补充内容
}
return RET_OK;
}
int pf_form_busvolt(int NOWISLAND)//根据调压设置的固定电压设置参考节点和PV节点的V的计算初始值
{
//根据调压设置的固定电压设置参考节点和PV节点的V的计算初始值
int I;
//平衡节点,以后调试
// printf("slack node: map_ACun[%d].vset==%f, vt_no==%d, vbase==%f \n\n",ref_bs_un,map_ACun[ref_bs_un].vset,map_ACun[ref_bs_un].vt_no,map_vt[map_ACun[ref_bs_un].vt_no-1].vbase);
// if(map_vt[map_ACun[ref_bs_un].vt_no-1].vbase>0.00001)
// {
// map_bus[lv_bus-1].v_pu=map_ACun[ref_bs_un].vset/map_vt[map_ACun[ref_bs_un].vt_no-1].vbase;
// printf("111:\n");
// }
if(map_bus[lv_bus-1].v_pu<0.9||map_bus[lv_bus-1].v_pu>1.1)
{
map_bus[lv_bus].v_pu=1.0;
}
//PV,对于此处，以后还要调试!!!
for(I=0;I<NOV;I++)
{
// if(lv_pv_rel_un[I]>0&&map_vt[map_ACun[pv_rel_un[I][0]].vt_no-1].vbase>0.00001)
// {
// map_bus[map_bs[REGS[I]].ibus].v_pu=map_ACun[pv_rel_un[I][0]].vset/map_vt[map_ACun[pv_rel_un[I][0]].vt_no-1].vbase;
// printf("22222222222222222,map_bus[%d].v_pu==%f,,vbase==%f,\n",map_bs[REGS[I]].ibus,map_bus[map_bs[REGS[I]].ibus].v_pu,map_vt[map_ACun[pv_rel_un[I][0]].vt_no-1].vbase);
// }
//
if(map_bus[map_bs[REGS[I]].ibus].v_pu<0.9||map_bus[map_bs[REGS[I]].ibus].v_pu>1.1)
{
map_bus[map_bs[REGS[I]].ibus].v_pu=1.0;
}
}
return RET_OK;
}
void pf_clear_iter(int NOWISLAND)
{
int I;
IT=0;
for(I=0;I<MAX_ITER;I++)
{
MXP_ITER[I]=0.0;//有功
MXQ_ITER[I]=0.0;//无功
MXPBUS_ITER[I]=0;
MXQBUS_ITER[I]=0;
}
}
int cjpq_newton(int NOWISLAND)
{
//基本原理见电力系统计算P542
//C[]:功率偏差,dPi和dQi
int I,J,JSQ;
float EI,VI,A8,B8,VJ,AN,C8,D8,AI,BI,CI,DI;
for(I=0;I<lv_bus;I++)
{
C[2*I]=0.0;
C[2*I+1]=0.0;
B9[2*I]=0.0;
B9[2*I+1]=0.0;
}
//PQ==>PV==>VO顺序排列节点
for(I=0;I<lv_bus;I++)
{
// printf("*********BUS_NO==%d**********************\n",I);
EI=map_bus[I].v_pu;//vi
VI=EI*EI;//vi*vi
A8=Ym[I].g.val*VI;//A8=Gii*vi*vi
B8=Ym[I].b.val*VI;//B8=Bii*vi*vi
// printf("G%d,%d==%f,B%d,%d==%f\n",I,I,Ym[I].g.val,I,I,Ym[I].b.val);
// printf("A8:Gii*vi*vi==%f, B8:Bii*vi*vi==%f\n\n",A8,B8);
C[2*I]=C[2*I]+A8;//对角元产生的Pi:Gii*vi*vi,节点功率
C[2*I+1]=C[2*I+1]-B8;//对角元产生的Qi:-Bii*vi*vi
// printf("C[%d]=%f, C[%d]==%f \n",2*I,C[2*I],2*I+1,C[2*I+1]);
matrix_jacb_newton[I].DJ_CJPQ[0]=B8;//Bii*vi*vi,DJ_CJPQ:放对角元值
matrix_jacb_newton[I].DJ_CJPQ[1]=B8;//Bii*vi*vi
matrix_jacb_newton[I].DJF_CJPQ[0]=-A8;//-Gii*vi*vi,DJF_CJPQ:放对角元值
matrix_jacb_newton[I].DJF_CJPQ[1]=A8;//Gii*vi*vi
// printf("matrix_jacb_newton[%d].DJ_CJPQ[0]==%f, matrix_jacb_newton[%d].DJ_CJPQ[1]=%f\n",I,matrix_jacb_newton[I].DJ_CJPQ[0],I,matrix_jacb_newton[I].DJ_CJPQ[1]);
// printf("matrix_jacb_newton[%d].DJF_CJPQ[0]==%f, matrix_jacb_newton[%d].DJF_CJPQ[1]=%f\n",I,matrix_jacb_newton[I].DJF_CJPQ[0],I,matrix_jacb_newton[I].DJF_CJPQ[1]);
for(JSQ=0;JSQ<Ym[I].lv_nz;JSQ++)
{
J=Ym[I].nz[JSQ].j;//列指
// printf("I==%d,J==%d\n",I,J);
// printf("G%d,%d=%f, B%d,%d=%f \n",I,J,Ym[I].nz[JSQ].g.val,I,J,Ym[I].nz[JSQ].b.val);
VJ=EI*map_bus[J].v_pu;//VJ=vi*vj
AN=map_bus[I].a_pu-map_bus[J].a_pu;//节点i和j的角度差: Oij
C8=VJ*Ym[I].nz[JSQ].g.val;//C8=vi*vj*Gij,
D8=VJ*Ym[I].nz[JSQ].b.val;//D8=vi*vj*Bij
A8=AN*(1.0-AN*AN/6.0);//sinOij
B8=1.0-AN*AN*(1.0-AN*AN/12.0)/2.0;//cosOij
AI=B8*C8;//AI=cosOij*(vi*vj*Gij)
BI=A8*D8;//BI=sinOij*(vi*vj*Bij)
CI=A8*C8;//CI=sinOij*(vi*vj*Gij)
DI=B8*D8;//DI=cosOij*(vi*vj*Bij)
// printf("AI=cosOij*(vi*vj*Gij)==%f, BI=sinOij*(vi*vj*Bij)==%f, CI=sinOij*(vi*vj*Gij)==%f, cosOij*(vi*vj*Bij)==%f \n",AI,BI,CI,DI);
// printf("C[%d]=%f, C[%d]==%f, C[%d]==%f, C[%d]==%f \n",2*I,C[2*I],2*I+1,C[2*I+1],2*J,C[2*J],2*J+1,C[2*J+1]);
C[2*I]=C[2*I]+AI+BI;//Pi,本端节点有功:Oij
C[2*I+1]=C[2*I+1]+CI-DI;//Qi,本端节点无功:Oij
C[2*J]=C[2*J]+AI-BI;//Pj,对端节点有功:Oij=-Oji
C[2*J+1]=C[2*J+1]-CI-DI;//Qj,对端节点无功:Oij=-Oji
// printf("C[%d]=%f, C[%d]==%f, C[%d]==%f, C[%d]==%f \n",2*I,C[2*I],2*I+1,C[2*I+1],2*J,C[2*J],2*J+1,C[2*J+1]);
matrix_jacb_newton[I].EJU_CJPQ[JSQ*2]=-CI+DI;//Hij, -sinOij*(vi*vj*Gij)+cosOij*(vi*vj*Bij)，有了Hij可以推导出Lij
matrix_jacb_newton[I].EJU_CJPQ[JSQ*2+1]=-AI-BI;//Nij, dQi/dVj,-cosOij*(vi*vj*Gij)-sinOij*(vi*vj*Bij),有了Nij可以得到Jij
// printf("H%d,%d=-sinOij*(vi*vj*Gij)+cosOij*(vi*vj*Bij)==%f, N%d,%d=-cosOij*(vi*vj*Gij)-sinOij*(vi*vj*Bij)==%f\n",I,J,matrix_jacb_newton[I].EJU_CJPQ[JSQ*2],I,J,matrix_jacb_newton[I].EJU_CJPQ[JSQ*2+1]);
matrix_jacb_newton[I].EJL_CJPQ[JSQ*2]=CI+DI;//Hji,取-Oij,表示对角转置的值(下三角)，sinOij*(vi*vj*Gij)+cosOij*(vi*vj*Bij)=-(sinOji*(vi*vj*Gij))+cosOji*(vi*vj*Bij)
matrix_jacb_newton[I].EJL_CJPQ[JSQ*2+1]=AI-BI;//Jji,取-Oij,表示对角转置的值(下三角)，cosOij*(vi*vj*Gij)-sinOij*(vi*vj*Bij)=cosOji*(vi*vj*Gij)+sinOji*(vi*vj*Bij)
// printf("H%d,%d=sinOij*(vi*vj*Gij)+cosOij*(vi*vj*Bij)==%f, J%d,%d=cosOij*(vi*vj*Gij)-sinOij*(vi*vj*Bij)==%f\n",J,I,matrix_jacb_newton[I].EJL_CJPQ[JSQ*2],J,I,matrix_jacb_newton[I].EJL_CJPQ[JSQ*2+1]);
}
// printf("C[%d]=%f, C[%d]=%f\n",2*I,C[2*I],2*I+1,C[2*I+1]);
matrix_jacb_newton[I].DJ_CJPQ[0]=matrix_jacb_newton[I].DJ_CJPQ[0]+C[2*I+1];//Hii,Bii*vi*vi+Qi
matrix_jacb_newton[I].DJ_CJPQ[1]=matrix_jacb_newton[I].DJ_CJPQ[1]-C[2*I+1];//Lii,Bii*vi*vi-Qi
matrix_jacb_newton[I].DJF_CJPQ[0]=matrix_jacb_newton[I].DJF_CJPQ[0]-C[2*I];//Nii,-Gii*vi*vi-Pi
matrix_jacb_newton[I].DJF_CJPQ[1]=matrix_jacb_newton[I].DJF_CJPQ[1]-C[2*I];//Jii,Gii*vi*vi-Pi
// printf("matrix_jacb_newton[%d].DJ_CJPQ[0]==%f, matrix_jacb_newton[%d].DJ_CJPQ[1]=%f\n",I,matrix_jacb_newton[I].DJ_CJPQ[0],I,matrix_jacb_newton[I].DJ_CJPQ[1]);
// printf("matrix_jacb_newton[%d].DJF_CJPQ[0]==%f, matrix_jacb_newton[%d].DJF_CJPQ[1]=%f\n",I,matrix_jacb_newton[I].DJF_CJPQ[0],I,matrix_jacb_newton[I].DJF_CJPQ[1]);
// printf("*******************************\n\n");
}
for(I=0;I<lv_bus;I++)
{
B9[2*I]=C[2*I];//Pi,节点结算功率
B9[2*I+1]=C[2*I+1];//Qi,节点计算功率
// printf("map_bus[%d].Pi=%f, P%d=%f\n",I,map_bus[I].Pi,I,C[2*I]);
// printf("map_bus[%d].Qi=%f, Q%d=%f\n",I,map_bus[I].Qi,I,C[2*I+1]);
C[2*I]=map_bus[I].Pi-C[2*I];//功率偏差,DPi
C[2*I+1]=map_bus[I].Qi-C[2*I+1];//功率偏差,DQi
//DEBUG
// if(fabs(B9[2*I])>0.00001&&fabs(B9[2*I+1])>0.00001&&(C[2*I]/B9[2*I]>1||C[2*I+1]/B9[2*I+1]>1))
// {
// printf("111:node_dp[%d]=%f, node_dq[%d]=%f\n",I,C[2*I]/B9[2*I],I,C[2*I+1]/B9[2*I+1]);
// printf("111:dp[%d]=%f, dq[%d]=%f\n",I,C[2*I],I,C[2*I+1]);
// printf("\n");
// }
//DEBUG
// printf("\n");
}
find_eps_max(NOWISLAND);
return RET_OK;
}
int find_eps_max(int NOWISLAND)
{
int I;
// printf("找最大功率偏差节点!\n");
if(DYN_FLAG==0)
{//非动态潮流
MXP_ITER[IT]=0.0;
for(I=0;I<lv_bus-1;I++)
{//PQ和PV节点计算P偏差,不考虑平衡节点
if(fabs(C[2*I])>fabs(MXP_ITER[IT]))
{
MXP_ITER[IT]=C[2*I];//MXW_ITER:计算的最大有功偏差
MXPBUS_ITER[IT]=I+1;//最大有功偏差的节点
}
}
// printf("max_dp==%f,bus_no==%d\n",MXP_ITER[IT],MXPBUS_ITER[IT]);
MXQ_ITER[IT]=0.0;
for(I=0;I<lv_bus-NOV-1;I++)
{//PQ节点考虑Q偏差,只考虑PQ节点,而PV和VO节点不计算dq
if(fabs(C[2*I+1])>fabs(MXQ_ITER[IT]))
{
MXQ_ITER[IT]=C[2*I+1];//MXW_ITER:最大Q偏差
MXQBUS_ITER[IT]=I+1;
}
}
// printf("max_dq==%f,bus_no==%d\n",MXQ_ITER[IT],MXQBUS_ITER[IT]);
}
return RET_OK;
}
int fact_newton(int NOWISLAN)
{
int i,MX,MA,J,K,I,JSQ,J1;
int flag;
float AI,BI;
float B1[2000],B2[2000];
for(I=0;I<2*lv_bus;I++)
{
matrix_fact_newton[I].UD_FACT=0.0;
matrix_fact_newton[I].lv_nz=0;
//DEBUG
for(J=0;J<2*lv_bus;J++)
{
matrix_fact_newton[I].nz[J].j=-1;
matrix_fact_newton[I].nz[J].UU_FACT=0.0;
matrix_fact_newton[I].nz[J].UL_FACT=0.0;
}
//DEBUG
}
//节点号由小到大：按PQ==>PV==>VO排序
//MA:初始设置MA=1
MA=1;//P/Q计算标志
// MX=0;//表示节点号(行)
for(K=0;K<2*lv_bus;K++)
{//循环整个雅可比矩阵
flag=0;
// printf("\n\n");
// printf("****************K=%d************\n",K);
MX=K/2;
if((fabs(Ym[MX].g.val)>0.00001)||(fabs(Ym[MX].b.val)>0.00001))
{
if(MX==lv_bus-1)
{ //平衡VO节点
if(MA==1)
{//P
MA=0;
}
else if(MA!=1)
{//Q
// if(typ_bs[map_bus[MX].ibs]!=BS_TYPE_PQ)
if(map_bs[map_bus[MX].ibs].bs_type!=BS_TYPE_PQ)
{
// MX=MX+1;
MA=1;
}
else
{
MA=1;
flag=1;
}
}
}
else if(MX!=lv_bus-1)
{
//非平衡VO节点
flag=1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -