liti.cpp

华中科技大学出版社出版《电力系统分析》（3版）（下册）求解11—5例题的源代码

#include<iostream>
#include<cmath>
using namespace std;

typedef double type;

//重载运算符
class complex
{
public:
	complex(){real=imag=0;}
	complex(type r,type i){real=r;imag=i;}
	void set(type r,type i);
	complex operator +(const complex &c);
	complex operator -(const complex &c);
	complex operator *(const complex &c);
	complex operator /(const complex &c);
	void operator =(const complex &c);
	bool operator !=(const complex &c);
	friend type getreal(const complex &c);
	friend type getimag(const complex &c);
	friend void print(const complex &c);
	friend complex gonge(const complex &c);
private:
	type real;
	type imag;
};
void complex::set (type r,type i)
{real=r;imag=i;}
inline complex complex::operator + (const complex &c)
{return complex(real+c.real,imag+c.imag );}
inline complex complex::operator - (const complex &c)
{return complex(real-c.real,imag-c.imag);}
inline complex complex::operator * (const complex &c)
{return complex(real*c.real -imag*c.imag ,real*c.imag +imag*c.real );}
inline complex complex::operator / (const complex &c)
{return complex((real*c.real+imag*c.imag)/(c.real*c.real+c.imag*c.imag),
				(imag*c.real-real*c.imag)/(c.real*c.real+c.imag*c.imag));}
inline void complex::operator = (const complex &c)
{real=c.real;imag=c.imag;} 
inline bool complex::operator !=(const complex &c)
{if(real !=c.real || imag != c.imag )return true;
else return false;}

type getreal (const complex &c)
{return c.real;}
type getimag (const complex &c)
{return c.imag ;} 

void print (const complex &c)
{
	if(c.imag<0)cout<<c.real <<c.imag <<'j';
	else cout<<c.real<<'+'<<c.imag<<'j';
}
complex gonge(const complex &c)
{
	return complex(c.real ,-c.imag);
}





