指派問題使用 python + lp_solve 解決

指派問題乃線性規劃的一種特例，它的特性是不須強調解為 0-1 變數或整數值，但最後算出來它卻一定就是 0 或 1 ，只是這種說法是學理上的，當我們使用程式來計算時，往往因為這些工具在計算過程中使用數值分析的方法，造成解的結果只會接近 0 或是 接近 1 ，而不是純正的 0 或 1。

這也就是我在第 73 行中，使用 v > 0 而不是 v == 1 的原因，如果是寫 v == 1 的話，有些 v 值是 0.999999 的，就不會顯現了。事實上，使用 v > 0.5 會更好。不過，我最後檢查時，發現 > 0 就可以秀出這 50 個 x 的值，也就算了。

lp_solve 函式庫安裝方法請見舊文。

 1 # -*- coding: utf8 -*-
 2 """ 問題： 
 3     
 4     指定 0, 1, 2, ..., 49 等 50 個不可重複的數字給 x0 ~ x49，例如 x0 = 12, x1 = 33, ...
 5 
 6     y = sin(1*x0) + sin(2*x1) + sin(3*x2) + ... + sin(50*x49)
 7 
 8     求解 y 之最大值？
 9 
10     解法：
11 
12     此問題可視為一種指派問題，也就是說有 0 ~ 49 等員工，要放到 x0 ~ x49 的職位去，
13     這樣決策變數就會變成 p00(值為 1 代表 x0=0), p01(值為 1 代表 x1=0),
14     p02 , ..., p49_49 等 2500 個決策變數，且其值必為 0 或 1 。
15 
16     雖然目標函式看起來是非線性的，但其實是線性的， y 函式的係數應該長得如下：
17 
18             x0          x1          x2          ...
19     0       0(C00)      0(C01)      0(C02)      ...
20     1       0.84(C10)   0.91(C11)   0.14(C12)   ...
21     2       0.91(C20)   -0.76(C21)  -0.28(C22)  ...
22     ...     ...         ...         ...         ...
23 
24     所以如果決策變數是 p20 = p01 = p12 = 1，其餘為 0 ，則代表 x0 = 2，x1 = 0，x2 = 1，
25     這樣 y = 0.91 + 0 + 0.14 = 1.05 。
26 
27     所以目標式可以寫成 y = C00 * p00 + C01 * p01 + ... + C49_49 * p49_49 。
28 
29     最後再加上限制式
30     
31     p00 + p01 + ... + p0_49 = 1
32     p10 + p11 + ... + p1_49 = 1
33     ...
34     p49_0 + p49_1 + ... + p49_49 = 1
35 
36     p00 + p10 + ... + p49_0 = 1
37     p01 + p11 + ... + p49_1 = 1
38     ...
39     p0_49 + p1_49 + ... + p49_49 = 1
40 
41     等 100 條限制式後，就能求 y 的最佳解。
42 
43 """
44 from math import sin
45 import lpsolve55 as L
46 
47 LENGTH = 50
48 C = []
49 
50 for i in xrange(LENGTH):
51     for j in xrange(LENGTH):
52         C.append(-1*sin((j+1)*i)) # lp_solve 預設解極小值問題，所以我把目標函數係數全乘以 -1
53 
54 lp = L.lpsolve('make_lp', 0, LENGTH**2)
55 L.lpsolve('set_verbose', lp, L.IMPORTANT)
56 ret = L.lpsolve('set_obj_fn', lp, C)
57 
58 for i in xrange(LENGTH):
59     p = [0,] * (LENGTH ** 2)
60     for j in xrange(i*LENGTH, i*LENGTH+LENGTH): p[j] = 1
61     ret = L.lpsolve('add_constraint', lp, p, L.EQ, 1)
62 
63     p = [0,] * (LENGTH ** 2)
64     for j in xrange(0, LENGTH):
65         p[j*LENGTH+i] = 1
66     ret = L.lpsolve('add_constraint', lp, p, L.EQ, 1)
67 
68 L.lpsolve('solve', lp)
69 print u'目標值： %s' % (L.lpsolve('get_objective', lp) * -1) #要乘以 -1 來還原目標值。
70 vars = L.lpsolve('get_variables', lp)[0]
71 print u'決策變數： %s' % vars
72 for (ij, v) in enumerate(vars):
73     if v > 0:
74         i = ij / LENGTH
75         j = ij % LENGTH
76         print 'x%s = %s, ' % (j, i),
77         if i % 5 + 1 == 5: print

目標值最佳解為 47.8620523191 。

