ai , b
IRm , i=1, ... ,n , c
IRn, требуется найти вектор x*
IRn такой, что1. Введение. Настоящая статья есть развитие работы [1] и посвящена вырожденным задачам линейного программирования [2]--[9]. При решении таких задач симплексным методом часто возникают большие последовательности итераций, на которых целевая функция практически не изменяется. Вероятность появления таких вырожденных итераций резко возрастает с увеличением размерности задачи и зачастую делает применение симплексного алгоритма бесполезным. Большинство классических методов теории вырожденного линейного программирования было посвящено лишь борьбе с зацикливанием, избежать вырожденных итераций при использовании таких методов не удавалось. Эти способы сводились либо к специальному выбору выводимого из базиса вектора (лексикографическое правило и правило случайного выбора [3]), либо к выбору вводимого в базис вектора (правило Данцига [4]), либо к одновременному выбору обоих этих векторов (правило Блэнда [2]).
По-видимому, наиболее эффективным методом борьбы с вырожденностью является метод Вольфа, предложенный впервые в [8] и модифицированный в работе [7]. Этот метод требует решения на каждой вырожденной итерации вспомогательной задачи линейного программирования. Ее решение позволяет либо сделать вывод об оптимальности текущего базиса основной задачи, либо заменить в нем сразу несколько векторов, что приводит к уменьшению целевой функции. Таким образом, процесс поиска оптимального решения исходной задачи становится строго монотонным. Однако метод Вольфа имеет два недостатка. Во-первых,размерность вспомогательной задачи совпадает с размерностью исходной. Кроме того, при решении вспомогательной задачи также могут возникнуть вырожденные итерации. В этом случае в процедуре преодоления вырожденности необходимо применять рекурсию.
В работе [1] предложен алгоритм, предусматривающий решение на каждой вырожденной итерации вспомогательной задачи линейного программирования, которая невырождена с вероятностью 1, и имеет меньше строк и столбцов в матрице ограничений, чем исходная. В настоящей работе дано оригинальное и дополненное изложение нового алгоритма, устраняющее некоторые неточности статьи [1], а также впервые представлены результаты решения с его помощью задач большой размерности. Кроме того, получены полезные модификации алгоритма и оценки для оптимума.
Сравнение предлагаемого алгоритма с алгоритмом Вольфа в модификации [7] показывает, что во многих случаях новый алгоритм окзывается более эффективным, хотя и требует дополнительных вычислений при построении вспомогательной задачи.
Будем рассматривать задачу линейного программирования,в которой
по заданным векторам
ai , b
IRm , i=1, ... ,n , c
IRn, требуется найти вектор x*
IRn такой, что
cTx*=
{ c Tx :
xi
ai = b, x = (x1 , ... , xn )T
0 } . (1)
Пусть x -- допустимое базисное решение задачи (1),
характеризуемое тем, что линейно независимы все k векторов-условий
ai , отвечающих положительным компонентам вектора x.
Очевидно, k 
m. В случае k < m текущее
допустимое базисное решение называется вырожденным, а число
m-k -- степенью вырожденности. Матрицу U, состоящую из
столбцов
ai , i
+ 
{ i : xi > 0, 1 i n },
будем называть матрицей строгого базиса[1] в отличие от составной базисной матрицы B=(U,V), где V -- такая произвольная m x (m-k) -матрица (не обязательно составленная из векторов ai ), что матрица B невырождена. Тогда пространство IRm разлагается в прямую сумму
IRm=
(U)
(V)
подпространств (U) и
(V) размерности k и m-k,
базисами которых являются столбцы матриц U и V соответственно.
Иными словами,
| ai= |  |
ui , ui+vi , |
i+ , i 0
{1,...,n} \ + , |
ui(U),
vi(V). |
(2) |
Пусть 
0
IRm -- любое решение неоднозначно разрешимой
относительно вектора двойственных переменных
системы уравнений
ci - Tai=0,
i+
T
U=cuT , (3)
где cu состоит из компонент ci
при i+.
Рассмотрим величины
i = ci
- 0T ai ,
1 i n , (4)
которые носят название относительных оценок [6]. Согласно (3)
, i = 0 при
i+ .
Достаточным условием оптимальности текущего допустимого базисного решения
является выполнение неравенства [3]
min
i
0 ,
(5)
где множество 0 определено
в (2). Условие (5) является необходимым для невырожденного допустимого
базисного решения [6]. Найдем необходимые и достаточные условия
оптимальности вырожденного допустимого базисного решения.
2. Основные теоремы. Рассмотрим вспомогательную задачу линейного программирования
{
i yi :
yi vi = d , yi0 ,
i0 } ,
(6)
где d (V)-- любой
вектор, для которого совместны ограничения в (6).
Замечание 1. Задача (6) имеет матрицу условий размером ( m - k ) x ( n - k ), так как в базисе, определяемом столбцами матрицы V (далее -- в базисе V), ее условия-равенства эквивалентны системе m - k скалярных уравнений. Ниже мы не будем пока переходить к базису V, чтобы не вводить новых координат.
Можно доказать следующую теорему 1 (критерий оптимальности). Допустимое базисное решение задачи (1) оптимально тогда и только тогда, когда конечно значение задачи (6).
Алгоритм уменьшения целевой функции в задаче (1) при неоптимальном допустимом базисном решении базируется на следующей теореме.
Теорема 2 (критерий неоптимальности). Текущее допустимое
базисное решение задачи (1) неоптимально тогда и только тогда, когда
найдется множество индексов S 
0 , | S | m - k + 1
такое, что выполняются следующие условия:
1. векторы ai , i S линейно
независимы;
2. существуют числа 
i > 0, i
S такие, что
i
vi = 0 , i
i < 0 ;
i > 0 , i S . (7)
Доказательство. Пусть 
B
0 -- множество,
состоящее из m - k индексов, соответствующих векторам текущего
базиса в задаче (6). По теореме 1 текущее допустимое базисное решение
задачи (1) неоптимально тогда и только тогда, когда целевая функция
задачи (6) неограничена на множестве своих допустимых решений. Согласно
теории линейного программирования [3], это эквивалентно тому, что на
некотором шаге симплексного метода решения задачи (6) находится небазисный
вектор vp , p
0 \ B
, такой, что коэффициенты 
j его
разложения по базису, определяемые однозначно из уравнений
vp=
j vj ,
(8)
удовлетворяют неравенствам
p
- j
j < 0 ; j
0 , j
B . (9)
Рассмотрим множество
S = { j : j
B , j < 0 } 
{ p },
составленное из индексов ненулевых компонент j
и индекса p. Примем
i =
- t i , i S ,
i 
p , p = t ,
(10)
где t -- произвольное положительное число. Тогда соотношения (8) и (9) записываются в виде (7). Теорема доказана.
Таким образом, при неоптимальности текущего допустимого базисного решения задачи (1) множество S , указанное в теореме 2, находится в процессе решения задачи (6). Покажем, как уменьшить при этом целевую функцию задачи (1). В соответствии с равенством в (7) и разложением (2) имеем
j vj = 0
j
aj(U)
j aj =
i
ai , (11)
где коэффициенты i определяются
однозначно, так как система уравнений в (11) относительно
i представляет собой k
линейно независимых уравнений с k неизвестными.
Замечание 2. Линейная комбинация
j aj является аналогом
вектора, вводимого в базис в обычном симплекс-методе.
Введем множество индексов
= { j : j
+ , j
> 0 } (12)
и, если оно не пусто, положим
 = { |
xjj |
: j } > 0 |
(13) |
Множество индексов, на котором достигается минимум в (13), обозначим
| R = {r : | xrr |
= }.
(14) |
Можно доказать следующую теорему 3 (метод уменьшения целевой функции). Пусть текущее
допустимое базисное решение x неоптимально. Тогда
1. если 
,
то множество индексов
+ =
+ S \ R
(15)
соответствует строгому базису задачи (1), отвечающему меньшему значению целевой функции
ci
xi = ci
xi +
j j ,
(16)
где числа j , j
S определены по формуле (10), а новые
переменные, соответствующие строгому базису, определяемому множеством
+ , вычисляются по
формулам
| xi= | |
xi - i , i , |
i+
\ R i S |
, | (17) |
2. если = , то значение
задачи (1) равно -
.
Замечание 3. Как следует из (11) и (13), величина
i при любых
i не зависит от множителя t в (10). Это означает, что его
выбор не влияет и на значение целевой функции в (16).
Замечание 4. Число столбцов в строгом базисе, определяемом
множеством + согласно
(15), равно
k = | + |
= k + | S | - | R |.
При этом, очевидно, | S | 2, т.е. в базис
всегда вводится не менее двух векторов. Рассмотрим некоторые частные случаи.
a) В общем случае, учитывая, что | R | k ,
| S | 2, получаем
k k + 2 - k = 2 ,
т.е. степень вырождения (m-1) после решения вспомогательной задачи не может реализоваться.
б) В случае, когда множество R содержит только один индекс, из базиса выводится один вектор. Тогда
k k + | S | - 1
k + 1,
т.е. при | R | = 1 степень вырождения уменьшается.
в) Пусть | S | = m - k + 1, т.е. в разложении (2) нет нулевых коэффициентов. Тогда k = m, и текущий базис становится невырожденным.
3. Описание алгоритма. Алгоритм метода, обоснованного теоремами 1 -- 3, зависит от выбора матрицы V, дополняющей матрицу строгого базиса U до базисной матрицы B. Приведем этот алгоритм сначала для случая, когда B найдена в результате обычных симплексных итераций, приведших к текущему вырожденному решению. Другой способ выбора базисной матрицы описан ниже в замечании 6.
Пусть giIRk,
hiIRm-k
-- векторы коэффициентов разложения вектора ai , не
входящего в строгий базис U, по столбцам матриц U и V
соответственно, т.е. ui =U gi , vi =
V hi , или
 |
