CC BY
17
Доклады БГУИР
2013 № 5 (75)
УДК 656.2-50: 519.8
БЫСТРАЯ ОЦЕНКА ИНТЕРВАЛОВ УСТОЙЧИВОСТИ РЕШЕНИЯ ЛИНЕЙНЫХ ЗАДАЧ О НАЗНАЧЕНИИ
МП. РЕВОТЮК, М.К. КАРОЛИ, П.М. БАТУРА
Белорусский государственный университет информатики и радиоэлектроники П. Бровки, 6, Минск, 220013, Беларусь
Поступила в редакцию 20 апреля 2013
Предложен эффективный алгоритм оценки интервалов устойчивости решений открытых и закрытых линейных задач о назначении, основанный на итерациях пересмотра результатов оптимизации. Экономный одношаговый переход к ближайшей вершине политопа задачи для всех ребер оптимального совершенного паросочетания позволяет практически на порядок снизить вычислительную сложность оценки устойчивости ее решения.
Ключевые слова: задача о назначении, интервал устойчивости решения, разностная схема, вычислительная сложность.
Постановка задачи
Решение классических закрытых линейных задач о назначении (ЛЗН) [1,2], формально записываемых в виде
I п п
I г=1 j=1
I Ху = 1 = 1; x1j > 0; г, j = 1, п\. (1)
¿=1 у=1 I
обычно есть вектор назначений строк матрицы коэффициентов ее столбцам: Я = {гу = г\х^ = 1, г,] = 1П }. (2)
Во многих случаях требуется оценка устойчивости назначения. При этом для каждого элемента е^ в (1) необходимо вычислить интервал (а^, Ъ), в котором значения таких элементов могут быть изменены без нарушения структуры существующего оптимального назначения (2). Очевидно, что е^ е (а^, Ъ^ ), г, у = 1, п, однако значения границ таких интервалов
не являются непосредственным результатом оптимизации назначения. Симплекс-метод -формальная основа лучших из известных алгоритмов решения ЛЗН, принципиально не гарантирует прохождения всех вершин политопа задачи (1).
Базовый алгоритм решения задачи
Рассматриваемая задача - специальный случай оценки чувствительности решений задач линейного программирования. Алгоритм определения интервалов устойчивости решения закрытых ЛЗН вида (1) предложен в [3] на основе модификации процедуры известного
венгерского метода. При этом для каждого элемента матрицы (е^, г, у = 1, п) интервал
устойчивости определятся путем оценки изменения функции критерия оптимизации после инверсии вхождения такого элемента в оптимальное совершенное паросочетание. Формально
ключевая процедура оценки интервалов устойчивости - реоптимизация решения ЛЗН после целенаправленного возмущения значения каждого элемента матрицы.
В настоящей работе рассмотрен способ улучшения предложенного в [3] алгоритма в направлении снижения его вычислительной сложности и распространения на случай открытых ЛЗН. Для этого предварительно проведем ревизию хорошо известных схем венгерского метода с целью эффективного выполнения этапа реоптимизации решений ЛЗН [1,2].
Открытые (асимметричные) ЛЗН с прямоугольными матрицами
ЧИ <
¿=1 j=1
шщ у у с у • Х у
УХу = УХу =1 Ху - 0г = 1т у = 1п \ (3)
¿=1 у=1 I
часто решают как задачи вида (1) после дополнения матрицы строками или столбцами с нулевыми элементами. Сканирование нулевых строк или столбцов является бесполезным, но необходимым этапом алгоритма венгерского метода [1].
Далее ЛЗН будут рассмотрены в постановке (3), когда т < п . При этом холостые итерации сканирования нулевых строк будут исключены. Случай т > п легко сводится к
решению (2) после транспонирования матрицы (с^, г, у = 1, п) и соответствующей замены
содержательной интерпретации индексов строк и столбцов.
Предварительно отметим, что для открытой ЛЗН вида (3), используемое далее понятие оптимального паросочетания соответствует паросочетанию степени т < п . С целью удобства выделения оптимального паросочетания вектор назначения строк столбцам определим так:
Я = |г] = г • (Гу = 1) + (т +1)• (Гу = 0), г = ^Ху = 1, к = 1т} Гу , у = ^ .
Так как по определению |Я| = п, то такой вектор формируется путем отображения оптимального решения на элементах, первоначально имеющих значения г ■ = т +1, у = 1, п.
Известно, что наиболее эффективные для решения задачи (1) алгоритмы венгерского метода строятся с учетом особенностей двойственной задачи
{т п . ___I
У и +У V- с^ - и - ^ - 0, г = 1, т, у = 1, п>. (4)
|=1 у=1 I
Здесь неизвестными являются потенциалы строк и столбцов. Значения потенциалов особого интереса не представляют, но определяют решение задачи (2). Схемы известных версий алгоритмов венгерского метода решения ЛЗН включают быстрый этап формирования начального назначения строк и итерационного дополнения решения для оставшихся строк. На начальном этапе пытаются выполнить назначение строк, используя операцию приведения матрицы задачи. Приведение состоит в вычитании из элементов столбцов минимальных элементов столбцов. Различные версии известных алгоритмов решения ЛЗН отличаются лишь приемами начального частичного назначения строк [1, 2].
Однако при реализации алгоритма решения ЛЗН с ориентацией на систематическую реоптимизацию решений этап формирования исходного решения не является узким местом. Инициализация векторов потенциалов строк и столбцов любой ЛЗН принципиально может быть реализована итерациями формирования оптимального паросочетания инкрементной версией алгоритма венгерского метода [4]. При этом удается учесть чувствительность целевой функции к изменениям переменных ЛЗН (см. рис. 1).
Предлагаемый шаблон класса решения ЛЗН аккумулирует ряд приемов сокращения холостых шагов при обработке последствий изменения размера и элементов матрицы задачи. Например, можно заметить, что изменение элементов матрицы влечет необходимость пересмотра результата оптимизации, если только меняется позиция нулевого элемента после операции приведения. Действительно, увеличение элемента с^ ^ с^ + 0 матрицы задачи (2),
когда Ху = 0 , не меняет решения для любых i е 1, m, j е 1, и . Аналогично, уменьшение
элемента су ^ с у — 0, когда Ху = 1, не нарушает соотношения с у = u + v. в задаче (2).
В других случаях требуется повтор итераций назначения строк. Любое изменение элементов матрицы должно отражаться значениями потенциалов. Если изменены элементы
строки i, то ее потенциал должен соответственно меняться: u = min {су — v ., j = 1, и} . В случае
изменения элементов столбца j его потенциал также должен меняться:
v = min {су — u, i = 1, m} .
template <class T> struct item { T a,b; }; // Интервал устойчивости
template <class T> struct to_min: public item<T> { // Класс стратегии поиска минимума static T inf() { return numeric_limits<T>::max()/2; } static bool lt(T &a, const T b) { if (a<=b) return false; a=b;
return true; }
void hide(T x)
void show(T x) };
a=-inf(); b=x; } a=x, b=inf(); }
template <typename E> struct Solution { // Результат решения ЛЗН vector<E> u,v; // вектор потенциалов строк и столбцов vector<int> q,r; // вектор текущего решения Solution(int M,int N): u(M),v(N),q(M),r(N) {}
void operator=(Solution &obj) { u=obj.u; v=obj.v; q=obj.q; r=obj.r; } };
template <typename E, template <typename T> class D> class LAP: public Solution<E> { // Решение и реоптимизация ЛЗН vector<int> prec, mark; vector<E> s; E maxs;
int stp; // номер итерации решения задачи vector<E> b; // локальная копия матрицы стоимости void lap(int i); public:
int m, n; // размерность задачи
vector<E *> c; // указатели строк матрицы стоимости void lap();
E setx(int i,int j, E x) { // изменение элемента матрицы
E cij=c[i][j]; c[i][j]=x; r[q[i]]=n;
lap(i); // назначение модифицированной строки
return cij; }
LAP(int M,int N,E *C): Solution(N,N),m(M),n(N),c(N), prec(N),s(N),mark(N), stp(N*N*N),b(N*M), maxs(D<E>::inf()) { for (int i=0, k=0; i<m; i++) { c[i]=&b[i*n];
for (int j=0; j<n; j++) c[i][j]=C[k++]; }
lap(); }
}
Рис. 1. Шаблон базового класса решения задач о назначении
Выявление отмеченных и других особых ситуаций предлагается совместить с расширением назначения и повтором назначения модифицированных строк (рис. 1).
Очевидно, что в случае, когда т < п, требуется организация итерации по строкам, для
которых столбцы не назначены. Список таких строк можно получить проверкой условия Г < т, у = 1, п . Нетрудно заметить, что подобная проверка не потребуется, если использовать вектор отображения назначения строк столбцам
О = [я, = Гу = г, у = 1 ^ 1 = 1 т}-
Векторы О и Я удобно формировать на каждой итерации алгоритма венгерского метода, совмещая операции обработки последствий подвергшихся изменению строк и столбцов согласно (5) или (6). Это допускает проведение реоптимизации решения ЛЗН в любой момент синхронно с коррекцией элементов матрицы (рис. 2).
template <typename E, template <typename T> class D> void LAP<E,D>::lap() { // назначение строк
for (int j=0; j<n; j++) { v[j]=0; r[j]=n; mark[j]=prec[j]=stp; } for (int i=0; i<m; i++) lap(i); // назначение реальных строк if (m<n) for (int i=m, j=0; j<n; j++) if (r[j]==n) {
u[i]=-v[j]; q[r[j]=i]=j; i++; // назначение виртуальной строки }
}
template <typename E, template <typename T> class D> void LAP<E,D>::lap(int i) { // назначение строки int j=0, k,row,col; E *b=c[i], h=maxs;
for (k=0; k<n; k++) if (D<E>::lt(h,s[k]=b[k]-v[k])) j=k; if (u[i]=h) for (k=0; k<n; k++) s[k]-=h; for (stp-=n; (row=r[j])<n; ) {
mark[j]=stp; h=maxs; b=c[row]; E uk=u[row]; col=j+stp; for (k=0; k<n; k++) if (mark[k]!=stp) {
if (D<E>::lt(s[k],((row<m)? b[k]:0)-uk-v[k])) prec[k]=col;
if (D<E>::lt(h,s[k])) j=k; }
if (h) for (u[i]+=h, k=0; k<n; k++) if (mark[k]==stp) { u[r[k]]+=h; v[k]-=h;
} else s[k]-=h; }
while ((k=prec[j]-stp)<n) { q[r[j]=r[k]]=j; j=k; } q[r[j]=i]=j;
}
Рис. 2. Основные функции класса решения задачи о назначении
Оценка интервалов устойчивости решения
Элементы матрицы ЛЗН можно разделить на два вида, отражая их принадлежность текущему оптимальному решению. Элементы такого решения соответствуют назначенным ребрам графа оптимального паросочетания
{/}, j) (r, < m), j = 1, n }= J;xt, yl), i = 1, m}.
Е. = ^
Оставшиеся элементы матрицы дополняют граф задачи (3):
Еи = ((г; j), г =1 m, у = 1,п }\ Ет .
Очевидно, что |ЕОТ | = т, а |Еи | = т • (п — 1), что приходится учитывать при оценке вычислительной сложности алгоритма. Обозначим 8т (Х, у) и 5и (Х, у) - допустимое изменение веса любого ребра х ^ у, где ( х, у) е Ет ^ Еи .
Интервал значений веса назначенного ребра, в котором назначение ребра остается неизменным, для задачи минимизации вида (3) легко записать из элементарных соображений:
Сп(Х,у) = (-»,сху +5т(Х,у)], (Х,у) еЕт; /:1П(Х,у) = [с —5и(х,у), + <ю), (х,у) еЕи.
Последнее означает, что существующий вес можно увеличить для назначенных ребер или уменьшить для остальных ребер без нарушения структуры текущего решения задачи минимизации.
При решении задачи максимизации границы интервалов зеркально отражаются относительно начала координат:
С*(х,У) = К-5т(х,у), + да),(х,у) еЕт ; Сх(х,у) = (-да, с^ +5М(х,у)],(х,у) еЕи.
Будем далее рассматривать случаи задачи минимизации. В таких случаях говорят, что ребро графа оптимального паросочетания скрыто, если его вес увеличен так, чтобы ребро больше не являлось частью существующего решения. Скрытие формально реализуется назначением веса из интервала значений (схй + 5т (х, у), + да) . Как известно, для элементов
оптимального паросочетания справедливо условие с = иж + V , (х, у) е Ет .
Пусть оценка оптимального решения ЛЗН есть
т п т п
20 =Е• х„ . (5)
¿=1 ]=1 ¿=1 ]=1
Рассмотрим ребра, принадлежащие оптимальному решению Ет . Можно показать, что если ребро х ^ у скрыто, то реоптимизация решения может быть проведена относительно строки х (при условии наличия других связей вершины х). Легко заметить, что процесс реоптимизации, начинающийся в вершине х, завершится в вершине у, потенциал которой тоже не изменится. Меняется только потенциал строки их (рис. 2). Из выражения (5) следует, что изменение оценок решения оценивается выражением
и1 -и0 = 21 -20, х е 1т. (6)
Здесь нулевой верхний индекс использован для пометки исходного, а единичный -нового решения. Далее покажем, что искомое значение 5т (х, у) = 21 — 20 . Очевидно, что назначение нового веса с ^ с + Ът (х, у) приведет к изменению оценки решения - 21 ^ 20.
Новое решение без ребра х ^ у имеет оценку = 20 +5т (х, у). Если после этого вес ребра х ^ у уменьшить так, что 5т (х, у) ^ 5т (х, у) — 0 , то последствия такого шага очевидны -
условие 21 > 20 останется справедливым, но ребро станет частью исходного оптимального паросочетания. Последнее означает, что с учетом (5)
С, (х,у) = (-да,с^ + 21 -20] = (-да,с^ + и\ -и0], (х,у) еЕт . (7)
Отсюда следует алгоритм построения интервала безопасного изменения значений веса ребра. Установим гарантированно приводящее к скрытию ребра значение с = +да. Фиксируя
в соответствии с (6) потенциалы строки х до и после проведения реоптимизации, получим 5т (х, у) = и\ - и° . Использование разности значений потенциалов исключает необходимость наивного прямолинейного вычисления оценок решений задачи (1), требующего т + п шагов.
Следуя [3], рассмотрим ребра, не принадлежащие оптимальному решению, когда с > иж + V , (х, у) е Ем . Если вес таких ребер менять в интервале (их + V , + да), то структура
решения остается неизменной.
Обозначим далее ем (х, у) = с -5м (х, у) . Будем считать, что ребро графа паросочетания
открыто, если его вес уменьшен так, что это ребро становится частью нового оптимального решения. Открытие ребра наступает в случае, когда с ^ еи (х, у) .
С целью определения еи (х, у), (х, у) е Ем, построим вспомогательный граф ^ , образуемый из графа оптимального паросочетания путем удаления всех дуг, инцидентных вершинам х и у. В таком графе будет (т-2) ребра графа оптимального паросочетания. Если
далее выполнить уменьшение веса с (х, у), то получим для графа ^, дополненного
ребром (х, у), оптимальное паросочетание степени т с оценкой 2а. Если для рассматриваемого ребра выполняется условие сху < 20 — 2а, то оценка решения, содержащее такое ребро, будет удовлетворять условию с + 2а < 20. Это противоречит предположению о том, что исходное оптимальное паросочетание было совершенным. Таким образом, значение еи (х, у) = 20 — 2а есть нижняя граница интервала устойчивости, поэтому
Сп(х, у) = [20 — 2а, + да), (х,у) еЕи. (8)
Схема ускорения процедуры оценки устойчивости
Однако реализация алгоритма построения интервала значений веса ребер на основе (8) не является эффективной. Для каждого из т • (п — 1) ребер придется решать ЛЗН, размер которых (т — 1)(п — 1). Учитывая дискретный характер процесса перемещения по вершинам симплекса при реализации алгоритма решения ЛЗН, предлагается вместо выражения (8) воспользоваться выражением (7), инвертируя направление шагов процесса построения интервала. Конечная граница интервала (х, у) после этого станет начальной границей интервала (х, у). Нулевой шаг в (7) становится решением ЛЗН для гарантированно приводящего к открытию ребра значения с = —да, а единичный шаг соответствует решению ЛЗН с исходной матрицей. В результате получаем
С, (х,У) = (—» + 20 — 2\+да] = (—да + < — и1х,+да], (х,у) е Еи. (9)
Таким образом, оценка выражения (9) требует лишь одношаговой реоптимизации исходной задачи (см. рис. 3). Отметим, что интервал бесконечных значений (—да, + да) для удобства программной реализации выражений (7) и (9) может безопасно отображаться на интервал (стп — 1, стах +1), где стп = , г = 1 т у =1 п} , стах = тах{с#, г =1, m, у = 1 п} .
template <class E, template <class T> class D>
class ILAP: public LAP<E,D> { // Оценка интервалов устойчивости vector<D<E>> ab; Solution<E> os;
void back(int i,int j, E x) { c[i][j]=x; *(Solution *)this=os; } void hide(int i,int j) { // окрытие ребра оптимального паросочетания E cij=setx(i,j,maxs); ab[i*n+j].hide(cij+u[i]-os.u[i]);
back(i,j,cij);
}
void show(int i,int j) { // открытие ребра оптимального паросочетания E cij=setx(i,j,-maxs); ab[i*n+j].show(-maxs+os.u[i]-u[i]);
back(i,j,cij);
}
public:
ILAP(int M,int N,E *C):LAP(M,N,C),ab(m*n),os(n,n) { os=*this;
for (int i=0; i<m; i++) { // построение матрицы интервалов int k=os.q[i];
for (int j=0; j<k; j++) show(i,j); hide(i,k);
for (int j=k+1; j<n; j++) show(i,j); }
} };
Рис. 3. Шаблон класса оценки интервалов устойчивости решения задач о назначении
Пример решения задачи о назначении с оценкой его устойчивости
Исходные данные и результаты решения примера задачи о назначении с оценкой интервалов его устойчивости представлены в таблице. Элементы таблицы - тройки
(ау, с^, Ь ), где интервал устойчивости задан парами (а у, Ъ ), i = 1, m, ] = 1, п . Оптимальное
назначение строк столбцам соответствует выделенным подчеркиванием элементам таблицы.
Исходные данные и результат решения открытой задачи о назначении (т = 4, п = 6)
1 2 3 4 5 6
1 33, 87, +œ 46, 62, +œ 57, 64, +œ 39, 43, +œ 57, 72, +œ -œ, 57, 61
2 -œ, 52, 76 41, 68, +œ 52, 94, +œ 36, 63, +œ 52, 76, +œ 50, 83, +œ
3 32, 77, +œ -œ, 48, 52 55, 64, +œ 34, 54, +œ 55, 59, +œ 55, 68, +œ
4 32, 58, +œ 45, 49, +œ 56, 71, +œ -œ, 42, 46 56, 66, +œ 56, 89, +œ
Заключение
Таким образом, определение интервалов устойчивости открытых и закрытых ЛЗН может проводиться посредством реоптимизации текущего оптимального решения, если инвертировать принадлежность дуг графа задачи соответствующему совершенному паросочетанию и учесть эту принадлежность направлением нумерации состояний. Инверсия реализуется заменой веса дуги на максимальное значение весов для дуг оптимального паросочетания или на минимальное значение для остальных дуг. Время получения оценок устойчивости на основе разности потенциалов изменяемых строк в первом приближении сокращается не менее, чем в п раз, по сравнению с известной версией [3]. Решение вспомогательных задач о назначении не требуется. Дополнительная память для хранения наследуемых значений потенциалов строк и столбцов не превышает объема 0(т + п) .
QUICK EVALUATION OF THE INTERVAL STABILITY OF THE LINEAR ASSIGNMENT PROBLEM SOLUTIONS
M.P. REVOTJUK, M.K. QARALEH, P.M. BATURA
Abstract
An efficient algorithm for evaluating the interval of solutions stability of opened and closed linear assignment problems based on a review of the results of optimization iterations is proposed. Economical one-step transition to the nearest vertex of the polytope of problem one for all the edges of the optimal perfect matching is almost an order of magnitude reduces the computational complexity of evaluating the stability of the current optimal solution.
Список литературы
1. David W. Pentico. // European Journal of Operational Research. 2007. № 176 (2). P. 774-793.
2. Bijsterbosch J., VolgenantA. // Annals of Operations Research. 2010. № 181 (1). P. 443-462.
3. Lantao Liu, Dylan A. Shell. // The International Journal of Robotics Research. 2011. № 30 (7). P. 936-953.
4. РевотюкМ.П., Батура П.М., Полоневич А.М. // Докл. БГУИР. 2011. № 1 (55). C. 55-62.
