Этот пример показывает, как использовать бинарное целочисленное программирование, чтобы решить классическую проблему коммивояжера. Эта проблема включает нахождение самого короткого закрытого тура (путь) через набор остановок (города). В этом случае существует 200 остановок, но можно легко заменить переменную nStops, чтобы получить различный проблемный размер. Вы будете решать начальную проблему и видеть, что решение имеет подтуры. Это означает, что найденное оптимальное решение не дает один непрерывный путь через все точки, но вместо этого имеет несколько разъединенных циклов. Вы будете затем использовать итеративный процесс определения подтуров, добавления ограничений и повторного выполнения оптимизации, пока подтуры не будут устранены.
Для основанного на проблеме подхода смотрите проблему Коммивояжера: основанный на проблеме.
Сгенерируйте случайные остановки в грубом многоугольном представлении континентальных США.
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')
Сформулируйте проблему коммивояжера для целочисленного линейного программирования можно следующим образом:
Сгенерируйте все возможные прохождения, имея в виду все отличные пары остановок.
Вычислите расстояние для каждого прохождения.
Функция стоимости, чтобы минимизировать является суммой расстояний прохождения для каждого прохождения на туре.
Переменные решения являются двоичным файлом, и сопоставленный с каждым прохождением, где каждый 1 представляет прохождение, которое существует в туре, и каждый 0 представляет прохождение, которое не находится в туре.
Чтобы гарантировать, что тур включает каждую остановку, включайте линейное ограничение, что каждая остановка включена точно два прохождения. Это означает одно прибытие и одно отклонение от остановки.
Поскольку существует 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 является бинарным вектором решения. Это - расстояние тура, который вы пытаетесь минимизировать.
Проблема имеет два типа ограничений равенства. Первое осуществляет это должно быть 200 общих количеств прохождений. Второе осуществляет ту каждую остановку, должен иметь два прохождения, присоединенные к нему (должно быть прохождение в каждую остановку и прохождение, отбыв из каждой остановки).
Задайте первый тип ограничения равенства, что у вас должны быть прохождения nStops в форме Aeq*x_tsp = beq.
Aeq = spones(1:length(idxs)); % Adds up the number of trips
beq = nStops;Чтобы задать второй тип ограничения равенства, что должно быть два прохождения, присоединенные к каждой остановке, расширяют матрицу Aeq как разреженную.
Aeq = [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+1,:) = whichIdxs'; % include in the constraint matrix end beq = [beq; 2*ones(nStops,1)];
Все переменные решения являются двоичным файлом. Теперь, установите аргумент intcon на количество переменных решения, поместите нижнюю границу 0 на каждом и верхней границе 1.
intcon = 1:lendist; lb = zeros(lendist,1); ub = ones(lendist,1);
Проблема готова к решению. Чтобы подавить итеративный вывод, выключите отображение по умолчанию.
opts = optimoptions('intlinprog','Display','off'); [x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);
segments = find(x_tsp); % Get indices of lines on optimal path lh = zeros(nStops,1); % Use to store handles to lines on plot lh = updateSalesmanPlot(lh,x_tsp,idxs,stopsLon,stopsLat); title('Solution with Subtours');
Как видно на карте, решение имеет несколько подтуров. Ограничения, заданные до сих пор, не предотвращают эти подтуры. В порядке предотвратить любой возможный подтур, вам было бы нужно невероятно большое количество ограничений неравенства.
Поскольку вы не можете добавить все подтуристические ограничения, проявите итерационный подход. Обнаружьте подтуры в текущем решении, затем добавьте ограничения неравенства, чтобы предотвратить те конкретные подтуры. Путем выполнения этого вы находите подходящий тур в нескольких итерациях.
Устраните подтуры по-разному ограничения. Пример того, как это работы - то, если у вас есть пять точек на подтуре, затем у вас есть пять строк, соединяющих те точки, чтобы создать подтур. Устраните этот подтур путем реализации ограничения неравенства, чтобы сказать, что там должно быть меньше чем или равно четырем строкам между этими пятью точками.
Еще больше найдите все строки между этими пятью точками и ограничьте решение не иметь больше чем четыре из этих существующих строк. Это - правильное ограничение потому что, если бы пять или больше из строк существовали в решении, то решение имело бы подтур (график с узлы и ребра всегда содержат цикл).
Функция detectSubtours анализирует решение и возвращает массив ячеек векторов. Каждый вектор в массиве ячеек содержит остановки, вовлеченные в тот конкретный подтур.
tours = detectSubtours(x_tsp,idxs); numtours = length(tours); % 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 = 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); 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); % Visualize result lh = updateSalesmanPlot(lh,x_tsp,idxs,stopsLon,stopsLat); % How many subtours this time? tours = detectSubtours(x_tsp,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
Решение представляет выполнимый тур, потому что это - один замкнутый цикл. Но действительно ли это - тур минимальной стоимости? Один способ узнать состоит в том, чтобы исследовать выходную структуру.
disp(output.absolutegap)
0
Малость абсолютного разрыва подразумевает, что решение или оптимально или имеет общую длину, которая является близко к оптимальному.