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

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

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

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

load('usborder.mat','x','y','xx','yy');
rng(3,'twister') % Makes 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) % tTest 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'*trips

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

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

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

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

Создать переменные и проблему

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

tsp = optimproblem;
trips = optimvar('trips',lendist,1,'Type','integer','LowerBound',0,'UpperBound',1);

Включите в задачу целевую функцию.

tsp.Objective = dist'*trips;

Ограничения

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

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

constr2trips = optimconstr(nStops,1);
for stop = 1:nStops
    whichIdxs = outedges(G,stop); % Identify trips associated with the stop
    constr2trips(stop) = sum(trips(whichIdxs)) == 2;
end
tsp.Constraints.constr2trips = constr2trips;

Решение начальной проблемы

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

opts = optimoptions('intlinprog','Display','off');
tspsol = solve(tsp,'options',opts)

tspsol = struct with fields:
    trips: [19900×1 double]

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

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

tspsol.trips = logical(round(tspsol.trips));
Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2),[],numnodes(G));
% Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,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

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

% Index of added constraints for subtours
k = 1;
while numtours > 1 % Repeat until there is just one subtour
    % Add the subtour constraints
    for ii = 1:numtours
        inSubTour = (tourIdxs == ii); % Edges in current subtour
        a = all(inSubTour(idxs),2); % Complete graph indices with both ends in subtour
        constrname = "subtourconstr" + num2str(k);
        tsp.Constraints.(constrname) = sum(trips(a)) <= (nnz(inSubTour) - 1);
        k = k + 1;        
    end
    
    % Try to optimize again
    [tspsol,fval,exitflag,output] = solve(tsp,'options',opts);
    tspsol.trips = logical(round(tspsol.trips));
    Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2),[],numnodes(G));
    % Gsol = graph(idxs(tspsol.trips,1),idxs(tspsol.trips,2)); % Also works in most cases
    
    % Plot new solution
    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

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