欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

c++ IPOPT 库 与 python非线性规(优)划（化）工具

程序员文章站 2022-06-18 17:15:34
...

本篇针对非线性规（优）划（化）问题，介绍两个好用的模型构建不求解工具：CppAD与pyomo
。这两个工具都可以只构建目标函数与约束函数的形式，而不需要求其雅可比、海森矩阵。两个工具都可以接受非线性模型形式，并且有与各种优化求解器的接口。
首先声明，这篇文章是转载的，先上链接 https://blog.csdn.net/u013468614/article/details/103624851
另外，本文的安装IPOPT部分，详情见我的另一篇文章

1. 背景

系统的真实模型都是非线性，在实际模型构建或使用过程中，我们会对模型进行线性化。例如，我们根据已有的模型求解最优控制问题：让机器人根据当前的状态输出最优的动作，从而使成本函数取最小值。

常见的优化库都需要我们求目标函数与约速函数的导数、雅可比矩阵以及海森矩阵。如果模型简单，姑且可以手动求解，模型一旦复杂些（例如大型的非线性模型，约束函数为微分方程or差分方程形式），非常容易出错。虽然也有一些库能够用于自动微分，例如python的：sympy、autograd、tangent等，c++的:Automatic differentiation（AD）。这些库仅仅用来微分求导的话是非常合适的，如果要跟最优化求解器（也即优化库）结合，需要相当多的结口函数。同样的情况，当问题变得复杂，头就大了，非常容易出错，且不好调试。

针对上面的问题，本篇介绍c++与python中将自动微分与优化求解器融合一体的库：CppAD (for c++), pyomo (for python)。我们只需要确定模型的目标函数、约束函数的形式，这些工具就可以帮我们算出在约束条件下，使目标函数取最值的解。
重点：这两个工具都不需要我们对模型进行线性化，不需要手动求雅可比、海森矩阵等。

2. CppAD [for C++]

2.1 CppAD简介

CppAD 是Bradley M. Bell针对运筹学 COIN-OR 的计算基础部分进行的一个开源项目。CppAD特点如下：

  • 支持正向、反向的优化求解模式;
  • 任何阶数的导数;
  • 使用自动微分（AD）的功能；
  • 稀疏模式；
  • 数值库模版（可以与其它AD库一起用）；
  • 很多使用例程；
  • 支持windows 与 Linux系统

CppAD源码：https://github.com/coin-or/CppAD
说明文档：https://coin-or.github.io/CppAD/doc/cppad.htm#Features.Base%20Type
应用例程集锦：https://coin-or.github.io/CppAD/doc/listallexamples.htm

2.2 CppAD安装

Ubuntu ：

sudo apt-get install cppad

CppAD需要调用系统中已经安装好的优化求解器，如非线性优化器 ipopt、Gurobi、GLPK。由于本篇介绍使用CppAD的例子中用到ipopt，此处先介绍一下ipopt的安装。

Ubuntu下:

sudo apt-get install gfortran
sudo apt-get install unzip
wget https://www.coin-or.org/download/source/Ipopt/Ipopt-3.12.7.zip
unzip Ipopt-3.12.7.zip && rm Ipopt-3.12.7.zip
sudo bash install_ipopt.sh Ipopt-3.12.7

2.3 一个例子

这个例子利用CppAD与Ipopt求解以下非线性规划问题：
minimizex1x4(x1+x2+x3)+x3s.t.x1x2x3x4≥25x12+x22+x32+x42=401≤x1,x2,x3,x4≤5(1) \begin{matrix} minimize & x_1x_4(x_1+x_2+x_3)+x_3 \\ s.t. & x_1x_2x_3x_4 \geq 25 \\ &x_{1}^{2} + x_{2}^{2}+x_{3}^{2}+x_{4}^{2}=40 \\ & 1 \leq x_{1}, x_{2}, x_{3}, x_{4} \leq 5 \end{matrix} \tag{1}
完整的c++代码如下

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

