yalmip + lpsolve + matlab 求解混合整数线性规划问题（MIP/MILP）

　　最近建立了一个网络流模型，是一个混合整数线性规划问题（模型中既有连续变量，又有整型变量）。当要求解此模型的时候，发现matlab优化工具箱竟没有自带的可以求解这类问题的算法（只有bintprog求解器，但是只能求解不含连续变量的二值线性规划问题）。于是在网上找了一些解决问题的途径，下面说说我试过的几种可能的解决方案，包括cplex、GLPK、lpsolve 和 yalmip。

cplex

　　首先想到的是IBM公司大名鼎鼎的cplex。cplex是IBM公司一款高性能的数学规划问题求解器，可以快速、稳定地求解线性规划、混合整数规划、二次规划等一系列规划问题。CPLEX 的速度非常快，可以解决现实世界中许多大规模的问题，它能够处理有数百万个约束 (constraint) 和变量 (variable) 的问题，而且一直刷新数学规划的最高性能记录。他的标准版本是一个windows下的IDE应用软件，但是开发人员能通过组件库从其他程序语言调用 CPLEX 算法。随标准版本一起发布的文件中包含一个名为matlab文件夹，将此文件夹添加到matlab的搜索路径下就可以在matlab下调用cplex高效地求解数学规划问题。

　　CPLEX Optimizer中文介绍：http://www.sstc.org.cn/components/detailview.aspx?id=ce16c50e-0059-417b-9806-c8b1d3224084

　　官方网址：http://www.sstc.org.cn/components/detailview.aspx?id=ce16c50e-0059-417b-9806-c8b1d3224084

　　遗憾的是，cplex是一款商业软件，可以从以上官方网址上下载免费试用版，使用时限是90天，而且试用版对问题规模有限制（我的问题有300个变量，370个约束，结果因为问题规模限制无法用试用版求解）。如果你要用cplex解决问题的话，可能还需要学习特定于cplex的建模语言。

　　值得一提的是，IBM公司一直对学术界有或多或少的支持，要想使用完整版的cplex，你可以参与IBM的学院计划，前提条件是你是大学/研究机构的老师/研究员，或者IBM公式的职员，通过这个网址：http://www-03.ibm.com/ibm/university/academic/pub/page/ban_ilog_programming? ，填写一个申请表格，通过审核之后你就有权限使用cplex的完整版，没有任何限制，和商业版完全一样的功能。

　　由于没钱买软件，试用版有规模限制，又是个学生不能参与学院计划，只好放弃这一途径==。

GLPK

　　在放弃了cplex之后搜寻其他解决方案的时候，我想起了GLPK。GLPK (GNU Linear Programming Kit，GNU线性编程工具)是GNU下的一个项目，用于建立大规模线性规划LP和混合型整数规划MIP问题，并对模型进行最优化求解。由于是GNU下的项目，因此没有商业非商业的版本限制，可以自由使用。

　　GLPK实现了对windows的支持，但是为此，你同样需要学习它的建模语言，并且所有的操作都在 glpsol.exe 提共的命令行下完成，比较不方便，且耗时长。如果要在matlab下使用，还需要下载额外的驱动文件。

　　GLPK英文介绍：http://www.gnu.org/software/glpk/

　　GLPK for windows：http://winglpk.sourceforge.net/

lpsolve（详细介绍，最好结合后面介绍的“yalmip”工具箱一起看）

　　在弄了一阵GLPK无果之后，我又转投lpsolve了。lpsolve是sourceforge下的一个开源项目，它的介绍如下：

　　Mixed Integer Linear Programming (MILP) solver lp_solve solves pure linear, (mixed) integer/binary, semi-cont and special ordered sets (SOS) models.lp_solve is written in ANSI C and can be compiled on many different platforms like Linux and WINDOWS

　　它是一个混合整数线性规划求解器，可以求解纯线性、（混合）整数/二值、半连续和特殊有序集模型。并且经过实际验证，有极高的求解效率。

