Об одном способе построения начального допустимого базиса в задачах оптимизации Текст научной статьи по специальности «Математика»

ОБ ОДНОМ СПОСОБЕ ПОСТРОЕНИЯ НАЧАЛЬНОГО ДОПУСТИМОГО БАЗИСА В ЗАДАЧАХ ОПТИМИЗАЦИИ

Е. А. Котельников

Институт вычислительной математики и математической геофизики СО РАН,

630090, Новосибирск, Россия

УДК 519.853.32

Предложен алгоритм построения начального допустимого базиса в задачах математического программирования с линейными ограничениями, заданными разреженными матрицами большой размерности. В данном случае тип целевой функции (линейная, квадратичная или нелинейная функция) несуществен.

Ключевые слова: оптимизация, линейные ограничения, начальный базис. In paper the algorithm of construction of initial admissible basis in problems of mathematical programming with the linear restrictions set by rarefied matrixes of the big dimension is offered. The criterion function type can be any. It can be linear function, quadratic or nonlinear function of a general view.

Key words: optimization, linear restrictions, initial basis.

Введение. Пусть система ограничений исходной задачи оптимизации представлена в форме

Ax = b; (1)

а < x < в- (2)

Здесь A — матрица размером m x n (m < n); a, x, в G Rn; b G Rm. Если в исходной задаче среди условий (1) имеются ограничения типа неравенств, то, введя дополнительные неотрицательные переменные, такие ограничения можно преобразовать к равенствам. В расширенной таким образом матрице A дополнительным переменным соответствуют единичные орты (с тем или иным знаком) пространства Rm. Исходная задача может содержать неограниченные (свободные) переменные Xj, что соответствует равенствам aj = — ж, (3j = Наличие таких переменных упрощает построение начального допустимого базиса.

Для построения начального допустимого базиса существует ряд методов. Во многих коммерческих пакетах (см. [1]) в фазе I (фазе поиска начального базиса) решается следующая

m

вспомогательная задача линейного программирования: минимизировать xn+j при усло-

j=i

виях

aT x ± Xn+i = bi, i =1, 2,...,m, a < x < в, Xn+i > 0, i = 1, 2,...,m. (3)

Здесь aT — i-я строка матрицы A.

Работа рекомендована к публикации Программным комитетом VIII Международной азиатской школы-семинара "Проблемы оптимизации сложных систем".

Знак в ограничениях (3) определяется следующим образом. В точке x с координатами

(ai, если ai > —ж,

ßi, если ai = —ж, ßi < +ж, (4)

0 в остальных случаях

вычисляется вектор b = b — Ax. Если для некоторых i bi < 0, то в i-м ограничении (3) должен быть знак "минус", в противном случае — знак "плюс". Это позволяет выбрать в качестве начальной допустимой точки x Е Rn+m вспомогательной задачи точку с координатами xi = xi (i = 1, 2,... ,n); xn+i = \bi \ (i = 1, 2,... , m), в качестве начальной базисной матрицы — диагональную матрицу размером m х m с элементами ±1 на диагонали.

Очевидно, что если допустимое множество исходной задачи не пусто, то в оптималь-

m

ном решении вспомогательной задачи xn+j = 0 и ее оптимальный базис может служить

j=i

начальным для фазы II (фазы решения исходной задачи).

Все модификации данного метода сводятся к попыткам уменьшить число искусственных переменных xn+i путем включения в начальный базис фазы I дополнительных или структурных переменных, на которые накладывается требование допустимости их значений. Будем считать, что дополнительная переменная xj имеет границы допустимости aj = 0, ßj = — ж, а искусственная — aj = ßj = 0.

Алгоритм. В настоящей работе предлагается способ построения начального допустимого базиса, позволяющий использовать на начальном этапе фазы I большее количество структурных и дополнительных переменных, уделяя основное внимание численной сбалансированности начальной базисной матрицы. Из всего набора столбцов матрицы A выбираем (если это возможно) m линейно независимых столбцов и формируем из них базисную матрицу B. Если это не удается сделать, то недостающую часть дополняем единичной подматрицей с соответствующим набором искусственных переменных. Проблема построения начальной базисной матрицы B рассматривается ниже.

Столбцы матрицы A, не вошедшие в матрицу B, объединим в матрицу небазисных столбцов N. Любому небазисному столбцу aj соответствует переменная x^j, значение которой всегда равно одному из ее граничных значений (x^j = a^j или x^j = ßNj). Далее каждой небазисной переменной xNj каким-либо образом (например, по правилу, заданному в (4)) присваивается значение xNj. Тогда, если все значения базисных переменных хв = B-1 (b — NXn) имеют допустимые значения, построение начального допустимого базиса заканчивается. В противном случае переопределим значения базисных переменных следующим образом: если xBi < aBi, то xBi = aBi; если aBi < xBi < ßBi, то xBi = xBi; если xBi > ßBi, то xBi = ßBi. Здесь ав, ßB Е Rm — векторы, составленные из компонент векторов а, ß, которые (компоненты) соответствуют базисным переменным.

Пусть r = xB — XB и b' = Br, тогда верно равенство BXB + NXn + b' = b. Поставим в соответствие столбцу b' искусственную переменную xn+1 и запишем систему ограничений (1) в виде системы

B xb + Nxn + b'xn+i = b,

решением которой является точка xT = (xB, xN ,xn+1) при xn+1 = 1. Допустимое значение переменной xn+1 есть только нулевое значение, поэтому на следующем этапе устраним недопустимость переменной xn+1, используя точку xo в качестве начальной.

Если столбец b' плохо масштабирован по сравнению со столбцами матрицы A, например если \\b'\\^ < min \\ai\\^ или \\b'\\^ > max \\ai\\^(ai — столбец матрицы A), то можно

нормировать вектор Ъ' : Ъ' = Ъ'/q, сбалансировав значения ||Ь' и i = 1, 2,...,n,

где q > 0 — некоторое число. Тогда новое начальное значение переменной xn+1 будет равно q. Переменную xn+1 нельзя оставлять небазисной, так как она не равна своему граничному значению. Поэтому, прежде чем устранить недопустимость переменной xn+1, введем ее в базис, заменив ею одну из базисных переменных. Для этого необходимо заменить вектор Ъ' на один из базисных столбцов, так чтобы не была нарушена невырожденность базисной матрицы и заменяемая базисная переменная Xßi имела значение aßi или ßßi.

Известно, что матрица B останется невырожденной при замене базисного столбца с номером i в базисе на столбец Ъ' если i-я компонента вектора B_1b' не равна нулю. Рассмотрим вектор f = B-1Ъ' = B-1(Br/q) = r/q и определим множество J номеров базисных переменных, таких что для любого j Е J xßj = aßj или xßj = ßßj. Множество J не пусто, так как в противном случае вектор Ъ' был бы равен нулю. Пусть I — множество ведущих

элементов в столбцах множества J, i0 = arg max и j0 Е J — номер столбца с номером

i&I

ведущего элемента i0. Тогда столбец aj0 можно вывести из матрицы B и перевести в разряд небазисных, а столбец Ъ' ввести в базис с номером ведущего элемента i0. Новая базисная

матрица невырождена, так как | = max — Xß^/q.

i

Последний этап фазы I состоит в исключении из базиса искусственной переменной xn+1 при условии допустимости всех остальных базисных переменных. С этой целью используем алгоритм, предложенный в [1]. В данном алгоритме на каждой итерации величина сдвига вдоль выбранного направления вычисляется таким образом, чтобы суммарная величина отклонений базисных переменных от их отрезков допустимости была минимальной. На каждой итерации фазы I коэффициенты целевой функции, подлежащей минимизации, задаются в виде

— 1, если xi < ai, + 1, если xi > ßi, 0 в остальных случаях.

Напомним, что дополнительная переменная xi имеет границы допустимости ai = 0, ßi = +ж, а искусственная — ai = ßi = 0. Коэффициенты ci используются для вычисления относительных оценок dj небазисных столбцов aj в модифицированном симплекс-методе [1]. На основе значений dj принимается решение о целесообразности ввода в базис столбца aj. Если на каком-либо шаге ни один небазисный столбец не может быть введен в базис (выполнен признак оптимальности), то полученное решение является оптимальным. Будем считать, что при xj = 0 j-я свободная переменная (aj = — ж, ßj = +ж) может быть небазисной, а знак относительной оценки dj при выборе j-го столбца для ввода в базис несуществен. В оптимальном решении dj свободной небазисной переменной xj должно быть равно нулю.

Пусть на текущей итерации для ввода в базис выбран небазисный столбец aj. Тогда, как известно, вектор сдвига pß базисных переменных можно вычислить по формуле pß = —B-1aj. Величину сдвига Л вдоль направления pß будем выбирать из условия минимума

суммы отклонений Si (Iß — список базисных переменных) вектора Xß = Xß + Лрв от ieiß

параллелепипеда aß < Xß < ßß, где

aßi — xßi, если xßi < aßi,

xßi — ßßi, если xßi > ßßi, (5)

0 в остальных случаях.

Будем считать, что небазисная переменная хм], которая вводится в базис, не нарушает границы допустимости, т. е. выполняется условие Л < вм] — ам].

Задача минимизации 8г сводится к задаче минимизации кусочно-линейной выпуклой

Шв

функции р(х) одной переменной при условии х > 0. Функция р является суммой кусочно-линейных функций рг, определенных для базисных переменных хг (г £ 1в):

1. Если хвг < авг и рвг < 0, то рг(х) = 0.

2. Если хвг < авг и рвг > 0, то

Рвг(ъ — х), если х < >Уг, Рг(х) = <( рвг(х — Пг), если х > Пг,

0, если тг < х < Пг,

где 7г = (авг — хвг)/pвi; Пг = (ввг — Хвi)/pвi.

3. Если хвг > ввг и рвг > 0, то рг(х) = 0.

4. Если хвг > ввг и рвг < 0, то

—'Рвг(ъ — х), если х <^г,

Рг(х) = { —Рвг(х — пг), если х > пг,

0, если < х < пг,

где ^г = —(хг — ввг)/рвг; Пг = —(хвг — авг)/рвг.

5. Если авг < хвг < ввг, то

I \ \ \рвг\(х — Пг), если х > пг, Рг(х) = < ~ '

[ 0, если х < пг,

где

= ( (хвг — авг)/рвг, если рвг < 0,

Пг \ (ввг — хвг)/рвг, если рвг > 0.

6. Если рвг = 0, то рг(х) = 0.

Упорядочим величины Пг (г = 1, 2,... ,т) по возрастанию их значений, запомним в соответствующем порядке номера переменных и обозначим через аг элементы отсортированного массива. Тогда на каждом отрезке [аг, ог+1] функция р линейна, а в точках аг имеет

изломы: р(0) = ^ рвг — ^ рвг. Здесь Д = {г : хвг < авг}; Ь = {г : хвг > Рвг}. Нетрудно г€/2 г^Ь

проверить, что р(0) равна относительной оценке небазисного столбца, вводимого в базис: ] = ОМ] — сТв В-1а] = св Рв = р(0).

Целью минимизации функции р является нахождение номера переменной, которая покидает список базисных столбцов, и величины сдвига вдоль рв. Для этого последовательно по

г

возрастанию индекса перебираются элементы аг и вычисляются величины тг = р(0)— ^ рвгк,

к=1

где %к — номер базисной переменной, стоящей на к-м месте после сортировки. Перебор осуществляется до тех пор, пока либо величина тг перестанет быть положительной, либо величина аг станет больше вм] — ам]. В первом случае номер переменной %к, при котором тгк < 0, выбирается в качестве номера столбца, выводимого из базисной матрицы, а в качестве искомой величины сдвига Л* вдоль вектора рв выбирается агк. Если после перебора всех аг величина тг остается положительной, то в качестве Л* выбирается вм] — ам] и состав базиса не меняется, а небазисная переменная хм] меняет свое граничное значение.

Фаза I считается законченной, если на некоторой итерации значения всех базисных переменных находятся в пределах своих границ допустимости, что означает принадлежность найденной точки множеству, определенному условиями (1) и (2), или если существуют недопустимые базисные переменные, при этом все относительные оценки небазисных переменных удовлетворяют условию оптимальности, что соответствует отсутствию допустимого решения в исходной задаче.

Вернемся к задаче построения начальной базисной матрицы B. Суть задачи состоит в следующем: из всего набора столбцов матрицы A следует выбрать m линейно независимых, так чтобы степень заполненности матрицы B-1 незначительно превышала степень заполненности матрицы B, а численная устойчивость разложения матрицы B-1 = U-1L-1 была достаточно высокой. Для решения данной задачи используем процедуру перепостроения обратных базисных матриц P3 (preassigned pivot procedure) [2], модифицировав ее с целью увеличения численной устойчивости разложения B-1 = U-1 L-1. Процедура P3 представляет собой эвристический алгоритм приведения с помощью перестановок строк и столбцов матрицы B к виду, по возможности близкому к нижней треугольной матрице. При этом предпринимается попытка минимизировать как число столбцов, имеющих ненулевые элементы над главной диагональю, так и количество ненулей в этих столбцах над главной диагональю. Столбцы, не имеющие ненулевых элементов над главной диагональю, называются нормальными, а имеющие — столбцами-выступами. Если матрица B состоит из нормальных столбцов и столбцов-выступов, причем последние расположены в порядке возрастания высоты столбцов, то новые ненулевые элементы в матрице B-1 могут возникнуть только в столбцах, которые соответствуют столбцам-выступам матрицы B, причем на месте нулевых элементов над главной диагональю. В данном случае k — высота столбца aj, если a^j = 0 и aij = 0 при i < k.

После упорядочивания строк и столбцов ведущие элементы выбираются на главной диагонали сверху вниз и слева направо. Однако при таком правиле выбора ведущих элементов обеспечить численную устойчивость разложения матрицы B-1 = U-1L-1 достаточно трудно. Поэтому при получении очередного ведущего элемента нужно проконтролировать его величину и в случае нарушения допустимого уровня перевести столбец в резерв, чтобы иметь возможность использовать его в более подходящих условиях. Для проверки величины ведущего элемента можно использовать стандартный контроль ведущих элементов на устойчивость и отложить использование столбца, если ведущий элемент gi не удовлетворяет условию \gi\ > umax \gk|, где 0 < u < 1; K — список индексов столбцов, для которых

k£K

еще не назначены ведущие элементы. В конце процесса разложения ведущие элементы в отложенных столбцах определяются путем частичного выбора в столбцах.

Описанная модификация процедуры P3 использовалась в качестве процедуры перепостроения обратной базисной матрицы в первых версиях пакетов программ для решения задач линейного и квадратичного программирования и их целочисленных вариантов [3, 4].

Применение модифицированной процедуры P3 для формирования начальной матрицы B из столбцов матрицы A упрощается за счет того, что на каждом шаге список столбцов-кандидатов на ввод в базис шире текущего списка кандидатов при перепостроении обратной базисной матрицы. При этом в рассматриваемом случае существует больше возможностей при наборе столбцов в так называемый нижний треугольник. Кандидатом на включение в "нижний треугольник" в текущий момент является столбец, у которого для всех номеров ненулевых элементов (кроме одного) уже определены ведущие элементы.

Процессы формирования начальной базисной матрицы B и вычисления вектора Ъ' можно модифицировать, выполняя их параллельно. На каждом шаге по правилам процедуры P3 набираем множество J — список номеров небазисных столбцов-кандидатов на ввод в базис, прошедших стандартный контроль устойчивости. Затем для каждого j Е J вычисляем значение j-й переменной Xj по текущим значениям остальных переменных. В качестве первоначального значения всем небазисным переменным Xj присваиваем значения Xj, определенные в (4). Для ввода в базис выбираем столбец с номером j0 = argmin{Sj||ajгде

Sj — величина отклонения, определенная в (5); вектор Ъ' пересчитывается с учетом 5j0, а в качестве нового значения Xj0 выбирается величина

если Xj0 < aj0, если Xj0 > ßj0, если aj0 < Xj0 < j.

Данная процедура позволяет ожидать, что в конечном счете может быть получено наименьшее возможное значение || Ъ' || .

Последняя модификация процессов построения начальной базисной матрицы B и вычисления вектора b', а также процедура устранения недопустимости, описанная выше, были использованы для одного из двух вариантов построения начального допустимого базиса во всех пакетах программ при решении задач оптимизации с линейными ограничениями [3-5]. Результаты расчетов показывают высокую численную устойчивость вычислительного процесса построения начального базиса.

Список литературы

1. Муртаф Б. Современное линейное программирование. Теория и практика. М.: Мир, 1984.

2. Hellerman E., Rarick D. C. Reinversion with the pressigned pivot procedure // Math. Program. 1971. V. 1, N 2. P. 195-216.

3. Забиняко Г. И. Пакет программ целочисленного линейного программирования // Дискретный анализ и исследование операций. 1999. Сер. 2. Т. 6, № 2. С. 32-41.

4. Забиняко Г. И., Котельников Е. А. Параллельный алгоритм целочисленного квадратичного программирования // Вычисл. технологии. 2004. Т. 9, № 1. С. 34-41.

5. Забиняко Г. И., Котельников Е. А., Рожин В. Е. Программы минимизации нелинейных функций при линейных ограничениях: Отчет / ВЦ СО РАН. № ГР 01.9.30 001317; Инв. № 02.9.70 004793. Новосибирск, 1997.

Котельников Евгений Алексеевич — ст. науч. сотр. Института вычислительной математики и математической геофизики СО РАН; тел. (383) 330-60-66

Дата поступления — 2.04.12

X

Jo

j0 >

Jo >

X

jo