using namespace std;

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);
        // variables
        AD<double> x1 = x[0];
        AD<double> x2 = x[1];
        AD<double> x3 = x[2];
        AD<double> x4 = x[3];
        // f(x) objective function
        fg[0] = x1 * x4 * (x1 + x2 + x3) + x3;
        // constraints
        fg[1] = x1 * x2 * x3 * x4;
        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; // number of varibles
    size_t ng = 2; // number of constraints
    Dvector x0(nx); // initial condition of varibles
    x0[0] = 1.0;
    x0[1] = 5.0;
    x0[2] = 5.0;
    x0[3] = 1.0;

    // lower and upper bounds for varibles
    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;
    // object that computes objective and constraints
    FG_eval fg_eval;

    // options
    string options;
    // turn off any printing
    options += "Integer print_level  0\n";
    options += "String sb            yes\n";
    // maximum iterations
    options += "Integer max_iter     10\n";
    //approximate accuracy in first order necessary conditions;
    // see Mathematical Programming, Volume 106, Number 1,
    // Pages 25-57, Equation (6)
    options += "Numeric tol          1e-6\n";
    //derivative tesing
    options += "String derivative_test   second-order\n";
    // maximum amount of random pertubation; e.g.,
    // when evaluation finite diff
    options += "Numeric point_perturbation_radius   0.\n";


    CppAD::ipopt::solve_result<Dvector> solution; // solution
    // solve the problem
    CppAD::ipopt::solve<Dvector, FG_eval>(options, x0, xl, xu, gl, gu, fg_eval, solution); 

    cout<<"solution: "<<solution.x<<endl;

    //
    //check some of the solution values
    //
    ok &= solution.status == CppAD::ipopt::solve_result<Dvector>::success;
    //
    double check_x[]  = {1.000000, 4.743000, 3.82115, 1.379408};
    double check_zl[] = {1.087871, 0.,       0.,       0.      };
    double check_zu[] = {0.,       0.,       0.,       0.      };
    double rel_tol    = 1e-6; // relative tolerance
    double abs_tol    = 1e-6; // absolute tolerance
    for(i = 0; i < nx; i++)
    {
        ok &= CppAD::NearEqual(
                    check_x[i], solution.x[i], rel_tol, abs_tol);
        ok &= CppAD::NearEqual(
                    check_zl[i], solution.zl[i], rel_tol, abs_tol);
        ok &= CppAD::NearEqual(
                    check_zu[i], solution.zu[i], rel_tol, abs_tol);
    }

    return ok;
}

int main()
{
    cout << "CppAD : Hello World Demo!" << endl;
    get_started();
    return 0;
}

输出如下：

CppAD : Hello World Demo!
solution: { 1, 4.743, 3.82115, 1.37941 }

我们只需要设置初始参数、初始状态、状态约束值、等式与不等式约束值、优化目标方程。其它的就交给CppAD吧！

3. pyomo [for Python]

3.1 pyomo简介

pyomo与CppAD的功能一致，最终我们仅仅需要把待求解问题构建好（可以是非线性，不需要手动求微分），然后扔给它就行了。但是，pyomo比CppAD好的一点是：pyomo具有深刻的模型概念，它把整个问题就当作是一个模型。CppAD虽然功能与pyomo相似，但是却没有这样一个比较让容易把握整体的概念。

3.2 pyomo安装

强列建议使用anaconda环境！
本人环境：Ubuntu 18.04 + Anaconda2.7 （python 2.7）
本人开始是尝试pip install来安装配置pyomo，出现很多问题，例如：能装上包，但是运行出错（可能是版本问题）。因此，本人干脆安装pyomo安装教程使用Anaconda环境。然后一切非常顺利。
Anaconda的一个比较全的安装教程。
在Ubuntu + Anaconda环境下安装pyomo与ipopt：

conda update conda
conda update anaconda
conda install -c conda-forge pyomo
conda install -c conda-forge pyomo.extras
conda install -c conda-forge coincbc
conda install -c conda-forge ipopt

当然，也可以安装基它可能会用到的优化求解库：

conda install -c conda-forge glpk

3.3 一个例子

此处仍采用式(1)作为求解对象来给出pyomo的一个应用例子。

import numpy as np
from pyomo.environ import *
from pyomo.dae import *

model = ConcreteModel()
model.x = Var(RangeSet(1, 4), bounds=(1, 25))
model.cons1 = Constraint(rule=lambda model: 40==model.x[1]**2+model.x[2]**2+model.x[3]**2+model.x[4]**2)
model.cons2 = Constraint(rule=lambda model: 25<=model.x[1]*model.x[2]*model.x[3]*model.x[4])
model.obj = Objective(expr = model.x[1]*model.x[4]*(model.x[1] + model.x[2] + model.x[3]) + model.x[3], sense=minimize) 
SolverFactory('ipopt').solve(model)
solutions = [model.x[i]() for i in range(1, 5)]
print("solutins is :"+str(solutions))

程序输出：

pyomo：Hello World!
solution: { 1.0, 4.742999644013376, 3.821149978907638, 1.3794082901545972 }

对比CppAD与pyomo求解的结果(如下)，两者是一致的（但是，pyomo的代码更为简洁）：

CppAD : Hello World Demo!
solution: { 1, 4.743, 3.82115, 1.37941 }

CppAD与pyomo时间成本比较：

pyomo：time cost is:0.0659720897675
CppAD: time cost is:0.004138

总结

为了简化求解最优问题的过程，本篇介绍了两种非线性规划工具：CppAD for C++；pyomo for python。
两者都能省去求雅可比、海森矩阵的过程。我们只需要定义好求解的模型即可。两个工具在同一问题上求解结果一致，但是CppAD比pyomo时间成本更小。但是，pyomo在实现代码量上明显小于CppAD,这是pyomo的优点。如果系统对实时性要求高，建议采用C++编程+CppAD，如果只是平时做简单算法测试，建议采用更简洁的python+pyomo组合。

————————————————
版权声明：本文为CSDN博主「风起了」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/u013468614/java/article/details/103624851

相关标签： 常用工具

上一篇： 非关系型数据库redis

下一篇： 如何在word中添加空白页(最前面无页眉页脚的空白页)

推荐阅读