ieee30pq分解法.cpp

电力系统潮流计算 采取PQ分解法 用VC++完成编程。

#include<iostream>
#include<fstream>
#include<cmath>
#include<iomanip>
#include"mymat.h"

using namespace std;
#define pi 3.1415926


void main()
{


	int i,j;
	double max;
    
	double G[31][31],B[31][31];
	
	ifstream GDATA("GDATA.txt",ios::in);
    if(!GDATA){
		cout<<"file can't be opened!"<<endl;
		exit(1);
	}
	else{
		for(i=1;i<=30;i++)
			for(j=1;j<=30;j++)
				GDATA>>G[i][j];
	}
	ifstream BDATA("BDATA.txt",ios::in);
    if(!BDATA){
		cout<<"file can't be opened!"<<endl;
		exit(1);
	}
	else{
		for(i=1;i<=30;i++)
			for(j=1;j<=30;j++)
				BDATA>>B[i][j];
	}

	ifstream PGDATA("PGDATA.txt",ios::in);
	if(!PGDATA){
		cout<<"error!"<<endl;
		exit(1);
	}
	double PG[31],QG[31],PL[31],QL[31];

	for(i=1;i<=30;i++)
		PGDATA>>PG[i];
	
	
	ifstream QGDATA("QGDATA.txt",ios::in);
	if(!QGDATA){
		cout<<"error!"<<endl;
		exit(1);
	}
	for(i=1;i<=30;i++)
		QGDATA>>QG[i];
		
	
	ifstream PLDATA("PLDATA.txt",ios::in);
	if(!PLDATA){
		cout<<"error!"<<endl;
		exit(1);
	}
	for(i=1;i<=30;i++)
		PLDATA>>PL[i];
		
	
	ifstream QLDATA("QLDATA.txt",ios::in);
	if(!QLDATA){
		cout<<"error!"<<endl;
		exit(1);
	}
	for(i=1;i<=30;i++)
		QLDATA>>QL[i];
		
	
	double Psp[31],Qsp[31];
	for(i=1;i<=30;i++){
		Psp[i]=PG[i]-PL[i];
        Qsp[i]=QG[i]-QL[i];
	}
	double U[31],theta[31];
	U[1]=1.0500;U[2]=1.0338;U[3]=1.0230;U[4]=1.0913;U[5]=1.0058,U[6]=1.0883;
	for(i=7;i<=30;i++)
		U[i]=1.000;
	for(i=1;i<=30;i++)
		theta[i]=0.000;



    ofstream IEEE14output("IEEE14output.txt",ios::out);

	for(i=1;i<=30;i++)
	IEEE14output<<'U'<<left<<setw(9)<<i;



	
	for(i=1;i<=30;i++)
		IEEE14output<<"theta"<<left<<setw(5)<<i;
	IEEE14output<<endl;

	double P[31],Q[31],dP[31],dQ[31],dU[31],dtheta[31];
    MyMat<double> D_P(29,1),D_Q(24,1),D_THETA(29,1),D_U(24,1),JP(29,29),JQ(24,24);
	for(i=0;i<29;i++)
		for(j=0;j<29;j++)
			JP[i][j]=B[i+2][j+2];

	for(i=0;i<24;i++)
		for(j=0;j<24;j++)
			JQ[i][j]=B[i+7][j+7];

	

	

	do{
		for(i=1;i<=30;i++)
		IEEE14output<<left<<setw(10)<<U[i];
		
		for(i=1;i<=30;i++)
			IEEE14output<<left<<setw(10)<<theta[i]/(2*pi)*360;
		IEEE14output<<endl;
		for(i=2;i<=30;i++)
		{ 
			P[i]=0;
			for(j=1;j<=30;j++)
				P[i]+=U[i]*U[j]*(G[i][j]*cos(theta[i]-theta[j])+B[i][j]*sin(theta[i]-theta[j]));
        }
		for(i=7;i<=30;i++){
			Q[i]=0;
			for(j=1;j<=30;j++)
				Q[i]+=U[i]*U[j]*(G[i][j]*sin(theta[i]-theta[j])-B[i][j]*cos(theta[i]-theta[j]));
		}

		for(i=2;i<=30;i++)
			dP[i]=Psp[i]-P[i];
		for(i=7;i<=30;i++)
			dQ[i]=Qsp[i]-Q[i];

		
		max=fabs(dP[2]);
		for(i=2;i<=30;i++){
			if(max<fabs(dP[i]))
				max=fabs(dP[i]);
		}
		for(i=7;i<=30;i++){
			if(max<fabs(dQ[i]))
				max=fabs(dQ[i]);
		}


		for(i=0;i<29;i++)
			D_P[i][0]=dP[i+2]/U[i+2];
        D_THETA=JP.inv()*D_P;
		for(i=0;i<29;i++)
			dtheta[i+2]=D_THETA[i][0]/U[i+2];

		for(i=0;i<24;i++)
			D_Q[i][0]=dQ[i+7]/U[i+7];
		D_U=JQ.inv()*D_Q;
		for(i=0;i<24;i++)
			dU[i+7]=D_U[i][0];
		for(i=2;i<=30;i++)
			theta[i]-=dtheta[i];
		for(i=7;i<=30;i++)
			U[i]-=dU[i];


		


	}while(max>pow(10,-8));



}