void main()
{
	complex yi,ling;
	yi.set(1.0,0.0);
	ling.set(0.0,0.0);

    //输入数据
	complex z[5][5];
	complex y[5][5];	
	z[1][2].set(0.10,0.40);
	y[1][2].set(0.0,0.01528);
	y[2][1]=y[1][2];	
	z[1][3].set(0.0,0.2727);//直接计算得出来
	complex z130,z310;
	z130.set(0.0,-2.9967);
	z310.set(0.0,2.724);
	y[3][1]=yi/z130;
	y[1][3]=yi/z310;//??	
	z[1][4].set(0.12,0.50);
	y[1][4].set(0.0,0.01920);
	y[4][1].set(0.0,0.01920);
	z[2][4].set(0.08,0.40);
	y[2][4].set(0.0,0.01413);y[4][2].set(0.0,0.01413);

	//赋重复值
	for(int i=1;i<5;i++)
	{
		for (int j=1;j<5;j++)
		{
			if(z[i][j]!=ling)
				z[j][i]=z[i][j];
		}
	}


	complex Y[5][5];//节点导纳矩阵
	//构成自导
	for(i=1;i<5;i++)
	{
		int j=1;
		while(j<5)
		{if(z[i][j]!=ling)Y[i][i]=Y[i][i]+(yi/z[i][j]);
		Y[i][i]=Y[i][i]+y[i][j];
		j++;}
	}
	//构成互导
	for(i=1;i<5;i++)
	{
		for(int j=1;j<5;j++)
		{
			if(i==j)continue;
			else if(z[i][j]!=ling)
			{Y[i][j]=ling-(yi/z[i][j]);}
		}
	}
	
	//输出节点导纳矩阵
	/*for(i=1;i<5;i++)
	{
		for(int j=1;j<5;j++)
		{
			print(Y[i][j]);
			cout<<'\t';
		}
		cout<<endl;
	}*/

	
	//获取实部导纳矩阵G与虚部导纳矩阵B
	type G[5][5];
	type B[5][5];	
	for(i=1;i<5;i++)
	{
		for(int j=1;j<5;j++)
		{
			G[i][j]=getreal(Y[i][j]);
			B[i][j]=getimag(Y[i][j]);
		}
	}
	
	//给定节点电压初值
	//计算deltP,deltQ,deltV2
	type e[5];type f[5];
	type sP[4],sQ[3],sV2[5];

	e[1]=e[2]=1.0;e[3]=1.1;e[4]=1.05;
	f[1]=f[2]=f[3]=f[4]=0.0;
	sP[1]=-0.3;sP[2]=-0.55;sP[3]=0.5;
	sQ[1]=-0.18;sQ[2]=-0.13;
	sV2[3]=1.1*1.1;
	
	type deltP[5],deltQ[5],deltV2[5];
	for(i=1;i<5;i++)
	{
		deltP[i]=0.0;deltQ[i]=0.0;deltV2[i]=0.0;
	}

L1:	for(i=1;i<=3;i++)// 求解deltP pq节点和 pv节点
	{
		int j=1;
		type temp1=0.0,temp2=0.0;
		while(j<5)
		{
			temp1=temp1+(G[i][j]*e[j]-B[i][j]*f[j]);
			j++;
		}
		j=1;
		while(j<5)
		{
			temp2+=G[i][j]*f[j]+B[i][j]*e[j];
			j++;
		}
		deltP[i]=sP[i]-e[i]*temp1-f[i]*temp2;
	}
	for(i=1;i<3;i++)// 求解deltQ pq节点
	{
		int j=1;
		type temp1=0.0,temp2=0.0;
		while(j<5)
		{
			temp1=temp1+(G[i][j]*e[j]-B[i][j]*f[j]);
			j++;
		}
		j=1;
		while(j<5)
		{
			temp2+=G[i][j]*f[j]+B[i][j]*e[j];
			j++;
		}
		deltQ[i]=sQ[i]-f[i]*temp1+e[i]*temp2;
	}
	i=3;
	while(i==3)//求解deltV2pv接点
	{
		deltV2[i]=sV2[i]-e[i]*e[i]-f[i]*f[i];
		i++;
	}
//	定义收敛条件，检验是否收敛
	type max=0.0;
	for(i=1;i<5;i++)
	{
		//cout<<deltP[i]<<'\t'<<deltQ[i]<<'\t'<<deltV2[i]<<endl;
		if(fabs(deltP[i])>fabs(max))max=deltP[i];
		if(fabs(deltQ[i])>fabs(max))max=deltQ[i];
		if(fabs(deltV2[i])>fabs(max))max=deltV2[i];
	}
//	cout<<max<<'\n'<<'\n'<<endl;
	if(fabs(max)<1e-5)//满足收敛条件计算平衡节点功率以及全部线路功率
	{
		complex V[5];
		for(int i=1;i<=4;i++)
		{
			V[i].set(e[i],f[i]);			
		}
	
		//计算功率：
		complex temp;
		complex Sp[5];		
			temp=ling;
			for(int j=1;j<=4;j++)
			{				
			   temp=temp+gonge(Y[4][j])*gonge(V[j]);			   
			}
			Sp[4]=V[4]*temp;
		cout<<"例题11-5"<<endl;
		cout<<"已知 z12=0.1.+j0.40"<<endl;
		cout<<"y120=y210=j0.01528"<<endl;
		cout<<"z13=j0.3,k=1.1"<<endl;
		cout<<"z14=0.12+j0.50"<<endl;
		cout<<"y140=y410=j0.01920"<<endl;
		cout<<"z24=0.08+j0.40"<<endl;
		cout<<"y240=y420=j0.01413"<<endl;
		cout<<"给定1、2为PQ节点，3为PV节点，4为平衡节点"<<endl;
		cout<<"已给定P1S+jQ1S=-0.30-j0.18"<<endl;
		cout<<"P2S+jQ2S=-0.55-j0.13"<<endl;
		cout<<"P3S=0.5,V3S=1.10,V4S=1.05<0,容许误差10^-5"<<endl;
		cout<<"求解平衡节点功率:"<<endl;
		print(Sp[4]);
		cout<<endl;
	}
	else//否则计算雅各比矩阵并修正e[i]\f[i]
	{
		type J[7][7];//雅各比矩阵
		for(int m=1;m<7;m++)
		{
			for(int n=1;n<7;n++)
			{
				J[m][n]=0.0;
			}
		}	
		//主对角线矩阵块，即当i=j		
		for(int i=1;i<4;i++)
		{
			int j=1;
			type temp=0.0;
			while(j<=4)
			{
				temp+=(G[i][j]*e[j]-B[i][j]*f[j]);
				j++;
			}
			J[(2*i-1)][(2*i-1)]=-temp-G[i][i]*e[i]-B[i][i]*f[i];
			temp=0.0;j=1;
			while(j<=4)
			{
				temp+=(G[i][j]*f[j]+B[i][j]*e[j]);
				j++;
			}			
			J[(2*i-1)][2*i]=-temp+B[i][i]*e[i]-G[i][i]*f[i];
			temp=0.0;j=1;			
			if(i<3)
			{
			while(j<=4)
			{
				temp+=(G[i][j]*f[j]+B[i][j]*e[j]);
				j++;
			}
			J[(2*i)][(2*i-1)]=temp+B[i][i]*e[i]-G[i][i]*f[i];
			temp=0.0;j=1;
			while(j<=4)
			{
				temp+=(G[i][j]*e[j]-B[i][j]*f[j]);
				j++;
			}
			J[2*i][2*i]=-temp+G[i][i]*e[i]+B[i][i]*f[i];
			}
			if(i==3)
			{
				J[(2*i)][(2*i-1)]=-2*e[i];
				J[2*i][2*i]=-2*f[i];
			}
		}//end of主对角线矩阵块，即当i=j
			
		//非主对角矩阵块，即当j!=i
		for(i=1;i<4;i++)
		{
			for(int j=1;j<4;j++)
			{
				if(i==j){continue;}
				if(i<3)
				{
					J[(2*i-1)][(2*j-1)]=-G[i][j]*e[i]-B[i][j]*f[i];
					J[2*i][2*j]=-J[(2*i-1)][(2*j-1)];
					J[2*i-1][2*j]=B[i][j]*e[i]-G[i][j]*f[i];
					J[2*i][2*j-1]=J[2*i-1][2*j];
				}
				else if(i==3)
				{
					J[(2*i-1)][(2*j-1)]=-G[i][j]*e[i]-B[i][j]*f[i];
					J[(2*i-1)][(2*j)]=B[i][j]*e[i]-G[i][j]*f[i];
					J[(2*i)][(2*j)]=0.0;
					J[(2*i)][(2*j-1)]=0.0;
				}

			}
		}//end of非主对角矩阵块，即当j!=i

	
		/*for(i=1;i<7;i++)
		{
		for(int j=1;j<7;j++)
		{
			cout<<J[i][j];
			cout<<'\t';
		}
		cout<<endl;
		}*///输出雅各比矩阵


		//定义增广矩阵并求解delte[i]\deltf[i]
		type JEx[7][8];//负雅各比矩阵增广矩阵
		for(m=1;m<7;m++)
		{
			for(int n=1;n<7;n++)
			{
				JEx[m][n]=0-J[m][n];
			}
		}
		JEx[1][7]=deltP[1];
		JEx[2][7]=deltQ[1];
		JEx[3][7]=deltP[2];
		JEx[4][7]=deltQ[2];
		JEx[5][7]=deltP[3];
		JEx[6][7]=deltV2[3];

		//输出增广矩阵
		/*for(i=1;i<7;i++)
		{
		for(int j=1;j<8;j++)
		{
			cout<<JEx[i][j];
			cout<<'\t';
		}
		cout<<endl;
		}*/


		//转化三角矩阵
		for(int c1=1;c1<=6;c1++)
		{
			for(int c3=1;c3<=6-c1;c3++)//定行
			{
				if(JEx[c1+c3][c1]==0.0)continue;
				type temp;
				temp=JEx[c1+c3][c1];						
				for(int c2=c1;c2<=7;c2++)//定列
				{				
				JEx[c1+c3][c2]=JEx[c1+c3][c2]/temp*JEx[c1][c1];
				JEx[c1+c3][c2]-=JEx[c1][c2];
				}			
			}
		}
		//输出三角矩阵
	/*	for(i=1;i<7;i++)
		{
		for(int j=1;j<8;j++)
		{
			cout<<JEx[i][j];
			cout<<'\t';
		}
		cout<<endl;
		}*/

		//求解矩阵方程
		type x[7];
		for(i=6;i>=1;i--)
		{
			x[i]=JEx[i][7]/JEx[i][i];
			for(int j=1;j<i;j++)
			{
				JEx[j][7]-=x[i]*JEx[j][i];
			}
		}
	/*	for(i=1;i<=6;i++)
			cout<<x[i]<<endl;*/
	
		type delte[4];
		type deltf[4];
		delte[1]=x[1];
		delte[2]=x[3];
		delte[3]=x[5];
		deltf[1]=x[2];
		deltf[2]=x[4];
		deltf[3]=x[6];
		for(i=1;i<=3;i++)
		{
			e[i]+=delte[i];
			f[i]+=deltf[i];
		}
		goto L1;		
	}//end of else
}