備忘:MIP演算法
來源:https://docs.python-mip.com/en/latest/quickstart.html#creating-models
快速開始
本章介紹了使用 Python-MIP 構建和優化模型所需的主要元件。方法及其引數的完整描述可在第 4 章中找到。
在 Python 程式碼中啟用 Python-MIP 的第一步是新增:
from mip import *
載入後,Python-MIP 將顯示其安裝的版本:
Using Python-MIP package version 1.6.2
建立模型
模型類代表優化模型。下面的程式碼使用預設設定建立一個空的混合整數線性規劃問題。
m = Model()
預設情況下,優化意義設定為最小化,選定的求解器設定為 CBC。如果 Gurobi 已安裝並配置,則將使用它。您可以使用建構函式的附加引數來更改模型的客觀意義或強制選擇特定的求解器引擎:
m = Model(sense=MAXIMIZE, solver_name=CBC) # use GRB for Gurobi
建立模型後,您應該包括決策變數、目標函式和約束。這些任務將在下一節中討論。
變數
使用該方法將決策變數新增到模型中add_var()
。沒有引數,一個帶域的變數R+R+被建立並返回它的引用:
x = m.add_var()
通過使用 Python 列表初始化語法,您可以輕鬆建立變數向量。y[0]
……,y[n-1]
並將它們的引用儲存在一個列表中。
n = 10
y = [ m.add_var(var_type=BINARY) for i in range(n) ]
其他變數型別是CONTINUOUS
(預設)和INTEGER
. 可以為變數指定的一些附加屬性是它們的下限和上限(分別為 和 )以及名稱(屬性lb
)。命名變數是可選的,如果您計劃以.LP 或 .MPS 檔案格式儲存模型(請參閱儲存、載入和檢查模型屬性),這將特別有用。以下程式碼建立了一個名為的整數變數ub
name
zCost
{−10,…,10}{−10,…,10}. 請注意,變數的引用儲存在名為z
.
z = m.add_var(name='zCost', var_type=INTEGER, lb=-10, ub=10)
您不需要儲存變數的引用,儘管編寫約束通常更容易這樣做。如果您不儲存這些引用,您可以在之後使用 Model 函式獲取它們var_by_name()
。以下程式碼檢索命名變數的引用zCost
並將其上限設定為 5:
vz = m.var_by_name('zCost')
vz.ub = 5
約束
約束是線性表示式,分別涉及變數、等於或等於、小於或等於和大於或等於的含義,以及一個常數==
。約束<=
>=
x+y≤10x+y≤10可以很容易地包含在模型中m
:
m += x + y <= 10
求和表示式可以用函式來實現xsum()
。如果對於揹包問題nn物品,每件都有重量wiwi,我們想包括一個約束來選擇具有二進位制變數的專案xixi尊重揹包容量cc,則可以使用以下程式碼在模型中包含此約束m
:
m += xsum(w[i]*x[i] for i in range(n)) <= c
在求和中條件包含變數也很容易。假設只有偶數索引項受到容量限制:
m += xsum(w[i]*x[i] for i in range(n) if i%2 == 0) <= c
最後,命名約束可能很有用。這樣做很簡單:線上性表示式之後包含約束的名稱,用逗號分隔。下面給出一個例子:
m += xsum(w[i]*x[i] for i in range(n) if i%2 == 0) <= c, 'even_sum'
與變數一樣,約束的引用可以通過它們的名稱來檢索。模型函式constr_by_name()
對此負責:
constraint = m.constr_by_name('even_sum')
目標函式
預設情況下,使用最小化意義建立模型。以下程式碼將目標函式更改為∑n−1i=0cixi∑i=0n−1cixi通過設定objective
我們示例模型的屬性m
:
m.objective = xsum(c[i]*x[i] for i in range(n))
為了指定目標是最小化還是最大化目標函式,包括了兩個有用的函式:minimize()
和maximize()
。下面是兩個使用示例:
m.objective = minimize(xsum(c[i]*x[i] for i in range(n)))
m.objective = maximize(xsum(c[i]*x[i] for i in range(n)))
您還可以通過將sense
模型屬性設定為MINIMIZE
或來更改優化方向MAXIMIZE
。
儲存、載入和檢查模型屬性
模型方法write()
和read()
可用於分別儲存和載入 MIP 模型。模型支援的檔案格式是LP 檔案格式,它更易讀,更適合除錯,而MPS 檔案格式,推薦用於擴充套件相容性,因為它是一種較舊且更廣泛採用的格式。呼叫該write()
方法時,副檔名(.lp 或 .mps)用於定義檔案格式。因此,要將m
使用 lp 檔案格式的模型儲存到檔案 model.lp 我們可以使用:
m.write('model.lp')
同樣,我們可以讀取模型,從而從讀取的 LP 或 MPS 檔案中建立變數和約束。讀取模型後,其所有屬性都可用,例如約束矩陣中的變數、約束和非零的數量:
m.read('model.lp')
print('model has {} vars, {} constraints and {} nzs'.format(m.num_cols, m.num_rows, m.num_nz))
優化和查詢優化結果
MIP 求解器執行分支切割 (BC) 演算法,該演算法將在有限時間內提供最佳解決方案。在許多情況下,這個時間可能對於您的需要來說太大了。幸運的是,即使完整的樹搜尋過於昂貴,通常在搜尋開始時就可以獲得結果。有時在處理第一個樹節點時會產生一個可行的解決方案,並花費大量額外的精力來改進對偶邊界,這是對最優解決方案成本的有效估計。當這個估計(最小化的下限)與找到的最佳解決方案(上限)的成本完全匹配時,搜尋就結束了。對於實際應用,通常執行截斷搜尋。這optimize()
執行配方優化的方法,接受可選的處理限制作為引數。以下程式碼執行分支和切割演算法以求解模型m
長達 300 秒。
1 2 3 4 5 6 7 8 9 10 11 12 13 |
m.max_gap = 0.05
status = m.optimize(max_seconds=300)
if status == OptimizationStatus.OPTIMAL:
print('optimal solution cost {} found'.format(m.objective_value))
elif status == OptimizationStatus.FEASIBLE:
print('sol.cost {} found, best possible: {}'.format(m.objective_value, m.objective_bound))
elif status == OptimizationStatus.NO_SOLUTION_FOUND:
print('no feasible solution found, lower bound is: {}'.format(m.objective_bound))
if status == OptimizationStatus.OPTIMAL or status == OptimizationStatus.FEASIBLE:
print('solution:')
for v in m.vars:
if abs(v.x) > 1e-6: # only printing non-zeros
print('{} : {}'.format(v.name, v.x))
|
可以使用額外的處理限制:max_nodes
限制搜尋樹中探索節點的最大數量,並max_solutions
在獲得多個可行解後完成BC演算法。指定結束搜尋的邊界應該有多緊也是明智的。模型屬性max_mip_gap
指定結束搜尋的上限與下限的允許百分比偏差。在我們的示例中,只要上下界的距離小於或等於 5%(參見第 1 行),就可以完成搜尋。
該optimize()
方法返回OptimizationStatus
BC搜尋的狀態( ): OPTIMAL
如果搜尋結束並找到最優解;FEASIBLE
如果找到了可行的解決方案,但沒有時間證明該解決方案是否最優; NO_SOLUTION_FOUND
如果在截斷搜尋中沒有找到解決方案;INFEASIBLE
或者INT_INFEASIBLE
如果模型不存在可行的解決方案; UNBOUNDED
如果缺少約束或ERROR
如果在優化過程中發生了一些錯誤。在上面的示例中,如果可行的解決方案可用(第 8 行),則列印值不為零的變數。還要注意,即使沒有可行的解決方案可用,雙重界限(最小化情況下的下限)也可用(第 8 行):如果執行了截斷執行,即求解器由於時間限制而停止,您可以檢查這個雙重界限是對檢查gap
屬性的解決方案質量的估計。
在樹搜尋過程中,通常會找到許多不同的可行解。求解器引擎將此解決方案儲存在解決方案池中。以下程式碼列印優化旅行商問題時找到的所有路線。
for k in range(model.num_solutions):
print('route {} with length {}'.format(k, model.objective_values[k]))
for (i, j) in product(range(n), range(n)):
if x[i][j].xi(k) >= 0.98:
print('\tarc ({},{})'.format(i,j))
效能調優
MIP 求解器的樹搜尋演算法提供了一組改進的可行解和下界。根據您的應用程式,您將對快速生成可行的解決方案更感興趣,而不是對可能需要昂貴計算的改進下限更感興趣,即使從長遠來看,這些計算證明值得證明找到的解決方案的最優性。模型屬性 emphasis
提供了三種不同的設定:
- 預設設定:
-
試圖在尋找改進的可行解決方案和改進的下限之間取得平衡;
- 可行性:
-
專注於在搜尋過程的最初時刻尋找改進的可行解決方案,啟用啟發式;
- 最優性:
-
啟用產生改進下界的過程,即使第一個可行解決方案的產生被延遲,也專注於修剪搜尋樹。
將此設定更改為 1 或 2 會觸發在搜尋樹的每個節點處處理的幾種演算法的啟用/停用,這些演算法會影響求解器的效能。即使平均而言,這些設定會如前所述改變求解器的效能,但根據您的公式,這些更改的影響可能會非常不同,通常值得在您的應用程式中使用這些不同的設定檢查求解器的行為。
另一個可能值得調整的引數是cuts
屬性,它控制在生成切割平面時應該花費多少計算工作量。