Salesman: Метод ветвей и границ


Дерево перебора

В задаче коммивояжёра перебор возможных путей ведёт к очень "ветвистому" дереву. Например, пусть коммивояжёр начинает движение из города 0 и есть ещё 4-е города. Вторым он может посетить один из городов 1,2,3,4. Если выбран путь 0 → 1, далее можно попасть в 2,3 или 4 и т.д.:

Такой перебор, вообще говоря, задваивает число путей, так как не учитывается возможность движения по замкнутому пути в любую сторону (0,1,2,3,4 эквивалентен 0,4,3,2,1). В результате, вместо (N-1)!/2 вариантов, их просматривается в 2 раза больше. Впрочем, проблема не в этом, а в общей ветвистости дерева, которое при больших N не позволяет использовать поиск в ширину (для N городов на глубине d имеем (N-1)(N-2)...(N-d) узлов). При направленном поиске и не слишком "агрессивной" эвристике движения в глубину, мы также рискуем получить очередь open не помещающуюся в памяти компьютера.


Построение бинарного дерева

Существенно менее ветвистым является бинарное дерево. Пусть на данном этапе известно, что город j ещё не пройден, а в город i коммивояжёр ещё не вошёл. Построим две ветви. На всех потомках левой ветви будем перебирать пути в которых присутствует переход j → i, а на правой, в которых переход j → i запрещен.

Рассмотрим, например, матрицу расстояний для 5 городов, приведенную ниже в первой таблице. Нижняя граница длины оптимального пути для этой таблицы равна 10 + 1 = 11. Для её получения, необходимо вычесть минимальные значение из каждой строки (2+3+3+1+1=10, вторая таблица), а затем из каждого столбца (только 1 в первом столбике, см. вторую таблицу):

0 1 2 3 4
0 - 5 8 9 2
1 5 - 3 3 4
2 8 3 - 7 6
3 9 3 7 - 1
4 2 4 6 1 -
0 1 2 3 4
0 - 3 6 7 0
1 2 - 0 0 1
2 5 0 - 4 3
3 8 2 6 - 0
4 1 3 5 0 -
0 1 2 3 4
0 - 3 6 7 0
1 1 - 0 0 1
2 4 0 - 4 3
3 7 2 6 - 0
4 0 3 5 0 -

Любой ноль в последней таблице является потенциально "хорошим перемещением". Выберем, например, переход 1 → 2. Тогда пространство решений можно разбить на два класса - содержащих переход 1 → 2 (левая ветвь) и не содержащих его (правая ветвь):

0 1 3 4
0 - 3 7 0
2 4 - 4 3
3 7 2 - 0
4 0 3 0 -
0 1 3 4
0 - 1 7 0
2 1 - 1 0
3 7 0 - 0
4 0 1 0 -
0 1 2 3 4
0 - 3 1 7 0
1 1 - - 0 1
2 4 0 - 4 3
3 7 2 1 - 0
4 0 3 0 0 -
0 1 2 3 4
0 - 3 6 7 0
1 1 - - 0 1
2 4 0 - 4 3
3 7 2 6 - 0
4 0 3 5 0 -

Для левой ветви выкинем из таблицы расстояний строку 1 (из этого города уже переходить нельзя) и столбец 2 (в этот город уже зашли). Кроме этого, заблокируем обратный переход 2 → 1, который теперь запрещён. Получится первая таблица выше. В ней снова вычтем минимальные значения из строк и столбцов, увеличив нижнюю границу на 3+2=5 (вторая таблица). Для правой ветви ставим прочерк в переходе 1 → 2 (четвёртая таблица) и, вычитая минимумы, также увеличим нижнюю границу на 5 (третья таблица).


Выбор перехода

Какой лучше выбрать ноль в таблице расстояний для создания очередных двух ветвей? Поиск оптимального решения будет происходить в основном в глубину, сначала по левым ветвям с уменьшающимися матрицами (удлиняющимся путём). Поэтому при выборе перехода j → i берётся такой ноль, для которого правая ветвь (содержащая большую матрицу) получает наибольший прирост в нижней границе расстояния. Для этого временно блокируется по очереди каждый ноль (заменяя его на прочерк) и производится минимизация строк и столбцов.

