nr_powerflow.cpp

牛啦潮流算法源代码

// NR_PowerFlow.cpp : Defines the initialization routines for the DLL.
//

#include "stdafx.h"
#include "NR_PowerFlow.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

//
//	Note!
//
//		If this DLL is dynamically linked against the MFC
//		DLLs, any functions exported from this DLL which
//		call into MFC must have the AFX_MANAGE_STATE macro
//		added at the very beginning of the function.
//
//		For example:
//
//		extern "C" BOOL PASCAL EXPORT ExportedFunction()
//		{
//			AFX_MANAGE_STATE(AfxGetStaticModuleState());
//			// normal function body here
//		}
//
//		It is very important that this macro appear in each
//		function, prior to any calls into MFC.  This means that
//		it must appear as the first statement within the 
//		function, even before any object variable declarations
//		as their constructors may generate calls into the MFC
//		DLL.
//
//		Please see MFC Technical Notes 33 and 58 for additional
//		details.
//

/////////////////////////////////////////////////////////////////////////////
// CNR_PowerFlowApp

BEGIN_MESSAGE_MAP(CNR_PowerFlowApp, CWinApp)
	//{{AFX_MSG_MAP(CNR_PowerFlowApp)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//    DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CNR_PowerFlowApp construction

CNR_PowerFlowApp::CNR_PowerFlowApp()
{
	// TODO: add construction code here,
	// Place all significant initialization in InitInstance
}

/////////////////////////////////////////////////////////////////////////////
// The one and only CNR_PowerFlowApp object

CNR_PowerFlowApp theApp;

/*
编写者：侯学勇
编写日期:2003-10-22
功能：计算潮流，DLL导出函数
参数：调用自己的接口结构体参数

*/
extern "C" _declspec(dllexport) bool CalculatePFOtherParam(struct tagGeneticForPQ *pParam,double *pfVoltage,double *pfPhase, double *pInpourP, double *pInpourQ)
{
		AFX_MANAGE_STATE(AfxGetStaticModuleState());
	
		//统一接口
		NodeParam   bpNode;
		BranchParam  bpBranch ;
		NRParam   bpNR ;
		bpNode.nCount = pParam->G_NodeNumber;
		bpNode.pTpIndex = pParam->G_NodeNum;
		bpNode.PQCount = 0;
		bpNode.PVCount = 0;
		int *pType = new int[bpNode.nCount];
		double *pPN = new double[bpNode.nCount];
		double *pQN = new double[bpNode.nCount];
		double *pFn =new double[bpNode.nCount];
		ZeroMemory(pType,sizeof(int)*bpNode.nCount);
		ZeroMemory(pPN,sizeof(double)*bpNode.nCount);
		ZeroMemory(pQN,sizeof(double)*bpNode.nCount);
		ZeroMemory(pFn,sizeof(double)*bpNode.nCount);

		double *pPQ_Uef = new double[2*(bpNode.nCount-1)];
		ZeroMemory(pPQ_Uef,sizeof(double)*2*(bpNode.nCount-1));
                bpNR.pPQ_Uef = pPQ_Uef;

		CStdioFile fp;
		CString szTemp;
		//fp.Open("c:\\test.txt", CFile::modeCreate | CFile::modeWrite);

		for (int i = 0; i < bpNode.nCount; i++)
		{
			if (pParam->G_Nodetype[i] == 0)
			{
				pType[i] = 1;
				bpNode.PQCount++;
			}
			else if (pParam->G_Nodetype[i] == 1)
			{
				pType[i] = 2;
				bpNode.PVCount++;
			}
			else
			{
				pType[i] = 0;
			}
			pPN[i] = pParam->G_J_gp[i]-pParam->G_J_lp[i];
			pQN[i] = pParam->G_J_gq[i]-pParam->G_J_lq[i];			
                         
			//szTemp.Format("%4d\t%8.4lf\t%8.4lf\n", pParam->G_NodeNum[i],pParam->G_J_gp[i],pParam->G_J_gq[i]);
			//fp.WriteString(szTemp);
		}
		
		//fp.Close();
		bpNode.pType = pType;
        bpNode.pPN = pPN;
        bpNode.pQN = pQN;
        bpNode.pFn = pFn;

		bpNode.pEn = pParam->G_Voltage;
		bpNode.pUN = pParam->G_Voltage;

		bpNR.nGBCount = bpNode.nCount*bpNode.nCount;
		bpNR.nJacobiCount = 2*(bpNode.PQCount+bpNode.PVCount);
		bpNR.pG = pParam->G_Conductance;
		bpNR.pB = pParam->G_B;  

		double *pTempJacobi =  new double[bpNR.nJacobiCount*bpNR.nJacobiCount];
		bpNR.pJacobi = pTempJacobi;

		/*高斯赛德尔迭代三次，给牛拉计算的节点电压初值*/
	   //	theApp.GAUSS_SEIDEL( &bpNode, &bpNR);
 
		theApp.Didainum = 50;
		theApp.RemainTimes = 50;

		bpNR.duEps = 0.00001;//计算精度

	    /*牛拉法迭代计算过程,如果计算收敛的话,将结果输出，*/
 		if(theApp.didai( &bpNode, &bpBranch, &bpNR) == 0)
 		{	
			//"计算成功!请查看计算结果!"
		//	AfxMessageBox("计算成功!请查看计算结果!");
			for (i = 0; i < bpNode.nCount; i++)
			{
				pfVoltage[i] = sqrt(bpNode.pEn[i]*bpNode.pEn[i]+bpNode.pFn[i]*bpNode.pFn[i]);
				pfPhase[i] = atan2(bpNode.pFn[i], bpNode.pEn[i]);
				pInpourP[i] = bpNode.pPN[i];
				pInpourQ[i] = bpNode.pQN[i];
			}
			delete [] pType;
			pType = NULL;
			delete [] pPN;
			pPN = NULL;
			delete [] pQN;
			pQN = NULL;
			delete [] pFn;
			pFn = NULL;
			delete [] pPQ_Uef;
			pPQ_Uef = NULL;
			delete []pTempJacobi;
			pTempJacobi = NULL;
			return TRUE;
		}
		else
		{
			//"计算不成功,请检查原始数据正确性!"
		//	AfxMessageBox("计算不成功,请检查原始数据正确性!");
			delete [] pType;
			pType = NULL;
			delete [] pPN;
			pPN = NULL;
			delete [] pQN;
			pQN = NULL;
			delete [] pFn;
			pFn = NULL;
			delete [] pPQ_Uef;
			pPQ_Uef = NULL;
			delete []pTempJacobi;
			pTempJacobi = NULL;
			return FALSE;
		} 		
}