各變數值如下：
x21 = 0, x32 = 1, x47 = 2, x33 = 3, x1 = 4,
x37 = 5, x16 = 6, x45 = 7, x11 = 8, x25 = 9,
x18 = 10, x30 = 11, x7 = 12, x17 = 13, x0 = 14,
x41 = 15, x36 = 16, x22 = 17, x49 = 18, x9 = 19,
x44 = 20, x26 = 21, x43 = 22, x13 = 23, x42 = 24,
x35 = 25, x8 = 26, x20 = 27, x39 = 28, x40 = 29,
x29 = 30, x10 = 31, x34 = 32, x4 = 33, x2 = 34,
x38 = 35, x24 = 36, x6 = 37, x46 = 38, x5 = 39,
x27 = 40, x28 = 41, x14 = 42, x23 = 43, x48 = 44,
x19 = 45, x31 = 46, x12 = 47, x15 = 48, x3 = 49,

=== 後記 ===

經老師指導後，使用

ret = L.lpsolve('set_binary', lp, [1,]*(LENGTH**2)) #大約加在第 59 行後

令決策變數為 0-1 二元變數後，計算時間馬上減少了 60% 。

x的，所有的糟糕事都擠在一起出現

今天早上，心血來潮，把線性規劃套在 51 個城市旅行問題上，演算法應該沒寫錯，但這個題目應該會跑很久。所以，我在程式寫完的那一刻要把它 commit 到我的版本控制器上，但發現網路不通。原來是，學弟們超流了，被計中斷網。

索性直接跑吧! 從上午 11 點一直跑到下午 4 點，才覺得這可能無止盡，想說還是把程式放到家裡的閘道器去 run 吧!反正它本來就不關機的。下面是我下的指令：

# /usr/bin/python2.4 TSP-51_byLP.py > TSP-51_byLP.log
# cp -rf TSP-51_byLP.*

原本在 TSP-51_byLP.* 後面，應該還要跟著一個資料夾的，但我太快按下 [Enter] 了。所以…。

x的，寫了一個早上的程式就變成它 run 出來的 log 。

x的，還是趁現在沒忘掉前，再寫一個吧!

Lp_solve 安裝

lp_solve 是以 C 語言寫成的 Mixed Integer Linear Programming (MILP) solver ，所以基本上，你是用 C 語言的人是可以直接把它的函式庫包在所寫的程式中的。

但我們是用 Python 的人，那要如何能在 python 語言把 lp_solve 載入呢!要作到這個，你首先要知道 lpsolve55.dll (windows)/liblpsolve55.so(unix) 就是 lp_solve 的函式庫，你只要有這個檔案，你就可以使用 C 語言來寫 MILP 應用了。

而用 Python 的人，則是要多安裝一個 lpsolve55.pyd(windows)/lpsolve55.so(unix) ，這個檔是 Python 與 lp_solve 函式庫的溝通者/包裝者，你的 Python 程式都是先載入 lpsolve55.pyd(windows)/lpsolve55.so(unix) ，然後它會幫你處理與 lpsolve55.dll (windows)/liblpsolve55.so(unix) 的資料傳輸。

Ubuntu 7.04 版： 請到 sourceforge.net 下載 5.5.0.10 的檔案。

1. lp_solve_5.5.0.10_dev.tar.gz - 裡面放的是 liblpsolve55.so
2. lp_solve_5.5.0.10_Python_source.tar.gz - 裡面放的是 lpsolve55.so 的程式碼

