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