はじめに

列生成法により配送計画の最適化を実行してみました。
理論的な説明は以下のサイトを参照してください。
動的計画法による厳密解やOR-Toolsとの比較も実施しました。

データの作成

import math
import pulp
import random

sd = 1
random.seed(sd)

# 店舗数
N = 10
# XY座標 先頭は倉庫
XY_list = []
for i in range(N+1):
    XY_list.append(random.sample(range(100),k=2))
# 移動距離
C_list = [[0]*(N+1) for _ in range(N+1)]
for i in range(N+1):
    for j in range(N+1):
        xi,yi = XY_list[i]
        xj,yj = XY_list[j]
        dx,dy = xi-xj,yi-yj
        C_list[i][j] = int(math.sqrt(dx*dx+dy*dy))
# 荷物量
W_list = random.sample(range(1,100),k=N+1)
W_list[0] = 0 # 倉庫はゼロ
# 最大積載量
W = max(W_list)*2//10*10

主問題

与えられた配送ルートの集合から移動距離を最小化する組合せを選択します。

# d_list[i]: ルートiの移動距離
# a_list[i][j]: ルートiに店舗jが含まれる場合は1
def solve_main(d_list,a_list):
    prob = pulp.LpProblem('delivery_main',sense=pulp.LpMinimize)
    # 変数 ルートを選択するかどうか
    M = len(a_list)
    x_list = [pulp.LpVariable(f'v_{i}',lowBound=0,cat='Binary') for i in range(M)]
    # 目的関数 移動距離の合計
    obj = pulp.lpDot(d_list, x_list)
    prob.setObjective(obj)
    # 制約条件 全ての店舗をカバーする
    for j in range(N):
        cf = [a_list[i][j] for i in range(M)]
        ct = pulp.LpConstraint(pulp.lpDot(cf,x_list),sense=pulp.LpConstraintGE,name=f'c_{j}',rhs=1)
        prob.addConstraint(ct)
    # 最適化実行
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    xv_list = [0]*M
    for i in range(M):
        if x_list[i].value() > 0:
            xv_list[i] = 1
    return (prob,xv_list)

双対問題

主問題を線形緩和した双対問題を解きます。

def solve_dual(d_list,a_list):
    prob = pulp.LpProblem('delivery_dual',sense=pulp.LpMaximize)
    # 変数
    y_list = [pulp.LpVariable(f'v_{i}',lowBound=0,cat='Continuous') for i in range(N)]
    # 目的関数
    obj = pulp.lpSum(y_list)
    prob.setObjective(obj)
    # 制約条件
    M = len(a_list)
    for i in range(M):
        ct = pulp.LpConstraint(pulp.lpDot(a_list[i],y_list),sense=pulp.LpConstraintLE,name=f'c_{i}',rhs=d_list[i])
        prob.addConstraint(ct)
    # 最適化実行
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    yv_list = [y.value() for y in y_list]
    return (prob,yv_list)

配送ルート問題

双対問題の制約条件を破る配送ルートを求めます。

