mathyw.cpp

C++最优化方法源程序

			direction[i] = -grad[i];
		}
		do
		{
			length = line_search2( pf, n, point, direction );	
			for( i=0; i<n; i++ )
			{
				point[i] += length* direction[i];
				dg[i] = grad[i];
			}
			value = 0;
			compute_grad( pf, n, point, grad );
			for( i=0; i<n; i++ )
			{
				value += grad[i];
			}
			value = sqrt( value );

			if( value >= precision )
			{
				for( i=0; i<n; i++ )
				{
					dp[i] = length* direction[i];   //计算point(k+1)- point(k)
					dg[i] = grad[i]- dg[i];          //计算grad(k+1)- grad(k)  
				}
				xandg = g_h_g = 0;
				for( int j=0; j<n; j++ )
				{
					xandg += dp[j]* dg[j];             //计算第一个分子
					g_h_g += dg[j]* dg[j]* f_hesan[j* n+ j];  //计算第二个分子
					for( i=j+1; i<n; i++ )
					{
						g_h_g += 2* dg[i]* dg[j]* f_hesan[j* n+ i];
					}						
					direction[j] = 0;
					for( i=0; i<n; i++ )
					{
						direction[j] += f_hesan[j* n+ i]* dg[i];  //计算 H* dg
					}
				}
				memset( matrix, 0, n* n* sizeof( double ) );
				for( i=0; i<n; i++ )
				{
					for( j=i; j<n; j++ )
					{
						for( int k=0; k<n; k++ )
						{
							matrix[i* n+ j] += direction[i]* dg[k]* f_hesan[k* n+ j]; //计算第二个分母
						}
					}
				}
				
				for( i=0; i<n; i++ )
				{
					int t = i* n+ i;
					f_hesan[t] += dp[i]* dp[i]/ xandg - matrix[t]/ g_h_g;
					for( j=i+1; j<n; j++ )
					{
						t = i* n+ j;
						f_hesan[t] += dp[i]* dp[j]/ xandg - matrix[t]/g_h_g; //计算新的拟hesan矩阵
						f_hesan[j* n+ i] = f_hesan[t];            //利用矩阵的对称性进行计算
					}

					direction[i] = 0;
					for( j=0; j<n; j++ )
					{
						direction[i] += f_hesan[i* n+ j]* grad[j];  //计算新的收索方向
					}
				}
			}
		}while( value > precision );
	}

	delete[] f_hesan;
	delete[] grad;
	delete[] matrix;
	delete[] direction;
	delete[] dp;
	delete[] dg;
}
/////////////////////////////////////////////////////////////////////////////////////
//void length_accelerate( double (* pf)( double * x ), int n, double * start, 
//                                                                     double * point )
//         步长加速的直接法求极小点。
//参数值：
//        pf: 被求的函数，
//        x: 所求函数的自变量，
//        n: 自变量的个数，
//        start: 提供的初始点， 
//        point: 记录收索的结果。
/////////////////////////////////////////////////////////////////////////////////////
void length_accelerate( double (* pf)( double * x ), int n, double * start, double * point )
{
	double precision = 1E-12;
	double length = 0.4;
	double proportion = 0.5;
	double * point_m = new double[n];
	double * point_zero = new double[n];
	double value_start, value_end;

	memcpy( point_zero, start, n* sizeof( double ) );
	memcpy( point, start, n* sizeof( double ) );
	do
	{
		for( int j=0; j<n; j++ )      //探测移动
		{
			value_start = pf( point );
			point[j] += length;
			if( pf( point ) >= value_start )
			{
				point[j] -= 2* length;
				if( pf( point ) >= value_start )
				{
					point[j] += length;
				}
			}
		}
		value_end = pf( point);
		if( value_end- value_start < 1E-15 )  //防止两值相等
		{
			do
			{
				for( int j=0; j<n; j++ )
				{
					point_m[j] = 2* point[j]- point_zero[j];
				}
				for( j=0; j<n; j++ )      //探测移动
				{
					value_start = pf( point_m );
					point_m[j] += length;
					if( pf( point_m ) >= value_start )
					{
						point_m[j] -= 2* length;
						if( pf( point_m ) >= value_start )
						{
							point_m[j] += length;
						}
					}
				}
				value_start = value_end;
				value_end = pf( point_m );
				memcpy( point_zero, point, n* sizeof( double ) );
				if( value_start > value_end )
				{
					memcpy( point, point_m, n* sizeof( double ) );
				}
			}while( (value_start- value_end) > 1E-15 );  //防止两值相等
		}
		length *= proportion;
	}while( length >= precision );

	delete[] point_m;
	delete[] point_zero;
}
/////////////////////////////////////////////////////////////////////////////////////
//void Powell_search( double (* pf)( double * x ), int n, double * start, double * point )
//         方向加速的Powell法求极小点。
//参数值：
//        pf: 被求的函数，
//        x: 所求函数的自变量，
//        n: 自变量的个数，
//        start: 提供的初始点， 
//        point: 记录收索的结果。
//函数调用：
//  一维收索函数double line_search1( double (*pf)( double * x ), int n, double * start,
//                                                                   double * direction )
/////////////////////////////////////////////////////////////////////////////////////
void Powell_search( double (* pf)( double * x ), int n, double * start, double * point )
{
	double * * direction = new double*[n];
	double * orientation = new double[n];
	double * dot         = new double[n];
	double * tempdot     = new double[n];
	double precision = 1E-5;
	double length, value, result, maxval;
	int numble;

	for( int i=0; i<n; i++ )
	{
		direction[i] = new double[n];
		memset( direction[i], 0, n* sizeof( double ) );
		direction[i][i] = 1;
	}
	memcpy( dot, start, n* sizeof( double ) );

	do
	{
		length = line_search1( pf, n, dot, direction[n-1] );
		for( i=0; i<n; i++ )
		{
			tempdot[i] = dot[i]+ length* direction[n-1][i];
		}
		memcpy( point, tempdot, n* sizeof( double ) );
		numble = 0;
		maxval = 0;
		result = pf( point );
		for( int j=0; j<n; j++ )
		{	
			length = line_search1( pf, n, point, direction[j] );
			for( i=0; i<n; i++ )
			{
				point[i] += length* direction[j][i];
			}
			value = pf( point );
			result -= value;
			if( result- maxval > precision )
			{
				maxval = result;
				numble = j;
			}
			result = value;
		}
		result = 0;
		for( i=0; i<n; i++ )
		{
			value = point[i]- tempdot[i];
			result += value* value;
		}
		result = sqrt( result );
		if( result- precision > precision )
		{
			for( i=0; i<n; i++ )
			{
				orientation[i] = point[i]- tempdot[i];   //???
			}
			value = pf( point );
			length = line_search1( pf, n, point, orientation );
			for( i=0; i<n; i++ )
			{
				tempdot[i] = point[i]+ length* orientation[i];
			}
			value = sqrt( ( value- pf( tempdot ) )/ maxval );
			if( fabs( length )- value > precision )
			{ 
				for( j=numble; j<n-1; j++ )
				{
					for( i=0; i<n; i++ )
					{
						direction[j][i] = direction[j+1][i];
					}
				}
				for( i=0; i<n; i++ )
				{
					direction[n-1][i] = orientation[i];
				}
				memcpy( dot, point, n* sizeof( double ) );
			}
			else
			{
				memcpy( dot, tempdot, n* sizeof( double ) );
			}
		}
	}while( result- precision > precision );
		
	for( i=0; i<n; i++ )
	{
		delete[] direction[i];
	}
	delete[] direction;
	delete[] dot;
	delete[] tempdot;
	delete[] orientation;
}
/////////////////////////////////////////////////////////////////////////////////////
//void Commen_condition( FUN Obj, FUN * pG_group, FUN * pH_group, int var_num, int g_num, 
//                                            int h_num, double * start, double * point );
//         一般约束的乘子法求极小点。
//参数值：
//        Obj:           被求的目标函数，
//        pG_group:      不等式约束第一个约束函数的指针，
//        pH_group:      等式约束第一个约束函数的指针，
//        var_num:       自变量的个数，
//        g_num:         不等式约束的个数，
//        h_num:		 等式约束的个数，
//        start:         提供的初始点， 
//        point:         记录收索的结果。
//函数调用：
//    方向加速法Powell_search( double (* pf)( double * x ), int n, double * start, double * point )
/////////////////////////////////////////////////////////////////////////////////////
typedef double ( * FUN )( double * );  //自定义函数指针类型为 FUN 

