optimize_without_constraint.c

求解无约束最优化问题的各种算法集合!

	}
	return FAIL;
}

//功能说明:dfp法求最优解
//入参说明:e--控制误差范围
//		   get_ak--取迭代步长的函数
//		   get_gx--取梯度函数
//		   get_G--取G矩阵函数
//		   get_funcval--计算函数值的函数
//		   para--参数列表
//		   paranum--参数个数
//返 回 值:成功:SUCCESS 并打印最优值,失败FAIL
int dfp(double e,//待完
	   int (* get_ak)(double *ak,
		 			  double u,
					  double sigma,
					  double *pk,
					  double *para,
					  int (* get_gk)(double *para,int paranum,double *gk),
					  double (* func)(double *para,int paranum),
					  int paranum),
	   int (* get_gk)(double *para,int paranum,double *gk),
	   double (* get_funcval)(double *para,int paranum),
	   double *para,int paranum)
{
	int steepcount=0,paracount,stepcount=1,paracount2;
	double ak=0,*gk1,*gk2,*pk,*tempgk,b=0,*h1,*h2,*sk,*yk,tempval;
	double *tempvector,*para1,*para2,*tempmatrix;
	char *position="fr";

	if (!(gk1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(gk2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(pk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(tempgk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(sk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(yk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(tempvector=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(h1=(double *)_malloc(position,sizeof(double)*paranum*paranum)))//h矩阵
		return FAIL;
	if (!(h2=(double *)_malloc(position,sizeof(double)*paranum*paranum)))//h矩阵
		return FAIL;
	if (!(tempmatrix=(double *)_malloc(position,sizeof(double)*paranum*paranum)))//h矩阵
		return FAIL;
	
	memcpy(para1,para,sizeof(double)*paranum);//复制参数
	//初始化h为单位阵
	for (paracount=0;paracount<paranum;paracount++)
	{
		for (paracount2=0;paracount2<paranum;paracount2++)
		{
			if (paracount==paracount2)
				h1[paracount*paranum+paracount2]=1.0;
			else
				h1[paracount*paranum+paracount2]=0.0;
		}
	}
	get_gk(para1,paranum,gk1);//取梯度值gk
	//循环迭代,直到取到符合条件最优解
	while(1)
	{
		stepcount++;
		//计算pk=-Hk*gk
		matrix_multiply(h1,paranum,paranum,gk1,paranum,1,pk);
		for (paracount=0;paracount<paranum;paracount++)
			pk[paracount]=-pk[paracount];
		get_ak(&ak,U,SIGMA,pk,para1,get_gk,get_funcval,paranum);//计算步长,step3
		for (paracount=0;paracount<paranum;paracount++)//计算xk+1,step4
			para2[paracount]=para1[paracount]+pk[paracount]*ak;
		get_gk(para2,paranum,gk2);//取梯度值gk+1
		if (vector_val(gk2,paranum)<=e)//当前已到精度要求
		{
			printf("最优解:");
			output(para2,paranum,get_funcval);//输出最优结果
			printf("其中乘除法 %d 次,加减法 %d 次,迭代%d次!\n",muticount,addcount,stepcount);			
 			free(gk1);
			free(gk2);
			free(pk);
			free(h1);
			free(h2);
			free(sk);
			free(yk);
			free(tempmatrix);
			free(para1);
			free(para2);
			return SUCCESS;
		}
		printf("       ");
		output(para2,paranum,get_funcval);//输出中间结果
		//计算sk,yk,step5
		for (paracount=0;paracount<paranum;paracount++)
		{
			sk[paracount]=para2[paracount]-para1[paracount];
			yk[paracount]=gk2[paracount]-gk1[paracount];
		}
		matrix_multiply(h1,paranum,paranum,yk,paranum,1,tempvector);//Hkyk
		matrix_multiply(tempvector,paranum,1,yk,1,paranum,tempmatrix);//Hkykyk'
		matrix_multiply(tempmatrix,paranum,paranum,h1,paranum,paranum,h2);//Hkykyk'Hk
		matrix_multiply(yk,1,paranum,h1,paranum,paranum,tempvector);//yk'Hk
		tempval=vector_inter_multiply(tempvector,yk,paranum);//yk'Hkyk
		matrix_mul_div(h2,paranum,paranum,h2,tempval,DIV);//Hkykyk'Hk/yk'Hkyk
		matrix_add_sub(h1,paranum,paranum,h2,paranum,paranum,h2,SUB);//Hk-Hkykyk'Hk/yk'Hkyk
		tempval=vector_inter_multiply(yk,sk,paranum);//yk'sk
		matrix_multiply(sk,paranum,1,sk,1,paranum,tempmatrix);//sksk'
		matrix_mul_div(tempmatrix,paranum,paranum,tempmatrix,tempval,DIV);//sksk'/yk'sk
		matrix_add_sub(h2,paranum,paranum,tempmatrix,paranum,paranum,h2,ADD);//Hk+1
		memcpy(h1,h2,sizeof(double)*paranum*paranum);
		memset(h2,0,sizeof(double)*paranum*paranum);
		memcpy(gk1,gk2,sizeof(double)*paranum);
		memset(gk2,0,sizeof(double)*paranum);
		memcpy(para1,para2,sizeof(double)*paranum);
	}
	return FAIL;
}

//功能说明:bfgs法求最优解
//入参说明:e--控制误差范围
//		   get_ak--取迭代步长的函数
//		   get_gx--取梯度函数
//		   get_G--取G矩阵函数
//		   get_funcval--计算函数值的函数
//		   para--参数列表
//		   paranum--参数个数
//返 回 值:成功:SUCCESS 并打印最优值,失败FAIL
int bfgs(double e,//待完
	   int (* get_ak)(double *ak,
		 			  double u,
					  double sigma,
					  double *pk,
					  double *para,
					  int (* get_gk)(double *para,int paranum,double *gk),
					  double (* func)(double *para,int paranum),
					  int paranum),
	   int (* get_gk)(double *para,int paranum,double *gk),
	   double (* get_funcval)(double *para,int paranum),
	   double *para,int paranum)
{
	int steepcount=0,paracount,stepcount=1,paracount2;
	double ak=0,*gk1,*gk2,*pk,*tempgk,b=0,*h1,*h2,*sk,*yk,tempval,*wk,wtemp1,wtemp2,*wtempvector;
	double *tempvector,*para1,*para2,*tempmatrix;
	char *position="bfgs";

	if (!(gk1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(gk2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(pk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(tempgk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(sk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(yk=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para1=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(para2=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(tempvector=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	if (!(wtempvector=(double *)_malloc(position,sizeof(double)*paranum)))
		return FAIL;
	
	if (!(h1=(double *)_malloc(position,sizeof(double)*paranum*paranum)))//h矩阵
		return FAIL;
	if (!(h2=(double *)_malloc(position,sizeof(double)*paranum*paranum)))//h矩阵
		return FAIL;
	if (!(tempmatrix=(double *)_malloc(position,sizeof(double)*paranum*paranum)))//h矩阵
		return FAIL;
	if (!(wk=(double *)_malloc(position,sizeof(double)*paranum)))//wk
		return FAIL;
	
	memcpy(para1,para,sizeof(double)*paranum);//复制参数
	//初始化h为单位阵
	for (paracount=0;paracount<paranum;paracount++)
	{
		for (paracount2=0;paracount2<paranum;paracount2++)
		{
			if (paracount==paracount2)
				h1[paracount*paranum+paracount2]=1.0;
			else
				h1[paracount*paranum+paracount2]=0.0;
		}
	}
	get_gk(para1,paranum,gk1);//取梯度值gk
	//循环迭代,直到取到符合条件最优解
	while(1)
	{
		stepcount++;
		//计算pk=-Hk*gk
		matrix_multiply(h1,paranum,paranum,gk1,paranum,1,pk);
		for (paracount=0;paracount<paranum;paracount++)
			pk[paracount]=-pk[paracount];
		get_ak(&ak,U,SIGMA,pk,para1,get_gk,get_funcval,paranum);//计算步长,step3
		for (paracount=0;paracount<paranum;paracount++)//计算xk+1,step4
			para2[paracount]=para1[paracount]+pk[paracount]*ak;
		get_gk(para2,paranum,gk2);//取梯度值gk+1
		if (vector_val(gk2,paranum)<=e)//当前已到精度要求
		{
			printf("最优解:");
			output(para2,paranum,get_funcval);//输出最优结果
			printf("其中乘除法 %d 次,加减法 %d 次,迭代%d次!\n",muticount,addcount,stepcount);			
			free(gk1);
			free(gk2);
			free(pk);
			free(h1);
			free(h2);
			free(sk);
			free(yk);
			free(tempmatrix);
			free(para1);
			free(para2);
			free(wk);
			free(wtempvector);
			return SUCCESS;
		}
		printf("       ");
		output(para2,paranum,get_funcval);//输出中间结果
		//计算sk,yk,step5
		for (paracount=0;paracount<paranum;paracount++)
		{
			sk[paracount]=para2[paracount]-para1[paracount];
			yk[paracount]=gk2[paracount]-gk1[paracount];
		}
		matrix_multiply(h1,paranum,paranum,yk,paranum,1,tempvector);//Hkyk
		memcpy(wtempvector,tempvector,sizeof(double)*paranum);
		matrix_multiply(tempvector,paranum,1,yk,1,paranum,tempmatrix);//Hkykyk'
		matrix_multiply(tempmatrix,paranum,paranum,h1,paranum,paranum,h2);//Hkykyk'Hk
		matrix_multiply(yk,1,paranum,h1,paranum,paranum,tempvector);//yk'Hk
		wtemp1=tempval=vector_inter_multiply(tempvector,yk,paranum);//yk'Hkyk
		matrix_mul_div(h2,paranum,paranum,h2,tempval,DIV);//Hkykyk'Hk/yk'Hkyk
		matrix_add_sub(h1,paranum,paranum,h2,paranum,paranum,h2,SUB);//Hk-Hkykyk'Hk/yk'Hkyk
		wtemp2=tempval=vector_inter_multiply(yk,sk,paranum);//yk'sk
		matrix_multiply(sk,paranum,1,sk,1,paranum,tempmatrix);//sksk'
		matrix_mul_div(tempmatrix,paranum,paranum,tempmatrix,tempval,DIV);//sksk'/yk'sk
		matrix_mul_div(sk,paranum,1,wk,wtemp2,DIV);
		matrix_mul_div(wtempvector,paranum,1,wtempvector,wtemp1,DIV);
		matrix_add_sub(wk,paranum,1,wtempvector,paranum,1,wtempvector,SUB);
		matrix_mul_div(wk,paranum,1,wk,pow(wtemp1,0.5),MUL);
		matrix_add_sub(h2,paranum,paranum,tempmatrix,paranum,paranum,h2,ADD);//Hk+1
		matrix_multiply(wk,paranum,1,wk,1,paranum,tempmatrix);
		matrix_add_sub(h2,paranum,paranum,tempmatrix,paranum,paranum,h2,ADD);
		memcpy(h1,h2,sizeof(double)*paranum*paranum);
		memset(h2,0,sizeof(double)*paranum*paranum);
		memcpy(gk1,gk2,sizeof(double)*paranum);
		memset(gk2,0,sizeof(double)*paranum);
		memcpy(para1,para2,sizeof(double)*paranum);
	}
	return FAIL;
}

double func(double *para,int paranum)
{
	Enter("func");
	return 50*pow(para[1]-para[0]*para[0],2)+pow(1-para[0],2);
}

int get_gk(double *para,int paranum,double *gk)
{
	Enter("get_gk");
	gk[0]=-200*para[0]*(para[1]-para[0]*para[0])-2*(1-para[0]);
	gk[1]=100*(para[1]-para[0]*para[0]);
	muticount+=6;
	addcount+=4;
	return SUCCESS;
}

 int get_G(double *para,int paranum,double *G)
 {
 	Enter("get_G");
	G[0]=-200*para[1]+600*para[0]*para[0]+2;
	G[1]=-200*para[0];
	G[2]=-200*para[0];
	G[3]=100.0;
	muticount+=5;
	addcount+=2;
	return SUCCESS;
 }

#define initcount	muticount=0;\
			     addcount=0
	
int example()
{
	double e=1E-12;

	double para[2]={0,0};
	initcount;
	printf("FR(不精确一维搜索):\n");
	fr(e,wolfe,get_gk,func,para,2,0);
	initcount;
	printf("阻尼Newtom法(不精确一维搜索):\n");
	newton(e,wolfe,get_gk,get_G,func,para,2);
	initcount;
	printf("BFGS(不精确一维搜索):\n");
	bfgs(e,wolfe,get_gk,func,para,2);	
	initcount;
	printf("最速下降法(不精确一维搜索):\n");
	fasterdown(e,wolfe,get_gk,func,para,2);
	initcount;
	printf("FR(N步重新开始,不精确一维搜索):\n");
	fr(e,wolfe,get_gk,func,para,2,1);
	initcount;
	printf("PRP:\n");
	prp(e,wolfe,get_gk,func,para,2,0);
	printf("PRP(N步重新开始,不精确一维搜索):\n");
	prp(e,wolfe,get_gk,func,para,2,1);
	initcount;
	printf("DFP(不精确一维搜索):\n");
	dfp(e,wolfe,get_gk,func,para,2);
	initcount;
	FirstStep=0.1;
	printf("阻尼Newtom法(精确一维搜索):\n");
	newton(e,forwardandback,get_gk,get_G,func,para,2);

	FirstStep=0.001;
	initcount;
	printf("最速下降法(精确一维搜索):\n");
	fasterdown(e,forwardandback,get_gk,func,para,2);

	printf("FR(精确一维搜索):\n");
	initcount;
	fr(e,forwardandback,get_gk,func,para,2,0);
	printf("FR(N步重新开始,精确一维搜索):\n");
	initcount;
	fr(e,forwardandback,get_gk,func,para,2,1);

	initcount;
	printf("PRP(精确一维搜索):\n");
	prp(e,forwardandback,get_gk,func,para,2,0);
	printf("PRP(N步重新开始,精确一维搜索):\n");
	initcount;
	prp(e,forwardandback,get_gk,func,para,2,1);
	initcount;
	printf("BFGS(精确一维搜索):\n");
	bfgs(e,forwardandback,get_gk,func,para,2);	

	initcount;
	printf("DFP(精确一维搜索):\n");
	dfp(e,forwardandback,get_gk,func,para,2);
	return SUCCESS;
}

int main()
{
	example();
	return SUCCESS;
}