容量制約付き配送計画問題

メモ

インストールが必要: この機能には、MATLAB Support Package for Quantum Computing が必要です。

容量制約付き配送計画は、ナップサック問題と巡回セールスマン問題を組み合わせた問題です。この問題では、車両 (または複数の車両) が地理的に分散した顧客群を訪問します。車両には容量制約があり、その容量とは車両がそれぞれの顧客に配送できる数量を指します。この問題には中央のデポがあり、車両は顧客群への訪問、すなわちルートの完了後にデポへ戻らなければなりません。この問題の目的は、顧客群を訪問するルートの総距離をコストとし、そのコストを最小化しながら顧客を訪問することです。

次の図は、単一の地点であるデポから出発する 4 つのルートを示しています。これらのルートは最適解ではありません。少なくともノード 2 とノード 3 は逆の順序で訪問すべきだからです。ノード 2 とノード 3 を含むルートには自己交差があり、このような交差は最適な巡回経路では発生しません。

Four tours starting from a central point, the depot

容量制約付き配送計画問題を解くには、Feld および共著者が示す手順に従います[1]。Feld は複数の解法を示していますが、この例ではそのうちの 1 つのみを使用します。

  • 1 回のルートで車両が訪問する顧客群を表すクラスターを作成します。このステップはナップサック問題です。

  • 各クラスターに対して巡回セールスマン問題を解きます。

クラスター作成問題は古典的アルゴリズムを使用して解きます。巡回セールスマン問題は QUBO 問題として定式化して解きます。QUBO 問題は、内部でタブー探索アルゴリズムを使用する solve 関数を使用して解きます。

問題データと変数

24 人の顧客についてランダムな整数地点を作成し、各地点の需要量を 100 ～ 2500 の範囲で 100 刻みのランダムな整数値とします。追加の 1 地点であるデポは、座標 [0,0] に配置します。

rng(1) % For reproducibility
numCustomers = 24; % Depot at [0 0] makes 25 locations
depot = [0 0]; % Depot at the origin
loc = [depot; randi([-50,50],numCustomers,2)]; % Integers from -50 to 50, 24-by-2
demand = 100*randi([1,25],numCustomers,1);
capacity = 6000;

地点と需要量をプロットします。

% Plot the locations with demands overlaid
figure;
scatter(loc(:,1),loc(:,2),'filled','SizeData',25);
hold on
text(loc(:,1),loc(:,2),["Depot"; num2str(demand)]);
title("Customer Locations and Demands");

24 random integer locations from -50 to 50 in each coordinate, with associated demands from 100 to 2500

座標ベクトルから顧客間の距離の行列を作成します。ピタゴラスの定理とユークリッド距離を使用します。

numLocations = numCustomers + 1;
[X,Y] = meshgrid(1:numLocations);
dist = hypot(loc(X(:),1)- loc(Y(:),1),loc(X(:),2) - loc(Y(:),2));
dist = reshape(dist,numLocations,numLocations);

問題データからサイト 1 のデポを除外して、この例の残りで使用する変数を設定します。

customerCoords = loc(2:end,:);
costMatrix = dist;
vehicleCapacity = capacity;

クラスターの作成

Feld [1]では、クラスターを作成するための次の 2 ステップの方法が示されています。

  1. 問題の座標と容量制約に基づいて初期クラスターを作成します。

  2. 初期クラスターを改善して、より短い巡回セールスマン ルートを得ます。

この例では、次の反復的な方法を使用して初期クラスターを作成します。

  1. デポから最も離れた顧客から開始します。

  2. クラスター内の顧客の現在の平均位置に最も近い顧客を選びながら、顧客を 1 人ずつ反復的に追加します。つまり、クラスターに追加する 2 人目の顧客は 1 人目に最も近い顧客とし、3 人目は最初の 2 人の平均位置に最も近い顧客とし、その後も同様に追加します。

  3. 次の顧客を追加すると容量制約を超える場合は、その時点でクラスターへの追加を停止します。この時点で最初のステップに戻り、新しいクラスターの作成を開始します。

クラスターは、customers、center、および demand のフィールドをもつ構造体として作成します。addToCluster 補助関数を使用して顧客をクラスターに追加します。 

clusters = [];

% CLUSTER GENERATION

% Select the first node of the first cluster
[~,idx] = max(costMatrix(1,:));

currentCluster.Customers = idx;
currentCluster.Center = customerCoords(idx,:);
currentCluster.Demand = demand(idx);