static FUN    Multiplication_Obj;
static FUN    * Multiplication_G_group;
static FUN    * Multiplication_H_group;
static int    Multiplication_var_num;
static int    Multiplication_g_num;
static int    Multiplication_h_num;
static double Multiplication_M ;
static double * Multiplication_g_multi;
static double * Multiplication_h_multi;

void Commen_condition( FUN Obj, FUN * pG_group, FUN * pH_group, int var_num, int g_num, 
                                              int h_num, double * start, double * point )
{
	Multiplication_Obj      = Obj;
	Multiplication_G_group  = pG_group;
	Multiplication_H_group  = pH_group;
	Multiplication_var_num  = var_num;
	Multiplication_g_num    = g_num;
	Multiplication_h_num    = h_num;
	Multiplication_M        = 1;

	double precision        = 1E-12;
	Multiplication_g_multi  = new double[g_num];
	Multiplication_h_multi  = new double[h_num];
	double * temp           = new double[var_num];
	double factor           = 2;
	double value, answer, sum1, sum2 = 1;
	double function( double * x );

	memcpy( point, start, sizeof(double)* var_num );

	int i;
	do
	{
		memset( Multiplication_g_multi, 0, sizeof(double)* g_num );
		memset( Multiplication_h_multi, 0, sizeof(double)* h_num );
		memcpy(  temp, point, sizeof(double)* var_num );
	    Powell_search( function, var_num, temp, point );

		sum2 = sum1;
		sum1 = 0;
		for( i=0; i<h_num; i++ )
		{
			value = (* (pH_group+ i))( point );
			sum1  += value* value;
		}
		for( i=0; i<g_num; i++ )
		{
			value  = (* (pG_group+ i))( point );
			answer = -Multiplication_g_multi[i]/ Multiplication_M;
			if( value < answer )
				value = answer;
			sum1 += value* value;
		}

//		sum2 = 0;
//		for( i=0; i<var_num; i++ )
//		{
//			sum2 += point[i]* point[i];
//		}
		if( sum1- precision > precision  )
		{
			for( i=0; i<g_num; i++ )
			{
				value = Multiplication_g_multi[i]+ Multiplication_M*( (*(pG_group+ i))( point ) );
				if( value < 0 )
					value = 0;
				Multiplication_g_multi[i] = value;
			}
			for( i=0; i<h_num; i++ )
			{
				Multiplication_h_multi[i] += Multiplication_M*( (*(pH_group+ i))( point ) );
			}
			Multiplication_M *= factor;
		}
	}while(  sum1- precision > precision  );

	delete[] Multiplication_g_multi;
	delete[] Multiplication_h_multi;
	delete[] temp;
}
double function( double * x )
{
	int i;
	double result;
	double value = Multiplication_Obj(x);
	for( i=0; i<Multiplication_h_num; i++ )
	{
		result = (* (Multiplication_H_group+ i))( x );
		value += Multiplication_h_multi[i]* result+ result* result *Multiplication_M/ 2;
	}
	for( i=0; i<Multiplication_g_num; i++ )
	{
		result = (* (Multiplication_G_group+ i))( x ) + Multiplication_g_multi[i]/ Multiplication_M;
		if( result < 0 )
			result = 0;
		result = result* result- Multiplication_g_multi[i]* Multiplication_g_multi[i]/ Multiplication_M/ Multiplication_M;
		result *= Multiplication_M/ 2;
		value += result;
	}
	return value;
}

//*/