ieee14g_s.cpp

IEEE14节点的牛顿拉夫逊算法进行潮流计算

#include<iostream>
#include<complex>
#include<fstream>
using namespace std;
#include<cmath>
#include<iomanip>
#define pi 3.1415926

complex<float> conj(complex<float>);
float Q(complex<float> [15][1],complex<float> [15][15],int);
float AM(complex<float>);


void main()
{   
	float realY,imagY;
	float PL[15][1],PG[15][1],QL[15][1],QG[15][1],value;
	complex<float> Y[15][15];
	int i,j;
	ifstream YDATA1("YDATAC1_10.txt",ios::in);
	if(!YDATA1){
		cout<<"Adata cannot be opened!"<<endl;
		exit(1);
	}
	else{
		for(i=1;i<15;i++)
			for(j=1;j<=10;j++){
				YDATA1>>realY>>imagY;
				Y[i][j].real(realY);
				Y[i][j].imag(imagY);
				}
	}


		ifstream YDATA2("YDATAC11_14.txt",ios::in);
	    if(!YDATA2){

		    cout<<"Adata cannot be opened!"<<endl;
		    exit(1);
		}
	    else{
		
			for(i=1;i<=14;i++)
			for(j=11;j<=14;j++){
				YDATA2>>realY>>imagY;
				Y[i][j].real(realY);
				Y[i][j].imag(imagY);
				}
	}
       

		
		
			
			ifstream SDATA("SDATA.txt",ios::in);
	    if(!SDATA){
				cout<<"Sdata cannot be opened!"<<endl;
		        exit(1);
		}
	    else{
	
	
			for(i=1;i<=14;i++){
			SDATA>>value;
		    PG[i][0]=value;
			}
		
			for(i=1;i<=14;i++){
			SDATA>>value;
		    QG[i][0]=value;
			}
		
			for(i=1;i<=14;i++){
			SDATA>>value;
			PL[i][0]=value;
			}
		
			for(i=1;i<=14;i++){
			SDATA>>value;
			QL[i][0]=value;
			}
		}



	float Psp[15][1],Qsp[15][1];
	for(i=1;i<=14;i++){
		Psp[i][0]=PG[i][0]-PL[i][0];
		Qsp[i][0]=QG[i][0]-QL[i][0];
	}

	complex<float> U[15][1];
	U[1][0].real(1.06);
	U[1][0].imag(0);

	for(i=2;i<=14;i++){
		U[i][0].real(1);
		U[i][0].imag(0);
	}
	float dU=0;
	complex<float> temp;

	ofstream IEEE14("IEEE14.txt",ios::out);
    int k=0;
			
	do{

		IEEE14<<left<<setw(10)<<k;
		for(i=1;i<=14;i++){
			IEEE14<<left<<setw(10)<<AM(U[i][0])<<left<<setw(10)<<atan(imag(U[i][0])/real(U[i][0]))/(2*pi)*360;
		}
		IEEE14<<endl;
		
		

	
		
		dU=0;
		for(i=2;i<=14;i++){
			if(i==2){
				    temp=U[i][0];
					complex<float> Ssp(Psp[i][0],-Qsp[i][0]);
							U[i][0]=Ssp/conj(U[i][0])/Y[i][i];
							for(j=1;j<=14;j++){
								if(j==i)
										continue;
									else
										U[i][0]-=Y[i][j]/Y[i][i]*U[j][0];
								}
                U[i][0]*=1.045/AM(U[i][0]);
				temp=U[i][0]-temp;
				dU+=AM(temp);
				if(Q(U,Y,i)>0.5)
					Qsp[i][0]=0.5;
				else
					if (Q(U,Y,i)<-0.4)
						Qsp[i][0]=-0.4;
					else
						Qsp[i][0]=Q(U,Y,i);
			}
			else
				if(i==3){
					    temp=U[i][0];
						complex<float> Ssp(Psp[i][0],-Qsp[i][0]);
							U[i][0]=Ssp/conj(U[i][0])/Y[i][i];
							for(j=1;j<=14;j++){
								if(j==i)
										continue;
									else
										U[i][0]-=Y[i][j]/Y[i][i]*U[j][0];
								}
					U[i][0]*=1.01/AM(U[i][0]);
					temp=U[i][0]-temp;
				dU+=AM(temp);
				if(Q(U,Y,i)>0.4)
					Qsp[i][0]=0.4;
				else
					if (Q(U,Y,i)<0)
						Qsp[i][0]=0;
					else
						Qsp[i][0]=Q(U,Y,i);
					
				}
				else
					if(i==6){
						    temp=U[i][0];
							complex<float> Ssp(Psp[i][0],-Qsp[i][0]);
							U[i][0]=Ssp/conj(U[i][0])/Y[i][i];
							for(j=1;j<=14;j++){
								if(j==i)
										continue;
									else
										U[i][0]-=Y[i][j]/Y[i][i]*U[j][0];
								}
						U[i][0]*=1.07/AM(U[i][0]);
						temp=U[i][0]-temp;
				dU+=AM(temp);
				if(Q(U,Y,i)>0.24)
					Qsp[i][0]=0.24;
				else
					if (Q(U,Y,i)<-0.06)
						Qsp[i][0]=-0.06;
					else
                        Qsp[i][0]=Q(U,Y,i);
					}
					else
						if(i==8){
							temp=U[i][0];
								complex<float> Ssp(Psp[i][0],-Qsp[i][0]);
							U[i][0]=Ssp/conj(U[i][0])/Y[i][i];
							for(j=1;j<=14;j++){
								if(j==i)
										continue;
									else
										U[i][0]-=Y[i][j]/Y[i][i]*U[j][0];
								}
							U[i][0]*=1.09/AM(U[i][0]);
							temp=U[i][0]-temp;
				dU+=AM(temp);
				if(Q(U,Y,i)>0.24)
					Qsp[i][0]=0.24;
				else
					if (Q(U,Y,i)<-0.06)
						Qsp[i][0]=-0.06;
					else
                            Qsp[i][0]=Q(U,Y,i);
						}
						else{
								temp=U[i][0];

							complex<float> Ssp(Psp[i][0],-Qsp[i][0]);
							U[i][0]=Ssp/conj(U[i][0])/Y[i][i];
							for(j=1;j<=14;j++){
								if(j==i)
										continue;
									else
										U[i][0]-=Y[i][j]/Y[i][i]*U[j][0];
								}
							temp=U[i][0]-temp;
				            dU+=AM(temp);

							}
		}



		
	
k++;
	}while(dU>pow(10,-5));




	






	
	















}


complex<float> conj(complex<float> x)
{
	x.real(real(x));
	x.imag(-imag(x));
	return x;
}


float Q(complex<float> U[15][1],complex<float> Y[15][15],int p){
	complex<float> S(0,0);
	S=conj(U[p][0])*Y[p][p]*U[p][0];
	for(int J=1;J<=14;J++){
		if(J==p)
			continue;
		else
			S+=conj(U[p][0])*Y[p][J]*U[J][0];
	}
	return -imag(S);

    
}

float AM(complex<float> x){
	return sqrt(pow(real(x),2)+pow(imag(x),2));
}