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

Question

Я пытаюсь реализовать симплексный алгоритм, следуя правилам, которые мне дали на курсе по оптимизации. Проблема в

min c'*x    s.t.
    Ax = b
    x >= 0

Все векторы предполагаются столбцами, ' обозначает транспонирование. Алгоритм также должен возвращать решение двойного LP. Соблюдайте следующие правила:

Здесь 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, Но здесь начинается цикл. Я предполагаю, что этого не должно было случиться, но, следуя правилам, я не вижу, как это предотвратить.

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

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

Answer 1

Проблема была решена. Оказалось, что пункт 7 в примечаниях был неправильным. Вместо этого пункт 7 должен быть