四月 27, 2011

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 指派問題使用 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% 。

三月 17, 2008

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 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的,還是趁現在沒忘掉前,再寫一個吧!

十月 9, 2007

hoamon's sandbox
hoamon
hoamon's sandbox is about »

tag cloud

» 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 下。確定你有安裝 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, ...)

biggo.com.tw

A Django site.