Симплекс-метод решения задачи ЛП, Простой симплексный метод Опции

[trminator]

Задача:
Требуется найти максимум функции:

Код

c[1]*x[1] + c[2]*x[2] + ... + c[n]*x[n] -> max

при ограничениях:

Код

a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] <= b[1]
...
a[m,1]*x[1] + a[m,2]*x[2] + ... + a[m,n]*x[n] <= b[m]

x[i] >= 0, x = 1, ..., n

Другой вариант записи задачи:

Код

C*X -> max
A[M,N]*X[N] <= B[M]
X >= 0

Решение задачи можно искать с помощью простого симплексного метода, например, как в такой программе (программа написана для компилятора Дельфи или 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        │   │
   │                 │   1      │   │
   │      A[M,N]     │     1    │ B │
   │                 │       1  │   │
   │                 │         1│   │
   ├─────────────────┼──────────┼───┤
   │     -С          │    0     │ 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 - правые части ограничений

Правильный ответ для этого примера:

Код

x[1] = 8.000
x[2] = 0.000
x[3] = 64.000
min f(x) = 344.000