　　sourceforge主页：http://sourceforge.net/projects/lpsolve/?source=directory

　　从以上主页上可以下载lpsolve的IDE版本，界面比较简陋，类似于如下的样子：

　　以上是用IDE工具建模求解，如果要在matlab下使用lpsolve，需要在网址http://sourceforge.net/projects/lpsolve/files/lpsolve/5.5.2.0/ 提供的文件列表中下载类似lp_solve_5.5.2.0_MATLAB_exe_win32（只针对windows 32位操作系统，其他操作系统请选择对应版本下载）的zip文件。

　　由于我的问题就是用的lpsolve解决的，在这里详细介绍一下，以lp_solve_5.5.2.0_MATLAB_exe_win32为例，过程如下：

　　　　1. 下载。将下载的zip解压后，得到以下文件结构：

　　　　　　bin目录下有matlab插件所必须的.mexw32文件和函数库API（.dll）。ex开头的一系列文件是自带的一些demo，教你如何在matlab下建模和求解。mxlpsove.m 是建模的核心函数，一个线性规划模型的所有配置和求解都是通过这个函数完成的。lp_maker.m 和 lp_solve.m 是对mxlpsolve.m的高层包装，简化了模型建立和求解的过程（后面会详细介绍）。

　　　　2. 准备驱动文件。在解压的bin目录下找到mxlpsolve.mexw32和mxlpsolve.dll两个文件，拷贝到解压根目录下（这两个文件就是matlab调用lpsolve的驱动文件），然后将此根目录添加到matlab搜索路径下（试试 pathtool 命令）。

　　　　3. 准备dll库文件。到这里不够，还需要lpsolve55.dll文件，真正求解问题的算法在这个函数库中。在lpsolve项目sourceforge首页下载安装一个IDE版本的程序，在安装目录下可以找到此dll文件，然后将此文件放到系统文件夹C:\Windows\System32下。也可以从我分享的这个链接下载到：http://vdisk.weibo.com/s/z4URgeDRTzBPd

　　　　4. 代码、求解。至此就可以在matlab下进尽情使用lpsolve了。以一个具体的例子说明用lpsolve求解数学规划问题的方法。

　　　　　假设我们要用matlab解决如下线性规划问题：

max 4x1 + 2x2 + x3
s. t. 2x1 + x2 <= 1
x1 + 2x3 <= 2
x1 + x2 + x3 = 1
x1 >= 0
x1 <= 1
x2 >= 0
x2 <= 1
x3 >= 0
x3 <= 2

　　　　　　matlab语句如下：

　　　　　　>> f = [4 2 1];
　　　　　　>> A = [2 1 0; 1 0 2; 1 1 1];
　　　　　　>> b = [1; 2; 1];
　　　　　　>> l = [ 0 0 0];
　　　　　　>> u = [ 1 1 2];