unclustered.Coords = customerCoords;
unclustered.Customers = 1:numCustomers;
unclustered.Coords(idx,:) = [];
unclustered.Customers(idx) = [];

% While not all the nodes have beed added
while ~isempty(unclustered.Customers)
    % Find the next closest node to the cluster
    distance = sqrt(sum((currentCluster.Center - unclustered.Coords).^2,2));
    [~,nextIdx] = min(distance);
    next = unclustered.Customers(nextIdx);
    % Remove the node from the unclustered list
    unclustered.Coords(nextIdx,:) = [];
    unclustered.Customers(nextIdx) = [];

    % If the customer can be added without violating capacity constraints,
    % add the customer to the current cluster; otherwise, start a new cluster
    if currentCluster.Demand + demand(next) <= vehicleCapacity
        currentCluster = addToCluster(currentCluster,next,customerCoords,demand);
    else
        clusters = [clusters, currentCluster];
        currentCluster.Customers = next;
        currentCluster.Center = customerCoords(next,:);
        currentCluster.Demand = demand(next);
    end
end
clusters = [clusters, currentCluster];

初期クラスターを作成した後、それらを改善できます。そのためには、新しいクラスターの容量制約を超えない範囲で、顧客を現在のクラスターよりも新しいクラスターの平均位置に近づけられる場合、その顧客を別のクラスターに再割り当てします。可能なステップがなくなるまで、または改善ステップを 10 回実行するまで、クラスター改善ステップの実行を続けます。

クラスターの改善には addToCluster 補助関数と removeFromCluster 補助関数を使用します。

% CLUSTER IMPROVEMENT
iterations = 0;
% For each cluster
while iterations < 10
    for i = 1:numel(clusters)
        % And each customer within the cluster
        for customer = clusters(i).Customers
            % Calculate the customer's distance from the center
            d_i = sqrt(sum((clusters(i).Center - customerCoords(customer,:)).^2));
            % For each alternative cluster
            for j = [1:i-1, i+1:numel(clusters)]
                % Calculate the customer's distance to the center of the
                % alternative cluster
                d_j = sqrt(sum((clusters(j).Center - customerCoords(customer,:)).^2));
                % Move the customer to the alternative cluster if it is closer and
                % capacity constraints are met
                if d_j < d_i && clusters(j).Demand + demand(customer) < vehicleCapacity
                    % Remove the customer from the original cluster
                    % and add the customer to the new cluster
                    clusters(i) = removeFromCluster(clusters(i), ...
                        customer, customerCoords, demand);
                    clusters(j) = addToCluster(clusters(j), ...
                        customer, customerCoords, demand);                    
                    break
                end
            end
        end
    end
    iterations = iterations + 1;
end

デポをサイト 1、最初の顧客をサイト 2 とするように、顧客ラベルを元の問題に合わせて調整します。

nRoutes = numel(clusters);
for i = 1:nRoutes
    clusters(i).Customers = clusters(i).Customers + 1;
end

巡回セールスマン問題の作成と求解

各クラスターは、デポから出発してデポへ戻る 1 つの車両ルートで訪問する顧客群を表します。各ルートについて、車両は可能な限り短い距離を移動する必要があります。この問題は、各クラスターに対する巡回セールスマン問題 (TSP) です。

各クラスターについて、そのクラスター内の顧客の距離行列を作成します。複数の TSP とその解を構造体にまとめます。各構造体には、その TSP に対応する顧客、顧客間の距離、および最短距離のルートである TSP の解が含まれます。

QUBO 問題へ変換して解くための TSP を作成します。

TSPsolutions = cell(nRoutes,1);
Routes = cell(nRoutes,1);
customerCoords = loc; % Return depot to list

各顧客クラスターについてペアワイズ距離を計算し、対応する TSP を定式化します。convertTSPtoQUBO 補助関数を使用して、各 TSP を QUBO 問題へ変換します。solvemyTSP 補助関数を使用して TSP を解きます。 

for rt = 1:nRoutes
    % Compute pairwise distance
    cluster = clusters(rt);
    coords = [depot; customerCoords(cluster.Customers,:)];
    M = height(coords);
    [X,Y] = meshgrid(1:M);
    dist = hypot(coords(X(:),1)- coords(Y(:),1),coords(X(:),2) - coords(Y(:),2));
    
    % Collect TSP data
    tsp.CostMatrix = reshape(dist,M,M);
    tsp.Customers = cluster.Customers;
    % Solve this TSP
    TSPsolutions{rt} = solvemyTSP(tsp);
    Routes{rt} = TSPsolutions{rt}.Route;