# y_list[j]: 双対問題による店舗jの解
def solve_child(y_list):
    prob = pulp.LpProblem('delivery_child',sense=pulp.LpMinimize)
    # 変数 i,j間を移動する場合は1
    z_list = [[None]*(N+1) for _ in range(N+1)]
    for i in range(N+1):
        for j in range(N+1):
            if i == j: continue
            z_list[i][j] = pulp.LpVariable(f'z_{i}_{j}',cat='Binary')
    u_list = [pulp.LpVariable(f'u_{i}',lowBound=1,upBound=N,cat='Integer') for i in range(N)]
    # 目的関数
    e_list = []
    for i in range(N+1):
        for j in range(N+1):
            if i == j: continue
            if j > 0:
                e_list.append((C_list[i][j]-y_list[j-1]) * z_list[i][j])
            else:
                e_list.append((C_list[i][j]) * z_list[i][j])
    prob.setObjective(pulp.lpSum(e_list))
    # 最大積載量制約
    e_list = []
    for i in range(N+1):
        for j in range(N+1):
            if i == j: continue
            e_list.append(W_list[j] * z_list[i][j])
    ct = pulp.LpConstraint(pulp.lpSum(e_list),sense=pulp.LpConstraintLE,name=f'c_w',rhs=W)
    prob.addConstraint(ct)
    # 店舗制約
    from_list = [None]*(N+1) # 流出
    to_list = [None]*(N+1) # 流入
    for i in range(N+1):
        a_list = []
        b_list = []
        for j in range(N+1):
            if j == i: continue
            a_list.append(z_list[i][j])
            b_list.append(z_list[j][i])
        from_list[i] = pulp.lpSum(a_list)
        to_list[i] = pulp.lpSum(b_list)
    for i in range(N+1):
        if i == 0:
            ct = pulp.LpConstraint(from_list[0],sense=pulp.LpConstraintEQ,name=f'c_from_0',rhs=1)
            prob.addConstraint(ct)
            ct = pulp.LpConstraint(to_list[0],sense=pulp.LpConstraintEQ,name=f'c_to_0',rhs=1)
            prob.addConstraint(ct)
        else:
            ct = pulp.LpConstraint(from_list[i],sense=pulp.LpConstraintLE,name=f'c_from_{i}',rhs=1)
            prob.addConstraint(ct)
            ct = pulp.LpConstraint(to_list[i],sense=pulp.LpConstraintLE,name=f'c_to_{i}',rhs=1)
            prob.addConstraint(ct)
            ct = pulp.LpConstraint(from_list[i]-to_list[i],sense=pulp.LpConstraintEQ,name=f'c_eq_{i}',rhs=0)
            prob.addConstraint(ct)
    # 経路制約
    for i in range(N):
        for j in range(N):
            if i == j: continue
            e = u_list[i] + 1 - 1000 * (1-z_list[i+1][j+1]) - u_list[j]
            ct = pulp.LpConstraint(e,sense=pulp.LpConstraintLE,name=f'c_{i}_{j}',rhs=0)
            prob.addConstraint(ct)
    # 最適化実行
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    # 新しい移動経路aと移動距離d
    a = [0]*(N)
    d = 0
    for i in range(N+1):
        for j in range(N+1):
            if i == j: continue
            if z_list[i][j].value() > 0.5:
                if j > 0:
                    a[j-1] = 1
                d += C_list[i][j]
    return (prob,d,a)

列生成法による最適化

上で定義した関数を使用して、列生成法による最適化を実行します。

# 初期ルート
a_list = [] # 行が一つのルートに対応
d_list = [] # ルートの移動距離
for i in range(N):
    a = [0]*(N)
    a[i] = 1
    a_list.append(a)
    d_list.append(C_list[0][i+1]*2) # 往復

for k in range(100):
    prob,y_list = solve_dual(d_list,a_list)
    prob,d,a = solve_child(y_list)
    v = prob.objective.value()
    print('child',k,v)
    if v >= -0.1: break
    d_list.append(d)
    a_list.append(a)

prob,x_list = solve_main(d_list,a_list)
print('main',pulp.LpStatus[prob.status],prob.objective.value())
d = 0
for i in range(len(a_list)):
    if x_list[i] == 1:
        print(a_list[i],d_list[i])
        d += d_list[i]
print(d)

以下の結果が得られました。

main Optimal 571.0
[0, 1, 0, 0, 1, 1, 1, 0, 0, 0] 217
[1, 0, 0, 1, 0, 0, 0, 1, 0, 1] 223
[0, 0, 1, 0, 0, 0, 0, 0, 1, 0] 131
571

移動距離の合計は571で、3台のトラックを使用しています。
[]内の0,1は訪問する店舗を示しており、全ての店舗を一回ずつ回っていることが分かります。
新しい経路は29個追加されました。

動的計画法による最適化

店舗数Nが小さければ動的計画法で最適解を得ることができます。
指定された店舗をすべて回る移動距離の最小値を求める関数を定義しておいて、
店舗の部分集合を列挙します。

