對偶單純形法演算法精要

郝hai發表於2024-09-01

單純形法是線性規劃中最經典且廣泛應用的求解方法,透過在可行解的邊界上移動,逐步逼近最優解。它從一個初始基本可行解開始,不斷最佳化目標函式值,直到找到最優解。對偶單純形法則是單純形法的一種變形,尤其適用於特定型別的線性規劃問題。不同於標準的單純形法,對偶單純形法從一個對偶可行但原始不可行的初始解出發,透過逐步改善解的可行性和最優性,最終找到最優解。對偶單純形法在快速調整和重新最佳化方面表現出色,特別是在處理新的約束或變數被新增到現有模型中的情況時,這使得它在實際應用中具有獨特的優勢。

單純形法比較 對偶單純形流程圖

一、對偶單純形演算法

‌對偶單純形法是一種用於求解線性規劃問題的數學方法,由美國數學家‌C.萊姆基在1954年提出。這種方法基於‌對偶理論,透過迭代過程逐步搜尋原始問題的最優解。與單純形法不同,對偶單純形法從滿足對偶可行性條件出發,透過迭代逐步逼近原始問題的最優解。在迭代過程中,始終保持基解的對偶可行性,使不可行性逐步消失。‌
對偶單純形法的基本思想是:從原規劃的一個對偶可行解出發;然後檢驗原規劃的基本解是否可行,即是否有負的分量,如果有小於零的分量,則進行迭代,求另一個基本解,此基本解對應著另一個對偶可行解(檢驗數非正)。如果得到的基本解的分量皆非負,則該基本解為最優解。也就是說,對偶單純形法在迭代過程中始終保持對偶解的可行性(即檢驗數非正),使原規劃的基本解由不可行逐步變為可行,當同時得到對偶規劃與原規劃的可行解時,便得到原規劃的最優解。

1.1 對偶問題演算法步驟

a)根據線性規劃典式形式,建立初始對偶單純形表。此表對應原規劃的一個基本解。要求檢驗數行各元素一定非正,原規劃的基本解允許有小於零的分量;
b)若基本解的所有分量皆非負,則得到原規劃的最優解,停止計算;若基本解中有小於零的分量\(b_j < 0\),並且\(b_j\)所在行各系數\(a_{ij} \geq 0\),則原規劃沒有可行解,停止計算;若\(b_j < 0\),並且存在\(a_{ij} < 0\),則確定\(x_i\)為出基變數,並計算

\[\theta = \min \left\{ \frac{\sigma_j}{a_{ij}} \middle| a_{ij} < 0 \right\} = \frac{\sigma_k}{a_{ik}} \]

確定\(x_k\)為進基變數。若有多個\(b_j < 0\),則選擇最小的進行分析計算;
c)以\(a_{ik}\)為中心元素,按照與單純形類似的方法,在表中進行迭代計算,返回步驟 b)。

1.2 練習

例1:用對偶單純形法求解下面線性規劃

\[\begin{aligned} \text{min}\quad w= 2x_1 + 3x_2 + 4x_3 \\ \text{s.t. } \begin{cases} x_1 + 2x_2 + x_3 \geq 3 \\ 2x_1 - x_2 + 3x_3 \geq 4 \\ x_1, x_2, x_3 \geq 0 \end{cases} \end{aligned} \]

\[\begin{aligned} \text{max } \quad w= -2x_1 - 3x_2 - 4x_3 \\ \text{s.t. } \begin{cases} -x_1 - 2x_2 - x_3 + x_4 = -3 \\ -2x_1 + x_2 - 3x_3 + x_5 = -4 \\ x_1, \ldots, x_5 \geq 0 \end{cases} \end{aligned} \]

對偶單純形法演算法精要

得原問題的最優解為$x_1 = 11/5,x_2 = 2/5,x_3 = 3$;最優值為$w = 28/5$。

二、例題

利用對偶單純形法求解線性規劃問題:

\[\begin{cases} \min f(x) = 4x_{1} + 12x_{2} + 18x_{3} \\ \text{s.t. } x_{1} + 3x_{3} \geq 3 \\ \quad \quad 2x_{2} + 2x_{3} \geq 5 \\ \quad \quad x_{1}, x_{2}, x_{3} \geq 0 \end{cases} \]

解:先標準化:

\[\begin{cases} \max f(x) =- 4x_{1} - 12x_{2} - 18x_{3} \\ \text{s.t.}\quad x_{1} + 3x_{3} - x_{4} = 3 \\ \quad \quad 2x_{2} + 2x_{3} - x_{5} = 5 \\ \quad \quad x_{1}, x_{2}, x_{3}, x_{4}, x_{5} \geq 0 \end{cases} \]

