nljstate.cpp

在系统辨识中,NLJ的优化算法用于系统辨识

// NLJSTATE.cpp : Defines the entry point for the console application.
//

#include "stdio.h"
#include "math.h"
#include "fstream.h"
int nNljParaNum; 
double dTargetVal =1e5;
double dStep =0.005;
double GetFunVal(double x[], double a[], int i)
{
	//////////////////////////////////////////////////////
	/// 根据给定的状态向量和系数求取状态方程i的导数值  ///
	/// x: 状态向量								       ///
	/// a: 系数数据								       ///
	/// i：第i个方程								   ///
	//////////////////////////////////////////////////////
	double dF = a[4]+8.75;
	double dI =0.9;
	double val = 0;
	switch(i)
	{
	case 0:
		val = a[4] - dF * x[0] -a[0] * x[0] * x[1] - a[3] * x[0] * x[5]* pow(dI, 0.5);
		break;
	case 1:
		val = 7.0 - dF * x[1] - a[0] * x[0] * x[1] - 2 * a[1] * x[1] * x[2];
		break;
	case 2:
		val = 1.75 -dF * x[2] - a[1] * x[1] *x[2];
		break;
	case 3:
		val = -dF * x[3] + 2 * a[0] * x[0] * x[1] - a[2] * x[3] * x[4];
		break;
	case 4:
		val =  -dF * x[4] + 3 * a[1] * x[1] * x[2] - a[2] * x[3] * x[4];
		break;
	case 5:
		val = -dF * x[5] + 2 * a[2] * x[3] * x[4] - a[3] * x[0] * x[5] * pow(dI, 0.5);
		break;
	case 6:
		val = -dF * x[6] + 2 * a[3] * x[0] * x[5] * pow(dI, 0.5);

	}
	return val;
}

void GetStateOutput(double da[],double dArrX[],  double step, double y[], int nStepNum)
{
	//////////////////////////////////////////////////////
	/// 根据给定的状态向量和系数求取状态方程i的导数值  ///
	/// dArrX: 初始状态向量						       ///
	/// a: 系数数据								       ///
	/// step：四阶龙格库塔法步长					   ///
	/// y:输出状态向量								   ///
	/// nStepNum：输出数据组数						   ///
	//////////////////////////////////////////////////////
	const int MAXNUM = 7;
	const int L = 40;
	double dArrK[MAXNUM][4];
	double x[MAXNUM];
	for (int is = 0; is <= nStepNum ; is++)
	{

		//------------------getkij
		for (int n = 0; n < 4; n++)
			for(int i = 0; i < MAXNUM; i++)
			{ 
				if(n == 0)
					for(int j = 0; j < MAXNUM; j++)
						x[j] = dArrX[j];
				else if(n < 3)
					for(int j = 0; j < MAXNUM; j++)
						x[j] = dArrX[j] + step / 2 * dArrK[j][n-1]; 
				else
					for(int j = 0; j < MAXNUM; j++)
						x[j] = dArrX[j] + step * dArrK[j][n-1]; ;
		dArrK[i][n] = GetFunVal(x, da, i);
			}
		//------------------endkij
		for(int j = 0; j < MAXNUM; j++)
			dArrX[j] += (dArrK[j][0] + 2 * dArrK[j][1] + 
								2 * dArrK[j][2] + dArrK[j][3]) *
							step /6;

		for(j = 0; j < MAXNUM; j++)
				y[is * MAXNUM + j ] = dArrX[j];	
	}
}

bool NljIdentMathord(double dNljX[], int nNljParaNum, double dNljRange[], int nSearchDeep)
{
	//////////////////////////////////////////////////////
	/// 使用NLJ算法求取最优解                          ///
	/// dNLJX: 辨识参数的初始值及结果			       ///
	/// nNljParaNum:辨识参数个数				       ///
	/// dNljRange：搜索范围内   					   ///
	/// nSearchDeep：搜索深度						   ///
	//////////////////////////////////////////////////////
	const int MAXNUM = 7;
	const int L = 40;
	double temp[MAXNUM];
	double opt[MAXNUM];
	double dNljRand[101][20];
	bool outfind = false, infind;
	double y[L][MAXNUM];
	double dX0[] = {0.1883, 0.2507, 0.0467, 0.0899, 0.1804, 0.1394, 0.1046};
	double dY0[L][MAXNUM];
	double x[MAXNUM];
	double ise = 0;
	ifstream InfRand, InfTest;
	InfRand.open("rand.txt");
	int nSearchWidth = 101;
	for(int i = 0; i < nSearchWidth; i++)
		for(int j = 0; j < 20; j++)
			InfRand >> dNljRand[i][j];
	InfRand.close();
	InfTest.open("simu.txt");
	for(i = 0; i < L; i++)
		for(int j = 0; j < MAXNUM; j++)
			InfTest >>dY0[i][j];
	InfTest.close();

	for(i = 0; i< nSearchDeep; i++)
	{
		//  获得前一步取得的最优值opt
		for(int k = 0; k < nNljParaNum; k++)
			opt[k] = dNljX[k];
		// 获取本次的最优值
		infind = false;
		bool bSatisfied = true;
		for(k = 0; k < nSearchWidth; k++)
		{
			for(int m = 0; m < nNljParaNum; m++)
				temp[m] = dNljX[m] + dNljRand[k][(i*nNljParaNum+m)%20] * dNljRange[m];
			for( m = 0; m < nNljParaNum; m++)
				if(temp[m] < 0)
					bSatisfied = false;
			if(!bSatisfied)
				continue;
			for(m = 0; m < MAXNUM; m++)
				x[m] = dX0[m];
			GetStateOutput(temp, x, dStep, &y[0][0], L);
			ise = 0;
			for(m = 0; m < L; m++)
				for(int n = 0; n < MAXNUM; n++)
					ise += (y[m][n] - dY0[m][n]) * (y[m][n] - dY0[m][n]);//ITAEswitch(nIdentType)

			if (ise > dTargetVal)
				continue;
			dTargetVal = ise;
			for( m = 0; m<nNljParaNum; m++)
					opt[m] = temp[m];
		}
		if(infind)
		{
			for(int k = 0; k < nNljParaNum; k++)
				dNljX[k] = opt[k];
			outfind = true;
		}
		double fai = 0.97;
		fai = pow(0.97, pow(i, 1.34));
		for(int j = 0; j < nNljParaNum; j++)
			dNljRange[j] = fai * dNljRange[j];	
	}
	return outfind;	
}

int main(int argc, char* argv[])
{
	double a0[] = {50, 100, 100, 50, 10};
	const int nNljNum = 5;
	int nSearchDeep = 50;
	double r[nNljNum];
	for(int i = 0; i < nNljNum; i++)
		r[i] = a0[i] * 0.5;
	if(NljIdentMathord(a0, nNljNum, r, nSearchDeep))
	{
		for(i = 0; i< nNljNum;i++)
			cout << a0[i] << endl;
	}
	cout << "the NLJ target val is " << dTargetVal << endl;
	return 0;
}