几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得

时间:2022-03-01 14:09:09

            最近在做network flow方面的优化问题,归纳起来是求解线性规划问题,于是尝试了几种优化工具,下面把自己的使用心得写下来,因为自己在搜集资料的时候发现网上这方面的好资源不是非常多,比如对优化工具的探讨大多在一个比较浅的层次上,我就深刻感觉到在使用中遇到问题往往除了官方资料很难找到答案,但官方资料太过庞杂,初学者不可能完全掌握,所以到现在对于有些工具我还有一些疑问,写出来希望有精通这方面的人可以给我解答补充。

       首先当然是Matlab自带的linprog,下面通过一个Max Flow的例子来简单说明一下,请看下图:

       几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得

      这是一个很经典的问题,所有讲到网络流问题一般都会用到这个例子,下面解释一下它的求解过程,这也基本是用linprog求解lp问题的代码格式,这里有五个变量(e1,e2,e3,e4,e5),求解每条边上的流数,f是目标函数,最大流其实就是从S点出去的流,是e1+e2边流的和,由于linprog中第一个变量默认求得是最小值,所以目标函数要取反;

下面就是约束条件了,A <= b是不等式约束,Aeq = beq是等式约束,具体到网络问题中,每条边的流量不能超过它的capacity constraint,对于每个结点,进入的流量应当等于流出的流量(流量守恒),这两类约束根据图看代码一目了然,最后的lb是变量的下界约束,这里当然是零,还可以有一个ub,表示变量的上界约束,最后输入到linprog中的变量一般有7个,分别就是前面分析得目标函数:f;不等式约束:A <= b;等式约束:Aeq = beq,变量上下界约束:lb,ub.

       这个例子当然很简单,但是将其扩展到复杂的网络拓扑时写代码其实非常不方便;不过让我放弃linprog最主要的原因还是跟后来的两个工具(lpsolve,gurobi)相比,同样的问题它解出的值总是不一样,比如网络流问题,lpsolve与gurobi可以解出整数值,但用linprog则一般解出的都是小数值,这对我的问题不是一个合适的解。再尝试其它的一些LP问题(非网络流),lpslove与gurobi的求解结果基本相似,但linprog求解结果和这两者有极大的差异,这种情况让我不大敢放心使用这种工具了。当然出现这种情况也有可能是我某个调用参数设置得不对,但我找过很多资料均不能回答我的问题,所以暂且搁置。

       lpsolve,yalmip的安装及使用说明可以参考这一篇博客的内容http://www.cnblogs.com/kane1990/p/3428129.html, lpsolve如果直接在它的IDE环境中使用是非常方便的,基本上LP表达式怎样写就可以怎样运行,举一个例子说明,例子来源http://www.mathcs.emory.edu/~cheung/Courses/323/Syllabus/NetFlow/,顺便说一下,这个网站上NetFlow的

 课件生动易懂,个人认为非常适合初学者。如下的一个图求解最大流问题(这幅图已经标注出求解结果了):

        几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得

lpsolve IDE环境中只需输入如下的文本,是不是非常简单直观:

            几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得         几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得

然后按一个运行按钮(红框标注),出现以下的求解结果,红框标注的是结果,蓝框标注的是求解信息,包括耗费时间等等:

      几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得  

       关于lpsolve IDE环境的详细使用说明可以参考其使用手册(lpsolve下载主页上可*下载),看到这里也许有人会问,上面的例子很简单,如果对于复杂的网络拓扑,自己手动输入这些表达式显然是不现实的, 那该怎么办,好在lpsolve可以集成在别的开发环境中,它提供了一整套API,可供调用,具体也请参考使用手册,上面提到的博客里对于matlab调用lpsolve有简单的说明,这里补充说一下,mxlpsolve('write_lp',lp,'a.lp')这个语句可以生成IDE环境里可直接执行的脚本文件(C,JAVA等接口也有类似语句),这样复杂的问题可以用高级编程语言建模,然后生成LP脚本文件单独在IDE中运行,是不是很方便。

以上这个Max Flow的例子直接用matlab自带的linprog工具箱求解,写出代码是:

f = [ -1 -1 -1 0 0 0 0 0 0 0 0 0 ];
A = [];
b = [];
Aeq = [1 0 0 -1 -1 0 0 0 0 0 0 0
       0 1 0 0 0 -1 -1 -1 0 0 0 0
       0 0 1 0 0 0 0 0 -1 0 0 0
       0 0 0 1 0 1 0 0 0 -1 0 0
       0 0 0 0 1 0 1 0 1 0 -1 0
       0 0 0 0 0 0 0 1 0 0 0 -1 ];
beq = zeros(6,1);
lb = zeros(12,1);
ub = [3;2;2;5;1;1;3;1;1;4;2;4];
x = linprog(f,A,b,Aeq,beq,lb,ub);
求解结果如下,可以看到虽然最大流的值是正确的,但其它边上的流却与lpsolve求解结果有很大差别:

几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得

接触了yalmip后,可以用yalmip方便的建模,基本LP表达式怎么写,yalmip语句就怎么写,用yalmip为LP问题建模是一个非常方便的过程,至此我认为我的问题用MATLAB + Yalmip + lpsolve就解决了,直到我发现了以下两个问题,让我不得不放弃这一套工具。

为了学习用Yalmip建模,再调用一种合适的solver求解,我找到了湖北工程学院学报上一篇论文《基于Yalmip工具箱的整数规划模型求解方法》,主要看中它贴出了源代码,而且实例又非常简单,作为入门yalmip的材料个人认为很合适,我把它的代码中的solver换成了lpsolve,但没想到求解了一个下午都没有求出来,这个问题只有84个变量,按照lpsolve宣称的性能本不应如此,难道是代码的效率不高?但这是个很简单的程序,按理说不存在代码优化问题。然后我按照源代码的求解器换成了gurobi,结果如论文所述,1s不到就解好了,关于gurobi的安装,使用包括入门学习后文会详述。

到了这一步我打算放弃lpsolve,但是不打算放弃yalmip,就像前面一直强调的用yalmip建模代码书写实在非常方便易学,而且建模和求解可以轻松分开,比如前面的程序,yalmip建好模后,可以调用不同的solver,实在很高效。可是我用一些小程序测试yalmip时出现了一些问题。

比如lpsolve自带的测试程序example1.m,存放位置是../lp_solve_5.5.2.0_MATLAB_exe_win64,自动生成的LP求解式如下:

/* Objective function */
min: +COLONE +3 COLTWO +6.24 COLTHREE +0.1 COLFOUR;

/* Constraints */
THISROW: +78.26 COLTWO +2.9 COLFOUR >= 92.3;
THATROW: +0.24 COLONE +11.31 COLTHREE <= 14.8;
LASTROW: +12.68 COLONE +0.08 COLTHREE +0.9 COLFOUR >= 4;

/* Variable bounds */
COLONE >= 28.6;
18 <= COLFOUR <= 48.98;
            这样一个简单的程序用yalmip建模后在在matlab中居然无法求解,yalmip建模代码如下:

x = sdpvar(4,1);
f = x(1) + 3*x(2) + 6.24*x(3) + 0.1*x(4);
F = [];
F = [F 78.26*x(2) + 2.9*x(4) >= 92.3];
F = [F 0.24*x(1) + 11.31*x(3) <= 14.8];
F = [F 12.68*x(1) + 0.08*x(3) + 0.9*x(4) >= 4];
F = [F x(1) >= 28.6];
F = [F 18 <= x(4) <= 48.98];
ops = sdpsettings('solver','lpsolve');
solvesdp(F,-f,ops);
z = double(f);
x = double(x);

matlab中求解的输出语句如下:

几种优化工具(linprog,lpsolve,yamlip,gurobi)使用心得

无论换成哪一个求解器都无法求解,又经过其他代码的测试发现,yalmip建模后它似乎规定了求解变量必须是整型值,即使用它宣称的实型变量定义sdpvar,有的问题我即使用yalmip + lpsolve/gurobi能求解,但求解值均不优于直接在lpsolve或gurobi中求解结果,因为它隐含规定了每个变量都必须是整型,所以以上这类问题当求解变量值不可能是整型时它就无法得出结果,对于有些问题,变量也可以是实型解,也可以是整型解,则yalmip + lpsolve/gurobi能求解出结果,但当实型解是最优结果时,它的结果就不是最优的。解决这个问题我想需要知道yalmip调用解析器的结果,最好能够调试它具体向解析器输送的用该解析器表示的约束式是什么形式,还有它定义变量究竟有什么规定,需要设置什么参数才可以把变量定义为实型(单纯sdpvar显然是不行的),这两个问题我在网上找了很久都没有答案。

