Решение задачи можно искать с помощью простого симплексного метода, например, как в такой программе (программа написана для компилятора Дельфи или FPC, но должна пойти и на TP, если изменить пару строк):
program simple_sim;
{$APPTYPE CONSOLE}
uses SysUtils; const mm = 100; nn = 100;
var A : array[1..mm, 1..nn] of double; fun : array[1..nn] of integer; // Коэффициенты целевой функции m, n : integer; // m ограничений, n переменных.
basis : array[1..nn] of integer; // Здесь храним номера базисных переменных i, j : integer; x : array[1..nn] of double; // Здесь будут значения переменных при расшифровке плана
procedure solve; var i, j, i0, j0 : integer; tmp : double; opt : boolean; begin opt := false; repeat j0 := 1; i0 := 0; while (j0 < m+n+1) and (A[m+1, j0] >= 0) do inc(j0); if A[m+1, j0] >= 0 then opt := true;
if not opt then begin tmp := 10000; for i := 1 to m do if (A[i, j0] > 0) and (A[i, m+n+1] / A[i, j0] < tmp) then begin tmp := A[i, m+n+1] / A[i, j0]; i0 := i end; // i0 - выводим, j0 - добавляем basis[i0] := j0; // Ввод нового элемента в базис // [i0, j0] - ведущий эл-т в Гауссе: for i := 1 to m + 1 do if i <> i0 then begin tmp := A[i, j0]; for j := 1 to m + n + 1 do A[i,j] := A[i,j] - A[i0,j]*tmp/A[i0,j0]; end; tmp := A[i0, j0]; for j := 1 to m + n + 1 do A[i0, j] := A[i0, j] / tmp; end; until opt; end;
begin assign(input, 'input.txt'); reset(input); // -------Ввод данных--------------------------- read(n); read(m);
for i := 1 to n do read(fun[i]); //Читаем коэффициенты целевой функции
for i := 1 to m do for j := 1 to n do read(A[i, j]);
for i := 1 to m do read(A[i, n+m+1]); // Читаем правые части ограничений
for i := 1 to m do // Вводим дополнительные переменные A[i, n+i] := 1; fillchar(A[m+1], sizeof(A[m+1]), 0);
// базис из доп. переменных for i := 1 to m do basis[i] := n + i; for j := 1 to n do A[m+1,j] := -fun[j]; // Оценки для небазисных переменных = -fun[j], для базисных - 0
solve; // DO IT! +)
// -- вывод базиса -- for i := 1 to m do if basis[i] <= n then x[basis[i]] := A[i, m+n+1];
for i := 1 to n do writeLn('x[', i, '] = ', x[i]:0:3); writeLn('min f(x) = ', A[m+1, m+n+1]:0:3); end.
Пояснения: Симплекс-метод решает ЗЛП в стандартной форме:
Код
C*X -> max A*X = B X >= 0
А дана задача в канонической форме. Чтобы привести задачу из исходной (канонической) формы, вводим дополнительные переменные:
Код
C*X + 0*Y -> max A*X + E*Y = B X, Y >= 0
где Е - единичная матрица.
Все структуры данных, нужные для работы метода, хранятся в одной симплексной таблице, которая перед вызовом самого метода должна выглядеть так:
На каждом шаге алгоритма: 1) Выбираем первый "слева" столбец j0, в котором в последней строке - отрицательное число. Если такого столбца нет, то заканчиваем работу алгоритма 2) Выбираем строку, для которой отношение <элемент последнего столбца> / A[i, j0] минимально. 3) Выполняем со всей таблицей преобразование Жордана-Гаусса с ведущим столбцом j0 и ведущей строкой i0, то есть добиваемся, чтобы в столбце j0 все элементы стали равными 0, кроме i0-го элемента, который будет равен 1
После окончания работы алгоритма в последнем столбце будут храниться значения переменных, на которых достигается максимум функции, а само значение функции будет лежать в правом нижнем углу таблицы.
Замечание: в найденном решении будет m (а то и меньше) неотрицательных элементов, а не n. Такое решение называется базисным.
Как вариант, можно выбирать столбец j0 как столбец с максимальным по модулю отрицательным нижним элементом. Так алгоритм потребует меньше итераций, но может зациклиться (а то правило выбора пары i0, j0, которое используется в этом алгоритме, гарантирует, что зацикливаться программа не будет. Правило называется правилом Бленда).
P.S. Пример входного файла, которым питается программа:
Код
3 2 3 4 5 1 5 3 4 1 2 200 160
Формат входных данных: n m - размерность матрицы fun - коэффициенты целевой функции A - сама матрица b - правые части ограничений