Реализация бесконечного цикла симплексного метода

Здесь A_J обозначает столбцы из A с индексами в J и x_J, x_K обозначает элементы вектора x с индексами в J или K соответственно. Вектор a_s — это столбец s матрицы A.

Теперь я не понимаю, как этот алгоритм заботится о состоянии x >= 0, но я решил попробовать и выполнить его шаг за шагом. Я использовал для этого Matlab и получил следующий код.

X = zeros(n, 1);
Y = zeros(m, 1);

% i. Choose starting basis J and K = {1,2,...,n} \ J
J = [4 5 6]  % for our problem
K = setdiff(1:n, J)

% this while is for goto
while 1
    % ii. Solve system A_J*\bar{x}_J = b.
    xbar = A(:,J) \ b

    % iii. Calculate value of criterion function with respect to current x_J.
    fval = c(J)' * xbar

    % iv. Calculate dual solution y from A_J^T*y = c_J.
    y = A(:,J)' \ c(J)

    % v. Calculate \bar{c}^T = c_K^T - u^T A_K. If \bar{c}^T >= 0, we have
    % found the optimal solution. If not, select the smallest s \in K, such
    % that c_s < 0. Variable x_s enters basis.
    cbar = c(K)' - c(J)' * inv(A(:,J)) * A(:,K)
    cbar = cbar'

    tmp = findnegative(cbar)
    if tmp == -1  % we have found the optimal solution since cbar >= 0
        X(J) = xbar;
        Y = y;
        FVAL = fval;
        return
    end

    s = findnegative(c, K)  %x_s enters basis

    % vi. Solve system A_J*\bar{a} = a_s. If \bar{a} <= 0, then the problem is
    % unbounded.
    abar = A(:,J) \ A(:,s)

    if findpositive(abar) == -1  % we failed to find positive number
        disp('The problem is unbounded.')
        return;
    end

    % vii. Calculate v = \bar{x}_J / \bar{a} and find the smallest rho \in J,
    % such that v_rho > 0. Variable x_rho exits basis.
    v = xbar ./ abar
    rho = J(findpositive(v))

    % viii. Update J and K and goto ii.
    J = setdiff(J, rho)
    J = union(J, s)
    K = setdiff(K, s)
    K = union(K, rho)
end

Функции findpositive(x) и findnegative(x, S) возвращают первый индекс положительного или отрицательного значения в x. S — это набор индексов, на которые мы смотрим. Если S опущено, проверяется весь вектор. Точки с запятой опущены в целях отладки.

Проблема, на которой я тестировал этот код, заключается в

c = [-3 -1 -3 zeros(1,3)];
A = [2 1 1; 1 2 3; 2 2 1];
A = [A eye(3)];
b = [2; 5; 6];

Причина для zeros(1,3) и eye(3) в том, что проблема заключается в неравенствах и нам нужны переменные резерва. Я установил начальный базис равным [4 5 6], потому что в примечаниях говорится, что начальный базис должен быть установлен на слабых переменных.

Теперь во время выполнения происходит то, что при первом запуске while переменная с индексом 1 входит в базу (в Matlab индексы начинаются с 1), а 4 выходит из нее, и это разумно. При втором прогоне в базу входит 2 (так как это наименьший индекс такой, что из него выходят c(idx) < 0 и 1. Но теперь на следующей итерации снова входит 1 в базу и я понимаю почему он входит, ведь это самый маленький индекс, такой что c(idx) < 0 Но тут начинается зацикливание Я предполагаю, что этого не должно было случиться, но следуя правилам я не вижу, как это предотвратить.

Я предполагаю, что должно быть что-то не так с моей интерпретацией заметок, но я просто не вижу, где я не прав. Я также помню, что когда мы решали ЛП на бумаге, мы каждый раз обновляли нашу субъективную функцию, поскольку, когда переменная входила в базис, мы удаляли ее из субъективной функции и выражали эту переменную в subj. функцию с выражением из одного из равенств, но я предполагаю, что это другой алгоритм.

Любые замечания или помощь будут высоко оценены.