对偶单纯形法算法精要

单纯形法是线性规划中最经典且广泛应用的求解方法，通过在可行解的边界上移动，逐步逼近最优解。它从一个初始基本可行解开始，不断优化目标函数值，直到找到最优解。对偶单纯形法则是单纯形法的一种变形，尤其适用于特定类型的线性规划问题。不同于标准的单纯形法，对偶单纯形法从一个对偶可行但原始不可行的初始解出发，通过逐步改善解的可行性和最优性，最终找到最优解。对偶单纯形法在快速调整和重新优化方面表现出色，特别是在处理新的约束或变量被添加到现有模型中的情况时，这使得它在实际应用中具有独特的优势。

单纯形法比较	对偶单纯形流程图

一、对偶单纯形算法

‌对偶单纯形法是一种用于求解线性规划问题的数学方法，由美国数学家‌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=0$；最优值为$w = 28/5$。

二、例题

例2：利用对偶单纯形法求解线性规划问题：

$$\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求解

例3：假设食物的各种营养成分、价格如下表：

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. 对偶单纯形法