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;							//列足码暂存单元
	int k;							//临时记数单元