Реоптимизация кратчайших путей приращений при решении асимметричных задач коммивояжера Текст научной статьи по специальности «Математика»

Доклады БГУИР

2011 № 3 (57)

УДК 656.2-50: 519.8

РЕОПТИМИЗАЦИЯ КРАТЧАЙШИХ ПУТЕЙ ПРИРАЩЕНИЙ ПРИ РЕШЕНИИ АСИММЕТРИЧНЫХ ЗАДАЧ КОММИВОЯЖЕРА

М П. РЕВОТЮК, П.М. БАТУРА, А.М. ПОЛОНЕВИЧ

Белорусский государственный университет информатики и радиоэлектроники П. Бровки, 6, Минск, 220013, Беларусь

Поступила в редакцию 1 апреля 2011

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

Ключевые слова: задача коммивояжера, разностная схема, вычислительная сложность.

Постановка задачи

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

7_„ = min

Ё Ё cvxv

i=1 j=1

Zx.. = Ё x.. = 1; x.. > 0, i, j = 1,n;

ij ij ij

i =1 j =1

ut - v + nxj < n - 1, i = 2, n, j = 2,n, i ф j

(1)

Для решения так называемых асимметричных задач вида (1), когда су Ф ср, 7, ] е 1, п ,

наиболее эффективным из точных методов считается метод ветвей и границ [1]. Схема алгоритма метода ветвей и границ может использовать разные способы порождения дерева вариантов. Исторически первый способ, предложенный Литтлом, использует операцию приведения матриц. Современный подход базируется на решении линейных задач о назначении (ЛЗН), анализе получающихся замкнутых циклов и, если таких циклов более одного, последующем переборе вариантов разрыва циклов. Установлено, что размер дерева вариантов в этом случае оказывается меньше. Рекурсия обхода дерева ЛЗН строится на матрице расстояний, где разрывы циклов задаются назначением бесконечных значений длин запрещаемых дуг [2].

Основной проблемой использования метода ветвей и границ, как и других процедур исчерпывающего поиска, является повышенные требования к объему памяти. В частности, прямолинейная реализация метода ветвей и границ для асимметричной задачи коммивояжера размерностью п характеризуется потребностью в памяти порядка п ■ (п2 + 2п). Использование разностной схемы представления состояний процесса ветвления позволяет получить зависимость 0(п2). Например, реализация такой схемы для случая ветвления на двоичных деревьях методом Литтла характеризуется потребностью в памяти п2 + 2п ■ (п + 3) [3].

Предмет рассмотрения - способ построения эффективной как по памяти, так и быстродействию разностной схемы ветвления на множестве ЛЗН. Алгоритм решения ЛЗН будет явно

учитывать возможность частичного наследования результатов решения порождающей задачи. Потенциальный эффект реализации такого приема объясняется тем, что вычислительная сложность решения независимых ЛЗН - 0(п3), а пересчета после изменения строки матрицы ЛЗН - 0(п2) [3]. Для простоты изложения далее будут опущены приемы сокращения количества ветвей дерева вариантов [1, 2, 5], не влияющие на общность предлагаемых решений.

Алгоритм ветвления и его переменные состояния

При прямолинейной реализации метода ветвей и границ время решения (1) в основном определяется временем решения в каждом узле дерева ЛЗН фиксированной размерности

,=i j=i

mm w 7 cj • xy

Z Xj =Z хг] = 1; хг] > 0; i, j = 1, n\. (2)

i=i j=i I

С целью определения переменных состояния процесса поиска оптимального решения (1) рассмотрим процесс ветвления более детально. Ветвление дерева вариантов задачи коммивояжера удобно проводить, используя стратегию DFS (Depth First Search) на векторе решения

R =

{i\xp = 1, j = 1,n} = {r(j), j = 1,n} . (3)

Поиск оптимального решения начинается с создания глобальных объектов - вектора лучшего текущего назначения Rmin и его оценки Ymin с начальными значениями Rmin = 0 и Ymin = да . Далее следует вызов процедуры анализа корневого узла дерева для матрицы C .

Процедура анализа узла на уровне k , получив матрицу Ck, включает следующие шаги.

n

Шаг 1. Решение ЛЗН Ck ^ Rh и оценка целевой функции Y = Z cr(j) j .

j=1

Шаг 2. Если Ymin < Y , то выход.

