【CppAD】优化问题求解例程

示例来源：https://cppad.readthedocs.io/ipopt_solve_get_started.cpp.html

使用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 优化后的结果保存在这里

举个例子

定义优化问题的目标函数和约束表达形式：

假设一个4维优化变量组成的优化问题，待求解的目标函数为：
$$f(x)=x_1\cdot x_4\cdot (x_1+x_2+x_3)+x_3$$
同时这个优化问题有两个约束限制了变量的范围：
$$\begin{gathered} g_1(x)=x_1\cdot x_2\cdot x_3\cdot x_4\ge 25 \\ {g}_2(x)=x_1^2+x_1^2+x_1^2+x_1^2=40 \end{gathered}$$
根据这三个表达式，可以构建一个仿函数来描述这个优化问题的表达式：

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;
         // g_1 (x)
         fg[1] = x1 * x2 * x3 * x4;
         // g_2 (x)
         fg[2] = x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4;
         //
         return;
      }
   };
}

在上面的代码中，fg变量存储了三个表达式，fg[0]中存储了目标函数的表达式，fg[1] 和 fg[2] 分别存储了两个约束的表达式，约束表达式的上限与下限将会在之后给出。

定义优化变量的的初始值和上下限

除了两个优化变量之外，需要给独立变量设置初值xi，并限制独立变量的范围，每一个变量对应的上下限根据索引值对应到了相应的xu,xl数组元素中。
在这个例子中，假设每个变量都用相同的取值范围[1, 5]，即：
$$1\le x_1,x_2,x_3,x_4\le 5$$

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;
   }

定义约束表达式的上下限

之前在FG_eval中定义了约束的表达式，$g_1(\mathbf{x})>25$, $g_2(\mathbf{x})=40$
gl,gu中索引为0的量对应了fg[1]即$g_1(\mathbf{x})$的上下限，同理索引1对应了fg[2]即$g_2(\mathbf{x})$的上下限

 // 约束表达式的上下限
   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";

其他优化器的设置项可以在这里查询到：ipopt_options

优化结果

	// 在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]获取对应的优化结果。

最后，如果要系统学习CppAD，还是要去CppAD文档进行学习

完整的程序

#include <iostream>
#include <cppad/ipopt/solve.hpp>

// 这是CppAD官方网站的一个例子，具体计算了：
// object function: x_1 * x_4 * (x_1 + x_2 + x_3) + x_3
// 约束条件为：
// x_1 * x_2 * x_3 * x_4 >= 25
// x_1 ^ 2 + x_2 ^ 2 + x_3 ^ 2 + x_4 ^ 2 = 40
// x_1 , x_2 , x_3 , x_4 取值范围均为[1,5]

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;
}