nlflow_1.java

牛顿-拉夫逊法潮流计算程序

package Powerflow;
import BaseCalculate.*;

/**
 * 牛顿——拉夫逊法潮流计算程序
 * @author Administrator
 */
public class NLflow_1 {
	Plural[][] Y;//节点导纳矩阵
	Plural[] U  ;//各节点电压
	Plural[] v;//未设置平衡节点时，输出电压空值
	Mpqv[] PQU;//各节点的状态，有功，无功或电压值
	int n;//节点个数
	int s=0;//平衡节点号
	double[][] J;//雅克比矩阵
	double[] dPQU;//不平衡量
	double[] dU;//电压修正量
	/**
	 * 牛顿——拉夫逊法潮流计算主程序
	 * @param args 参数YY是一个Plural（复数）型的二维数组，存放节点导纳矩阵，
	               要求Y中的元素必须指向一个plural（复数）型的对象；
	 * @param args 参数UU是一个Plural（复数）型的一维数组，存放各节点电压
	 * @param args 参数PQU是一个Mpqv型的一维数组，
	               其中Mpqv型中包含各节点的状态，有功，无功或电压值
	 * @return 返回值是一个Plural（复数）型一维数组（各节点电压）
	 */
	public   Plural[] flow(Plural[][] YY,Plural[] UU,Mpqv[] PQU1) {
		Y=YY;
		PQU=PQU1;
		
		this.n=Y.length;
		U=new Plural[n];
		for(int i=0;i<n;i++){     
			U[i]=UU[i].Clone();   //电压赋值
			/*if(PQU[i].mark==0){   //寻找平衡节点号
				this.s=i;
			}*/
		}
		/*if(s==0){                //未设置平衡节点时，输出电压空值
			return v;
		}*/
		int isCal=0;
		for(int i=0;i<n;i++){     
			if(PQU[i].mark==1){//寻找平衡节点号
				this.s=i;
				isCal++;
			}
		}
		if(isCal!=1){                //未设置平衡节点时，输出电压空值
			return null;
		}
	
		for(int i=0;i<10;i++){ 
			J=J();//求雅克比矩阵
			dPQU=dPQU();//
			dU=Mat.multi(Mat.inv(J),dPQU);//最多循环次数
			int sl=convergence(0.0001); //判断电压是否收敛
			U=New_U();
			if(sl==1){            //若收敛，则强制跳出循环
				break;
			}
		}
		return U;
	}
	/**
	 * 求节点注入电流的实部
	 * @param 参数i是一个int型的数（节点号）
	 * @return 返回值是一个double型的数（节点注入电流的实部）
	 */
	public  double aii(int i){
		double aii=0;
		for(int k=0;k<n;k++){
		aii=aii+Y[i][k].a*U[k].a-Y[i][k].b*U[k].b;
		}
		return aii;
	}
	/**
	 * 节点注入电流的虚部
	 * @param 参数i是一个int型的数（节点号）
	 * @return 返回值是一个double型的数（节点注入电流的虚部）
	 */
    public  double bii(int i){
		double bii=0;
		for(int k=0;k<n;k++){
		bii=bii+Y[i][k].a*U[k].b+Y[i][k].b*U[k].a;
		}
		return bii;
    }
    /**
     * 求雅可比矩阵中的元素Hij
     * @param 参数i是一个int型的数（行号）
     * @param 参数j是一个int型的数（列号）
     * @return 返回值是一个double型的数（雅可比矩阵中的元素Hij）
     */
	public double Hij(int i,int j){
		double Hij=0;
		if(i!=j){
			Hij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b;
		}
		else{
			Hij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b+bii(i);
		}
		
		return Hij;
	}
	/**
	 * 求雅可比矩阵中的元素Nij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Nij）
	 */
	public double Nij(int i,int j){
		double Nij=0;
		if(i!=j){
			Nij=Y[i][j].a*U[i].a+Y[i][j].b*U[i].b;
		}
		else{
			Nij=Y[i][j].a*U[i].a+Y[i][j].b*U[i].b+aii(i);
		}
		return Nij;
	}
	/**
	 * 求雅可比矩阵中的元素JRij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Jij或Rij）
	 */
	public double JRij(int i,int j){
		double JRij=0;
		if(i!=j){
			if(PQU[i].mark==1){
			    JRij=-Y[i][j].a*U[i].a-Y[i][j].b*U[i].b;  //PQ节点时的计算公式
			}
			else{
				JRij=0;  //PV节点时的计算公式
			}
		}
		else{
			if(PQU[i].mark==1){
			    JRij=-Y[i][j].a*U[i].a-Y[i][j].b*U[i].b+aii(i);  //PQ节点时的计算公式
			} 
			else{
				JRij=2*U[i].b;   //PV节点时的计算公式
			}
		}
		return JRij;
	}
	/**
	 * 求雅可比矩阵中的元素LSij
	 * @param 参数i是一个int型的数（行号）
	 * @param 参数j是一个int型的数（列号）
	 * @return 返回值是一个double型的数（雅可比矩阵中的元素Lij或Sij）
	 */
	public double LSij(int i,int j){
		double LSij=0;
		if(i!=j){
			if(PQU[i].mark==1){
			    LSij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b;  //PQ节点时的计算公式
			}
			else{
				LSij=0;  //PV节点时的计算公式
			}
		}
		else{
			if(PQU[i].mark==1){
			    LSij=-Y[i][j].b*U[i].a+Y[i][j].a*U[i].b-bii(i);  //PQ节点时的计算公式
			}
			else{
				LSij=2*U[i].a;  //PV节点时的计算公式
			}
		}
		return LSij;
	}
	/**
	 * 求注入有功的不平衡量
	 * @param 参数i是一个int型的数（节点号）
	 * @return  返回值是一个double型的数（注入有功的不平衡量）
	 */
	public double dP(int i){
		double dP=0;
		for (int j=0;j<n;j++){
			dP+=U[i].a*(Y[i][j].a*U[j].a-Y[i][j].b*U[j].b)+U[i].b*(Y[i][j].a*U[j].b+Y[i][j].b*U[j].a);
		}
		dP=PQU[i].p-dP;
		return dP;
	}
	/**
	 * 求注入无功或节点电压平方的不平衡量
	 * @param 参数i是一个int型的数（节点号）
	 * @return  返回值是一个double型的数（注入无功或节点电压平方的不平衡量）
	 */
	public double dQU2(int i){
		double dQU2=0;
		if(PQU[i].mark==1){
			for (int j=0;j<n;j++){         //PQ节点时的计算公式
				dQU2+=U[i].b*(Y[i][j].a*U[j].a-Y[i][j].b*U[j].b)-U[i].a*(Y[i][j].a*U[j].b+Y[i][j].b*U[j].a);
			}
			dQU2=PQU[i].qv-dQU2;
		}
		else{
			dQU2=PQU[i].qv*PQU[i].qv-U[i].a*U[i].a-U[i].b*U[i].b;  //PV节点时的计算公式
		}
		return dQU2;
	}
	/**
	 * 求雅可比矩阵
	 * @return 返回值是一个double型的二维数组（雅可比矩阵）
	 */
	public double[][] J(){
		double[][] J=new double[2*n-2][2*n-2];
		for(int i=0;i<s;i++){     //求平衡节点号之前的各节点的2x2阶子阵
			for(int j=0;j<s;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2][j*2]=Hij(i,j);
				    J[i*2][j*2+1]=Nij(i,j);
					J[i*2+1][j*2]=JRij(i,j);
					J[i*2+1][j*2+1]=LSij(i,j);
				 }
				else{
					J[i*2][j*2]=0;
					J[i*2][j*2+1]=0;
					J[i*2+1][j*2]=0;
					J[i*2+1][j*2+1]=0;
				}
			}
			for(int j=s+1;j<n-1;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2][j*2-2]=Hij(i,j);
				    J[i*2][j*2-1]=Nij(i,j);
			        J[i*2+1][j*2-2]=JRij(i,j);
				    J[i*2+1][j*2-1]=LSij(i,j);
				}
				else{
					J[i*2][j*2-2]=0;
					J[i*2][j*2-1]=0;
					J[i*2+1][j*2-2]=0;
					J[i*2+1][j*2-1]=0;
				}
			}
		}
		for(int i=s+1;i<n;i++){     //求平衡节点号之后的各节点的2x2阶子阵
			for(int j=0;j<s;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2-2][j*2]=Hij(i,j);
				    J[i*2-2][j*2+1]=Nij(i,j);
				    J[i*2-1][j*2]=JRij(i,j);
				    J[i*2-1][j*2+1]=LSij(i,j);
				}
				else{
					J[i*2-2][j*2]=0;
					J[i*2-2][j*2+1]=0;
					J[i*2-1][j*2]=0;
					J[i*2-1][j*2+1]=0;
				}
			}
			for(int j=s+1;j<n;j++){
				if(Y[i][j].a!=0||Y[i][j].b!=0){
				    J[i*2-2][j*2-2]=Hij(i,j);
				    J[i*2-2][j*2-1]=Nij(i,j);
				    J[i*2-1][j*2-2]=JRij(i,j);
					J[i*2-1][j*2-1]=LSij(i,j);
				}
				else{
					J[i*2-2][j*2-2]=0;
					J[i*2-2][j*2-1]=0;
					J[i*2-1][j*2-2]=0;
					J[i*2-1][j*2-1]=0;

				}
			}
		}
		return J;
	}
	/**
	 * 求不平衡量矩阵
	 * @return 返回值是一个double型的一维数组（不平衡量）
	 */
	public double[] dPQU(){
		double[] dPQU=new double[2*n-2];
		for(int i=0;i<s;i++){
			dPQU[2*i]=dP(i);      //求注入有功的不平衡量   
			dPQU[2*i+1]=dQU2(i);  //求注入无功或节点电压平方的不平衡量
		}
		for(int i=s+1;i<n;i++){
			dPQU[2*i-2]=dP(i);
			dPQU[2*i-1]=dQU2(i);
		}
		return dPQU;
	}
	/**
	 * 求各节点电压修正量
	 * @return 返回值是一个double型的一维数组（各节点电压修正量）
	 */
	public double[] dU(){
		double[] dU=Mat.multi(Mat.inv(J),dPQU);
		return dU;
	}
	/**
	 * 求各节点电压新值
	 * @return 返回值是一个Plural型的一维数组（各节点电压新值量）
	 */
	public Plural[] New_U(){
		//double[] dU=dU();
		for(int i=0;i<s;i++){
			U[i].a+=dU[2*i+1];
			U[i].b+=dU[2*i];                  
		}
		for(int i=s+1;i<n;i++){
			U[i].a+=dU[2*i-1];
			U[i].b+=dU[2*i-2]; 
		}
		return U;
	}
	/**
	 * 判断收敛
	 * @param 参数aa是一个double型的数（收敛精度）
	 * @return 返回值是一个int型的数（即电压是否收敛，1表示收敛，0表示不收敛）
	 */
    public int convergence(double aa){
    	int result=1;
    	//double[] dU=dU();
    	for(int i=0;i<2*n-2;i++){
    		if(Math.abs(dU[i])>aa){    //求各节点电压修正量的绝对值
    			result=0;              //若某个节点电压修正量的绝对值大于收敛精度aa
    			break;                 //则置result为0，既不收敛，强制退出循环
    		}	
    	}
    	return result;
    }
}