Формулировка задачи

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

Сгенерируйте остановки

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

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 бинарных переменных (# переменные = 200 выбирают 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 возвращает вектор с количеством подтура к который каждое ребро belongs.

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

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

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

disp(output.absolutegap)

     0

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