Документация

Чертите карту и остановки

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

figure;

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
plot(x,y,'Color','red'); % draw the outside border
hold on
% Add the stops to the map
plot(stopsLon,stopsLat,'*b')
hold off

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

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

Вычислите расстояния между точками

Поскольку существует 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'*trips

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

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

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

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

Включайте целевую функцию в проблему.

tsp.Objective = dist'*trips;

Ограничения равенства

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

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

constrips = sum(trips) == nStops;
tsp.Constraints.constrips = constrips;

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

constr2trips = optimconstr(nStops,1);
for stops = 1:nStops
    whichIdxs = (idxs == stops);
    whichIdxs = any(whichIdxs,2); % start or end at stops
    constr2trips(stops) = 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]

Визуализируйте решение

hold on
segments = find(tspsol.trips); % Get indices of lines on optimal path
lh = zeros(nStops,1); % Use to store handles to lines on plot
lh = updateSalesmanPlot(lh,tspsol.trips,idxs,stopsLon,stopsLat);
title('Solution with Subtours');

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

Подсовершите поездку по ограничениям

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

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

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

detectSubtours функция анализирует решение и возвращает массив ячеек векторов. Каждый вектор в массиве ячеек содержит остановки, вовлеченные в тот конкретный подтур.

tours = detectSubtours(tspsol.trips,idxs);
numtours = length(tours); % 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
        subTourIdx = tours{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);
        a = false(length(idxs),1);
        for jj = 1:length(variations)
            whichVar = (sum(idxs==subTourIdx(variations(jj,1)),2)) & ...
                       (sum(idxs==subTourIdx(variations(jj,2)),2));
            a = a | whichVar;
        end
        tsp.Constraints.(sprintf('subtourconstr%i',k)) = sum(trips(a)) <= length(subTourIdx)-1;
        k = k + 1;
    end
    % Try to optimize again
    [tspsol,fval,exitflag,output] = solve(tsp,'options',opts);

    % Visualize result
    lh = updateSalesmanPlot(lh,tspsol.trips,idxs,stopsLon,stopsLat);

    % How many subtours this time?
    tours = detectSubtours(tspsol.trips,idxs);
    numtours = length(tours); % 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

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