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

最新推荐文章于 2024-09-15 07:39:36 发布

碎步の流年

最新推荐文章于 2024-09-15 07:39:36 发布

阅读量5.1k

点赞数 4

分类专栏： C++

原文链接：https://blog.csdn.net/u013468614/article/details/103624851

版权

C++ 专栏收录该内容

43 篇文章 10 订阅

订阅专栏

本篇针对非线性规（优）划（化）问题，介绍两个好用的模型构建不求解工具：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求解以下非线性规划问题：
$\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;
}

对应的CMakeLists.txt，为

cmake_minimum_required( VERSION 2.8 )
project(test)

set( CMAKE_BUILD_TYPE "Release" )
set( CMAKE_CXX_FLAGS "-O3" )
set( CMAKE_CXX_FLAGS "-std=c++11")

# 添加头文件
include_directories("/usr/local/include")
# 添加库
link_directories( /usr/local/lib )
# 添加需要编译的可执行文件
add_executable(test test.cpp)
# 链接到对应的库
target_link_libraries(${PROJECT_NAME} ipopt)

输出如下：

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