把 liblpsolve55.so 放到你的 /usr/lib 下(而其他解出來的檔案則拿到 lp_solve_5.5.0.10_Python_source 的資料夾中，與 extra 同一層)。接下來確定你有安裝 python-dev 套件，如果沒有，你應該知道如何安裝它吧!(# sudo apt-get install python-dev)，接下來解開 lp_solve_5.5.0.10_Python_source.tar.gz ，進入到 extra/Python 中，執行 # sudo python setup.py install 。

如果沒有意外出現的話，你可以在 python 直譯器中見到如下訊息：

Python 2.5.1 (r251:54863, May  2 2007, 16:56:35) 
[GCC 4.1.2 (Ubuntu 4.1.2-0ubuntu4)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from lpsolve55 import lpsolve
>>> lpsolve()
lpsolve  Python Interface version 5.5.0.5
using lpsolve version 5.5.0.10

Usage: [ret1, ret2, ...] = lpsolve('functionname', arg1, arg2, ...)

CMCLass: 派送問題(1)

受強烈地震之災，USA--A區有 12 軍事基地 毀損，被迫需搬遷至 B 區 12 軍事基地。搬遷成本與軍事基地間搬遷直線距離成正比。

兩區的平面座標如下：

A區	X 座標	Y 座標		B區	X 座標	Y 座標
1	1130	863		1	1031	1206
2	1705	1283		2	1046	1000
3	1326	1736		3	1803	1809
4	975	825		4	1588	1821
5	1286	807		5	1300	1482
6	909	1143		6	1939	1228
7	1579	1608		7	1147	1703
8	1162	1118		8	1319	1254
9	876	1573		9	1300	1182
10	1198	1748		10	1817	1269
11	801	1134		11	1945	1605
12	956	1226		12	1349	1488

非常典型的 lp 問題，方程式也非常好設定。首先我們要列出這 a[12] 到 b[12] 的所有相對距離，共有 12 x 12 = 144 個，這也表示決策變數 Xij 有144個，當 Xij = 1 時，表示由 i 基地搬遷至 j 基地，當 Xij = 0 時，表示「沒事發生」。

所以目標函數就是所有相對距離 x 決策變數的總和。

接下來設定 25 條方程式：

all Xij >= 0。

sum j=1~12 Xij = 1 對所有的 i 成立。(實在不知道該如何展現數學方程式，只好用口語的方式)

sum i=1~12 Xij = 1 對所有的 j 成立。(實在不知道該如何展現數學方程式，只好用口語的方式)

這樣就有 25 條方程式了。跑一下 lp 程式就可以得到答案。

那麼 Python 如何跑 lp 程式呢!首先請裝 GLPK 函數庫，這是由 gnu 組織開發的 linear programming toolkit ，功能還可以，但如果決策變數太多，有可能會跑不完。

理論上，你裝了 GLPK 後，就可以執行 lp 計算，但是直接用 GLPK 的話，必須用 mathprog 格式餵資料給 GLPK ，而這個格式太複雜了，所以我先去載入一個輸出 mathprog 格式的 Python 類別 pulp ，透過這個類別幫我控制 GLPK 函式庫。

以下是我的 Python 程式：

def caldest(From, To):
----from math import sqrt
----dest = []

----for start in From:
--------subdest = []
--------for end in To:
------------subdest.append(round(sqrt((start[0] - end[0])**2 + (start[1] - end[1])**2), 0))
--------dest.append(subdest)

----return dest

def optima(coef):
----from pulp import *
----(Dest, var, obj) = (coef, [], 0)
----prob = LpProblem("dest", LpMinimize)

----for i in xrange(len(Dest)):
--------var.append([])
--------for j in xrange(len(Dest[0])):
------------(str_i, str_j) = ("%03d" % (i + 1), "%03d" % (j + 1))
------------var[i].append(LpVariable("var"+str_i+'_'+str_j, 0, 1, LpInteger))
------------obj += Dest[i][j] * var[i][j]

----prob += obj, 'OBJ.'

----for i in xrange(len(Dest)):
--------st = 0
--------for j in xrange(len(Dest[0])):
------------st += var[i][j]
--------prob += st == 1

----for j in xrange(len(Dest[0])):
--------st = 0
--------for i in xrange(len(Dest)):
------------st += var[i][j]
--------prob += st == 1

----#prob.writeLP("dest.lp")
----prob.solve()

----print "Status:", LpStatus[prob.status]

----for v in prob.variables():
--------if v.varValue >= 1:
------------print v.name[3:6] + ' => ' + v.name[7:10]

----print "objective=", value(prob.objective)

if __name__ == '__main__':
----A = [[1130, 863],
--------[1705, 1283],
--------[1326, 1736],
--------[975, 825],
--------[1286, 807],
--------[909, 1143],
--------[1579, 1608],
--------[1162, 1118],
--------[876, 1573],
--------[1198, 1748],
--------[801, 1134],
--------[956, 1226]]

----B = [[1031, 1206],
--------[1046, 1000],
--------[1803, 1809],
--------[1588, 1821],
--------[1300, 1482],
--------[1939, 1228],
--------[1147, 1703],
--------[1319, 1254],
--------[1300, 1182],
--------[1817, 1269],
--------[1945, 1605],
--------[1349, 1488]]

----dest = caldest(A, B)
----optima(dest)

執行後，得到搬遷方式為：

001 => 009
002 => 010
003 => 003
004 => 002
005 => 006
006 => 001
007 => 011
008 => 008
009 => 007
010 => 004
011 => 005
012 => 012

最小路徑為 4412.0