ortools系列-整數規劃

  • 1. 整數規劃
  • 2. 混合整數規劃CBC
  • 3. 混合整數規劃CP-SAT
  • 4. 參考

1. 整數規劃

有一些問題，要求一些但不是全部變數約束為整數，這類問題可以用混合整數規劃(MIP)來解決，混合整數規劃也稱為混合整數線性規劃(MILP)。這裡有幾個例子:

  • 銷售各種消費品的公司需要決定每個月生產多少，以實現利潤最大化。在這類問題中，變數表示離散項的數量，如汽車或電視機，這些項必須是整數。
  • 包裹遞送公司需要將卡車分配到不同的路線，以滿足他們的遞送時間表，同時最小化成本。有時，解決此類問題的最佳方法是讓變數表示要將哪些卡車分配到哪些路線的二進位(0-1)決策。如果給定的卡車被分配到特定的路線，變數為1，否則為0。由於變數只能取0或1的值，這也是一個整數優化問題。（特別地，這是一個布爾優化問題，OR-Tools有專門的技術來解決這個問題。）

OR-Tools提供了幾種方法求解此類問題：

  • 混合整數規劃求解器(MIP)
  • 約束規劃求解器(CP)
  • 最小成本流求解器

其中前兩個是一般性的求解器，能解決任何整數規劃問題，最後一個最小成本流求解器用於求解那些能建模成網路流的問題，對於這類特殊問題，最小成本流求解器比一般的整數規劃和約束規劃要快的多。

2. 混合整數規劃CBC

ortools提供一個混合整數規劃的介面用於調用不同的求解器，ortools內置了Coin-or branch and cut (CBC)，當然你也可以使用第三方的求解器，比如SCIP,GLPK,Gurobi。CBC是一個開源的混合整數規劃求解器，用C++開發。

我們看下面這個問題：

max x + 10y

egin{array}{l} x + 7y le 17.5 \ x le 3.5 \ x ge 0 \ y ge 0 \ x,y {mathop{ m integers} olimits} end{array}

問題的解空間如下：

我們看ortools是如何解決這個問題的。

from ortools.linear_solver import pywraplp

def main():

