PageRank



轉移公告

計劃把 http://blog.hoamon.info/ 文章全部轉移至 http://www.hoamon.info/blog/ 這裡,而本 Blogger 站台的文章近 500 篇,我預計在 2014-12-31 前移轉完畢,完成後 http://blog.hoamon.info/ 將只作代轉服務,一律把舊連結如 http://blog.hoamon.info/index.html 轉成 http://www.hoamon.info/blog/index.html ,敬請舊雨新知互相走告。

新文章只發佈在 http://www.hoamon.info/blog/

何岳峰 敬上

2007年3月15日 星期四

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

2 則留言:

  1. Via your blog, I found out the existence of pulp. It may be easier than swiging glpk, which I have been struggling with for days.

    Thanks!

    回覆刪除
  2. but after days, i found the lp_solve lib is better. you can also use python to call the lp_solve. lp_solve 安裝

    回覆刪除

Related Posts Plugin for WordPress, Blogger...