gi hi |
 |
= B-1ai , | i
0 . |
(18) |
Тогда из соотношений (11) следует
=
j
gj ,
(19)
а задача (6) может быть записана в эквивалентной форме
{
i yi :
yi hi = f , yi 0 ,
i
0 } , (20)
где f = B-1 d . Задача (20) содержит m - k ограничений, причем столбцы hi , соответствующие столбцам матрицы V , образуют, очевидно, единичную матрицу. Это позволяет сразу указать начальное допустимое базисное решение задачи (20), а в качестве компонент вектора f выбрать произвольные положительные числа.
В соответствии с изложенным выше алгоритм борьбы с вырожденностью может быть записан в следующем виде.
1. Из условия (18) определяются векторы gi и
hi , i
0 .
2. Строится вектор f , компоненты которого -- положительные случайные числа.
3. Решается задача (20), в которой начальный базис -- единичная матрица.
4. Если задача (20) имеет конечное значение, то текущий базис B задачи (1) оптимален. Иначе -- переход на шаг 5.
5. Из условия (10) при t = 1 определяются коэффициенты
j , j S .
6. По формуле (19) находится вектор .
7. По формуле (13) определяется параметр .
8. В соответствии с условием теоремы 3 строится множество
+ индексов,
соответствующих новому строгому базису задачи (1).
Замечание 5. При случайном выборе компонент вектора f задача (20) с вероятностью 1 невырождена. В случае появления вырождения достаточно лишь заново задать компоненты вектора f в соответствии с шагом 2 алгоритма и пересчитать значения базисных переменных, после чего продолжить решение.
Замечание 6. Возможны и другие способы выбора базисной матрицы B. Например, можно принять
| U = | |
U1 U2 |
|
, | B = | |
U1 0 U2 Im-k |
|
, |
где U1 -- невырожденная k x k - матрица, Im-k -- единичная матрица порядка m - k . Тогда
| B-1 = | |
U1-1 -U2 U1-1 |
0 Im-k |
|
, | gi = U1-1 ui , | hi = vi - U2 U1-1 ui , | i0 |
, |
т.е. в данном случае достаточно обращать матрицу размерности k x k. Поэтому такой выбор базисной матрицы B может оказаться полезным при k << m.
4. Эквивалентный критерий оптимальности и дополнения к алгоритму.
Теорема 3 является обобщением обычного симплекс-метода. При этом роль
относительной оценки i < 0 играет
величина (см. теорему 3)
D() =
i
i ,
которую естественно называть обобщенной относительной оценкой.
Рассмотрим минимальную из величин D()
при условии нормировки на вектор , заданный в (10)
с точностью до множителя:
D* =
{ D() :
i
vi = 0,
i =1 , i
0 , i S } .
Последняя задача оптимизации фактически является задачей линейного программирования, которая в обычном виде имеет вид
D* =
{
ii :
i
vi = 0,
i =1 , i
0 , i
0 } . (21)
Величина D* является аналогом минимальной относительной
оценки min =
i
в обычном симплекс-методе, и для невырожденного допустимого базисного
решения, очевидно, совпадает с ней. Приведем необходимый для дальнейшего
другой критерий оптимальности задачи (1).
Теорема 4.
1. Текущее допустимое базисное решение оптимально
тогда и только тогда, когда D* 0 или
условия в (21) несовместны. В последнем случае текущее допустимое базисное
решение является единственным допустимым решением задачи (1).
2. Если D* < 0 , то множество S , введенное в теореме 2,
находится из решения задачи (21) и соответствует любому текущему базису
этой задачи, для которого D() < 0. При этом целевая
функция задачи (1) уменьшается по формуле (16).
Доказательство легко проводится от противного с использованием теорем 2 и 3.
Непосредственно из теоремы (2) вытекает следующая лемма.
Лемма 1. (связь между задачами (21) и (6)) Если задача (6) не
имеет конечного значения, то множество B
{ p } определяет допустимый базис задачи (6), соответствующий вектору , ненулевые компоненты
которого определяются из соотношения (10) при t = (1 -
i
)-1 . При этом D() < 0 .
На практике при решении задачи (1) обычным симплекс-методом достаточное
условие оптимальности (5) заменяется условием
-оптимальности
min
- ,
(22)
где - заданное положительное число. Такая замена
может быть обоснована следующей леммой.
Лемма 2. Если оптимальное решение задачи (1) удовлетворяет условию
xi*
M , (23)
где M > 0 - заданное число, то для любого неоптимального допустимого базисного решения x справедливо неравенство
cT x*
cT x + M
min .
(24)
Доказательство леммы 2 приведено в Приложении.
Аналогично доказывается, что в вырожденном случае имеет место более
точная оценка
cT x*
cT x + MD* .
(25)
Поэтому вместо эквивалентного критерия оптимальности D*
0 на практике может быть использован
критерий -оптимальности
D*
-.
(26)
Однако непосредственная проверка этого критерия может быть затруднена в силу вырожденности задачи (21). Поэтому на практике целесообразно использовать следующий результат.
Лемма 3. Значение D* задачи (21) удовлетворяет неравенству
D*
D() +
min ,
где min - минимальная
относительная оценка для текущего базиса в задаче (21).
Доказательство леммы 1 следует из формулы (24), примененной к задаче (21) при M = 1 . Другое доказательство этого утверждения приведено в [9].
С использованием приведенных лемм шаг 5 описанного выше алгоритма может быть модифицирован следующим образом.
5.1 Согласно лемме 1 строится допустимое базисное решение
задачи (21), компоненты которого вычисляются по
формуле
| i = |
1 1 -
i |
|
- i , i
B , 1 , i = p, 0 иначе. |
5.2 Ищется величина D() (согласно лемме 1,
D() < 0).
5.3 Для текущего базиса задачи (21), определяемого множеством индексов
B
{ p } , ищется минимальная относительная оценка
min.
5.4 Если выполняется условие D() +
min -
, то согласно лемме 3 текущий базис задачи (1) является -оптимальным , и вычисления завершаются. В противном случае выполняется один из двух следующих шагов.
5.5 Если min
0 , то решение задачи (21), построенное на шаге 5.1, является
оптимальным. В этом случае числа i
определяются из соотношения (10) при t = (1 -
i
)-1 . Переход к шагу 6.
5.6 Если min < 0 , то решение,
построенное на шаге 5.1, неоптимально. В этом случае можно попытаться
уменьшить найденное значение D() , решая
задачу (21) начиная с базиса, определяемого множеством индексов
B
{ p } . Критерий окончания -- выполнение достаточных условий оптимальности для
задачи (21) или появление при ее решении вырожденных итераций. В результате
определяется вектор , для которого
D( ) D() < 0 .
Тогда в качестве чисел i ,
i S , выбираются ненулевые компоненты
вектора .
Замечание 7. (оценивание целевой функции с положительными
коэффициентами) На практике неравенство (26) часто является лишь косвенным
признаком того, что решение близко к оптимальному, так как не всегда
известно, что справедливо условие (23). Однако в случае ci > 0 ,
1 i n имеет место
непосредственная оценка значения задачи (1) :
| cT x
cT x* |
cT x 1 - Dc* |
, | (27) |
где x -- текущее допустимое базисное решение,
Dc* =
{
ii :
i
hi = 0 , ci
i =1 , i
0 , i
0 } .
Оценка (27) следует из неравенства
cT x*
cT x +
cT x*D* ,
аналогичного неравенству (25). В невырожденном случае оценка (27) переходит в оценку
| cT x* |
cT x 1 - 
( i / ci ) |
, | |
полученную ранее в [11] другим способом.
5. Практические результаты. Рассмотрим использование предложенного алгоритма на примере решения описанной в [7] задачи составления расписания авиаперевозок.
Пусть m -- число заказов на авиаперевозки или заданий, n -- число возможных способов доставки грузов, которые условимся называть режимами работы. Будем предполагать, что в каждом режиме выполняется один или несколько заказов по доставке. Обозначим cj -- суммарные затраты на проведение j-го режима. Введем матрицу A с элементами
| aij = | |
1 , если в j-м режиме работы выполняется i-е задание, 0 в противном случае. |
Обозначим также через xj число использований при проведении работы j-го режима. Тогда задача минимизации расходов при выполнении заказов будет иметь вид
{ cTx :
xi ai= e ,
xi IN ,
1 i n } . (28)
где e = (1, ... ,1)T
IRm, и первое условие означает, что каждое задание
будет выполнено ровно один раз.
Задача (28) представляет собой задачу целочисленного программирования
и может быть решена при помощи соответствующего метода. При этом на
предварительном этапе требуется решить задачу линейного программирования,
полученную из задачи (28) заменой условия xi
IN на более слабое
xi 0 при
1 i n.
При решении такой задачи с использованием стандартного симплекс-метода
часто наблюдаются большие последовательности вырожденных итераций. Приведем
результаты применения описанного выше алгоритма, а также алгоритма,
описанного в [7].
При численных расчетах решалась задача для n=3135 и m=80. Коэффициенты матрицы A выбирались случайно из множества { 0,1 }. Коэффициенты cj также задавались случайно в интервале от 1000 до 20000. Результаты расчетов приведены в таблицах.
| Номер итерации | Размер строгого базиса |
Значение целевой функции | Время расчетов (сек.) |
|---|---|---|---|
| 1 | 80 | 2965.5 | |
| 624-631 | 69 | 1677.7 | 3.61 |
| 674-735 | 41 | 1674.5 | 29.48 |
| 898-927 | 52 | 1673.5 | 30.16 |
| 1043-1069 | 36 | 1673.1 | 13.15 |
| 1130-1153 | 32 | 1672.2 | 11.57 |
| 1243-1279 | 27 | 1671.6 | 18.29 |
| 1323 | 25 | 1670.7 | 0.59 |
Выполнено итераций: всего -- 1324, вырожденных -- 607.
Общее время расчетов -- 13 минут 36.96 секунд.
| Номер итерации | Размер строгого базиса |
Размерность подзадачи | Число итераций в подзадаче |
Значение целевой функции | Время расчетов |
|---|---|---|---|---|---|
| 1 | 80 | 2868.3 | |||
| 554 | 77 | 80 X 3135 | 1 | 1681.2 | 0.38 |
| 560 | 65 | 80 X 3135 | 13 | 1672.3 | 5.49 |
| 564 | 54 | 80 X 3135 | 21 | 1671.3 | 8.51 |
| 566 | 41 | 80 X 3135 | 36 | 1670.8 | 17.08 |
| 568 | 26 | 80 X 3135 | 170 | 1670.3 | 72.89 |
| 569 | 23 | 80 X 3135 | 150 | 1670.2 | 97.49 |
Выполнено итераций: всего -- 569, вырожденных -- 36.
Общее время расчетов -- 9 минут 12.61 секунд.
| Номер итерации | Размер строгого базиса |
Размерность подзадачи | Число итераций в подзадаче |
Значение целевой функции | Время расчетов (сек.) |
|---|---|---|---|---|---|
| 1 | 80 | 2887.3 | |||
| 554 | 78 | 2 X 3057 | 2 | 1675.1 | 1.10 |
| 560 | 65 | 15 X 3070 | 7 | 1672.1 | 5.22 |
| 564 | 41 | 39 X 3094 | 82 | 1671.1 | 24.34 |
| 565 | 32 | 48 X 3103 | 240 | 1670.8 | 59.59 |
| 566 | 28 | 52 X 3107 | 100 | 1670.6 | 34.88 |
Выполнено итераций: всего -- 566, вырожденных -- 22.
Общее время расчетов -- 7 минут 16.87 секунд.
Использование стандартного алгоритма приводит к последовательностям вырожденных итераций, число которых достигает 50. Целевая функция при этом практически остается постоянной. Всего было сделано 1324 итерации, и затем выполнились достаточные условия оптимальности.
Использование нового алгоритма и алгоритма Вольфа позволило сократить число итераций при решении более чем вдвое. При использовании алгоритма Вольфа на каждой вырожденной итерации решалась вспомогательная задача той же размерности, что и исходная, при использовании нового алгоритма размерность вспомогательной задачи зависела от числа ненулевых базисных переменных. Это позволило сократить время расчетов по сравнению с алгоритмом Вольфа примерно на 20 процентов.
Было проведено большое число экспериментов, при которых решались задачи различной размерности с различными условиями. В большом количестве случаев новый алгоритм оказывался более эффективным и позволял сократить время вычислений порой вдвое. В то же время в ряде случаев более эффективным оказывался алгоритм Вольфа. Это не позволяет сделать однозначный вывод о преимуществе нового метода борьбы с вырождением перед методом Вольфа, но говорит о его высокой эффективности и целесообразности проведения дальнейших исследований в этой области.