# 一台のトラックでkに属する店舗をすべて回って倉庫に戻る移動距離の最小値
INF = 1<<61
def solve_child(k):
    w = 0
    for i in range(N):
        if k&(1<<i) > 0:
            w += W_list[i+1]
    if w > W:
        return INF
    s_list = [k]
    s = (k-1)&k
    while s > 0:
        s_list.append(s)
        s = (s-1)&k
    s_list.reverse()
    dp = dict()
    for s in s_list:
        dp[s] = [INF]*(N)
        for i in range(N):
            if s&(1<<i) == 0: continue
            d = INF
            for j in range(N):
                if s&(1<<j) == 0: continue
                if j == i: continue
                d = min(d,dp[s-(1<<i)][j]+C_list[j+1][i+1])
            if d == INF:
                dp[s][i] = C_list[0][i+1]
            else:
                dp[s][i] = d
    ans = INF
    for i in range(N):
        ans = min(ans,dp[k][i]+C_list[i+1][0])
    return ans

# 複数のトラックを使用した場合の移動距離の最小値
N2 = 1<<N
dp = dict()
child = dict()
for k in range(N2):
    a = solve_child(k)
    child[k] = (k,a)
    s = (k-1)&k
    while s > 0:
        c = solve_child(s)
        b = c+dp[k-s]
        if a > b:
            a = b
            child[k] = (s,c)
        s = (s-1)&k
    dp[k] = a
k = N2-1
while k > 0:
    s,c = child[k]
    x = [0]*(N)
    for i in range(N):
        if s&(1<<i) > 0:
            x[i] = 1
    print(x,c)
    k -= s
print(dp[N2-1])

列生成法よりも1だけ改善した解が得られました。

[1, 0, 0, 1, 0, 0, 0, 0, 0, 1] 222
[0, 0, 1, 0, 0, 0, 0, 0, 1, 0] 131
[0, 1, 0, 0, 1, 1, 1, 1, 0, 0] 217
570

OR-Toolsによる最適化

GoogleのOR-Toolsでも実行してみました。

データ作成は省略。

"""Capacited Vehicles Routing Problem (CVRP)."""
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data['distance_matrix'] = C_list
    data['demands'] = W_list
    data['vehicle_capacities'] = [W]*N
    data['num_vehicles'] = N
    data['depot'] = 0
    return data

def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    total_distance = 0
    total_load = 0
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        route_distance = 0
        route_load = 0
        route_node = []
        while not routing.IsEnd(index):
            node_index = manager.IndexToNode(index)
            route_load += data['demands'][node_index]
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(previous_index, index, vehicle_id)
            route_node.append(node_index)
        if route_distance > 0:
            print(vehicle_id,route_distance,route_node)
        total_distance += route_distance
        total_load += route_load
    print(f'Total {total_distance}')

def main():
    """Solve the CVRP problem."""
    # Instantiate the data problem.
    data = create_data_model()

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),data['num_vehicles'], data['depot'])

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        # Convert from routing variable Index to distance matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['distance_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Capacity constraint.
    def demand_callback(from_index):
        """Returns the demand of the node."""
        # Convert from routing variable Index to demands NodeIndex.
        from_node = manager.IndexToNode(from_index)
        return data['demands'][from_node]

    demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # null capacity slack
        data['vehicle_capacities'],  # vehicle maximum capacities
        True,  # start cumul to zero
        'Capacity')

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.FromSeconds(1)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)

if __name__ == '__main__':
    main()

最適解と一致する結果が得られました。

0 217 [0, 8, 5, 7, 2, 6]
2 131 [0, 3, 9]
9 222 [0, 10, 1, 4]
Total 570

おわりに

店舗数N=10として、乱数シードを変えながら、列生成法を100回実行してみました。
列生成法により最適解が得られる割合は50%でした。
また最適解でない場合の、最適解との差異は平均16(3%)となりました。
同様にOR-Toolsでも100回実行したところ、98回、最適解が得られました。
最適解との差異は、1と5であり、非常に精度の良い解が得られました。
【追記】処理時間の上限値を5秒に変更したところ、全て最適解が得られました。