為了將 $$ B = (p_{4}, p_{5}) $$ 作為初始基,對約束條件兩邊同乘 \(-1\)

\[\begin{cases} -x_{1} - 3x_{3} + x_{4} = -3 \\ -2x_{2} - 2x_{3} + x_{5} = -5 \end{cases} \]

\[C^{T}-c_{B}^{T}BA = (-4, -12, -18, 0, 0)^{T} \]

對應的初始單純形表如下:

變數 $$ x_{1} $$ $$ x_{2} $$ $$ x_{3} $$ $$ x_{4} $$ $$ x_{5} $$
$$ f $$ 0 -4 -12 -18 0 0
$$ x_{4} $$ -3 -1 0 -3 1 0
$$ x_{5} $$ -5 0 -2 -2 0 1

當前解 $$ x = (0, 0, 0, -3, -5)^{T} $$ 為非可行解,但是它是對偶可行解。

因為:

\[\min\left\{ -3, -5 \right\} = -5, \quad r = 2, \quad x_{j2} = x_{5} \text{ 為離基變數} \]

\[\frac{b_{0s}}{b_{rs}} = \min\left\{ \frac{b_{0j}}{b_{rj}} \mid b_{rj} < 0 \right\} = \min\left\{ \frac{-12}{-2}, \frac{-18}{-2} \right\} = 6 = \frac{b_{02}}{b_{22}}, \quad s = 2 \]

因此,進基變數 $$ x_{s} = x_{2} $$,且 $$ b_{22} = -2 $$ 為主元。

經旋轉變換後,有:

變數 $$ x_{1} $$ $$ x_{2} $$ $$ x_{3} $$ $$ x_{4} $$ $$ x_{5} $$
$$ f $$ 30 -4 0 -6 0 -6
$$ x_{4} $$ -3 -1 0 -3 1 0
$$ x_{2} $$ 5/2 0 1 1 0 -1/2

此時 $$ b_{i0} \geq 0 , (i=1,...,m) $$ 不成立,仍為不可行解。

因為:

\[b_{r0} = \min\left\{b_{i0} \mid b_{i0} < 0, i=1,...,m \right\} = \min\left\{ -3\right\} = -3, \quad x_{4} \text{ 為離基變數}, \quad r = 1 \]

\[\frac{b_{0s}}{b_{rs}} = \min\left\{ \frac{b_{0j}}{b_{rj}} \mid b_{rj} < 0 \right\} = \min\left\{ \frac{-4}{-1}, \frac{-6}{-3} \right\} = 2 = \frac{b_{03}}{b_{r3}}, \quad s = 3 \]

因此,進基變數 $$ x_{s} = x_{3} $$,且 $$ b_{13} = -3 $$ 為主元。

經旋轉變換後,有:

變數 $$ x_{1} $$ $$ x_{2} $$ $$ x_{3} $$ $$ x_{4} $$ $$ x_{5} $$
$$ f $$ 36 -2 0 0 -2 -6
$$ x_{3} $$ 1 1/3 0 1 -1/3 0
$$ x_{2} $$ 3/2 -1/3 1 0 1/3 -1/2

此時 $$ b_{i0} \geq 0, (i=1,...,m) $$ 成立,是可行解,演算法終止。

最優解和最優值:

\[\bar{x} = \left(0, \frac{3}{2}, 1\right)^{T} \quad f(\bar{x}) = 36 \]

三、對偶單純形的Python求解

假設食物的各種營養成分、價格如下表:

Food Energy (能量) Protein (蛋白質) Calcium (鈣) Price
Oatmeal (燕麥) 110 4 2 3
Whole milk (全奶) 160 8 285 9
Cherry pie (櫻桃派) 420 4 22 20
Pork with beans (豬肉) 260 14 80 19

要求我們買的食物中,至少要有2000的能量,55的蛋白質,800的鈣,怎樣買最省錢?

設買燕麥、全奶、櫻桃派、豬肉為\(x_1, x_2, x_3, x_4\),於是可以寫出如下的不等式組:

\[\begin{aligned} & \min \quad 3x_1 + 9x_2 + 20x_3 + 19x_4 \quad \text{money} \\ & \text{s.t.} \\ & 110x_1 + 160x_2 + 420x_3 + 260x_4 \geq 2000 \quad \text{energy} \\ & 4x_1 + 8x_2 + 4x_3 + 14x_4 \geq 55 \quad \text{protein} \\ & 2x_1 + 285x_2 + 22x_3 + 80x_4 \geq 800 \quad \text{calcium} \\ & x_1, x_2, x_3, x_4 \geq 0 \end{aligned}\]