Шаг 3. Поиск в решении цикла обхода вершин минимальной длины гк с Rk .

Шаг 4. Если rk с Rk , т.е. цикл не гамильтонов, то переход к шагу 6.

Шаг 5. Сохранение результата - Ymin = Y , Rmin = Rk, а затем - выход.

Шаг 6. Для всех вершин цикла поочередно создать копию матрицы Ck ^ Ck+1, установить запрет посещения других вершин цикла и выполнить процедуру анализа узла, соответст-

Ck+1 i .

В настоящее время лучшим по времени решения ЛЗН (2) является метод кратчайшего пути приращений SAP (Shortest Augmenting Path) [4]. Развивая идею классического венгерского метода, метод SAP базируется на двойственной задаче линейного программирования

{n n _I

Z u +Z vj cv - u - vj > 0 i, j=1 n \.

В известных реализациях метода SAP чаще всего внимание уделяется эффективности решения независимых ЛЗН вида (2) c квадратной матрицей. Однако алгоритм метода SAP пригоден после модификации для реоптимизации задач с изменяющимися прямоугольными матрицами, требуя сохранения лишь вектора потенциалов.

Алгоритм решения ЛЗН методом SAP включает следующие шаги:

1) формирование начального назначения и потенциалов различающихся столбцов по результатам поиска минимальных элементов в строках;

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

Первый шаг можно опустить, тогда ЛЗН будет фактически решена инкрементной версией алгоритма венгерского метода [3]. Потенциалы столбцов при этом должны быть нулевыми.

Изменение элементов матрицы ЛЗН влечет необходимость пересмотра результата оптимизации, если только меняется позиция нулевого элемента после операции приведения. Действительно, увеличение элемента с^ ^ с* матрицы задачи (2), когда = 0, не меняет решения для любых I, j е 1, п . Аналогично, уменьшение элемента с^ ^ с*, когда = 1, не нарушает соотношения с у = ui + vj■ в задаче (2). В общем случае, для решения ЛЗН с измененной матрицей требуется повтор итераций назначения лишь для строк

М * =

{ ^ (с* > су ) Л (Ху = 0) V (с* < су ) Л (ХР = 1), j = 1 п) . (4)