end

plotClusters 補助関数を使用して得られたルートをプロットします。

plotClusters(clusters,loc,Routes)

Five colored line plots showing routes starting from and ending at depot.

各ルートで必要となる容量を表示します。

d = arrayfun(@(x)x.Demand,clusters)
d =

        5000        5800        5400        4600        5700

各ルートの需要量は、容量上限である 6000 未満です。

補助関数

次のコードは addToCluster 補助関数を作成します。この補助関数は updateCluster 補助関数を使用していることに注意してください。

function cluster = addToCluster(cluster,next,customerCoordinates,demandVector)
    cluster.Customers = [cluster.Customers, next];
    cluster = updateCluster(cluster,next,customerCoordinates,demandVector,1);
end

次のコードは removeFromCluster 補助関数を作成します。この補助関数は updateCluster 補助関数を使用していることに注意してください。

function cluster = removeFromCluster(cluster,next,customerCoordinates,demandVector)
    cluster.Customers(cluster.Customers == next) = [];   
    cluster = updateCluster(cluster,next,customerCoordinates,demandVector,-1);
end

次のコードは updateCluster 補助関数を作成します。

function cluster = updateCluster(cluster,next,customerCoordinates,demandVector,s)
% Update the cluster as well as the lists of available coordinates and
% customers

currentN = numel(cluster.Customers);
previousN = currentN - s;

% Update center
newX = (cluster.Center(1)*previousN + s*customerCoordinates(next,1))/(currentN);
newY = (cluster.Center(2)*previousN + s*customerCoordinates(next,2))/(currentN);
cluster.Center = [newX,newY];

% Update demand
cluster.Demand = cluster.Demand + s*demandVector(next);
end

次のコードは plotClusters 補助関数を作成します。

function plotClusters(clusters,customerCoords,Routes)
    % Plot the result of the clustering algorithm in 2D space. Plot the
    % routes as well if they are provided. 

    % Plot the depot
    f = figure;

    scatter(customerCoords(1,1), customerCoords(1,2),"filled");
    ax = f.CurrentAxes;
    text(customerCoords(1,1), customerCoords(1,2),"Depot");
    hold on

    % Plot each cluster and label customers
    for i = 1:numel(clusters)
        customer = clusters(i).Customers;
        scatter(customerCoords(customer,1),...
                customerCoords(customer,2),"filled",SizeData=25); % Locations
        text(customerCoords(customer,1),...
               customerCoords(customer,2),num2str(customer')); % Labels
    end

    % If the routes are provided, plot them
    if nargin > 2
        % Reset color index - 1 is for the depot, so start at 2
        colorIdx = 1;
        for k = 1:numel(Routes)
            route = Routes{k};
            % Advance color index for the next route
            colorIdx = colorIdx + 1;
            % Append depot as the start of the route
            cu = [1 clusters(k).Customers]; 
            tr = cu(route); % Route in original indices
            for customer = 1:numel(tr) - 1
                ax.ColorOrderIndex = colorIdx; % Keep current route color
                plot(customerCoords(tr(customer:customer+1),1),...
                    customerCoords(tr(customer:customer+1),2));
            end
            ax.ColorOrderIndex = colorIdx;
            plot(customerCoords([tr(end),tr(1)],1),...
                customerCoords([tr(end),tr(1)],2)); 
        end
    end
    drawnow
    hold off
end

次のコードは convertTSPtoQUBO 補助関数を作成します。

function QP = convertTSPtoQUBO(dist)
%   QP = CONVERTTSPTOQUBO(DIST) returns a QUBO problem from the traveling salesperson
%   problem specified by the distance matrix DIST. DIST is an N-by-N
%   nonnegative matrix where DIST(i,j) is the distance between locations
%   i and j.

% Copyright 2023 The MathWorks, Inc.

N = size(dist,1);
% Create constraints on routes
A = eye(N);
B = ones(N);
Q0 = kron(A,B);
Q1 = kron(B,A);
% Create upper diagonal matrices of distances
v = ones(N-1,1);
A2 = diag(v,1);
Q2 = kron(B,A2); % Q2 has a diagonal just above the main diagonal in each block
C = kron(dist,B);
Q2 = Q2.*C; % Q2 has an upper diagonal dist(i,j)
% Create dist(j,i) in the upper-right corner of each block
E = zeros(N);
E(1,N) = 1;
Q3 = kron(B,E); % Q3 has a 1 in the upper-right corner of each block
CP = kron(dist',B); % dist' for D(j,i)
Q3 = Q3.*CP; % Q3 has dist(j,i) in the upper-right corner of each block
% Add the multipliers
M = max(max(dist));
QN = sparse(M*(Q0 + Q1)*N^2 + Q2 + Q3);
% Symmetrize
QN = (QN + QN.')/2;

% Include the constant and linear terms
c = -4*ones(N^2,1)*M*N^2;
d = 2*N*M*N^2;

QP = qubo(QN,c,d);

end

次のコードは solvemyTSP 補助関数を作成します。この補助関数は solveTSPwithTabu 補助関数と convertSolutionToRoute 補助関数を使用していることに注意してください。 

function TSPsolution = solvemyTSP(tsp)
% solvemyTSP solves the TSP by first converting it to a QUBO problem 
% and then using tabu search to find a solution. 

% Inputs: 
%   tsp: Structure with fields
%       customers: the customers for the tsp
%       costMatrix: the cost matrix for the tsp
% Outputs: 
%   TSP_solution: Structure with fields
%       Route: The order of the customers in the best route found. For example,
%       if there are 5 customers in this tsp, the Route might be [2 3 1 5 4]
%       for customers [13 6 3 12 7]. 
%       customers: The customers in the current tsp, taken directly
%       from the tsp input
%       quboFval: The fval returned by the qubo algorithm

n = numel(tsp.Customers);

% If the number of customers is less than 3, the solution is trivial and
% there is no need to solve the QUBO. 
if n < 3
    TSPsolution.Route = 1:n;
    TSPsolution.Customers = tsp.Customers;
    TSPsolution.QuboFval = tsp.CostMatrix(1,2) + tsp.CostMatrix(2,1);
    return
end
    
% Solve with tabu search
TSPsol = solveTSPwithTabu(tsp.CostMatrix);

% Convert the QUBO solution to a route
TSPsolution.Route = convertSolutionToRoute(TSPsol);
TSPsolution.QuboFval = TSPsol.BestFunctionValue;
TSPsolution.Customers = tsp.Customers;
end

次のコードは solveTSPwithTabu 補助関数を作成します。この補助関数は convertTSPtoQUBO 補助関数と checkTSPConstraints 補助関数を使用していることに注意してください。 

function solution = solveTSPwithTabu(costMatrix)
% Solves the TSP with tabu search. Convert the TSP to a QUBO
% according to the penalty term and return the best feasible solution. 
% Convert the TSP problem to QUBO
Q = convertTSPtoQUBO(costMatrix);
x = solve(Q);

validSolutions = checkTSPConstraints(x);
% If no valid solution is found yet, try again up to 10 times
if validSolutions == false
   it = 1;
   while (validSolutions == false) && it < 10
       x = solve(Q);
       validSolutions = checkTSPConstraints(x);
       it = it + 1;
   end
end

if validSolutions == false
   solution = [];
else
   solution = x;
end
end

次のコードは checkTSPConstraints 補助関数を作成します。

function validSolutions = checkTSPConstraints(x)
% The solutions are encoded as [x11, x12, x13, ... x21, x23, ...]'
% where xij indicates that customer i is visited in the jth position of the
% route. To establish that the constraints are satisfied, verify
% that each customer is visited only once and each position in the route
% has only one customer assigned to it.
N = size(x.BestX,1);
n = sqrt(N);
y = reshape(x.BestX,n,n);

% Ensure each customer has only one position assigned
r = all(sum(y,1) == 1);

% Ensure each route position (row) on has one customer assigned
s = all(sum(y,2) == 1);
validSolutions = r && s;
end

次のコードは convertSolutionToRoute 補助関数を作成します。

function Route = convertSolutionToRoute(sol)
% The solutions are encoded as [x11, x12, x13, ... x21, x23, ...]'
% where xij indicates that customer i is visited in the jth position of the
% route. Decode the solution into a route.
solution = sol.BestX;
n = sqrt(numel(solution));
selected_indices = find(solution);
order = mod(selected_indices - 1, n) + 1;
[~, node_order] = sort(order);
nodes = 1:n;

Route = nodes(node_order);
end

参照

[1] Feld, Sebastian, Christoph Roch, Thomas Gabor, Christian Seidel, Florian Neukart, Isabella Galter, Wolfgang Mauerer, and Claudia Linnhoff-Popien. A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer. Available at https://arxiv.org/abs/1811.07403.

参考

関数

オブジェクト

トピック