На левой ветви, после вычёркивания строки и столбца, всегда блокируется обратный переход. При этом необходимо учесть, что по мере спуска по дереву могут возникать несвязные участки пути. Их следует склеивать, как только появляется такая возможность. Например, если были выбран переход 1 → 2, а в одном из дочерних узлов переход 2 → 3, то их необходимо объединить в один участок пути 1 → 2 → 3.

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

Перейдём к деталям реализации алгоритма ветвей и границ на JavaScript.


Узел бинарного дерева Salesway

Метод ветвей и границ в модуле salesman.js использует дополнительный класс Salesway, экземпляры которого являются узлами бинарного дерева:

function Salesway(num, cities)
{
   this.num   = num;                              // размерность матриц расстояний
   this.dists = new Array(num*num);               // оставшиеся варианты расстояний
   this.src   = new Array(num);                   // города из которых можно ещё перейти
   this.des   = new Array(num);                   // города в которые можно ещё перейти
   this.frw   = new Array(cities);                // участки пути в одну сторону
   this.rev   = new Array(cities);                // участки пути в обратную сторону
}
Заголовки строк и столбцов (номера городов) матрицы расстояний для оставшихся городов находятся в массивах src и des соответственно. Для экономии памяти, матрица расстояний dists хранится в одномерном массиве. Он содержит в себе кроме расстояний между оставшимися городами, информацию о заблокированных переходах (прочерки), расстояния которых считаются равными -1. Корневой узел бинарного дерева инициализируется по исходной матрице расстояний между городами dists следующим образом:
Salesway.prototype.init = function (dists)
{
   for(var j=0, k=0; j < dists.length; j++){
      for(var i=0; i < dists.length;  i++)
         this.dists[k++] = dists[j][i];           // квадратная матрица в линейном массиве
      this.src[j] = this.des[j] = j;              // в именах строк и колонок все города
      this.frw[j] = this.rev[j] = -1;             // пока нет участков пути
   }
   this.boundLen = this.value = 0;                // накопленная оценка снизу и её ценность
}

Работу с квадратной матрицей размерности num x num, хранящейся в линейном массиве dists, проиллюстрируем на примере функции, которая находит минимальное значение в i-той строке или i-той колонке (если есть параметр col и он не нулевой):

Salesway.prototype.getMin = function (i, col)
{
   var min = Infinity;                            // минимум разрешённых расстояний
   for(var j=0; j < this.num; j++){
      var m = col? this.dists[j*this.num+i]: this.dists[i*this.num+j];
      if(m===0) return 0;                         // меньше 0 быть не может
      if(m > 0 && m < min) min = m;               // только не отрицательные расстояния
   }
   return min;
}


Получение нижней границы оставшегося пути

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

Salesway.prototype.minimizeRows = function ()
{
   var sum = 0;
   for(var j=0; j < this.num; j++){
      var min = this.getMin(j);                  // ищем минимум в строке
      if(min < Infinity){
         sum += min;
         for(var i=0; i < this.num; i++)         // вычитаем его из элементов строки
            if( this.dists[j*this.num+i] >=0 )
               this.dists[j*this.num+i] -= min;
       }
   }
   return sum;
}
а колонок - функция minimizeCols:
Salesway.prototype.minimizeCols = function ()
{
   var sum = 0;
   for(var i=0; i < this.num; i++){
      var min = this.getMin(i, true);             // ищем минимум в колонке
      if(min < Infinity){
         sum += min;
         for(var j=0; j < this.num; j++)          // вычитаем его из элементов столбика
            if( this.dists[j*this.num+i] >=0 )
               this.dists[j*this.num+i] -= min;
      }
   }
   return sum;
}
Результат их работы используется в функции:
Salesway.prototype.calcBound = function ()
{
   this.boundLen += this.minimizeRows() + this.minimizeCols();
   this.value = this.boundLen*this.num;
   return this.boundLen;
}
Она не только вычисляет нижнюю границу пути boundLen, но и задаёт параметр value по которому будет сортировать узлы приоритетная очередь (см. ниже). Чтобы более активно опускаться вниз дерева, нижняя граница умножается на размерность матрицы (число городов в оставшемся пути).


Выбор очередного перехода

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