"""
對偶單純形法
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 = [-3, -9, -20, -19, 0, 0, 0]  # Coefficients of the objective function
constraints_matrix = [
    [-2000, -110, -160, -420, -260, 1, 0, 0],
    [-55, -4, -8, -4, -14, 0, 1, 0],
    [-800, -2, -285, -22, -80, 0, 0, 1]
]  # Constraints matrix (with RHS)
z = [0, -3, -9, -20, -19, 0, 0, 0]  # Z-row
Xb = [5, 6, 7]  # Basis variables (indices of slack variables)
Cb = [0, 0, 0]  # Coefficients of basis variables

# 為便於表示而生成
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)

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

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

    # Finding the index of the minimum value in the first column
    min_index = first_column.index(min(first_column))

    # 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)
    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)

    # 根據 Xb 提取列以形成矩陣 B
    B = [[row[xb] for xb in Xb] for row in constraints_matrix]

    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

    # 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()

# 輸出最優解和最優值
def print_optimal_solution():
    # 最優解:基變數的值
    optimal_solution = [0] * len(Cj)  # 初始所有變數的值為0
    for i, xb in enumerate(Xb):
        optimal_solution[xb - 1] = constraints_matrix[i][0]  # 設定基變數的值

    # 最優值:計算目標函式值
    optimal_value = sum(Cj[i] * optimal_solution[i] for i in range(len(Cj)))
    optimal_value = -optimal_value  # 取相反數

    # 列印最優解和最優值,保留兩位小數
    print(f"最優解(變數的取值):{[round(val, 2) for val in optimal_solution]}")
    print(f"最優值(目標函式值的相反數):{optimal_value:.2f}")

# 迭代執行對偶單純形法
while any(row[0] < 0 for row in constraints_matrix):
    fu()
    show()

# 輸出最後的最優解和最優值
print_optimal_solution()
#最終的單純形表
Updated basis variables Xb:
[1, 6, 2]
Updated coefficients of basis variables Cb:
[-3, 0, -9]
----------------------------------------------------------------------------------------------------
|               Cj               |  -3.00   |  -9.00   |  -20.00  |  -19.00  |   0.00   |   0.00   |   0.00   |
----------------------------------------------------------------------------------------------------
|    Cb    |    Xb    |    b     |    X1    |    X2    |    X3    |    X4    |    X5    |    X6    |    X7    |
----------------------------------------------------------------------------------------------------
|  -3.00   |    X1    |  14.24   |   1.00   |   0.00   |   3.74   |   1.98   |  -0.01   |   0.00   |   0.01   |
|   0.00   |    X6    |  23.63   |   0.00   |   0.00   |  11.38   |  -3.96   |  -0.04   |   1.00   |  -0.01   |
|  -9.00   |    X2    |   2.71   |   0.00   |   1.00   |   0.05   |   0.27   |   0.00   |   0.00   |  -0.00   |
----------------------------------------------------------------------------------------------------
|          Z          |   0.00   |   0.00   |   0.00   |  -8.31   |  -10.67  |  -0.03   |   0.00   |  -0.02   |
----------------------------------------------------------------------------------------------------

最優解(變數的取值):[14.24, 2.71, 0, 0, 0, 23.63, 0]
最優值(目標函式值的相反數):67.10

總結

單純形法和對偶單純形法是線性規劃中最常用的兩種演算法,它們在解決最佳化問題方面發揮著關鍵作用。單純形法透過從一個初始的基本可行解出發,沿著多面體的邊逐步移動,尋找使目標函式值不斷改善的鄰近頂點,直到找到最優解。相較之下,對偶單純形法則採用了不同的思路。它從一個對偶可行解開始,重點在於原問題的對偶問題的可行性,而非原問題的直接可行性。透過逐步調整,使得對偶條件和原始問題的可行性條件同時得到滿足,從而找到最優解。對偶單純形法特別適用於那些初始解不滿足原問題約束但滿足對偶約束的情況。這在需要快速適應變化的約束條件或在求解過程中頻繁調整約束的場景中具有明顯優勢。
兩種方法的有效結合和運用,可以大幅提升線性規劃問題的求解效率和靈活性。理解和掌握這兩種方法對於最佳化問題的解決至關重要,它不僅幫助最佳化計算速度和準確性,還能在面對不同型別的線性規劃問題時,選擇最適合的演算法。無論是在理論研究還是實際應用中,這兩種方法都提供了強大的工具,幫助我們應對各種複雜的決策問題。

對偶單純形法演算法精要

參考文獻

  1. 線性規劃-單純形演算法詳解
  2. 對偶單純形法

相關文章