问题描述
我正在尝试使用 Python 和 numpy 实现 revised simplex method (RSM) 算法。我坚持要么只处理最大化(在 2x4、5x5 等小矩阵上正确,而在较大矩阵上错误),要么在大多数最小化情况下进入无限循环。下面的代码演示了我实现最小化的尝试:
def __init__(self,A: np.ndarray,b: np.ndarray,c: np.ndarray):
base_size = A.shape[0] # Number of basis vectors
I = np.eye(base_size) # Identity matrix,an initial basis
self.extended_matrix = np.concatenate((A,I),axis=1) # Extended matrix of the system
self.free_size = A.shape[1] # Number of free vectors
self.b = b # Right parts of the constraints
self.base_inv = I # Initial inverse basis matrix
self.c = np.concatenate((c,np.zeros(base_size))) # Objective function quotients including those related to the basis variables
self.c_i = [i for i,coeff in enumerate(self.c)] # Indices for all variables
@property
def c_f_indices(self):
"""
Indices of the free variables.
"""
return self.c_i[:self.free_size]
@property
def c_T_B(self):
"""
Objective function quotients related to the basis variables.
"""
c_B_indices = self.c_i[self.free_size:] # Basis variables indices.
return self.c[c_B_indices]
@property
def c_T_f(self):
"""
Objective function quotients related to the free variables.
"""
return self.c[self.c_f_indices]
@property
def free(self):
"""
Free vectors.
"""
return self.extended_matrix[:,self.c_f_indices]
@property
def y_T(self):
"""
Lagrange multipliers.
"""
return self.c_T_B @ self.base_inv
@property
def deltas(self):
"""
Net evaluations.
"""
return (self.y_T @ self.free) - self.c_T_f
def _swap(self,exits: int,enters: int) -> None:
"""
In-/excluding respective vectors into/from the basis.
"""
self.c_i[enters],self.c_i[exits + self.free_size] = self.c_i[exits + self.free_size],self.c_i[enters]
def optimize(self):
while any([delta > 0 for delta in self.deltas]): # < 0 in case of maximization
x_B = self.base_inv @ self.b # Current basis variables
enters_base = np.argmax(self.deltas) # Vector to enter the basis; argmin in case of maximization
# Finding the vector to leave the basis:
alpha = self.base_inv @ self.free[:,enters_base]
try:
exits_base = np.argmin([b/a if a > 0 else np.inf for b,a in zip(x_B,alpha)])
assert alpha[exits_base] != 0
except (ValueError,AssertionError):
raise Exception("No solutions")
# Finding the E_r matrix,which current inverse basis will be left-multiplied with in order to achieve the new inverse basis:
zetas = [-alpha[j] / alpha[exits_base] if j != exits_base else 1/alpha[exits_base] for j,a_j in enumerate(alpha)]
E = np.eye(self.free.shape[0])
E[:,exits_base] = zetas
self.base_inv = E @ self.base_inv
# In-/excluding respective vectors into/from the basis:
self._swap(exits_base,enters_base)
return self.c_T_B @ self.base_inv @ self.b # Final objective function value
我也尝试过对 c_f_indices 进行排序,但仍然陷入无限循环。 similar RSM implementation 在较大的矩阵(例如 16x8)上也会产生错误的结果,但在较小的矩阵上效果很好。
问题的根源在哪里?
解决方法
正如 Erwin Kalvelagen 已经提到的,正常的 Dantzig 枢轴规则会导致循环和停滞,在这种情况下,目标值在很长一段时间内没有改善。
一般来说,这种现象被称为 LP 的退化。有限的数值精度和舍入误差会导致问题的退化。这就是为什么它通常在大型 LP 中更为普遍。
有几种方法可以解决这个问题。正如 Erwin 所提到的,最常用的是扰动。但是,如果您将此作为学习项目进行,我建议您使用一种解决方案,在该解决方案中您使用更精细的透视规则,例如 Zadeh's Rule 或 Cunningham's Rule,您只需保留一个表或变量以确保循环后,您可以选择一个不同的变量来输入基数。