本页对应的英文页面已更新,但尚未翻译。 若要查看最新内容,请点击此处访问英文页面。

推销员行程问题:基于求解器

此示例说明如何使用二元整数规划来求解经典的销售员差旅问题。此问题涉及找到一条历经一系列停留点(城市)的最短环程(路径)。在本例中有 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 个停留点,因此有 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'*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);

使用 intlinprog 进行优化

此问题已可以求解。要隐藏迭代输出,请关闭默认的迭代输出。

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

从地图上可以看出,解有几个子环程。到当前为止指定的约束没有阻止这些子环程的发生。要防止出现任何可能的子环程,您将需要极大数量的不等式约束。

子环程约束

由于无法添加所有子环程约束,需要采用迭代方法。检测当前解中的子环程,然后添加不等式约束以防止发生这些特定子环程。通过执行此操作,您只需几次迭代即可找到合适的环程。

使用不等式约束消除子环程。具体的工作原理如下:如果您在一个子环程中有五个点,则您有五条连接这些点的线来创建该子环程。可实现一个不等式约束来消除此子环程,即,这五个停留点之间存在的线必须小于或等于四条。

还要找到这五个停留点之间的所有线,并限制解中这些线不超过四条。这是正确的约束,因为如果解中存在五条或更多条线,则该解将有子环程(具有 n 个节点和 n 条边的图始终包含一个环路)。

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

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

相关主题