Salesway.prototype.selectPair = function ()
{
   var ii=0, jj=1, minI, minJ, max=-1;
   for(var j=0; j < this.num; j++)
      for(var i=0; i < this.num; i++){            // бежим по всем ячейкам матрицы,
         if(this.dists[j*this.num+i]) continue;   // которые равны нулю
         this.dists[j*this.num+i]=-1;             // временно блокируем этот переход
         var mJ = this.getMin(j);                 // получаем минимум в строке j
         var mI = this.getMin(i,true);            // и в колонке i
         this.dists[j*this.num+i]=0;              // убираем блокировку
         if(mI+mJ > max){ max=mI+mJ; ii=i; jj=j; minI=mI; minJ=mJ; }
      }
    return {max:max, minI:minI, minJ:minJ, ii:ii, jj:jj};
}
Функция возвращает объект, свойствами которого является прирост нижней границы max, индексы jj, ii строки и столбца в матрице расстояний, где стоит "лучший" ноль и минимумы minJ и minI в этой строке и столбце (после блокирования нуля). Этот объект res = way.selectPair() используется далее для блокировки перехода:
Salesway.prototype.blockWay = function (res)
{
   this.dists[res.jj*this.num+res.ii] = -1;     // блокируем путь jj-ii
   if(res.minJ)                                 // если минимум в строке не нулевой
      for(var i=0; i < this.num; i++) if(this.dists[res.jj*this.num+i] > 0)
         this.dists[res.jj*this.num+i] -= res.minJ;
   if(res.minI)                                 // если минимум в колонке не нулевой
      for(var j=0; j < this.num; j++) if(this.dists[j*this.num+res.ii] > 0)
         this.dists[j*this.num+res.ii] -= res.minI;
}


Сшивка участков пути

Для сшивки отдельных участков пути, прямые и обратные переходы хранятся в двух массивах: frw и rev. Значение frw[j] означает номер города в который идёт переход из города j (или -1, если перехода ещё нет), rev[i] - означает номер города который предшествует городу i (или -1, если его ещё нет). Например, пусть для восьми городов зафиксированы два участка пути: 0 → 5 → 2 и 3 → 4 → 7 → 6. Тогда массивы имеют следующие элементы:

  j:   0  1  2  3  4  5  6  7            i:     0  1  2  3  4  5  6  7
frw = [ 5,-1,-1, 4, 7, 2,-1, 6]          rev = [-1,-1, 5,-1, 3, 0, 7, 4]
Приведём пример работы с этими массивами в функции, возвращающей массив массивов с участками сформированного пути:
Salesway.prototype.getPath = function ()
{
   var res = [];                                  // результирующий массив массивов

   for(var k=0; k < this.rev.length; k++){        // по массиву с обратными путями
      if(this.rev[k] >= 0 || this.frw[k] < 0)
         continue;                                // ищем начало последовательности

      var ar = [k];
      var n = k;
      while( this.frw[n] >= 0 ){
         n = this.frw[n];
         ar.push(n);
      }
      res.push(ar);                               // добавляем цепочку городов
   }
   return res;
}
Когда путь сформирован полностью, результирующий массив будет состоять из одного элемента, являющегося массивом списка городов.


Создание нового участка пути

Левая ветвь бинарного дерева получается в результате добавления к уже частично построенному пути перехода wj -> wi, где wj,wi - положение в оставшейся матрице расстояний:

Salesway.prototype.setWay = function (way, wj, wi)
{
   this.copyWay(way, wj, wi);                     // копируем way, исключая строку,столбец

   var cityJ = way.src[wj], cityI = way.des[wi];  // города перехода  cityJ -> cityI
   var backJ, backI;                              // обратный переход backI -> backJ
   for(var k=0; k < this.num; k++)                // находим его индексы в массиве
      if(this.src[k]===cityI){ backJ=k; break; }  // src
   for(var k=0; k < this.num; k++)                // и в массиве
      if(this.des[k]===cityJ){ backI=k; break; }  // des
   this.dists[backJ*this.num+backI] = -1;         // блокируем переход cityI -> cityJ

   this.frw[cityJ] = cityI;                       // ставим переход из cityJ -> cityI
   var end = cityI;                               // и ищем конец этой цепочки переходов
   while( this.frw[end] >=0 )
      end = this.frw[end];

   this.rev[cityI] = cityJ;                       // запоминаем обратный переход
   var beg = cityJ;                               // и ищем начало этой цепочки переходов
   while( this.rev[beg] >=0 )
      beg = this.rev[beg];

   var kJ=0, kI=0;                                // блокируем end -> beg в матрице
   for(; kJ < this.num; kJ++)                     // для этого ищем индекс kJ конца end
      if(this.src[kJ]===end) break;               // в строках src
   for(; kI < this.num; kI++)                     // и индекс kI
      if(this.des[kI]===beg) break;               // в колонках des
   this.dists[kJ*this.num+kI] = -1;               // и блокируем его
}
Функция копирования массивов copyWay c исключением строки wj и колонки wi выглядит следующим образом:
Salesway.prototype.copyWay = function (way, wj, wi)
{
   for(var j=0, k=0; j < way.num; j++)            // копируем заголовки строчек
      if(j!==wj) this.src[k++] = way.src[j];

   for(var i=0, k=0; i < way.num; i++)            // копируем заголовки колонок
      if(i!==wi) this.des[k++] = way.des[i];

   for(var j=0, k=0; j < way.num; j++)            // копируем матрицу расстояний
      if(j!==wj)                                  // удаляем строку wj
         for(var i=0; i < way.num; i++)
           if(i!==wi)                             // удаляем столбец wi
              this.dists[k++] = way.dists[j*way.num+i];

   for(var k=0; k < way.frw.length; k++){         // копируем участки пути
      this.frw[k] = way.frw[k];
      this.rev[k] = way.rev[k];
   }
   this.boundLen = way.boundLen;
}


Подготовка к методу ветвей и границ

Salesman.prototype.boundInit = function()
{
   var root = new Salesway(this.N, this.N);       // корневой узел дерева
   root.init(this.dists);                         // помещаем в него исходную матрицу
   root.calcBound();                              // вычисляем нижнюю границу
   this.open = new Queue();                       // приоритетная очередь
   this.open.lt = function(a,b) { return a.value < b.value };
   this.open.unshift(root);                       // помещаем в неё корень
}

Таймерная часть метода

Salesman.prototype.boundTimer = function()
{
   var num = this.loops;                          // итераций за один "тик" таймера
   var time = window.performance.now();
   do{
      if( this.open.empty() ){                    // очередь пустая
         this.stop();                             // убиваем таймер
         break;
      }

      var rt = this.open.shift();                 // берём узел из начала списка
      if(rt.boundLen >= this.minLen){             // уже есть путь лучше
         this.cntCuts++;                          // число отсечений
         continue;
      }

      var lf = new Salesway(rt.num-1, this.N);    // левая ветвь дерева
      var res = rt.selectPair();                  // ищем максимальный прирост оценки
      lf.setWay(rt, res.jj, res.ii);              // копируем, исключая res.jj, res.ii
      lf.calcBound();                             // вычисляем нижнюю границу длины пути

      rt.boundLen += res.max;                     // уточняем границу снизу
      if(rt.boundLen < this.minLen){              // сначала правую ветку
         rt.blockWay(res);                        // блокируем путь res.jj - res.ii
         rt.value = rt.boundLen/(this.N+1-rt.num);// по value сортируем очередь
         this.open.unshift(rt);                   // запихиваем в неё правую ветвь
         this.cntNode++;                          // число сформированных узлов
      }
      if(lf.num===1){                             // дошли до конца
         this.cntWays++;
         if(lf.boundLen < this.minLen){           // и это путь лучший (? финал)
            this.minLen = lf.boundLen;            // запоминаем его
            this.minWay = lf.getPath()[0];        // получаем массив пути
            this.timeMin = this.timeTot + (window.performance.now() - time);
         }
      }
      else if(lf.boundLen < this.minLen){         // оценка снизу меньше минимальной
         lf.value = lf.boundLen/(this.N+1-lf.num);
         this.open.unshift(lf);                   // помещаем узел левой ветки
         this.cntNode++;                          // число сформированных узлов
      }
      else
         this.cntCuts++;                          // число отчечений
   } while(--num);                                // крутимся loops раз
   this.timeTot += window.performance.now() - time;

   this.outInfo(", open: "+this.open.length);     // выводим текущие результаты
}