?? pf30.cpp
字號:
{
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;*/
//從發電機節點數據中歸納出平衡節點、PV節點、PQ節點的新節點名（號）和對
//應的舊發電機序號，并對平衡節點和PV節點修正其節點類型標志
for(i=0;i<Num_Node;i++)Node_Flag[i]=1; //節點類型賦初值1(PQ節點)
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<<"PV Generators Number > PVMAX!"<<endl;exit(7);
}
}
}
}
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]; //下三角某行非零元所對
//應的按列壓縮存儲序列
//中的序號
//形成節點導納矩陣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 X,B; //線路參數工作單元
l=0;
for(i=0;i<Num_Node-Num_Swing;i++) //初始化
{
Y_Diag[i][1]=0.0;
Num_Y_UpTri[i]=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]; //對應的舊線路順序號
X=LLine[k_old].RXBK[1]; //取線路電抗值
B=-1.0/X; //不計線路電阻后的線路支路電納
if(j>=Num_Node-Num_Swing) //左為普通節點，右為平衡節點
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; //多回線
Y_Diag[i][1]=Y_Diag[i][1]+B;
if(i!=j) //非接地支路
{
Y_Diag[j][1]=Y_Diag[j][1]+B;
if(Flag==0) //第一回線
{
Y_UpTri[l][1]=-B;
Foot_Y_UpTri[l]=j;
Num_Y_UpTri[i]++;
l++;
if(l>NODEMAX*NODEFACTOR)
{
cout<<"Number of none-zero elements of "\
"up_triangle > NODEMAX*NODEFACTOR!"<<endl;
exit(8);
}
}
else //多回線
{
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-B;
}
}
}
}
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; */
}
//形成節點導納矩陣2(包括線路充電容納及非標準變比的影響)子程
void Y_Bus2(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;
l=0;
for(i=0;i<Num_Node-Num_Swing;i++) //初始化
{
Y_Diag[i][0]=0.0;
Y_Diag[i][1]=0.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][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B+BK;
}
else if(Line_Flag[k]==1) //非標準變比在左側節點
{
Y_Diag[i][0]=Y_Diag[i][0]+1.0/BK/BK*G;
Y_Diag[i][1]=Y_Diag[i][1]+1.0/BK/BK*B;
}
else if(Line_Flag[k]==2) //非標準變比在右側節點
{
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) //接地支路(變壓器支路不直接接地)
{
Y_Diag[i][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B+BK;
}
else //非接地支路
{
if(Line_Flag[k]==0) //普通支路
{
Y_Diag[i][0]=Y_Diag[i][0]+G;
Y_Diag[i][1]=Y_Diag[i][1]+B+BK;
Y_Diag[j][0]=Y_Diag[j][0]+G;
Y_Diag[j][1]=Y_Diag[j][1]+B+BK;
if(Flag==0) //第一回線
{
Y_UpTri[l][0]=-G;
Y_UpTri[l][1]=-B;
l++;
}
else //多回線
{
Y_UpTri[l][0]=Y_UpTri[l][0]-G;
Y_UpTri[l][1]=Y_UpTri[l][1]-B;
}
}
else if(Line_Flag[k]==1) //非標準變比在左側節點
{
Y_Diag[i][0]=Y_Diag[i][0]+1.0/BK/BK*G;
Y_Diag[i][1]=Y_Diag[i][1]+1.0/BK/BK*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]=-1.0/BK*G;
Y_UpTri[l][1]=-1.0/BK*B;
l++;
}
else //多回線
{
Y_UpTri[l-1][0]=Y_UpTri[l-1][0]-1.0/BK*G;
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-1.0/BK*B;
}
}
else //非標準變比在右側節點
{
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]+1.0/BK/BK*G;
Y_Diag[j][1]=Y_Diag[j][1]+1.0/BK/BK*B;
if(Flag==0) //第一回線
{
Y_UpTri[l][0]=-1.0/BK*G;
Y_UpTri[l][1]=-1.0/BK*B;
l++;
}
else //多回線
{
Y_UpTri[l-1][0]=Y_UpTri[l-1][0]-1.0/BK*G;
Y_UpTri[l-1][1]=Y_UpTri[l-1][1]-1.0/BK*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];
}
}
}
//稀疏導納矩陣上三角及對角元的結果輸出
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;
*/
}
//復數A和B相乘得C（直角坐標形式）
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（直角坐標形式）
void Comp_Div(double C[2],double A[2],double B[2])
{
double t[2],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]; //因子表對角元素:0-有功因
//子表;1-無功因子表。
double Fact_UpTri[NODEMAX*NODEFACTOR][2]; //因子表上三角非零元素:0-
//有功因子表;1-無功因子表。
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; //列足碼暫存單元
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -