最近建立了一個網絡流模型,是一個混合整數線性規划問題(模型中既有連續變量,又有整型變量)。當要求解此模型的時候,發現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中還未安裝):
>> yalmiptest
+++++++++++++++++++++++++++++++++++++++++++++++
| Searching for installed solvers |
+++++++++++++++++++++++++++++++++++++++++++++++
| Solver| Version/module| Status|
+++++++++++++++++++++++++++++++++++++++++++++++
| LPSOLVE| MXLPSOLVE| found|
| CPLEX| IBM| found|
| CPLEX| IBM| found|
| CPLEX| IBM| found|
| LINPROG| | found|
| QUADPROG| | found|
| LMILAB| | found|
| FMINCON| geometric| found|
| FMINCON| standard| found|
| FMINSEARCH| | found|
| BNB| | found|
| BINTPROG| | found|
| CUTSDP| | found|
| BMIBNB| | found|
| KKTQP| | found|
| NONE| | found|
| GUROBI| MEX| not found|
| CPLEX| CPLEXINT| not found|
| GLPK| GLPKMEX-CC| not found|
| GLPK| GLPKMEX| not found|
| CDD| CDDMEX| not found|
| NAG| e04mbf| not found|
| NAG| e04naf| not found|
| CLP| CLPMEX-LP| not found|
| XPRESS| MEXPRESS 1.1| not found|
| XPRESS| MEXPRESS 1.0| not found|
| XPRESS| FICO| not found|
| XPRESS| FICO| not found|
| QSOPT| MEXQSOPT| not found|
| OSL| OSLPROG| not found|
| MOSEK| LP/QP| not found|
| MOSEK| SOCP| not found|
| MOSEK| GEOMETRIC| not found|
| CPLEX| CPLEXMEX| not found|
| BPMPD| | not found|
| CLP| CLPMEX-QP| not found|
| OOQP| | not found|
| QPIP| | not found|
| QPAS| | not found|
| LINDO| MIQP| not found|
| SEDUMI| 1.1| not found|
| SEDUMI| 1.3| not found|
| SEDUMI| 1.05| not found|
| SEDUMI| 1.03| not found|
| SDPT3| 4| not found|
| SDPNAL| 0.1| not found|
| LOGDETPPA| 0.1| not found|
| SPARSECOLO| 0| not found|
| SDPT3| 3.1| not found|
| SDPT3| 3.02| not found|
| SDPT3| 3.0| not found|
| SDPA| M| not found|
| DSDP| 5| not found|
| DSDP| 4| not found|
| SDPLR| | not found|
| CSDP| | not found|
| MAXDET| | not found|
| PENSDP| PENOPT| not found|
| PENSDP| TOMLAB| not found|
| PENBMI| PENOPT| not found|
| PENBMI| TOMLAB| not found|
| SDPNAL| | not found|
| LMIRANK| | not found|
| VSDP| 0.1| not found|
| MPT| | not found|
| MPLCP| | not found|
| KYPD| | not found|
| STRUL| 1| not found|
| PENNON| standard| not found|
| SNOPT| geometric| not found|
| SNOPT| standard| not found|
| LINDO| NLP| not found|
| IPOPT| standard| not found|
| IPOPT| geometric| not found|
| GPPOSY| | not found|
| SPARSEPOP| | not found|
| POWERSOLVER| | not found|
+++++++++++++++++++++++++++++++++++++++++++++++
有了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 |
開放網絡 |
網絡開放獲取,建模語言簡單、統一,自帶基本求解算法,並支持集成幾乎所有其他求解器。推薦使用。 |
暫無。 |