При этом достигается снижение вычислительной сложности пересчета задачи (1) на величину (п - М* ' 0(п2). Из выражения (4) следует, что релаксация на множестве ЛЗН при решении задачи коммивояжера методом ветвей и границ соответствует условию |М | = 1. Последнее объясняет сокращение времени решения задачи коммивояжера приблизительно в п -1 раз.

Из (5) следует, что изменение элементов матриц в любом узле дерева вариантов производится только путем увеличения значения. Учитывая, что потенциал измененной строки будет уточняться на итерации пересчета, предварительное уточнение потенциалов не требуется.

Вектор ^ пригоден как для выявления циклов, не являющихся гамильтоновыми, так и исключает необходимость хранения решения в матричном виде.

Можно заметить, что если k - некоторая вершина цикла в решении задачи (2), то последовательность г1 (е) = {г(0) = k, (г^) = ^ |г(0 Ф ^|с ^ только тогда соответствует га-

мильтонову циклу, когда условием остановки является г(п -1) = k, k е 1,п. Если цикл не га-мильтонов, т.е. г1 (е) с R1, то необходимо породить множество задач уровня 1 +1. Для этого следует указать цикл минимальной длины, выбрав вершину входа в цикл (для наследования)

^ = а^тт {|г1 (k)|, k е 1, п|.

Правило порождения ЛЗН тривиально - для каждой вершины обнаруженного цикла необходимо запретить посещение других вершин этого цикла. При этом матрица очередной ЛЗН отличается от предыдущей одной строкой, в которой часть элементов заменяются значением, не меньшим значения

стах = таХ {сР : i Ф j, i = 1 n, j = 1, п| . (5)

Построчное изменение матриц проводится по следующему закону: сГ = 4 + стах ^ е Г1 (е ) Л j е г1 (^ )), и j = Щ . (6)

Так как изменения матрицы С выполняются построчно, то достаточно сохранить изменяемые элементы в буфере размером п . Однако запрет на использование элемента матрицы можно установить смещением его текущего значения на величину стах . В этом случае дополнительный буфер для сохранения элементов строки не требуется. Надежность реализации такого способа очевидна, если выполняется условие размещения в машинном слове, выделяемом для хранения элементов матрицы, значения стах • п .

Фильтрация бесперспективных вариантов

В случае решения задачи коммивояжера естественно использовать взаимосвязь матриц рекурсивно порождаемых ЛЗН. Например, можно учитывать возможности прерывания итераций решения задачи ЛЗН. При решении ЛЗН на основе (3) имеется возможность получения нижних оценок целевой функции. Это позволяет прервать анализ бесперспективного варианта матрицы, используя глобальное значение рекордной оценки среди просмотренных листьев дерева [3]. На основании теории двойственности, нижняя оценка целевой функции Zk на итерации k решения ЛЗН венгерским методом определяется выражением

¿k = &k +Zvk, (7)

i=1 j=1

где uk, vj - значения потенциалов строк и столбцов на этой итерации. Если выясняется, что

Zk > Zrec, где Zrec - значение целевой функции лучшего из просмотренных вариантов, то итерации анализа рассматриваемой ЛЗН можно прекратить. Нижняя оценка целевой функции в этот момент совпадает с оценкой оптимального решения, но выражение (7) эффективнее для вычисления по сравнению с прямолинейным использованием выражения (1).

В алгоритме метода SAP промежуточные итерации проводятся с использованием лишь

потенциалов столбцов vk, j = 1, n . Потенциалы строк можно определить после завершения итераций, используя соотношение ur( j) = cr( j)j - vj, j = 1, n. В этом случае r (j) e 1, n, j = 1, n, что

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

Zk = &(j),j • (r( j) > °).

j=1

Шаблон класса решения задачи

Рассмотренная схема поиска решения задачи коммивояжера представлена полной версией исходного текста шаблона класса на языке С++(рис. 1-3).

template <cla ss cost> struct task { // Дескриптор узла дере ва вариантов

int е, i; // Индексы головного и текущего элемента цикла

vector<int> г; // Вектор решения - назначенные строки для столбцов

ve ct or< со s t> v; Вектор потенциалов столбцов

taste (int : v(n),r(n) {>

template <cla ss cost> class TSP { Класс решения задачи комм ивояжера

int dots,n,* mark,*cols; // Вектор наилучшего текущего реше ния

cost **c, y. с max, inf;

cost values(Int { // Оценка целевой функции текущего н азначения

cost v=0;

for tint j= 0; j<n; j++) v4=c[r[j] ] [j];

return v;

ve ct or< со s t> dtmp;

vector<int> mark,rows,cols,cist,prnt;

coat laps(ta sk<cost> *xn); // Решение ЛЗН методом SAP

void lapi(ta sk<cost> *xn, int e); // Итерация решения ЛЗН

cost lapn(ta sk<cost> *хп, task<cost> +xo); // Реоптимизали я ЛЗН

void spit(ta sk<cost> *xn); // Порождение подзадач ЛЗН

void best(ta sk<cost> *хп, cost x); // Анализ структуры реш ения ЛЗН

Public:

TSP(int N, со st ^C): // размерность и указатель массива исходных данных

n (N) , с (С) , dots (0) ,mark (N) ,dtmp (N) ,rows (N) ,cols (N) ,clst (N) , prnt (N) {

c_max=0;

for (int i= 0; i<n; i++> for (int j = 0; j<n,- j++)

if ((i!=j)s&(c max<c[i][j])) с max=c[i] [j];

y-inf-(++C max)* n;

for (i=Q; i*cn / clst[i]=i, mark [i++] =dots) с [ i] [i] =c_raax;

f void store(c ost x, int *r> { // Сохранение назначения

y=x;

for (int j= 0; j<n; j++) cols[j]=r[j] ;

void operate r() 0 { it Старт решения в корне дерева задач

tas k<cost> аг(п);

best(&ar,laps(Ear));

int *solutio }; n() { return Ecols[0]; } // Выборка вектора ре шения задачи

Рис. 1. Шаблон класса решения задачи коммивояжера

Атомарной задачей, соответствующей узлу дерева, здесь выступает процесс решения ЛЗН (рис. 2) и анализ результата (рис. 3). При этом проверяется наличие гамильтонова цикла и порождаются новые узлы, если цикл не гамильтонов.

После этапа решения ЛЗН в случае необходимости ветвления имеем список порождаемых задач следующего уровня. Список задач полностью представлен элементами вектора решения (3) текущей ЛЗН.

template <class cost> // Решение ЛЗН метопам SAP

inline cost TSP<cost>: : lsps (tssk<cost> *лп) {

vector<int> unchoozen_rows (n}f

cost ж jhj *a, *v=4xn->v[0] ;

int iT j, krrar ^г=£жп->г [0] ;

for (i=0 ; i-in; i++) rows [i]=-l;

for (j=n; —j>=0; r[j]=k, v[j]=x|- { // Инициализация вектор а назначений

for (д=с [k=0] [ j] , 1=1; i<n; i++:<

if ( (h=c[i] [j]J<x) { x=h; k=i; }

if (rows[k]<QJ rows[k]=j; else k=-l;

s for (m=i=0; i<n ; 1++J // Инициализация вектора потенлиалов столбцов

if ( (k=rows [i] )<0J unchoozen rows[m++]=i;

else {

for (jt=inf, a=c[i], j=0; j<n; j++!(

if (t3I=WM((hFH[j]-v[j] }<*}> x=h;

v [ k] -=x ;

} for (k=m; —k>=0; } { // Окончательное фораачрозание ■зазначения строк

lapi ixn r unchoozen rows[k]};

for (x=0r j=0; j<n; // Оценка низшей границы целеаой функции

if ( ( (i=r [j]} 5=0} (x+=c[i] [j] )>у) } return inf;

} return values (r) ;

} template <class cost> // Итерация решетя ЛЗН метопом SAP

inline void TSP<cost>: : lapi { task<costs- ^xn, int s'/ {

cost x rh, fr *a=c [e ] r *4=fixn-Sir [0] ;

int i, j, k,*r=&Jtn->r [0] ;

for (j=0; j<n; prnt [j++]=e} dtnp[ j]=a [j]-v[ j];

for (int head=0r tail=0r last;;) {

if (tail=head) {

last=head; x=dt:rap [cist [tail++] ] ;

for (k=taili k<n; k++} if ((h=dtmp[j=clst[k]] {

if (h<at} { K=h ; tail=head; }

cist [k]=clst [tail]; cist [tail++]=j ;

} for (k=head; k<tail; k++) if (r [j=clst [k] ]<0J goto backtr ack;

} i=r [j=clst [head++]]; a=c[i]; h=a [j]-v[ j] -к;

for Сk=tail; k<n; k++; {

j=clst [k] ;

if ( (f=a[j]-v[j]-h}<dtrap[j]i {

pmt [ j]=i;

if (f=H) {

if (r[j]<0J goto backtrack;

cist[k]=c1st[tail]; cist[tail++]=j;

} dtmp[j]=f; }

i } backtrack: do { // Фиксация кратчайшего пути приращений

i=r [j]=pmt [ j] ; k=j; j=rows[i]; rows [i]=k;

} while(il=e};

for (k=0i k<last; v[j]+=dtrap [j]-«} j=clst [k++] ; }

Рис. 2. Шаблоны основных функций решения задач о назначении

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

Построенный шаблон класса пригоден для непосредственного использования в одно-поточных последовательных вычислениях. При распараллеливании шаблон применим для анализа отдельной ветви дерева вариантов. Основным достоинством рассмотренной здесь разностной схемы реализации ветвления на множестве ЛЗН является экономное использование памяти. Оценка потребности в памяти - п2 + (п -1) • (п +1), где первое слагаемое - исходная матрица задачи (1), второе - память стека вариантов ветвления.

template <class cost> // Анализ структуры решения ЛЗН inline void TSP<cost>::best(task<cost> cost x) {

dots++; xn->e=0; int k=0,m=0,+г=£хп->г[0];

do \ m++; mark[k]=dots; k=r[k]; } while (mark [k]I=dots); if (m<n) { // Поиск иикла минимальной длины for (int 1=1? ; ) { while ( (Kn) && (mark[i]==dots} } i++; if (i==n) break; int k=i, t=0;

do { t++; mark[k]=dots; k=r[k]; } while (mark[k]3=dots); if (m>t) itt=t, xn->e=i;

}

splt(xn); // Переход к ветвлению } eise störe(x,r); // Сохранение лучшего из просмотренных назначений

}

template Cclass cost> // порождение подзадач ЛЗН inline void TSP<cost>::splt(task<cost> + xo) { task<cost> an(n); xo—>i=xo—>e; do \

cost +b=c[xo->i]; int j=xo->r[xo->i];

do { b[j]+=c_max; j=xo->r[j]; } while (j 3 =xo->i) ;

cost x=lapn (&an, xo) ; // Реоптимивадия после изменения строки матрицы if (у>х) best(&an,x); // Анализ структуры решения ЛЗН j=xo->r[xo->i];

do { b [j ] —=с шах; j=xo—>r[j]; } while (j 3 =xo->i) ; } while ((xo —>i=xo—>r[xo->i])!=xo->e);

}

: : < с La 3 3 ссзт> // Еесптимхзацкя ЛЗН inline ссзт ZE?<cc зт> : : lapn (t аз '.-:<cc зт> rxn, таз!-кссзт> 'xc: { £c~ iinT j-0 ; j<n; xn->v[j]=xo->v[]] , j—) roT.vs |xn->:: [ j ] zxc->: [ j ]] =j ;

~_xc->L ] zxn->i | rc "3 [xq->± ] ] =-1; la^i (xn, xc->i :■ ; // : : : строка e д^: = ес : ' : j ra Turn vaIueа(£ xn- : 1 ])r _)_

Рис. 3. Функции организации ветвления задач о назначении

Экспериментальная оценка времени решения задач

В таблице приведены результаты экспериментов по оценке среднего времени решения задач коммивояжера процедурами реализации классического алгоритма венгерского метода [5] и алгоритма кратчайшего пути приращений (рис. 1-3). Результаты фиксировались для трех вариантов применения алгоритмов каждого вида: 1) независимое решение ЛЗН; 2) реоптимиза-ция порождаемых ЛЗН в узлах дерева вариантов; 3) реоптимизация порождаемых ЛЗН с прерыванием решения бесперспективных вариантов. Размеры матриц случайных исходных данных, генерируемых по равномерному закону распределения - п=50.. .150.

Оценка времени решения задач коммивояжера различными процедурами

Размерность задачи (n) Среднее время решения, сек (Celeron 1,7 Ггц, 512 Мбайт)

Венгерский метод Метод SAP

1 2 3 1 2 3

50 0,183 0,028 0,009 0,036 0,011 0,009

60 0,342 0,036 0,012 0,059 0,036 0,013

70 0,576 0,082 0,027 0,126 0,027 0,022

80 0,879 0,122 0,038 0,282 0,047 0,043

90 2,513 0,241 0,081 0,299 0,050 0,046

100 3,234 0,283 0,099 0,481 0,084 0,061

110 9,854 0,344 0,124 1,018 0,188 0,146

120 5,599 0,666 0,238 1,038 0,164 0,112

130 8,048 0,861 0,320 1,545 0,101 0,134

140 14,983 0,607 0,228 2,726 0,281 0,216

150 14,447 0,696 0,260 4,571 0,157 0,293

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

Заключение

Реализация метода ветвей и границ предложенным открытым для расширения шаблоном класса для решения асимметричных задач коммивояжера базируется на наследовании решений порождающих задач о назначении, при котором оценка вариантов порожденных задач проводится методом коррекции дерева кратчайших путей приращений. Это приводит к снижению вычислительной сложности задачи коммивояжера в первом приближении на порядок. Дополнительная память для хранения наследуемых значений потенциалов столбцов и вектора решения не превышает объема O(2n2).

REOPTIMIZATION OF THE SHORTEST AUGMENTING PATHS IN ASYMMETRIC TRAVELING SALESMAN PROBLEM

MP. REVOTJUK, P.M. BATURA, A.M. POLONEVICH

Abstract

The problem of the solution of asymmetric traveling salesman problem, based on branch and bound techniques with linear assignment problems relaxation, is considered. Inheritance of the result's data of previous problems and its reoptimization allows to decreasing time of reception of the new solution on branch's tree path. The algorithm of reoptimization, based on a Shortest Augmenting Path method, is offered.

Литература

1. Miller D, Pekny J. // Science. 1991. Vol. 251. P. 754-761.

2. MahshidA.F., Rosnah M. Y. // European Journal of Scientific Research. 2009. Vol. 29, №3, P. 349-359.

3. РевотюкМ.П., Батура П.М., Полоневич А.М// Докл. БГУИР. 2011. №1, С. 55-62.

4. Jonker R., VolgenantA. // Computing. 1987. Vol. 38. P. 325-340.

5. Zhang W. // Journal of Artificial Intelligence Research. 2004. Vol. 20. P. 471-497.