www.gusucode.com > optim 案例源码 matlab代码程序 > optim/TravellingSalesmanExample.m

    %% Travelling Salesman Problem
% This example shows how to use binary integer programming to solve the
% classic travelling salesman problem. This problem involves finding the
% shortest closed tour (path) through a set of stops (cities). In this case
% there are 200 stops, but you can easily change the |nStops| variable to get
% a different problem size.  You'll solve the initial problem and see that
% the solution has subtours. This means the optimal solution found doesn't
% give one continuous path through all the points, but instead has several
% disconnected loops. You'll then use an iterative process of determining
% the subtours, adding constraints, and rerunning the optimization until
% the subtours are eliminated.

%   Copyright 2014 The MathWorks, Inc. 

%% Draw the Map and Stops
% Generate random stops inside a crude polygonal representation of the
% continental U.S.

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,xx,yy) % 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

%% Problem Formulation
% Formulate the travelling salesman problem for integer linear
% programming as follows:
%
% * Generate all possible trips, meaning all distinct pairs of stops.
%
% * Calculate the distance for each trip.
%
% * The cost function to minimize is the sum of the trip distances for each
% trip in the tour.
%
% * The decision variables are binary, and associated with each trip, where
% each 1 represents a trip that exists on the tour, and each 0 represents a
% trip that is not on the tour.
%
% * To ensure that the tour includes every stop, include the linear
% constraint that each stop is on exactly two trips. This means one arrival
% and one departure from the stop.

%% Calculate Distances Between Points
% Because there are 200 stops, there are 19,900 trips, meaning 19,900
% binary variables (# variables = 200 choose 2).
%
% Generate all the trips, meaning all pairs of stops.

idxs = nchoosek(1:nStops,2);

%%
% Calculate all the trip distances, assuming that the earth is flat in
% order to use the Pythagorean rule.

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

%%
% With this definition of the |dist| vector, the length of a tour is
%
% |dist'*x_tsp|
%
% where |x_tsp| is the binary solution vector. This is the distance of a
% tour that you try to minimize.

%% Equality Constraints
% The problem has two types of equality constraints.  The first enforces
% that there must be 200 trips total. The second enforces that each stop
% must have two trips attached to it (there must be a trip to each stop and a
% trip departing each stop).
%
% Specify the first type of equality constraint, that you must have
% |nStops| trips, in the form |Aeq*x_tsp = beq|.


Aeq = spones(1:length(idxs)); % Adds up the number of trips
beq = nStops;

%%
% To specify the second type of equality constraint, that there needs to be
% two trips attached to each stop, extend the |Aeq| matrix as sparse.

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)];

%% Binary Bounds
% All decision variables are binary. Now, set the |intcon| argument to the
% number of decision variables, put a lower bound of 0 on each, and an
% upper bound of 1.

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

%% Optimize Using intlinprog
% The problem is ready to be solved. Call the solver.
opts = optimoptions('intlinprog','Display','off');
[x_tsp,costopt,exitflag,output] = intlinprog(dist,intcon,[],[],Aeq,beq,lb,ub,opts);

%% Visualize the Solution

hold on
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');

%%
% As can be seen on the map, the solution has several subtours.  The
% constraints specified so far do not prevent these subtours from
% happening.  In order to prevent any possible subtour from happening, you
% would need an incredibly large number of inequality constraints.

%% Subtour Constraints
% Because you can't add all of the subtour constraints, take an iterative
% approach. Detect the subtours in the current solution, then add
% inequality constraints to prevent those particular subtours from
% happening. By doing this, you find a suitable tour in a few iterations.
%
% Eliminate subtours with inequality constraints. An example of how this
% works is if you have five points in a subtour, then you have five lines
% connecting those points to create the subtour. Eliminate this subtour by
% implementing an inequality constraint to say there must be less than or
% equal to four lines between these five points.
%
% Even more, find all lines between these five points, and constrain the
% solution not to have more than four of these lines present. This is a
% correct constraint because if five or more of the lines existed in a
% solution, then the solution would have a subtour (a graph with $n$ nodes
% and $n$ edges always contains a cycle).
%
% The |detectSubtours| function analyzes the solution and returns a cell
% array of vectors.  Each vector in the cell array contains the stops
% involved in that particular subtour.

tours = detectSubtours(x_tsp,idxs);
numtours = length(tours); % number of subtours
fprintf('# of subtours: %d\n',numtours);

%%
% Include the linear inequality constraints to eliminate subtours, and
% repeatedly call the solver, until just one subtour remains.

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

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

%% Solution Quality
% The solution represents a feasible tour, because it is a single closed
% loop. But is it a minimal-cost tour? One way to find out is to examine
% the output structure.

disp(output.absolutegap)

%%
% The smallness of the absolute gap implies that the solution is either
% optimal or has a total length that is close to optimal.