　　　　　　>> lp=mxlpsolve('make_lp', 1, 3);
　　　　　　>> mxlpsolve('set_verbose', lp, 3);
　　　　　　>> mxlpsolve('set_obj_fn', lp, f);
　　　　　　>> mxlpsolve('add_constraint', lp, A(1, :), 1, b(1));
　　　　　　>> mxlpsolve('add_constraint', lp, A(2, :), 1, b(2));
　　　　　　>> mxlpsolve('add_constraint', lp, A(3, :), 0, b(3);
　　　　　　>> mxlpsolve('set_lowbo', lp, l);
　　　　　　>> mxlpsolve('set_upbo', lp, u);
　　　　　　>> mxlpsolve('write_lp', lp, 'a.lp');
　　　　　　>> mxlpsolve('get_mat', lp, 1, 2)

　　　　　　>> mxlpsolve('solve', lp)

　　　　　　>> mxlpsolve('get_objective', lp)

　　　　　　>> mxlpsolve('get_variables', lp)

　　　　　　>> mxlpsolve('get_constraints', lp)
　　　　最后不要忘了用

　　　　　　>> mxlpsolve('delete_lp', lp)

　　　　释放模型占用的内存。

　　　　如果需还要其他功能，请参考包含完整API文档的网址（重要，推荐看!!!）：http://web.mit.edu/lpsolve/doc/MATLAB.htm 

　　从以上的过程我们看到用 lpsolve 建立一个规划问题的代码还是够麻烦的，想必你刚开始看到上面那些语句的时候，也是一头雾水。不要着急，对于这类简单的问题，还有更简便的方法。好在 lpsolve 为我们提供了一种简化的途径，我们注意到以上文件列表中有一个lp_maker.m和lp_solve.m文件。lp_maker.m文件的功能是创建一个（混合整数）线性规划问题，调用格式类似于其他matlab自带的优化工具箱，你只需要为它提供f、A、b、l、u几个矩阵，它会自动为你实现创建模型、设置目标函数、添加约束的过程。help一下可以看到如下帮助：

　　>> help lp_maker % 嫌麻烦可以跳过此段帮助信息
　　LP_MAKER Makes mixed integer linear programming problems.

　　SYNOPSIS: lp_handle = lp_maker(f,a,b,e,vlb,vub,xint,scalemode,setminim)
　　　　make the MILP problem
　　　　　　max v = f'*x
　　　　　　　　a*x <> b
　　　　　　　　　　vlb <= x <= vub
　　　　　　　　　　x(int) are integer

　　ARGUMENTS: The first four arguments are required:
　　　　f: n vector of coefficients for a linear objective function.
　　　　a: m by n matrix representing linear constraints.
　　　　b: m vector of right sides for the inequality constraints.
　　　　e: m vector that determines the sense of the inequalities:
　　　　　　e(i) < 0 ==> Less Than
　　　　　　e(i) = 0 ==> Equals
　　　　　　e(i) > 0 ==> Greater Than
　　　　vlb: n vector of non-negative lower bounds. If empty or omitted,
　　　　　　then the lower bounds are set to zero.
　　　　vub: n vector of upper bounds. May be omitted or empty.
　　　　xint: vector of integer variables. May be omitted or empty.
　　　　scalemode: Autoscale flag. Off when 0 or omitted.
　　　　setminim: Set maximum lp when this flag equals 0 or omitted.

　　OUTPUT: lp_handle is an integer handle to the lp created.

　　而 lp_solve.m 的调用格式与lp_maker.m类似，唯一的不同是，lp_solve.m 在创建模型的同时还求解模型，求解结果直接返回给输出参数。help一下帮助文档如下：

　　>> help lp_solve % 嫌麻烦可以跳过此段帮助信息

　　LP_SOLVE Solves mixed integer linear programming problems.
　　SYNOPSIS: [obj,x,duals,stat] = lp_solve(f,a,b,e,vlb,vub,xint,scalemode,keep)
　　　　solves the MILP problem
　　　　　　max v = f'*x
　　　　　　　　a*x <> b
　　　　　　　　vlb <= x <= vub
　　　　　　　　x(int) are integer

　　ARGUMENTS: The first four arguments are required:
　　　　f: n vector of coefficients for a linear objective function.
　　　　a: m by n matrix representing linear constraints.
　　　　b: m vector of right sides for the inequality constraints.
　　　　e: m vector that determines the sense of the inequalities:
　　　　　　e(i) = -1 ==> Less Than
　　　　　　e(i) = 0 ==> Equals
　　　　　　e(i) = 1 ==> Greater Than
　　　　vlb: n vector of lower bounds. If empty or omitted,
　　　　　　then the lower bounds are set to zero.
　　　　vub: n vector of upper bounds. May be omitted or empty.
　　　　xint: vector of integer variables. May be omitted or empty.
　　　　scalemode: scale flag. Off when 0 or omitted.
　　　　keep: Flag for keeping the lp problem after it's been solved.
　　　　　　If omitted, the lp will be deleted when solved

　　OUTPUT: A nonempty output is returned if a solution is found:
　　　　　obj: Optimal value of the objective function.
　　　　　x: Optimal value of the decision variables.
　　　　　duals: solution of the dual problem.

　　例如，同样解决以上线性规划问题，可以用如下语句简化过程(lp_maker版)：

　　>> f = [4 2 1];
　　>> A = [2 1 0; 1 0 2; 1 1 1];
　　>> b = [1; 2; 1];
　　>> l = [ 0 0 0];
　　>> u = [ 1 1 2];

　　>> lp = lp_maker(f, A, b, [-1; -1; 0], l, u, [], 1, 0);
　　>> solvestat = mxlpsolve('solve', lp)
　　solvestat =
     　　0
　　>> obj = mxlpsolve('get_objective', lp)
　　obj =
    2.5000
　　>> x = mxlpsolve('get_variables', lp)
　　x =
    　　0.5000
         　　0
    　　0.5000
　　>> mxlpsolve('delete_lp', lp)
或者只需要一句（lp_solve版）：
　　>> [obj,x,duals,stat] = lp_solve(f, A, b, [-1; -1; 0], l, u, [], 1, 0);

　　高层次的包装带来简便的同时也会让我们失去对问题更精细化的控制。例如，要使用 lp_solve.m 和 lp_maker.m，你必须事先知道约束系数矩阵A，然而对于很多实际问题，由于问题规模太大或者其他限制，你不能事先知道A矩阵，而是要用嵌套的for循环一步步建立起约束条件的时候，这两个高层包装就显得力不从心了。

yalmip（重点推荐！！）

　　最后登场的是yalmip，本来问题到 lpsolve 似乎就已经解决了，为什么还要介绍yalmip呢？因为我在解决这个问题的时候，其实是先遇到 yalmip，之后才遇到 lpsolve 的；再者，对于我的问题，lp_maker.m和lp_solve.m两个封装也无能为力；而且yalmip有它独特的优点，在这里不得不介绍。

　　此网址有它的详细介绍和下载链接：http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Main.Download （我猜测yalmip是耶鲁大学出品的，但是竟然找不到支持的证据。）

　　可以说，yalmip是一位“集大成者”，它不仅自己包含基本的线性规划求解算法，比如linprog（线性规划）、bintprog（二值线性规划）、bnb（分支界定算法）等，他还提供了对cplex、GLPK、lpsolve等求解工具包更高层次的包装。更为可贵的是，yalmip真正实现了建模和算法二者的分离，它提供了一种统一的、简单的建模语言，针对所有的规划问题，都可以用这种统一的方式建模；至于用哪种求解算法，你只需要通过一次简单的参数配置指定就可以了，甚至不用你指定，yalmip会自动为你选择最适合的算法。

　　总而言之，你只需要知道在matlab下如何用yalmip的方式建模，而不需要单独针对每一种工具包学习新的建模语法；而且yalmip 的建模语法非常简单，简单到你只需要记住四个命令就可以了：

1. 创建决策变量：

　　>> x = sdpvar(m, n [, option])：创建m*n的连续型决策变量矩阵，option是对矩阵的一些参数指定。

　　相应的，如果要创建整型或二值型决策变量，matlab语句分别为：　　

　　>> x = intvar(m, n, [option])

　　>> x = binvar(m, n, [option])

2. 添加约束：

　　>> F = set(constraint [, tag])：创建一个以constraint指定的约束，可选参数tag可以给该约束指定一个字符串标记。重要的是constraint的表达也非常简单，例如如果有 x1 + x2 + x3 <= 3 的约束，直接写：

　　>> x = sdpvar(3, 1);

　　>> F = set(x(1) + x(2) + x(3) <= 3, 'cost bound1');

　　如果要继续添加约束也非常简单，支持用+直接相连：

　　>> F = F + set(constraint1 [, tag1]);

　　>> F = F + set(constraint2 [, tag2]);……

　　例如，如果继续限制 x 只能取[0, 1]之间的值，则：

　　>> F = F + set(0 <= x <= 1, ‘upper and lower bound’);

　　一句话就搞定了，是不是非常简单。！

3. 参数配置

　　这个比较简单，语句如下：

　　>> ops = sdpsettings(option1, value1, option2, value2, ……)

　　例如语句

　　>> ops = sdpsettings('solver', 'lpsolve', 'verbose', 2);

　　'solver' 参数指定程序用lpsolve求解器（如果已经安装，否则会报错），如果不指定 ‘solver’ 参数，他会根据决策变量类型自动挑选已安装的、最适合的求解器；'verbose' 指定显示冗余度（冗余度越大，你就可以看到越详细的求解过程信息）。

4. 求解

　　这个也只有一句话：

　　>> result = solvesdp(F, f, ops) 求解一个数学规划（最小化）问题，该问题的目标函数由 f 指定，约束由 F 指定，ops指定求解参数，最后的结果存储在result结构体中。

还是以前面那个问题作为例子，如果用yalmip的话，只需要如下简单几句：

　　>> x = sdpvar(3, 1);
　　>> f = [4 2 1] * x;
　　>> F = set(2*x(1) + x(2) <= 1);
　　>> F = F + set(x(1) + 2 * x(3) <= 2);
　　>> F = F + set(x(1) + x(2) + x(3) == 1);
　　>> F = F + set(0 <= x(1) <= 1) + set(0 <= x(2) <= 1) + set(0 <= x(3) <= 2);
　　>> ops = sdpsettings('solver', 'lpsolve', 'verbose', 2);
　　>> result = solvesdp(F, -f, ops);

　　如果你想用 cplex 求解器求解，只需要将以上的‘solver’参数的‘lpsolve’改成‘cplex’，其他任何地方都不需要做改动。

　　除此以外，yalmip还支持几乎所有其他的求解算法，在matlab下输入yalmiptest命令可以得到所有支持的算法以及它们的安装状态（其中cplex和lpsolve是我安装的，其他status为found的表示默认支持，not found表示支持但matlab中还未安装）：

　　有了yalmip，你不再需要针对每一种工具包去学习特定的建模语言（比如用cplex要专门学习cplex的建模语言，用lingo要专门学习lingo的建模语言，还有GLPK、lpsolve、Matlab自带的求解器等等，如果每一种求解器都要学习新的建模语言的话，这个工作量是可想而知的）。相反，如果你选择使用yalmip，那么你只需要学习yalmip一种建模语法，因为yalmip真正实现了建模和算法的分离，所有的问题都可以用统一的方法建模，如果需要使用不同的求解器，只需要一句简单的配置即可。因此，yalmip不仅仅是一个线性规划求解器，更强大的地方在于，它提供了一个统一的建模平台，支持现有的几乎所有的求解算法。有了yalmip，一切都变得简单起来。

总结

　　我的问题总共有300个变量，其中120个连续变量，120个0-1变量，60个整型变量；另外还有370个约束（不包括变量本身的上下限、整型约束）。由于yalmip自带的bnb算法求解出了错误的结果（这是个已经reported的bug，在最近的release版本中修复了），而且效率极差，一次求解要花费大概5分钟的时间，最后我在matlab下，用 yalmip 建模，求解器采用 lpsolve 把问题解决了，而且由于用了lpsolve，效率得到极大的提升，每一次求解只花费不到5秒钟中的时间。

　　将以上四个工具总结如下：

工具	来源	优点	缺点
cplex	IBM	高效，快速，稳定，功能全面，适合大型商业化解决方案；可以通过学院计划免费获取。	付费，试用版有规模限制，需要学习特定建模语言。
GLPK	GNU	开源，免费；适合linux用户。	windows下繁琐，需要学习特定建模语言。
lpsolve	sourceforge	开源，免费。	需要学习特定建模语言；建模语言较繁琐。
yalmip	开放网络	网络开放获取，建模语言简单、统一，自带基本求解算法，并支持集成几乎所有其他求解器。推荐使用。	暂无。

posted on 2013-11-17 19:07 balabala已被注册阅读(89595) 评论(74) 编辑收藏举报

刷新页面返回顶部

kane1990