另外我还有一个发现,当我用flow = sdpvar(n,n)定义一个方阵变量,再用F = [F flow >= 0 ];定义对此变量的约束时,Warning: Solver not applicable (lpsolve)这种形式的警告,换用任何解析器都不行,解决方案只能是这样定义变量flow = sdpvar(n,n,'full');或把flow矩阵中的每个变量的约束完整的列出来,例如flow(1,1) > 0; flow(1,2) > 0等等,总之感觉yalmip建模语言看似简单,其实其中还是有很多玄机的,想用它来解决实际问题,还是要花一番功夫去学习,对于想用它解决复杂问题的初学者来说未必是一个好的选择。

同时研究了一下lpsolve及yalmip编码,这两种工具都要求程序员显式地将变量组织成矩阵,对于具体的网络流问题(抽象成数据结构就是图结构),这种要求就会造成写代码时非常的不方便,而且还需要很多冗余的空间(尤其对于网络拓扑是稀疏图的情况),不过这当然也是跟我后来使用的工具比较得出的结论。
       下面隆重介绍一下gurobi以及在python环境中调用gurobi, 首先gurobi的求解能力宣称>=cplex,这个没有测试过,不辨真假,但是从我求解问题的测试结果看,它的求解能力是高于lpsolve的,而且还可以申请到一个免费使用一年的许可证(这点就优于cplex的90天的许可证,而且还有求解规模限制)。

       它的官方网站是http://www.gurobi.com/,在这里可以下到最新版本以及申请许可证,安装很简单,windows下下载对应的版本双击即可安装,安装好后第一次运行时会提示要许可证,这时输入申请的许可证号,如果生成的.lic文件不放在默认的C://users//user_name目录下,则需要设置一下环境变量,gurobi激活程序会提示详细步骤,非常简单。

        将其集成到Matlab安装环境中仅需在Matlab中运行其安装目录下的win64/matlab/gurobi_setup.m.

       我的问题是通过在python环境中调用gurobi得到解决的,下面叙述一下在python环境中集成gurobi的步骤,首先要选对python版本,我的gurobi版本是5.6.3,则对应的python版本是2.7.8(任意一个python2.7版本均可),同时python和gurobi必须同时是32bit或64bit,如果一个是32bit,一个是64bit则安装一定会失败,这两点注意了,一般后面的步骤就不会失败了。然后运行gurobi安装目录下win64/bin/pysetup.bat,没有报错则表明安装成功。在python环境下运行win64/examples/python目录下的任何一个实例代码,如果跑出了结果则表明gurobi环境已经正确集成到python中了

       顺便说一下gurobi安装目录下自带了一些资料非常好,首先是win64/docs下的quickstart.pdf,指导gurobi环境的安装使用及其与其它编程语言的集成方法和实例,refman.pdf介绍了它支持的语言调用它的API接口,而examples.pdf则是各个实例的代码,各个实例可直接运行代码存放于win64/examples目录下,例如解决网络流问题可以参考win64/examples/python目录下的netflow.py,我觉得用python调用gurobi一个很大的好处就是从程序员的角度来说变量可以用m.addVar()随意定义,约束也可以用m.addConstr()随意定义,不像用其他方式解LP问题最终都要程序员显式地把变量或约束条件定义成矩阵形式,在解决特定问题时这种要求往往使代码结构非常复杂且需要很多冗余的存储空间,另外python由于它的语言特点写代码也非常方便,我的问题尝试过用matlab需要几百行,还需要定义一些复杂的存储结构,但python建模过程(不包含打印求解结果)一百行都不到,而且python的字典结构,tuplelist结构等等对于表示网络拓扑图是非常方便的,因此python+gurobi则是我最终采用的方法。

        很简单的一个问题,我用了将近一个月的时间去学习实践,包括对LP,network flow问题的理解以及具体优化工具的使用,虽然时间有些长了,也走了不少弯路,但是所得到的还是值得的。