Cplex混合整數規劃求解(Python API)

碼頭牛牛發表於2023-10-01
PythonAPI

絕對的原創！罕見的Cplex-Python API混合整數規劃求解教程！這是我盯了一天的程式一條條寫註釋悟出來的•́‸ก

一、問題描述

求解有容量限制的的設施位置問題，使用Benders分解。模型如下：

\[min\quad\sum^{locations}_{j=1}fixedCost_j//open_j+\sum^{locations}_{j=1}\sum^{clients}_{i=1}cost_{ij}×supply_{ij} \]

\(s.t.\)

\[\sum^{locations}_{j=1}supply_{ij}=1\quad\quad\forall{i\in{clients}} \]

\[\sum^{clients}_{i=1}supply_{ij}\leq{capacity}_{ij}×open_{j}\quad\quad\forall{j\in{locations}} \]

\[0\leq{supply}_{ij}\leq{1}\quad\quad\forall{i\in{clients}};\forall{j\in{locations}} \]

\[open_i=0或1\quad\quad\forall{i\in{clients}} \]

二、程式

1. sys.exit函式：優雅地退出程式

摘自博文：Python中的sys.exit函式：優雅地退出程式_Python 筆記_設計學院

在Python程式設計中，程式在執行過程中可能會遇到需要停止程式的情況，如果不加處理，程式執行到中途就被強制停止的話，可能會導致資料丟失，甚至可能會讓程式異常崩潰。因此，對於Python程式退出的處理，我們可以使用Python的內建函式sys.exit()，進行優雅地退出程式。

（1）sys.exit函式的基本用法

  • Python內建函式\(sys.exit()\)，用於退出程式。如果該函式被呼叫時不帶任何引數，那麼Python直譯器將會以狀態碼0來退出程式，表示程式執行成功，並返回控制檯。下面是\(sys.exit()\)基本用法的示例程式碼：
import sys
sys.exit()
  • 如果該函式被呼叫時帶有整數引數n，那麼Python直譯器將會以狀態碼n來退出程式，並返回控制檯。下面是以狀態碼1退出程式的示例程式碼：
import sys
sys.exit(1)

注意，當狀態碼不等於0時，表示程式執行發生某種異常或錯誤，需要進一步處理。

（2）sys.exit程式終止時的清理工作

在程式結束之前，我們可能需要完成一些清理工作，比如關閉一些檔案、釋放一些記憶體等。如果程式直接使用強制停止的方式結束，就有可能使這些工作被忽略或未能完全完成。此時，我們可以利用try/finally語句來完善這部分清理工作：

import sys
import time

try:
    # code block here
    time.sleep(5)
finally:
    # closing file or releasing resource
    print('clean up resources')
    sys.exit()

2. 程式實現

（1）函式簡介

\(cpx = cplex.Cplex()\)

  • \(cpx.variables.add(obj,lb,ub,type)\)

    • 簡介：cpx.cariables.add()是用於向模型中新增變數的方法，在新版的Cplex中，為了節約記憶體空間，通常以range的格式儲存資料。所以為了方便檢視，一般還要將其轉換為list。

    • 引數詳解：

      • obj(list): 變數的目標函式係數列表；

      • lb(list): 變數的下界列表。預設為0；

      • ub(list): 變數的上界列表。預設為正無窮。

      • types(list): 變數的型別。可以是以下值：'C': 連續變數;'I': 整數變數;'B': 二進位制變數。預設值為 'C'。

    • 注意：引數只支援一維列表的輸入，如果是二維及以上的列表需要使用for迴圈進行輸入

  • \(cplex.SparsePair(ind,val)\)

    • 簡介：cplex.SparsePair()用於表示稀疏的線性表示式。它主要用於定義線性約束和目標函式，特別是在涉及大量變數且多數係數為0的情況下。使用稀疏表示可以顯著提高效率和節省記憶體。

    • 引數詳解：

      • ind(list): 一個包含變數索引的列表。這些索引指向模型中的特定變數（也就是用於輸入決策變數）。這個列表為int型資料列表時，代表的是決策變數對應索引列表；為字串型資料時，代表的是決策變數。如\(x_1+2x_3=0\)，ind=["x1","x3"]或ind=[0,2]

      • val(list): 與 ind 中的變數索引相對應的係數列表（與決策變數一一對應的係數）

    • cplex.SparsePair()有每個引數之間相加的意思

    • 注意：ind和val引數只支援一維列表的輸入

  • \(cpx.linear\_constraints.add(lin\_expr=[cplex.SparsePair(ind,val)],senses,rhs)\)

    • 簡介：方法用於向模型中新增線性約束

    • 引數詳解：

      • lin_expr: 這是線性表示式的列表，表示線性約束的左側。每個線性表示式由變數的索引和相應的係數組成，通常使用cplex.SparsePair來表示。
      • senses: 這是一個字元列表，表示每個約束的符號（'L'表示“<=”，'E'表示“=”，'G'表示“>=”）
      • rhs: 這是一個數字列表，表示每個約束的右側值。

    • 注意：上面所有引數只支援一維列表的輸入

    舉例：

    \[x+y=5 \]

    \[2x-y\leq{10} \]

  cpx.linear_constraints.add(
      lin_expr=[[["x", "y"], [1, 1]], [["x", "y"], [2, -1]]],
      senses=["E", "L"],
      rhs=[5, 10]
  )

  • \(cpx.long\_annotations.add(name, defval)\)

    • 簡介：cpx.long_annotations.add()一般用於新增長註解(long annotations)。長註解是Cplex用來為決策變數、約束等模型元素新增後設資料或“註解”的機制。這些註解可能會影響求解器如何解決模型，尤其是在高階策略和方法中。

    • 引數詳解：

      • name (string): 註解的名稱。Cplex預定義了一些註解名稱，如cpx.long_annotations.benders_annotation，用於Benders分解策略。

      • defval (int): 註解的預設值。這是當註解沒有明確為某個模型元素設定值時使用的值。

    • 返回值：這個函式返回新新增註解的索引

  • \(cpx.solve()\)

    • 簡介：用於求解在 CPLEX 物件中定義的最佳化問題。

    • cpx.solve() 函式會啟動相應的演演算法來求解該模型，具體有以下操作：

      • 選擇適當的演演算法。

      • 執行求解過程。

      • 儲存結果。

  • \(cpx.solution.get\_status\_string()\)

    • 簡介：使用cpx.solution.get_status_string()可獲得關於這個狀態的描述性字串。這個函式特別有用，因為它可以讓你快速瞭解模型解的狀態，並據此採取相應的決策。

    • 返回值：

      • "optimal": 表示找到了最優解。

      • "infeasible": 表示問題是不可行的。

      • "unbounded": 表示問題是無界的。

      • "feasible": 表示找到了一個可行解，但不一定是最優的。

      • "integer optimal solution": 表示在整數線性規劃或混合整數線性規劃問題中，已經找到了最優的整數解。

      • ... （還有其他可能的狀態）

  • \(cpx.solution.get\_objective\_value()\)

    • 簡介：用於從 CPLEX 求解器的當前解中獲取目標函式的值。

    • 返回值：目標函式的最優值

  • \(cpx.parameters.mip.tolerances.integrality.get()\)

    • 簡介：用於獲取整數容差 (integrality tolerance) 引數的當前值，返回的是Cplex預設的整數容差 (float資料)

    • 整數容差：整數容差定義了一個變數距離其最近的整數值可以有多遠，而仍然被認為是整數。例如，如果整數容差設定為0.1，那麼一個值為0.9或1.1的變數仍然會被認為滿足整數約束。但如果整數容差設定得更小，例如0.001，那麼這兩個值就不會被認為滿足整數約束。這個引數的目的是為了處理數值誤差和確保求解器在數值上是穩健的。

    • 自定義整數容差：如設定為0.001：cpx.parameters.mip.tolerances.integrality.set(0.001)

  • \(cpx.solution.get\_values()\)

    • 簡介：用於從當前解決方案中檢索變數 (決策變數) 的值。

    • 返回值：此函式返回一個列表，其中每個元素對應於模型中一個決策變數的值。

    • 輸出：

      • 輸出所有決策變數的值：values = cpx.solution.get_values()，如\(values=[1,2,3]\)，那麼代表\(x_1=1;x_2=2;x_3=3\)

      • 指定決策變數輸出：如果只想輸出\(x_1\)和\(x_3\)的值，那麼根據其索引，可以寫成：values = cpx.solution.get_values([0, 2])

（2）程式碼(Cplex官方程式碼修改的)

import sys
import cplex

# 用於求解模型的Benders分解型別
# 定義Benders分解的三種型別：無Benders分解、自動Benders分解和帶註釋的Benders分解
NO_BENDERS = 1
AUTO_BENDERS = 2
ANNO_BENDERS = 3

# 輸出如何使用該指令碼的說明，並退出程式。
def usage():
    print("""\
    Usage: facility.py [options] [inputfile]
     where
       inputfile describes a capacitated facility location instance as in
       ../../../../examples/data/facility.dat. If no input file
       is specified read the file in example/data directory.
       Options are:
       -a solve problem with Benders letting CPLEX do the decomposition
       -b solve problem with Benders specifying a decomposition
       -d solve problem without using decomposition (default)
     Exiting...
    """)
    sys.exit(2)

# 解決有容量限制的設施位置問題(模型輸入部分)
def facility(bendersopt):
    """輸入引數（已知量）"""
    fixedcost=[ 480, 200, 320, 340, 300]
    cost=[[24, 74, 31, 51, 84],
          [57, 54, 86, 61, 68],
          [57, 67, 29, 91, 71],
          [54, 54, 65, 82, 94],
          [98, 81, 16, 61, 27],
          [13, 92, 34, 94, 87],
          [54, 72, 41, 12, 78],
          [54, 64, 65, 89, 89]]
    capacity=[3, 1, 2, 4, 1]

    num_locations = len(fixedcost)  #計算區域數量locations
    num_clients = len(cost)  #計算顧客數量clients

    # 建立Cplex模型例項
    cpx = cplex.Cplex()

    """輸入目標函式"""
    #下面為輸入決策變數為Open部分的目標函式
    #obj：變數的目標函式係數列表；lb:變數的下界列表。預設為0；ub:變數的上界列表。預設為正無窮。
    #types:變數的型別。可以是以下值：'C': 連續變數;'I': 整數變數;'B': 二進位制變數。預設值為 'C'。
    #lb是元素為0，1×num_location的列表；ub是元素為1，1×num_location的列表；這兩個對應的是決策變數open的上限和下限約束
    open_ = list(cpx.variables.add(obj=fixedcost,
                                   lb=[0] * num_locations,
                                   ub=[1] * num_locations,
                                   types=["B"] * num_locations))
    # 在Cplex中，當你向模型中新增變數或約束時，它會返回一個表示新新增物件的索引範圍的序列，這個序列用list來儲存
    # 比如上面的open_,代表的是目標函式係數由 fixedcost 定義，且其值範圍在 0 和 1 之間，最後會生成open_=[0,1,2,3,4](就是5個open變數的索引)

    #輸入決策變數為supply部分的目標函式
    supply = [None] * num_clients  #初始化supply[i]的部分，用於代表顧客i的索引
    for i in range(num_clients):
        # 目標：最小化使用某一位置的固定成本，以及從特定位置為客戶提供服務的成本。
        # 因為cost是一個二維list，所以需要透過for去一維維輸入，這裡的type預設為‘C’(連續變數)
        supply[i] = list(cpx.variables.add(obj=cost[i],
                                           lb=[0.0] * num_locations,
                                           ub=[1.0] * num_locations))
        # 經過迴圈之後，cpx.variable.add()就會在open_的基礎上，繼續往後面加變數，所以supply索引從5開始，supply=[[5, 6, 7, 8, 9], [10, 11, 12, 13, 14],...]

    """輸入約束條件"""
    # 定義每個客戶必須被分配的約束
    #   sum(j in nbLocations) supply[i][j] == 1  for each i in nbClients
    for i in range(num_clients):
        cpx.linear_constraints.add(
            lin_expr=[cplex.SparsePair(
                ind=supply[i], val=[1.0] * num_locations)],
            senses=["E"],
            rhs=[1.0])

    # 定義每個位置的容量必須被遵守的約束
    #   sum(i in nbClients) supply[i][j] <= capacity[j] * open_[j]
    for j in range(num_locations):
        #ind: 先把supply中的每一列取出來組成一個列表，再在列表後面append(open_[j])
        ind = [supply[i][j] for i in range(num_clients)] + [open_[j]]
        val = [1.0] * num_clients + [-capacity[j]]
        cpx.linear_constraints.add(
            lin_expr=[cplex.SparsePair(ind=ind, val=val)],
            senses=["L"],
            rhs=[0.0])

    """設定Benders分解（如果需要的話）"""
    # 這個程式碼沒有使用Benders分解

    # 帶註釋的Benders分解
    if bendersopt == ANNO_BENDERS:
        # 我們透過指定結構來執行Benders分解，
        # 透過使用註解告訴CPLEX哪些變數在主問題中。
        # 預設情況下，變數被分配值CPX_BENDERS_MASTERVALUE+1，因此進入工作區。
        # 變數 open_[j] 應該進入主問題，所以
        # 我們為它們分配值 CPX_BENDERS_MASTER_VALUE。
        mastervalue = cpx.long_annotations.benders_mastervalue
        idx = cpx.long_annotations.add(
            name=cpx.long_annotations.benders_annotation,
            defval=mastervalue + 1)
        objtype = cpx.long_annotations.object_type.variable
        cpx.long_annotations.set_values(idx, objtype,
                                        [(open_[x], mastervalue)
                                         for x in range(num_locations)])
        print("Solving with explicit Benders decomposition.")
    # 自動Benders分解
    elif bendersopt == AUTO_BENDERS:
        # 讓CPLEX自動分解問題。在有容量的設施位置問題中，
        # 主問題的變數應該是整數變數。透過將Benders策略引數設定為Full，
        # CPLEX會將所有整數變數放入主問題，將所有連續變數放入一個子問題，
        # 並進一步分解那個子問題（如果可能的話）。
        cpx.parameters.benders.strategy.set(
            cpx.parameters.benders.strategy.values.full)
        print("Solving with automatic Benders decomposition.")
    # 無Benders分解
    elif bendersopt == NO_BENDERS:
        print("Solving without Benders decomposition.")
    # 否則直接報錯
    else:
        raise ValueError("invalid bendersopt argument")

    """解決模型並顯示解決方案"""
    cpx.solve()  #求解最佳化問題
    print("Solution status =", cpx.solution.get_status_string())  #返回解的狀態
    print("Optimal value:", cpx.solution.get_objective_value())  #返回求得的目標函式值
    tol = cpx.parameters.mip.tolerances.integrality.get()  #返回整數容差引數的當前值
    values = cpx.solution.get_values()  #輸出決策變數的值

    # 遍歷所有可能的設施位置
    for j in [x for x in range(num_locations) if values[open_[x]] >= 1.0 - tol]:
        # 首先，open[j]為0-1變數，要認定open[j]=1，需滿足open[j]∈[1-tol,1+tol]，所以這裡用了values[open_[x]] >= 1.0 - tol]
        # 這裡先找出open[j]=1的量，然後進行迴圈遍歷

        # 同樣的，這裡先找出supply[x][j]為1的數值，並生成針對該區域j的服務顧客編號列表
        # 佔位符{0}表示的是設施地點編號j,佔位符{1}表示的是對應j的服務顧客列表
        print("Facility {0} is open, it serves clients {1}".format(
            j, " ".join([str(x) for x in range(num_clients)
                         if values[supply[x][j]] >= 1.0 - tol])))

def main():
    """處理命令列引數。"""
    # filename = "../../../examples/data/facility.dat"   # 預設的資料檔案路徑

    # 解析命令列引數，設定Benders分解選項或資料檔案路徑
    # 初始化Benders分解選項為NO_BENDERS
    benders = NO_BENDERS

    # 遍歷命令列引數
    for arg in sys.argv[1:]:
        # 判斷引數是否為選項（開始於"-"）
        if arg.startswith("-"):
            # 根據選項設定Benders分解的方式
            if arg == "-a":
                benders = AUTO_BENDERS
            elif arg == "-b":
                benders = ANNO_BENDERS
            elif arg == "-d":
                benders = NO_BENDERS
            else:  # 如果是未知選項則呼叫usage函式顯示幫助資訊
                usage()
        else:  # 如果引數不是選項，則認為它是檔案路徑
            filename = arg

    # 呼叫facility函式來解決問題
    facility(bendersopt=benders)


if __name__ == "__main__":
    main()

相關文章