CppAD优化问题求解

参考文章：
（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$
约束：
$g_1(x)=x_1.x_2.x_3.x_4\ge 25$
$g_2(x)=x_1^2+x_2^2+x_3^2+x_4^2=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适用于需要计算导数的数值计算，比如优化、拟合、求解微分方程等。

---