對偶理論和對偶單純形法——Python實現

郝hai發表於2024-06-01
Python

對偶單純形法是從對偶可行性逐步搜尋出原始問題最優解的方法。由線性規劃問題的對偶理論，原始問題的檢驗數對應於對偶問題的一組基本可行解或最優解；原始問題的一組基本可行解或最優解對應於對偶問題的檢驗數；原始問題約束方程的係數矩陣的轉置是對偶問題約束條件方程的係數矩陣。所以，在求解常數項小於零的線性規劃問題時，可以把原始問題的常數項視為對偶問題的檢驗數，原始問題的檢驗數視為對偶問題的常數項。

對偶理論和對偶單純形法——Python實現

一、對偶理論

若原問題為 \(LP_1\)

\[\max z=C X\\ s.t. \left\{\begin{array}{l}A X \leq b \\ X \geq b \end{array}\right. \]

則其對偶問題定義為 \(LP_2\)

\[\begin{aligned} & \qquad \min w=Y^T b \\ & \text { s.t. }\left\{\begin{array}{l} A^T Y \geq C^T \\ Y \geq 0 \end{array}\right. \end{aligned} \]

  • 弱對偶性: 如果 \(X\) 是原問題的可行解， \(Y\) 是對偶問題的可行解，則 \(C X \leq Y^T b\)證明: \(Y^T \geq Y^T A X=\left(Y^T A X\right)^T=X^T A^T Y \geq X^T C^T=C X\)推論:
  • 最優性: 如果 \(X\) 是原問題的可行解， \(Y\) 是對偶問題的可行解, 且 \(C X=Y^T b\), 則 \(X\) 和 \(Y\) 分別為原問題和對偶問題的最優解
  • 無界性: 如果原問題（對偶問題）具有無界解，則其對偶問題（原問題）無可行解
  • 強對偶性: 如果原問題有最優解，則其對偶問題也一定具有最優解，且max \(z=\min w\)

證明：設 \(\mathrm{B}\) 為原問題標準形式的可行解基，且其基解為最優解，則由最優性條件應有檢驗數全部小於等於零：

\[-C_B B^{-1} \leq 0, \theta=C-C_B B^{-1} A \leq 0 \]

從而可得:

\[A^T B^{-T} C_B^T \geq C^T, B^{-T} C_B^T \geq 0 \]

所以 \(B^{-T} C_B^T\) 是對偶問題的一個可行解，且

\[\left(B^{-T} C_B^T\right)^T b=C_B B^{-1} b=\max z \]

由最優性可知， \(B^{-T} C_B^T\) 為對偶問題的最優解，且 \(\max z=\min w\)

  • 互補昖馳性: 線上性規劃問題的最優解中，如果某一約束條件的對偶變數值非零，則該約束條件嚴格取等，反之，如果約束條件取不等式，則對應的對偶變數一定為o，可以表示為 \(Y^T(b-A X)=0\)
    證明: 由強弱偶性可知 \(C X=Y^T b\) ，

\[\begin{aligned} & \text { 又 } \because C X=X^T C^T \leq X^T A^T Y=Y^T A X \leq Y^T b \\ & \therefore Y^T A X=Y^T b \Longrightarrow Y^T(b-A X)=0 \end{aligned} \]

  • 基解互補性: 原問題及其對偶問題之間存在一對互補的基解，其中原問題的松馳變數對應對偶問題的變數對偶問題的剩餘變數對應原問題的變數.這些互相對應的變數如果在一個問題的解中是基變數則在另一個問題的解中是非基變數; 將這對互補的集解分別帶入原問題和對偶問題的目標函式中有 \(z=w\)。

二、對偶單純形法迭代-python實現

流程圖 原理 關係

\[\begin{align*} \text{min } Z &= 15x_1 + 24x_2 + 5x_3 \\ & \text{subject to: } \ & \\ & 6x_2 + x_3 \geq 2 \\ & 5x_1 + 2x_2 + x_3 \geq 1 \\ & x_i \geq 0, \quad (i = 1,2,3) \end{align*} \]

"""
對偶單純形法
min z=15*x1+24*x2+5*x3
      6*x2+x3>=2
      5*x1+2*x2+x3>=1
      x1,x2>=0
標準化後格式
max z=-15*x1-24*x2-5*x3
      -6*x2-x3+x4=-2
      -5*x1-2*x2-x3+x5=-1
      x1,x2,x3,x4,x5>=0     
"""
import numpy as np

#根據線性規劃要改變的資料
Cj = [-15, -24, -5, 0, 0]  # Coefficients of the objective function
constraints_matrix = [[-2, 0, -6, -1, 1, 0], [-1, -5, -2, -1, 0, 1]]  # Constraints matrix (with RHS)
z = [0, -15, -24, -5, 0, 0]  # Z-row
Xb = [4, 5]  # Basis variables (indices of slack variables)
Cb = [0, 0]  # Coefficients of basis variables

#為便於表示而生成
C = [0] + Cj
print(C)

# Output simplex table
def show():
    print("--" * 50)
    print("|               Cj               ", end='|')
    for i in Cj:
        print(f"{i:.2f}".center(10), end='|')
    
    print("--" * 50)
    print("|    Cb    |    Xb    |    b     ", end='|')
    for i in range(len(Cj)):
        print(("X" + str(i + 1)).center(10), end='|')
    
    print("--" * 50)
    for i in range(len(constraints_matrix)):
        print("|", end="")
        print(f"{Cb[i]:.2f}".center(10), end='|')
        print(("X" + str(Xb[i])).center(10), end='|')
        for j in range(len(constraints_matrix[0])):
            print(f"{constraints_matrix[i][j]:.2f}".center(10), end="|")
        
        print("--" * 50)
    print("|          Z          ", end='|')
    for i in range(len(z)):
        print(f"{z[i]:.2f}".center(10), end='|')
    
    print("--" * 50)
    print("**" * 50)
    
show()


def fu():
    global constraints_matrix, z, Xb, Cb  # Declare global variables

    # Extracting the first column
    first_column = [row[0] for row in constraints_matrix]
    print(first_column)

    # Finding the index of the minimum value in the first column
    min_index = first_column.index(min(first_column))
    print("Index of the minimum value in the first column:", min_index)

    # Extract the row corresponding to min_index
    min_row = constraints_matrix[min_index]

    # Calculate the ratios z / min_row
    ratios = [z[i] / min_row[i] if min_row[i] != 0 else float('inf') for i in range(len(z))]

    # Filter out the non-positive values from z and corresponding min_row
    filtered_ratios = [ratios[i] for i in range(len(ratios)) if z[i] < 0 and min_row[i] < 0]

    # Find the minimum ratio and its index in z
    if filtered_ratios:
        min_ratio = min(filtered_ratios)
        min_ratio_index = ratios.index(min_ratio)
        print("Minimum ratio:", min_ratio)
        print("Index of the minimum ratio in z:", min_ratio_index)
    else:
        print("No valid ratio found")
        return

    # Updating Xb
    Xb[min_index] = min_ratio_index 
    print("Updated basis variables Xb:")
    print(Xb)

    # Updating Cb based on Xb
    Cb = [C[Xb[i]] for i in range(len(Xb))]
    print("Updated coefficients of basis variables Cb:")
    print(Cb)

    # Extract columns according to Xb to form matrix B
    B = [[row[Xb[0]], row[Xb[1]]] for row in constraints_matrix]
    print("Matrix B:")
    for row in B:
        print([f"{val:.2f}" for val in row])

    try:
        # Calculate the inverse of matrix B
        B_matrix = np.array(B)
        B_inv = np.linalg.inv(B_matrix)
    except np.linalg.LinAlgError:
        print("Matrix B is singular, skipping this iteration.")
        return

    print("Inverse of matrix B (B_inv):")
    print(B_inv)

    # Update z
    z_update = (np.array(C[1:]) - np.dot(np.array(Cb), np.dot(B_inv, np.array(constraints_matrix)[:, 1:]))).tolist()
    z = [z[0]] + z_update  # Add the first element back to z

    # Update constraints_matrix with the inverse of B
    constraints_matrix = np.dot(B_inv, np.array(constraints_matrix)).tolist()

# Iterate until all elements in the first column of constraints_matrix are non-negative
while any(row[0] < 0 for row in constraints_matrix):
    fu()
    show()

三、對偶單純形法迭代示例

\[\begin{align*} \text{min } Z &= 9x_1 + 12x_2 + 15x_3 \\ & \text{subject to: } \ & \\ & 2x_1 + 2x_2 + x_3 \geq 10 \\ & 2x_1 + 2x_2 + 4x_4 \geq 12 \\ & x_1 + x_2 + 5x_4 \geq 14 \\ & x_i \geq 0, \quad (i = 1,2,3) \end{align*} \]

import numpy as np

Cj = [-9, -12, -15, 0, 0, 0]  # 要求的函式
constraints_matrix = [
    [-10, -2, -2, -1, 1, 0, 0],
    [-12, -2, -3, -1, 0, 1, 0],
    [-14, -1, -1, -5, 0, 0, 1]
]  # 每個子列表第一位為b
z = [0, -9, -12, -15, 0, 0, 0]
Xb = [4, 5, 6]
Cb = [0, 0, 0]

# 為便於表示而生成
C = [0] + Cj

# Output simplex table
def show():
    print("--" * 50)
    print("|               Cj               ", end='|')
    for i in Cj:
        print(f"{i:.2f}".center(10), end='|')
    print()
    
    print("--" * 50)
    print("|    Cb    |    Xb    |    b     ", end='|')
    for i in range(len(Cj)):
        print(("X" + str(i + 1)).center(10), end='|')
    print()
    
    print("--" * 50)
    for i in range(len(constraints_matrix)):
        print("|", end="")
        print(f"{Cb[i]:.2f}".center(10), end='|')
        print(("X" + str(Xb[i])).center(10), end='|')
        for j in range(len(constraints_matrix[0])):
            print(f"{constraints_matrix[i][j]:.2f}".center(10), end='|')
        print()
        
    print("--" * 50)
    print("|          Z          ", end='|')
    for i in range(len(z)):
        print(f"{z[i]:.2f}".center(10), end='|')
    print()
    
    print("--" * 50)
    print("**" * 50)

show()

def fu():
    global constraints_matrix, z, Xb, Cb  # Declare global variables

    # Extracting the first column
    first_column = [row[0] for row in constraints_matrix]
    print("First column:", first_column)

    # Finding the index of the minimum value in the first column
    min_index = first_column.index(min(first_column))
    print("Index of the minimum value in the first column:", min_index)

    # Extract the row corresponding to min_index
    min_row = constraints_matrix[min_index]

    # Calculate the ratios z / min_row
    ratios = [z[i] / min_row[i] if min_row[i] != 0 else float('inf') for i in range(len(z))]
    print("Ratios:", ratios)

    # Filter out the non-positive values from z and corresponding min_row
    filtered_ratios = [ratios[i] for i in range(len(ratios)) if z[i] < 0 and min_row[i] < 0]

    # Find the minimum ratio and its index in z
    if filtered_ratios:
        min_ratio = min(filtered_ratios)
        min_ratio_index = ratios.index(min_ratio)
        print("Minimum ratio:", min_ratio)
        print("Index of the minimum ratio in z:", min_ratio_index)
    else:
        print("No valid ratio found")
        return

    # Updating Xb
    Xb[min_index] = min_ratio_index 
    print("Updated basis variables Xb:")
    print(Xb)

    # Updating Cb based on Xb
    Cb = [C[Xb[i]] for i in range(len(Xb))]
    print("Updated coefficients of basis variables Cb:")
    print(Cb)

    # Extract columns according to Xb to form matrix B
    B = [[row[Xb[0]], row[Xb[1]], row[Xb[2]]] for row in constraints_matrix]
    print("Matrix B:")
    for row in B:
        print([f"{val:.2f}" for val in row])

    try:
        # Calculate the inverse of matrix B
        B_matrix = np.array(B)
        B_inv = np.linalg.inv(B_matrix)
    except np.linalg.LinAlgError:
        print("Matrix B is singular, skipping this iteration.")
        return

    print("Inverse of matrix B (B_inv):")
    print(B_inv)

    # Update z
    z_update = (np.array(C[1:]) - np.dot(np.array(Cb), np.dot(B_inv, np.array(constraints_matrix)[:, 1:]))).tolist()
    z = [z[0]] + z_update  # Add the first element back to z

    # Update constraints_matrix with the inverse of B
    constraints_matrix = np.dot(B_inv, np.array(constraints_matrix)).tolist()

# Iterate until all elements in the first column of constraints_matrix are non-negative
while any(row[0] < 0 for row in constraints_matrix):
    fu()
    show()
    input("按Enter鍵繼續...")

----------------------------------------------------------------------------------------------------
|               Cj               |  -9.00   |  -12.00  |  -15.00  |   0.00   |   0.00   |   0.00   |
----------------------------------------------------------------------------------------------------
|    Cb    |    Xb    |    b     |    X1    |    X2    |    X3    |    X4    |    X5    |    X6    |
----------------------------------------------------------------------------------------------------
|   0.00   |    X4    |  -10.00  |  -2.00   |  -2.00   |  -1.00   |   1.00   |   0.00   |   0.00   |
|   0.00   |    X5    |  -12.00  |  -2.00   |  -3.00   |  -1.00   |   0.00   |   1.00   |   0.00   |
|   0.00   |    X6    |  -14.00  |  -1.00   |  -1.00   |  -5.00   |   0.00   |   0.00   |   1.00   |
----------------------------------------------------------------------------------------------------
|          Z          |   0.00   |  -9.00   |  -12.00  |  -15.00  |   0.00   |   0.00   |   0.00   |
----------------------------------------------------------------------------------------------------
****************************************************************************************************
----------------------------------------------------------------------------------------------------
|               Cj               |  -9.00   |  -12.00  |  -15.00  |   0.00   |   0.00   |   0.00   |
----------------------------------------------------------------------------------------------------
|    Cb    |    Xb    |    b     |    X1    |    X2    |    X3    |    X4    |    X5    |    X6    |
----------------------------------------------------------------------------------------------------
|   0.00   |    X4    |  -7.20   |  -1.80   |  -1.80   |   0.00   |   1.00   |   0.00   |  -0.20   |
|   0.00   |    X5    |  -9.20   |  -1.80   |  -2.80   |   0.00   |   0.00   |   1.00   |  -0.20   |
|  -15.00  |    X3    |   2.80   |   0.20   |   0.20   |   1.00   |   0.00   |   0.00   |  -0.20   |
----------------------------------------------------------------------------------------------------
|          Z          |   0.00   |  -6.00   |  -9.00   |   0.00   |   0.00   |   0.00   |  -3.00   |
----------------------------------------------------------------------------------------------------
****************************************************************************************************
----------------------------------------------------------------------------------------------------
|               Cj               |  -9.00   |  -12.00  |  -15.00  |   0.00   |   0.00   |   0.00   |
----------------------------------------------------------------------------------------------------
|    Cb    |    Xb    |    b     |    X1    |    X2    |    X3    |    X4    |    X5    |    X6    |
----------------------------------------------------------------------------------------------------
|   0.00   |    X4    |  -1.29   |  -0.64   |  -0.00   |   0.00   |   1.00   |  -0.64   |  -0.07   |
|  -12.00  |    X2    |   3.29   |   0.64   |   1.00   |   0.00   |   0.00   |  -0.36   |   0.07   |
|  -15.00  |    X3    |   2.14   |   0.07   |   0.00   |   1.00   |   0.00   |   0.07   |  -0.21   |
----------------------------------------------------------------------------------------------------
|          Z          |   0.00   |  -0.21   |   0.00   |   0.00   |   0.00   |  -3.21   |  -2.36   |
----------------------------------------------------------------------------------------------------
****************************************************************************************************
----------------------------------------------------------------------------------------------------
|               Cj               |  -9.00   |  -12.00  |  -15.00  |   0.00   |   0.00   |   0.00   |
----------------------------------------------------------------------------------------------------
|    Cb    |    Xb    |    b     |    X1    |    X2    |    X3    |    X4    |    X5    |    X6    |
----------------------------------------------------------------------------------------------------
|  -9.00   |    X1    |   2.00   |   1.00   |   0.00   |   0.00   |  -1.56   |   1.00   |   0.11   |
|  -12.00  |    X2    |   2.00   |   0.00   |   1.00   |   0.00   |   1.00   |  -1.00   |   0.00   |
|  -15.00  |    X3    |   2.00   |   0.00   |   0.00   |   1.00   |   0.11   |   0.00   |  -0.22   |
----------------------------------------------------------------------------------------------------
|          Z          |   0.00   |   0.00   |   0.00   |   0.00   |  -0.33   |  -3.00   |  -2.33   |
----------------------------------------------------------------------------------------------------
****************************************************************************************************

總結

對偶理論線上性規劃中主要涉及原問題（稱為"主問題"）和其對應的對偶問題。每個線性規劃問題都有一個對應的對偶問題，透過解決對偶問題，可以獲得關於原問題的重要資訊。對偶問題的解提供了對原問題最優解的界限和性質的洞察。對偶定理表明，如果原問題有最優解，那麼對偶問題也有最優解，並且兩者的目標函式值相等。這種關係在經濟解釋、靈敏度分析和複雜問題分解等方面有廣泛應用。
對偶單純形法是一種用於解決線性規劃問題的迭代演算法，尤其適用於初始解不滿足原問題的可行性條件但滿足對偶問題的可行性條件的情況。該方法從對偶問題的可行解出發，透過逐步調整，使得原問題的約束條件逐漸滿足，最終找到原問題的最優解。具體步驟包括選擇入基變數和出基變數，以保持對偶問題的可行性，並透過對偶單純形表調整基解。

參考文獻

  1. python 實現對偶單純形法
  2. python 實現對偶單純形法
  3. 對偶理論

對偶理論和對偶單純形法——Python實現

相關文章