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

郝hai發表於2024-06-01

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

一、對偶理論

若原問題為 \(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   |
----------------------------------------------------------------------------------------------------
****************************************************************************************************

總結

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

參考文獻

python 實現對偶單純形法
python 實現對偶單純形法
對偶理論

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

解析對偶理論與對偶單純性法
2021-09-11
【BZOJ3265】志願者招募加強版線性規劃單純形法對偶原理
2015-05-12
UA SIE545 優化理論基礎4 對偶理論簡介4 求解對偶問題的割平面演算法
2020-12-02
優化演算法
關於原始對偶演算法（拉格朗日對偶）
2017-11-04
演算法
線性規劃對偶原理
2018-01-09
線性規劃的對偶問題——由拉格朗日對偶問題匯出
2024-06-23
調整陣列順序使奇數位於偶數前面,偶數和偶數之間的相對位置不變
2016-04-25
陣列
深入理解SVM之對偶問題
2018-04-11
運籌學——matlab實現單純形法
2020-12-17
Matlab
線性規劃對偶學習筆記
2024-06-10
筆記
Andrew Ng機器學習筆記+Weka相關演算法實現（四）SVM和原始對偶問題
2017-08-08
機器學習筆記演算法
ICML 2018 | 模型層面的對偶學習
2018-07-10
模型
機器學習，詳解SVM軟間隔與對偶問題
2020-09-09
機器學習
深入理解SVM，軟間隔與對偶問題
2020-09-03
對偶學習_劉鐵巖_微軟亞洲研究院
2018-03-18
微軟
簡易解說拉格朗日對偶（Lagrange duality）
2015-04-21
一文理解拉格朗日對偶和KKT條件
2018-09-09
【Oracle】對RID為偶數的記錄進行操作
2017-09-29
Oracle
偶的Java之路
2005-11-09
Java
HNOI2016礦區（計算幾何+對偶圖）
2019-02-23
php的oci_connect偶爾出現OCIEnvNlsCreate fail和缺失LD_LIBRARY_PATH的處理辦法
2011-08-30
PHPAI
運籌優化（四）--線性規劃之對偶問題和靈敏度分析
2019-01-12
優化
純js實現 vue 元件與 vue 單檔案元件對比
2024-04-10
JSVue元件
Table tr 奇數偶數行漸變色Jquery實現
2019-05-11
jQuery
偶滴工作環境
2010-01-21
偶翻譯的小說
2007-09-24
求陣列內所有偶數的和
2020-12-15
陣列
Python使用TCP實現簡單對話
2020-10-06
PythonTCP
Python實現簡單的excel對比工具
2021-04-22
PythonExcel
從朋友簽名偶得
2016-01-14
MySQL:雙主單寫主庫偶爾出現大量延遲的原因
2020-06-21
MySql
線性規劃之單純形演算法矩陣描述與python實現
2021-11-17
演算法矩陣Python
怎麼判斷Python數字中的偶數
2021-09-11
Python
純Python和PyTorch對比實現迴圈神經網路RNNCell及反向傳播
2018-12-15
PythonPyTorch神經網路RNN反向傳播
計算1000以內所有偶數和.
2013-11-11
陣列進行奇數和偶數操作把奇數放在陣列前面偶數放在陣列後面
2017-04-25
陣列
偶的伺服器：nio,困惑！
2006-06-28
伺服器
單例模式實現對比
2018-05-22
單例模式

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

一、對偶理論

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

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

總結

參考文獻

相關文章