Main Content

本页翻译不是最新的。点击此处可查看最新英文版本。

推销员行程问题：基于问题

打开实时脚本

此示例说明如何使用二元整数规划来求解经典的销售员差旅问题。此问题涉及找到一条历经一系列停留点（城市）的最短环程（路径）。在本例中有 200 个停留点，但您可以很轻松地更改 nStops 变量以得到不同规模的问题。对最初的问题进行求解后得到的解会包含子环程。这意味着找到的最优解并没有给出一条穿过所有点的连续路径，而是有几个独立的环路。然后，您将使用迭代过程来确定子环程，添加约束，并重新运行优化，直到消除子环程。

要了解如何通过基于求解器的方法处理此问题，请参阅推销员行程问题：基于求解器。

问题表示

将要进行整数线性规划的销售员差旅问题表示如下：

  • 生成所有可能的行程，意味着经过所有不同停留点对组。

  • 计算每次行程的距离。

  • 要最小化的代价函数是行程中每次旅行的旅行距离之和。

  • 决策变量是与每个行程相关联的二元变量，其中每个 1 表示环程中存在的一次行程，每个 0 表示不在环程中的一次行程。

  • 为确保环程包括每个经停留点，应加入这样一个线性约束：每个停留点都正好涉及两段行程。这意味着一段行程到达该停留点，一段行程离开该停留点。

生成停留点

在美国大陆的粗糙多边形表示中生成随机停留点

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) % Test if inside the border
        stopsLon(n) = xp;
        stopsLat(n) = yp;
        n = n+1;
    end
end

计算停留点之间的距离

由于存在 200 个停留点，因此有 19900 次行程，这意味着有 19900 个二元变量（# 变量 = 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 是表示解采用的行程的二元向量。这是您尝试最小化的旅行距离。

创建图和绘制地图

将问题表示为图。创建一个图，其中停留点是节点，行程是边。

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;

约束

创建线性约束，使每个停留点都有两个关联行程，因为每个停留点上必须有一次到达该停留点的行程和一次离开该停留点的行程。

使用图表示法，通过查找连接到某个停留点的所有边来识别在该停留点开始或结束的所有行程。对于每个停留点，创建使该停留点的行程总和等于 2 的约束。

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

绝对间隙小意味着该解或者是最优的，或者总长度接近最优。

相关主题