clpex,GLPK,lpsolve,yalmip简介与教程

合集下载

相关主题

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

最近建立了一个网络流模型，是一个混合整数线性规划问题（模型中既有连续变量，又有整型变量）。

当要求解此模型的时候，发现matlab优化工具箱竟没有自带的可以求解这类问题的算法（只有bintprog求解器，但是只能求解不含连续变量的二值线性规划问题）。

于是在网上找了一些解决问题的途径，下面说说几种可能的解决方案。

cplex
首先想到的是IBM公司大名鼎鼎的cplex。

cplex是IBM公司一款高性能的数学规划问题求解器，可以快速、稳定地求解线性规划、混合整数规划、二次规划等一系列规划问题。

CPLEX 的速度非常快，可以解决现实世界中许多大规模的问题，它能够处理有数百万个约束(constraint) 和变量 (variable) 的问题，而且一直刷新数学规划的最高性能记录。

他的标准版本是一个windows 下的IDE应用软件，但是开发人员能通过组件库从其他程序语言调用CPLEX 算法。

随标准版本一起发布的文件中包含一个名为matlab文件夹，将此文件夹添加到matlab的搜索路径下就可以在matlab下调用cplex高效地求解数学规划问题。

cplex IDE主界面（是不是很熟悉的界面？没错，cplex也是基于eclipse插件机制开发的。

）：
CPLEX Optimizer中文介绍：/components/detailvi ew.aspx?id=ce16c50e-0059-417b-9806-c8b1d3224084
官方网址：/components/detailview.aspx?id=ce16c 50e-0059-417b-9806-c8b1d3224084
遗憾的是，cplex是一款商业软件，可以从以上官方网址上下载免费试用版，使用时限是90天，而且试用版对问题规模有限制（我的问题有300个变量，370个约束，结果因为问题规模限制无法用试用版求解）。

如果你要用cpl ex解决问题的话，可能还需要学习特定于cplex的建模语言。

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

GLPK
在放弃了cplex之后搜寻其他解决方案的时候，我想起了GLPK。

GLPK (GNU Linear Programming Kit:GNU线性编程工具)是GUN下的一个项目，用于建立大规模线性规划LP和混合型整数规划MIP问题，并对模型进行最优化求解。

由于是GNU下的项目，因此没有商业非商业的版本限制，可以自由使用。

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

如果要在matlab下使用，还需要下载额外的驱动文件。

GLPK英文介绍：/software/glpk/
GLPK for windows：/
lpsolve
在弄了一阵GLPK无果之后，我又转投lpsolve了。

lpsolve是sourceforge下的一个开源项目，它的介绍如下：
Mixed Integer Linear Programming (MILP) solver lp_solve solves pure lin ear, (mixed) integer/binary, semi-cont and special ordered sets (SOS) mo dels.lp_solve is written in ANSI C and can be compiled on many differen t platforms like Linux and WINDOWS
它是一个混合整数线性规划求解器，可以求解纯线性、（混合）整数/二值、半连续和特殊有序集模型。

并且经过实际验证，有极高的求解效率。

sourceforge主页：/projects/lpsolve/?source=directory 从以上主页上可以下载lpsolve的IDE版本，界面比较简陋，类似于如下的样子：
如果要在matlab下使用lpsolve，需要在网址/projects/lp solve/files/lpsolve/5.5.2.0/ 提供的文件列表中下载类似于lp_solve_5.5.2.0_MA TLAB_exe_win32的zip文件。

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

将下载的zip解压后，得到以下文件结构：
bin目录下有matlab插件所必须的.mexw32文件和函数库API（.dll）。

ex 开头的一系列文件是自带的一些demo，教你如何在matlab下建模和求解。

mx lpsove.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下。

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文档的网址（重要，一定要看）：/lpsolve/doc/MATLAB.htm
从以上的过程我们看到用 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,setmi nim)
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_maker
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.
此时，同样解决以上线性规划问题，可以用如下语句简化过程：
>> 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)
或者只需要一句：
>> [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/p mwiki.php?n=Main.Download
可以说，yalmip是一位“集大成者”，它不仅自己包含基本的线性规划求解算法，比如linprog（线性规划）、bintprog（二值线性规划）、bnb（分支界定算法）等，他还提供了对cplex、GLPK、lpsolve等求解工具包更高层次的包装。

你只需要知道在matlab下如何用yalmip的方法建模和求解，就可以轻松地将你的求解算法替换成cplex、GLPK、lpsolve等，而不需要单独针对每一种工具包学习建模语法；而且yalmip的建模语法非常简单，简单到你只需要记住四个命令就可以了：
1. 创建决策变量：
>> x = sdpvar(m, n [, option])：创建m*n的连续型决策变量矩阵，optio n是对矩阵的一些参数指定。

相应的，如果要创建整型或二值型决策变量，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');
如果要继续添加约束也非常简单，支持用+直接相连（类似于.net 里面的事件处理程序）：
>> 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，一切都变得简单起来。