非线性优化--NLopt算法使用及C++实例
看这篇之前建议先看这篇,里面讲了非线性优化的原理即相关名词的概念,然后介绍了NLopt的使用方法,这个方法是基于C语言的,本片介绍一个NLopt的实例,用的C++语言。
在实例之前,先介绍下NLopt支持的算法,以及算法使用的注意事项
NLopt 支持的算法
NLopt 包含很多种不同的优化算法。
在头文件里面算法名称的枚举类型为
enum algorithm {
GN_DIRECT = 0,
GN_DIRECT_L,
GN_DIRECT_L_RAND,
GN_DIRECT_NOSCAL,
GN_DIRECT_L_NOSCAL,
GN_DIRECT_L_RAND_NOSCAL,
GN_ORIG_DIRECT,
GN_ORIG_DIRECT_L,
GD_STOGO,
GD_STOGO_RAND,
LD_LBFGS_NOCEDAL,
LD_LBFGS,
LN_PRAXIS,
LD_VAR1,
LD_VAR2,
LD_TNEWTON,
LD_TNEWTON_RESTART,
LD_TNEWTON_PRECOND,
LD_TNEWTON_PRECOND_RESTART,
GN_CRS2_LM,
GN_MLSL,
GD_MLSL,
GN_MLSL_LDS,
GD_MLSL_LDS,
LD_MMA,
LN_COBYLA,
LN_NEWUOA,
LN_NEWUOA_BOUND,
LN_NELDERMEAD,
LN_SBPLX,
LN_AUGLAG,
LD_AUGLAG,
LN_AUGLAG_EQ,
LD_AUGLAG_EQ,
LN_BOBYQA,
GN_ISRES,
AUGLAG,
AUGLAG_EQ,
G_MLSL,
G_MLSL_LDS,
LD_SLSQP,
LD_CCSAQ,
GN_ESCH,
NUM_ALGORITHMS /*不是一种算法 只是算法的数量*/
};
命名规律:
G/L代表的就是 全局(global)或者局部(local)优化
N/D代表的就是 不需导数 或者 需要导数 的优化
例如 LN_COBYLA 就是用的 COBYLA 算法 ,然后该算法用于局部(L)无需导数(N)的优化
算法选择
那么我们在使用的时候用哪种算法呢?随意选择吗?NO
对于一个确定的数学模型的优化问题,比较好的方式是 通过比较几种可用的算法,再确定选哪种,因为没有最佳算法,不同优化问题对于最佳的解决方式不同。
在比较算法的时候也要注意,不同算法适配的 函数值容差和参数容差 也应该是不同的。
有一种比较好的比较两种算对于一个数学模型的优劣的方式:
先运行一个算法得到最小值,然后运行第二个算法的时候设定就收敛到那个值,然后比较运行的时间
选择全局优化要注意的问题
当前,所有全局优化算法都要求对所有优化参数指定边界约束。
并且一定要注意的是
在这些算法中 只有 ISRES, AGS, 和 ORIG_DIRECT 支持非线性不等式约束
只有 ISRES 支持非线性等式约束
有一种很好的做法就是通过全局算法找到一个最优点,然后以这个最优点作为局部优化的起点,使结果更准确
全局优化算法在找到最优点时 比局部优化算法要 花的精力多很多
Code
数学模型:
求 :max(lnx1+lnx2) 目标函数
约束: p1x1+p2x2=5 等式约束
x1<=x2 不等式约束
x1>=0;x2>=0 边界约束
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#include <nlopt.hpp> //nlopt的头文件
引用头文件
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//声明一个 是否使用梯度的 全局变量
bool grad_bool=1;
声明一个 是否使用梯度的 全局变量
这样方便在切换 基于梯度的算法和不需要导数的算法的时候可以 改下这个变量就可以
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
double myfunc(const std::vector<double>& x,std::vector<double>&grad,void* f_data )
{
double result;//声明结果
if (grad_bool)
{
grad[0] = 1 / x[0];
grad[1] = 1 / x[1];
}
result = log(x[0])+log(x[1]);//计算结果
return result;//返回结果
}
目标函数是 max(lnx1+x2)
定义目标函数 注意函数的参数的形式 是不可以变化的
参数x 是 要优化的参数向量
grad返回为此时最优化参数的梯度 就是数学模型里里 对x求偏导的结果 如果要求偏导那么 lnx1+lnx2 对x1求偏导,及1/x1 也就是 代码的 grad[0] = 1 / x[0]; 这个地方
f_data 是要传入的参数
返回值 为在一个x向量下的目标函数的值
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
double myconstraint(const std::vector<double>& x,std::vector<double>&grad,void* f_data)
{
double *p = (double *) f_data;
if (grad_bool)
{
grad[0] = p[0];
grad[1] = p[1];
}
double result;//声明结果
result = x[0] * p[0] + x[1] * p[1] - 5;//计算结果
return result;//返回结果
}
定义等式约束函数
数学模型的等式是:p1x1+p2x2=5
p1和p2以外部参数的方式传进来
最好的方式换种写法,这样什么样的参数格式都可以套用
typedef struct{
double p1,p2;
}my_constraint_data;
声明要传入参数的结构体
double myconstraint(const std::vector<double>& x,std::vector<double>&grad,void* f_data)
{
my_constraint_data *d = (my_constraint_data*)f_data;
//double *p = (double *) f_data;
double p1 = d->p1;
double p2 = d->p2;
if (grad_bool)
{
grad[0] = p1;
grad[1] = p2;
}
double result;//声明结果
result = x[0] * p1 + x[1] * p2 - 5;//计算结果
return result;//返回结果
}
然后取参数的时候这样就可以
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
double myinconstraint(const std::vector<double>& x,std::vector<double>&grad,void* f_data)
{
if (grad_bool)
{
grad[0] = 1;
grad[1] = -1;
}
double result;//声明结果
result = x[0] - x[1] ;//计算结果
return result;//返回结果
}
不等式的约束函数
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
int main(int argc, char** argv)
{
//初始化节点
ros::init(argc, argv, "lidar_align");
//声明两个句柄
ros::NodeHandle nh;
然后就是main函数的,我是在ros下做的,所以加个节点的初始化
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//数学模型里的参数
//double p[2]={1,2};
my_constraint_data data = {1,2};
//待优化 求取的参数
std::vector<double> x{1,1};
声明数学模型里面 等式约束的 参数 my_constraint_data这个结构体在上面定义了
声明 待优化 求取的参数 1,1相当于初始化了
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//声明一个优化器
nlopt::opt opter;
声明一个NLopt的优化器
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/*局部优化算法*/
//opter = nlopt::opt(nlopt::LD_SLSQP,2); //有导数
//opter = nlopt::opt(nlopt::LN_COBYLA,2); //无导数
/*全局优化算法*/
opter = nlopt::opt(nlopt::GN_ISRES,2); //无导数 是可以的 x1 = 1.66667 x2= 1.66667 fmax = 1.02165
给优化器设置使用什么优化算法,并设置因子个数,就是x的个数
优化算法在nlopt.hpp里面有, 可以从里面选择想用的
但是注意并不是所有算法都可以计算出来,例如这个数学模型有 不等式约束和等式约束的 全局优化算法里面只有ISRES可以使用
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/*优化参数的边界约束*/
std::vector<double> lb {1,1};//注意参数的个数要对应上
std::vector<double> ub {10000,10000};//注意参数的个数要对应上
//设置 参数 边界
opter.set_lower_bounds(lb);//设置参数下限
opter.set_upper_bounds(ub);//设置参数上限
优化参数的边界约束
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/*设置优化精度*/
double tol = 1e-8;//容差
opter.set_xtol_rel(tol);
opter.set_force_stop(tol);
设置优化精度
控制优化什么时候停止
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/*设置目标函数 其中目标函数是自己定义的 */
// myfunc 是自己定义的函数名字 ,第二个参数为 可传入数据 没有数据则为NULL
opter.set_max_objective(myfunc,NULL);
设置目标函数
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/*设置数学模型的 等式 约束 和 不等式 约束 */
//myconstraint 是自己定义的函数名字 ,第二个参数为 可传入数据,
opter.add_equality_constraint(myconstraint,&data,tol);
opter.add_inequality_constraint(myinconstraint, NULL,tol);
设置数学模型的 等式 约束 和 不等式 约束
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/* 获取结果 最大值储存在f_max内 对应的向量储存在x内*/
double f_max = -10000;
nlopt::result res = opter.optimize(x,f_max);
执行优化计算
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
std::cout<<"x1 = " <<x[0]<<" x2= "<<x[1]<<" fmax = "<<f_max<<std::endl;
打印计算结果
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~