参考文章:
(1)【CppAD】优化问题求解例程
(2)CppAD说明文档
(3)自动驾驶中的模型预测控制(Model Predictive Control)
0 前言
使用CppAD并使用ipopt求解一个优化问题,最核心的函数:
CppAD::ipopt::solve<Dvector, FG_eval>(options, xi, xl, xu, gl, gu, fg_eval, solution);
option
:求解器的配置项
xi
:待求解的优化变量(包含初值)
xl
:待优化变量的下限
xu
:待优化变量的上限
gl
:约束表达式的下限
gu
:约束表达式的上限
fg_eval
:存储目标函数和约束变量的表达式
solution
:优化后的结果保存在这里
1 使用CppAD
定义优化问题的目标函数和约束表达形式:
目标函数
:
f
(
x
)
=
x
1
.
x
4
.
(
x
1
+
x
2
+
x
3
)
+
x
3
f(x)=x_1. \ x_4.(x_1+x_2+x_3)+x_3
f(x)=x1. x4.(x1+x2+x3)+x3
约束
:
g
1
(
x
)
=
x
1
.
x
2
.
x
3
.
x
4
≥
25
g_1(x)=x_1.x_2.x_3.x_4≥25
g1(x)=x1.x2.x3.x4≥25
g
2
(
x
)
=
x
1
2
+
x
2
2
+
x
3
2
+
x
4
2
=
40
g_2(x)=x^2_1+x^2_2+x^2_3+x^2_4=40
g2(x)=x12+x22+x32+x42=40
构建一个仿函数来描述这个优化问题的表达式:
namespace {
using CppAD::AD;
class FG_eval {
public:
typedef CPPAD_TESTVECTOR( AD<double> ) ADvector;
void operator()(ADvector& fg, const ADvector& x)
{
assert(fg.size() == 3);
assert(x.size() == 4);
// Fortran style indexing
AD<double> x1 = x[0];
AD<double> x2 = x[1];
AD<double> x3 = x[2];
AD<double> x4 = x[3];
// f(x)
fg[0] = x1 * x4 * (x1 + x2 + x3) + x3;
fg[1] = x1 * x2 * x3 * x4;
fg[2] = x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4;
return;
}
};
}
定义优化变量的的初始值和上下限
typedef CPPAD_TESTVECTOR( double ) Dvector;
// 独立变量的个数,在这个例子中,为x1~x4共4个独立变量
size_t nx = 4;
// 给独立变量设置初始值
Dvector xi(nx);
xi[0] = 1.0;
xi[1] = 5.0;
xi[2] = 5.0;
xi[3] = 1.0;
// 设置独立变量的上下限
Dvector xl(nx), xu(nx);
for(i = 0; i < nx; i++)
{ xl[i] = 1.0;
xu[i] = 5.0;
}
定义约束表达式的上下限
Dvector gl(ng), gu(ng);
gl[0] = 25.0; gu[0] = 1.0e19;
gl[1] = 40.0; gu[1] = 40.0;
优化器的设置
std::string options;
//关闭任何优化器的输出
options += "Integer print_level 0\n";
options += "String sb yes\n";
//设置最大迭代次数
options += "Integer max_iter 10\n";
// 设置收敛条件
// see Mathematical Programming, Volume 106, Number 1,
// Pages 25-57, Equation (6)
options += "Numeric tol 1e-6\n";
优化结果
// 在solution中保存优化结果
CppAD::ipopt::solve_result<Dvector> solution;
// 执行优化
CppAD::ipopt::solve<Dvector, FG_eval>(options, xi, xl, xu, gl, gu, fg_eval, solution);
// 判断求解是否成功
bool ok = true;
ok &= solution.status == CppAD::ipopt::solve_result<Dvector>::success;
如果要查看/使用优化的结果,可以通过 solution.x[i] 获取对应的优化结果。
完整的程序:
#include <iostream>
#include <cppad/ipopt/solve.hpp>
namespace
{
using CppAD::AD;
class FG_eval
{
public:
typedef CPPAD_TESTVECTOR(AD<double>) ADvector;
void operator()(ADvector &_fg, const ADvector &_x)
{
assert(_fg.size() == 3);
assert(_x.size() == 4);
AD<double> x1 = _x[0];
AD<double> x2 = _x[1];
AD<double> x3 = _x[2];
AD<double> x4 = _x[3];
// f(x) 描述object function
_fg[0] = x1 * x4 * (x1 + x2 + x3) + x3;
// g_1 (x) 描述第一个不等式约束
_fg[1] = x1 * x2 * x3 * x4;
// define g_2 (x) 描述第二个等式约束
_fg[2] = x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4;
return;
}
};
bool get_started(void)
{
bool ok = true;
size_t i;
typedef CPPAD_TESTVECTOR(double) DVector;
/* 独立变量的个数 */
size_t nx = 4;
/* 约束的数量 */
size_t ng = 2;
/* 对独立变量进行初始化 */
DVector xi(nx);
xi[0] = 1.0;
xi[1] = 5.0;
xi[2] = 5.0;
xi[3] = 1.0;
/* 建立x的上下界约束 */
DVector xl(nx), xu(nx);
for (size_t i = 0; i < nx; i++)
{
xl[i] = 1.0;
xu[i] = 5.0;
}
/* 建立约束g的上下界 */
DVector gl(ng), gu(ng);
gl[0] = 25.0;
gu[0] = 1.0e19; // 不等式约束
gl[1] = 40.0;
gu[1] = 40.0; // 等式约束
FG_eval fg_eval;
/* 添加配置文件 */
std::string options;
// 关闭所有printing
options += "Integer print_level 0\n";
options += "String sb yes\n";
// 最大迭代次数
options += "Integer max_iter 10\n";
// 使用一阶近似条件判断收敛
options += "Numeric tol 1e-6\n";
/* 设置保存结果的对象 */
CppAD::ipopt::solve_result<DVector> solution;
/* 开始解决问题 */
CppAD::ipopt::solve<DVector, FG_eval>(
options, xi, xl, xu, gl, gu, fg_eval, solution);
/* 观察结果 */
ok &= solution.status == CppAD::ipopt::solve_result<DVector>::success;
std::cout << "solution: x1: " << solution.x[0]
<< " x2: " << solution.x[1]
<< " x3: " << solution.x[2]
<< " x4: " << solution.x[3] << std::endl;
std::cout << "object func value: " << solution.obj_value << std::endl;
return ok;
}
};
int main(int argc, char **argv)
{
get_started();
return 0;
}
CMakeList.txt
cmake_minimum_required(VERSION 3.0.2)
project(cppa)
find_package(catkin REQUIRED COMPONENTS
roscpp
)
set(CMAKE_CXX_STANDARD 14)
include_directories(/usr/local/include)
link_directories(/usr/local/lib)
add_executable(CppAD_demo src/CppAD_demo.cpp)
target_link_libraries(CppAD_demo ipopt)
🍎 在cppad库中,Dvector和ADvector的区别
- Dvector是一个普通的向量类,用于存储和操作双精度浮点数向量。它类似于标准库中的std::vector,提供了一系列成员函数和操作符重载,可以对向量进行元素访问、插入和删除等操作
- ADvector是一个自动微分向量类,用于存储和操作自动微分变量的向量。它是基于Dvector的扩展,支持自动微分功能。自动微分是一种计算导数的方法,在ADvector中,可以通过对向量中的某个元素求导来得到整个向量的导数。ADvector提供了一系列成员函数和操作符重载,用于实现自动微分计算
Dvector适用于普通的数值计算,而ADvector适用于需要计算导数的数值计算,比如优化、拟合、求解微分方程等。