Пытаюсь реализовать
модифицированный симплекс-метод в мультипликативной форме на Python с использованием numpy. Постоянно сталкиваюсь с проблемой того, что он у меня работает только на максимизацию, причем на маленьких матрицах - большие считает с ошибками, значения целевой функции сильно отличаются от тех, которыми они должны быть. На минимизации чаще всего зацикливается. Теория изложена, нарпимер,
здесь. Вот пример моего кода для минимизации:
class Simplex(object):
def __init__(self, A: np.ndarray, b: np.ndarray, c: np.ndarray):
base_size = A.shape[0] # Кол-во базисных векторов
I = np.eye(base_size) # Единичная матрица - начальный базис системы.
self.extended_matrix = np.concatenate((A, I), axis=1) # Расширенная матрица системы
self.free_size = A.shape[1] # Кол-во свободных векторов
self.b = b # Правые части ограничений
self.base_inv = I # Начальная обратная матрица базиса равна его прямой матрице, т.к. является единичной
self.c = np.concatenate((c, np.zeros(base_size))) # Коэффициенты целевой ф-ции с учетом базисных переменных
self.c_i = [i for i, coeff in enumerate(self.c)] # Индексы всех переменных
@property
def c_f_indices(self):
"""
Индексы свободных переменных.
"""
return self.c_i[:self.free_size]
@property
def c_T_B(self):
"""
Коэффициенты целевой функции при базисных переменных.
"""
c_B_indices = self.c_i[self.free_size:] # Индексы базисных переменных.
return self.c[c_B_indices]
@property
def c_T_f(self):
"""
Коэффициенты целевой функции при свободных переменных.
"""
return self.c[self.c_f_indices]
@property
def free(self):
"""
Свободные векторы.
"""
return self.extended_matrix[:, self.c_f_indices]
@property
def y_T(self):
"""
Вектор относительных оценок (разрешающие множители по Канторовичу).
"""
return self.c_T_B @ self.base_inv
@property
def deltas(self):
"""
Текущие оценки векторов.
"""
return (self.y_T @ self.free) - self.c_T_f
def _swap(self, exits: int, enters: int) -> None:
"""
Включение/исключение в/из базис(-а) соответствующих векторов
"""
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
x_B = self.base_inv @ self.b # Текущие значения базисных переменных
enters_base = np.argmax(self.deltas) # Нахождение вектора, включаемого в базис. В случае максимизации будет argmin
# Нахождение вектора, исключаемого из базиса:
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")
# Нахождение матрицы соотношения E_r и вычисление новой обратной матрицы базиса системы:
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
# Включение/исключение в/из базис(-а) соответствующих векторов:
self._swap(exits_base, enters_base)
return self.c_T_B @ self.base_inv @ self.b # Итоговое значение целевой функции
Я пробовал, в том числе, сортировать индексы свободных переменных, но все равно приходил к зацикливанию. И, кстати, нашел
чью-то похожую реализацию, которая тоже допускает грубые ошибки на больших матрицах. Две недели ломания головы ни к чему не привели. Где может крыться ошибка? Заранее благодарю.