# 首先，調用CBC求解器
# 還記得我們在LP問題中調用的求解器嗎，
# 是 pywraplp.Solver.GLOP_LINEAR_PROGRAMMING
solver = pywraplp.Solver(SolveIntegerProblem,
pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# 和LP問題一樣，定義變數 x,y 並指定其定義域
# 這裡需要注意，我們定義的是Int類型的變數
x = solver.IntVar(0.0, solver.infinity(), x)
y = solver.IntVar(0.0, solver.infinity(), y)

# 添加約束，方法和LP一樣，通過制定係數來添加約束
# 如果能像cp-sat那麼方便就好了
# x + 7 * y <= 17.5
constraint1 = solver.Constraint(-solver.infinity(), 17.5)
constraint1.SetCoefficient(x, 1)
constraint1.SetCoefficient(y, 7)

# 添加約束
# x <= 3.5
constraint2 = solver.Constraint(-solver.infinity(), 3.5)
constraint2.SetCoefficient(x, 1)
constraint2.SetCoefficient(y, 0)

# 定義目標函數
# Maximize x + 10 * y.
objective = solver.Objective()
objective.SetCoefficient(x, 1)
objective.SetCoefficient(y, 10)
objective.SetMaximization()

# 求解問題，並列印結果，後面的代碼就很簡單了
result_status = solver.Solve()
# The problem has an optimal solution.
assert result_status == pywraplp.Solver.OPTIMAL

# The solution looks legit (when using solvers other than
# GLOP_LINEAR_PROGRAMMING, verifying the solution is highly recommended!).
assert solver.VerifySolution(1e-7, True)

print(Number of variables =, solver.NumVariables())
print(Number of constraints =, solver.NumConstraints())

# The objective value of the solution.
print(Optimal objective value = %d % solver.Objective().Value())
print()
# The value of each variable in the solution.
variable_list = [x, y]

for variable in variable_list:
print(%s = %d % (variable.name(), variable.solution_value()))

if __name__ == __main__:
main()

# 結果
Number of variables = 2
Number of constraints = 2
Optimal objective value = 23

x = 3
y = 2

通過代碼可以看到，其實和LP問題基本一樣，思路是很清晰，就是囉嗦呢。

3. 混合整數規劃CP-SAT

前面我們講了CBC解MIP的問題。約束規劃是一種有別於經典優化理論的優化方法。CP的基礎是可行性(尋找問題的可行解)，側重於約束和變數，而不是目標函數。對於許多類型的問題，CP可以比MIP求解器更快地找到最優解。

那我們應該如何選擇CP還是MIP呢。

對於一個可行點必須滿足所有約束的標準整數規劃問題，MIP求解器速度更快。在這種情況下，可行集是凸的:對於集合中的任意兩點，連接它們的線段完全位於集合中。

另一方面，對於高度非凸可行集的問題，CP-SAT求解器通常比MIP求解器快。當可行集由許多由or連接的約束定義時，就會出現這樣的問題，or意味著一個點只需要滿足其中一個約束就可以是可行的。

為了提高計算速度，CP-sat求解器和原始CP求解器都對整數進行運算。要解決某些約束具有非整數項的問題，必須首先將這些約束轉換為一個足夠大的整數。下面的示例說明瞭這一點。

max 2x + 2y + 3z \ s.t. \ egin{array}{l} x + frac{7}{2}y + frac{3}{2}{ m{z}} le 25 \ 3x - 5y + 7z le 45 \ 5x + 2y - 6z le 37 \ x,y,z ge 0 \ x,y,z {mathop{ m int}} egers end{array}

我們來看代碼：

from ortools.sat.python import cp_model

def main():

# 首先是創建模型
model = cp_model.CpModel()

# 定義變數，指定變數的範圍
# 根據題目可以知道，x,y,z 的最大值肯定不超過50
# var_upper_bound 相當於減少了搜索空間範圍
# 注意，我們這裡定義的是整數變數
var_upper_bound = max(50, 45, 37)
x = model.NewIntVar(0, var_upper_bound, x)
y = model.NewIntVar(0, var_upper_bound, y)
z = model.NewIntVar(0, var_upper_bound, z)

# 添加約束
# 注意到，我們第一個約束條件是 x + 7?2 y + 3?2 z ≤ 25
# 但是 CP-SAT 只能處理整數問題，所以需要對第一個約束處理，使得係數全部是整數
# 方法很簡單，左邊右邊都乘以2就可以了
model.Add(2 * x + 7 * y + 3 * z <= 50)
# 添加第二個約束
model.Add(3 * x - 5 * y + 7 * z <= 45)
# 添加第三個約束
model.Add(5 * x + 2 * y - 6 * z <= 37)

# 定義目標函數
model.Maximize(2 * x + 2 * y + 3 * z)

# 求解並列印結果
solver = cp_model.CpSolver()
status = solver.Solve(model)

if status == cp_model.OPTIMAL:
print(Maximum of objective function: %i % solver.ObjectiveValue())
print()
print(x value: , solver.Value(x))
print(y value: , solver.Value(y))
print(z value: , solver.Value(z))

if __name__ == __main__:
main()

# 結果
Maximum of objective function: 35

x value: 7
y value: 3
z value: 5

看起來也不複雜呢，因為我們的例子也不複雜啊，現實生活中的大部分問題都是大規模整數規劃問題，說不定跑半天都沒出結構。

ortools的文檔中還是一個例子，是用 original-cp 來求解上面的問題，和之前一樣的套路，感興趣的同學可以看看。

4. 參考

  • ortools-Original-CP-Solver
  • Coin-or-branch-and-cut

大家看完記得關注點贊.

master蘇.


推薦閱讀：
相關文章