/*
编写者：侯学勇
编写日期:2003-07-22
功能：迭代计算后修正电压及形成功率误差量
输入参数：
	nCount：系统总节点数；
	pType：指向支路数据首节点号数组的首地址；
	pEn：指向支路数据尾节点号数组的首地址；
	pEn：节点电压实部；
	pFn：节点电压虚部；
	pUN：节点电压模值；
	pPN：节点注入有功功率；
	pQN：节点注入无功功率；	
	pG：指向导纳矩阵实部的首地址；
	pB：指向导纳矩阵虚部的首地址；
	pPQ_Uef：修正量的首地址；
返回参数：无
备注:
     迭代计算后电压修正及功率误差量形成
*/ 
void CNR_PowerFlowApp::xiuzhen(NodeParam *bpNode, NRParam *bpNR)
{
	int  nCount = bpNode->nCount;
	int *pType = bpNode->pType;
	int nPVCount = bpNode->PVCount;
	double *pEn = bpNode->pEn;
	double *pFn = bpNode->pFn;
	double *pUN = bpNode->pUN;
	double *pPN = bpNode->pPN;
	double *pQN = bpNode->pQN;
	double *pG = bpNR->pG;
	double *pB = bpNR->pB;
	double *pPQ_Uef = bpNR->pPQ_Uef;


	int i=0,k=0;
	//电压修正(第一次也修正，但是等于原值（－0）,相当于未修正)
	for(i=0,k=0;i<nCount;i++)
	{	
		if(pType[i]==0)
			continue;//平衡节点则进入下一次循环
		pEn[i]+=pPQ_Uef[2*k+1];
		pFn[i]+=pPQ_Uef[2*k+0];
		k++;
	}
	//求解功率不平衡量
	double sigmaP,sigmaQ;
	for(i=0,k=0;i<nCount;i++)
	{	
		if(pType[i]==0) 
		{
			continue;//平衡节点则进入下一次循环
		}
		if(pType[i]==1)//PQ节点
		{
			sigmaP=0;sigmaQ=0;
			for(int j=0;j<nCount;j++)
			{
				sigmaP+=pEn[i]*(pG[i*nCount+j]*pEn[j]-pB[i*nCount+j]*pFn[j])+pFn[i]*(pG[i*nCount+j]*pFn[j]+pB[i*nCount+j]*pEn[j]);
				sigmaQ+=pFn[i]*(pG[i*nCount+j]*pEn[j]-pB[i*nCount+j]*pFn[j])-pEn[i]*(pG[i*nCount+j]*pFn[j]+pB[i*nCount+j]*pEn[j]);
			}
			pPQ_Uef[2*k+0]=pPN[i]-sigmaP;
			pPQ_Uef[2*k+1]=pQN[i]-sigmaQ;
			k++;
		}
		else if(pType[i]==2)//PV节点
		{			
			sigmaP=0;
			for(int j=0;j<nCount;j++)
				sigmaP+=(pEn[i]*(pG[i*nCount+j]*pEn[j]-pB[i*nCount+j]*pFn[j])+pFn[i]*(pG[i*nCount+j]*pFn[j]+pB[i*nCount+j]*pEn[j]));
			
			pPQ_Uef[k*2+0]=pPN[i]-sigmaP;
			pPQ_Uef[k*2+1]=pUN[i]*pUN[i]-((pEn[i]*pEn[i])+(pFn[i]*pFn[i]));
			k++;
		}			
	}
}

/*
编写者：侯学勇
编写日期:2003-10-22
功能：形成雅克比矩阵
输入参数：
	pG：指向导纳矩阵实部的首地址；
	pB：指向导纳矩阵虚部的首地址；
	pType：指向节点类型数组的首地址；
	pEn：指向节点电压实部的首地址；
	pFn：指向节点电压虚部的首地址；
	nCount1：系统总节点数；
	pJacobi：指向雅克比矩阵（一维数组）的首地址；
	nCount2：系统非平衡节点总数；
返回参数：无
备注:
根据修正后的电压及导纳矩阵形成雅克比矩阵
*/ 
void CNR_PowerFlowApp::JACOBI(NodeParam *bpNode, NRParam *bpNR)
{
	int nCount1 = bpNode->nCount;
	int *pType = bpNode->pType;
	double *pEn = bpNode->pEn;
	double *pFn = bpNode->pFn;

	double *pG = bpNR->pG;
	double *pB = bpNR->pB;