Постановка проблемы

Сформулируйте задачу коммивояжера для целочисленного линейного программирования следующим образом:

Генерировать остановки

Генерировать случайные остановки внутри сырого полигонального представления континентального США.

load('usborder.mat','x','y','xx','yy');
rng(3,'twister') % Makes a plot with stops in Maine & Florida, and is reproducible
nStops = 200; % You can use any number, but the problem size scales as N^2
stopsLon = zeros(nStops,1); % Allocate x-coordinates of nStops
stopsLat = stopsLon; % Allocate y-coordinates
n = 1;
while (n <= nStops)
    xp = rand*1.5;
    yp = rand;
    if inpolygon(xp,yp,x,y) % Test if inside the border
        stopsLon(n) = xp;
        stopsLat(n) = yp;
        n = n+1;
    end
end

Расчет расстояний между точками

Поскольку существует 200 остановок, существует 19 900 поездок, что означает 19 900 двоичных переменных (# variables = 200 choose 2).

Создайте все командировки, то есть все пары остановок.

idxs = nchoosek(1:nStops,2);

Рассчитайте все расстояния поездки, предполагая, что земля плоская, чтобы использовать правило пифагора.

dist = hypot(stopsLat(idxs(:,1)) - stopsLat(idxs(:,2)), ...
             stopsLon(idxs(:,1)) - stopsLon(idxs(:,2)));
lendist = length(dist);

С этим определением dist вектор, длина тура

dist'*x_tsp

где x_tsp - двоичный вектор решения. Это расстояние тура, которое вы пытаетесь минимизировать.

Создание графика и рисование карты

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

G = graph(idxs(:,1),idxs(:,2));

Отображение остановок с помощью графика. Постройте график узлов без рёбер графика.

figure
hGraph = plot(G,'XData',stopsLon,'YData',stopsLat,'LineStyle','none','NodeLabel',{});
hold on
% Draw the outside border
plot(x,y,'r-')
hold off

Ограничения

Создайте линейные ограничения, согласно которым каждая остановка имеет две связанные поездки, поскольку должна быть поездка на каждую остановку и поездка на каждую остановку.

Aeq = spalloc(nStops,length(idxs),nStops*(nStops-1)); % Allocate a sparse matrix
for ii = 1:nStops
    whichIdxs = (idxs == ii); % Find the trips that include stop ii
    whichIdxs = sparse(sum(whichIdxs,2)); % Include trips where ii is at either end
    Aeq(ii,:) = whichIdxs'; % Include in the constraint matrix
end
beq = 2*ones(nStops,1);

Двоичные границы

Все переменные принятия решений являются двоичными. Теперь установите intcon аргумент к числу решающих переменных, поставьте нижнюю границу 0 для каждой и верхнюю границу 1.

intcon = 1:lendist;
lb = zeros(lendist,1);
ub = ones(lendist,1);

Оптимизация с помощью intlinprog

Проблема готова к решению. Для подавления итеративного вывода отключите отображение по умолчанию.

opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

Создайте новый график с перемещениями решения в виде кромок. Для этого округляйте решение, если некоторые значения не являются целыми числами, и преобразуйте полученные значения в logical.

x_tsp = logical(round(x_tsp));
Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
% Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2)); % Also works in most cases

Визуализация решения

hold on
highlight(hGraph,Gsol,'LineStyle','-')
title('Solution with Subtours')

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

Подчиненные ограничения

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

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

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

Определите подтексты, определив подключенные компоненты в Gsol, граф, построенный с рёбрами в текущем решении. conncomp возвращает вектор с номером подчиненного элемента, которому принадлежит каждое ребро.

tourIdxs = conncomp(Gsol);
numtours = max(tourIdxs); % number of subtours
fprintf('# of subtours: %d\n',numtours);

# of subtours: 27

Включите линейные ограничения неравенства, чтобы устранить подтексты, и многократно вызовите решатель, пока не останется только один подтекст.

A = spalloc(0,lendist,0); % Allocate a sparse linear inequality constraint matrix
b = [];
while numtours > 1 % Repeat until there is just one subtour
    % Add the subtour constraints
    b = [b;zeros(numtours,1)]; % allocate b
    A = [A;spalloc(numtours,lendist,nStops)]; % A guess at how many nonzeros to allocate
    for ii = 1:numtours
        rowIdx = size(A,1) + 1; % Counter for indexing
        subTourIdx = find(tourIdxs == ii); % Extract the current subtour
%         The next lines find all of the variables associated with the
%         particular subtour, then add an inequality constraint to prohibit
%         that subtour and all subtours that use those stops.
        variations = nchoosek(1:length(subTourIdx),2);
        for jj = 1:length(variations)
            whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...
                       (sum(idxs==subTourIdx(variations(jj,2)),2));
            A(rowIdx,whichVar) = 1;
        end
        b(rowIdx) = length(subTourIdx) - 1; % One less trip than subtour stops
    end

    % Try to optimize again
    [x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,A,b,Aeq,beq,lb,ub,opts);
    x_tsp = logical(round(x_tsp));
    Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2),[],numnodes(G));
    % Gsol = graph(idxs(x_tsp,1),idxs(x_tsp,2)); % Also works in most cases
    
    % Visualize result
    hGraph.LineStyle = 'none'; % Remove the previous highlighted path
    highlight(hGraph,Gsol,'LineStyle','-')
    drawnow
    
    % How many subtours this time?
    tourIdxs = conncomp(Gsol);
    numtours = max(tourIdxs); % number of subtours
    fprintf('# of subtours: %d\n',numtours)
end

# of subtours: 20

# of subtours: 7

# of subtours: 9

# of subtours: 9

# of subtours: 3

# of subtours: 2

# of subtours: 7

# of subtours: 2

# of subtours: 1

title('Solution with Subtours Eliminated');
hold off

Качество решения

Решение представляет собой выполнимый обход, поскольку представляет собой единый замкнутый контур. Но стоит ли это тур с минимальными затратами? Один из способов выяснить это - изучить структуру вывода.

disp(output.absolutegap)

     0

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