对偶理论和对偶单纯形迭代——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)
# 根据 Xb 提取列以形成矩阵 B
B = [[row[xb] for xb in Xb] for row in constraints_matrix]
print("矩阵 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 + 3x_2 + x_3 \geq 12 \\
& x_1 + x_2 + 5x_3\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)
# 根据 Xb 提取列以形成矩阵 B
B = [[row[xb] for xb in Xb] for row in constraints_matrix]
print("矩阵 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("按回车键继续...")
----------------------------------------------------------------------------------------------------
| 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 |
----------------------------------------------------------------------------------------------------
****************************************************************************************************
总结
对偶理论在线性规划中主要涉及原问题（称为"主问题"）和其对应的对偶问题。每个线性规划问题都有一个对应的对偶问题，通过解决对偶问题，可以获得关于原问题的重要信息。对偶问题的解提供了对原问题最优解的界限和性质的洞察。对偶定理表明，如果原问题有最优解，那么对偶问题也有最优解，并且两者的目标函数值相等。这种关系在经济解释、灵敏度分析和复杂问题分解等方面有广泛应用。
对偶单纯形法是一种用于解决线性规划问题的迭代算法，尤其适用于初始解不满足原问题的可行性条件但满足对偶问题的可行性条件的情况。该方法从对偶问题的可行解出发，通过逐步调整，使得原问题的约束条件逐渐满足，最终找到原问题的最优解。具体步骤包括选择入基变量和出基变量，以保持对偶问题的可行性，并通过对偶单纯形表调整基解。
参考文献
