nolinear1.cpp

运筹学中的解非线性规划问题的一种方法,zoutendijk可行方向法

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

//#include "stdafx.h"
#include <math.h>
#include "stdio.h"
#include "lineSearch_sm.cpp"

//#include "define.h"

#define LIMIT_GENERATION 100
#define EPSLON 0.00000001



extern void sm();
extern double sm_min;
extern double sm_x[NUM_variable];


void zou_init();
void zou_getRestriction();
//double vector_Multiply(double *a,double *x,int dim);

//static double A[NUM_restriction][NUM_variable]={1.0,1.0,1.0,5.0,-1.0,0.0,0.0,-1.0};
//static double B[NUM_restriction] = {2.0,5.0,0.0,0.0};
//static double current_result[LIMIT_GENERATION] ;//当前解
static bool isStrict[NUM_restriction];//if the restriction is strict
static double A1[NUM_restriction][NUM_variable];//A1: matrix of sm
static double A2[NUM_restriction][NUM_variable];//A2: to decide lamata_max
static double B1[NUM_restriction];//B1: b of sm
static double B2[NUM_restriction];//B2: to decide lamata_max
static int dim_B1;//dimension of B1 :紧约束的个数
static int dim_B2;//dimension of B2 :非紧约束的个数
static double grads[NUM_variable];//梯度
/*=================================================================*/

double f(double g, double h)
{
return (2.0*g*g +2.0*h*h -2.0*g*h -4.0*g -6.0*h);
}


double fx(double g, double h)
{
return (2.0*g - 2.0*h - 4);
}

double fy(double g, double h)
{
	return (-2.0*g + 2.0*h - 6);
}

//double fz(double g, double h,double i)


/*

*/ 
double lineSearch(double XK[2],double dk[2],double max)
{
double x[100];
//double y[100];
//double z[100];

int i=0;

double p;//=0.01;//步长
if(max<0.0)
	p = 0.01;
else if (max >1.0)
    p = 0.01;
else
    p = max * 0.001;


//初始值
x[0] = 0;
//y[0] = 0;
//z[0] = z0;

for(i=0;i<99;i++)//迭代100次
{
	//x[i + 1] = x[i] - p * fx(x[i], y[i]);
	//y[i + 1] = y[i] - p * fy(x[i], y[i]);
	//z[i + 1] = z[i] - p * fz(x[i], y[i],z[i]);
	x[i+1] = x[i] - p * ( fx(XK[0]+x[i]*dk[0],XK[1]+x[i]*dk[1])*dk[0]+ fy(XK[0]+x[i]*dk[0],XK[1]+x[i]*dk[1])*dk[1]);
	
	
	if(max >0.0)
	{
		if(x[i+1] > max)
		{
			if(f(XK[0]+x[i]*dk[0],XK[1]+x[i]*dk[1]) < f(XK[0]+max*dk[0],XK[1]+max*dk[1]))
				x[99] = x[i];
			else
				x[99] = max;
			break;
		}
	}
	
}


//printf("%f,%f,%f",x[99],y[99],f(x[99],y[99]));
printf("lamata:%f,min(f(XK+lamata*dk)):%f",x[99],f(XK[0]+x[99]*dk[0],XK[1]+x[99]*dk[1]));
return f(XK[0]+x[99]*dk[0],XK[1]+x[99]*dk[1]);
}

/*=================================================================*/
int main()
{
	int i,step;
	bool flag = true;
	double b_[NUM_variable] ,d_[NUM_variable],lamata_max;//用于求线搜索的步长限制
	double lamata_k;//步长
	zou_init();
	for(step=0;step<LIMIT_GENERATION;step++)
	{
		/*已知xk A B 求grads A1 A2 B1 B2 isStrict dim_B1 dim_B2*/
		zou_getRestriction();
		/*调用sm求线性规划问题，得到可行改进方向dk*/
		sm(NUM_variable,dim_B1+1,B1,A1,xk);//得到结果sm_min ,sm_x[](相当于dk)
		/*判断是否结束*/
		if(vector_Multiply(grads,sm_x,NUM_variable) < EPSLON)
		{
			//打印结果xk fun(xk)，结束
			printf("\n最优解：(");
			for(i = 0 ; i<NUM_variable ; i++)
			{
				printf(" %f ",sm_x[i]);
			}
			printf(")");
			printf("\n最优值：%f \n",sm_min);
			return 1;
		}
		/*求线搜索问题的步长限制lamata_max*/
		for(i=0;i<dim_B2;i++)
		{
			b_[i] = B2[i] - vector_Multiply(A2[i], xk,NUM_variable);
			d_[i] = vector_Multiply(A2[i], sm_x,NUM_variable);
			if (d_[i] >0.0)
				flag = false;
		}
		if(flag)
			lamata_max = -1.0;//此时lamata_max无界，用一个负数表示
		else
		{
			lamata_max = b_[0]/d_[0];
			for(i=1;i<dim_B2;i++)
			{
				if(lamata_max > b_[i]/d_[i])
					lamata_max = b_[i]/d_[i];
			}
		}

		/*线搜索求步长lamata*/
		lamata_k = lineSearch(xk, sm_x, lamata_max);		
		/*根据方向dk和步长lamata_k得到下一个可行解xk*/
		for(i = 0 ; i<NUM_variable ; i++)
		{
			xk[i] = xk[i]+lamata_k * sm_x[i];
		}		
	}

	return 0;
}

void zou_init()//赋值：x0  A ,B, f
{
	int i;

	for( i=0 ; i<NUM_variable ; i++)
	{
      xk[i] = x0[i];	
	}
}

void zou_getRestriction()//已知xk A B 求grads A1 A2 B1 B2 isStrict dim_B1 dim_B2
{
	int i,j;
	
	fun_grads(xk,grads);//求出grads,即 需要输入到sm中的目标函数

	dim_B1 =0;
	dim_B2 =0;
	for(i=0 ; i<NUM_restriction ; i++)
	{	
		if(vector_Multiply(A[i],xk,NUM_variable)==B[i])
		{
			isStrict[i] = true;
			for(j=0;j<NUM_variable;j++)
			{
				A1[dim_B1][j] = A[i][j];
			}
			B1[dim_B1] = B[i];
			dim_B1 ++;
		}
		else
		{
			isStrict[i] = false;
			for(j=0;j<NUM_variable;j++)
			{
				A2[dim_B2][j] = A[i][j];
			}			
			B2[dim_B2] = B[i];
			dim_B2 ++;
		}
	}
}