Пересмотренный симплексный метод входит в бесконечный цикл

Я пытаюсь реализовать алгоритм пересмотренный симплекс-метод (RSM) с использованием Python и numpy. Я застрял либо в том, что он работает только с максимизацией (правильно для крошечных матриц, таких как 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, но все равно получаю бесконечный цикл. похожая реализация RSM также дает неверные результаты на больших матрицах (например, 16x8) и работает на малюсеньких нормально.

Где корень проблемы?