OR-Toolsの概要とサンプルコードの実行

OR-Toolsの概要

OR-Toolsは線形計画問題(LP; Linear Programming)、混合整数計画問題(MIP; Mixed-Integer Programming)、制約プログラミング(CP; Constraint Programming)、配送計画問題(VRP; Vehicle Routing Programming)を解く目的で開発されたオープンソースのソフトウェアです。

OR-Toolsのインストール

PythonからOR-Toolsを用いる場合、下記のコマンドを実行することでPyPIからダウンロードすることができます。

$ pip install ortools

OR-Toolsのサンプルコードの実行

\begin{align}
\mathrm{Objective} \quad &: 3x+y \quad \longrightarrow \quad \mathrm{Maximize} \\
\mathrm{Constraint} \quad &: x+y \leq 2, \, 0 \leq x \leq 1, \, 0 \leq y \leq 2
\end{align}

以下、当項では上記の問題に対しortoolsのサンプルコードを実行します。

from ortools.init.python import init
from ortools.linear_solver import pywraplp


def main():
    print("Google OR-Tools version:", init.OrToolsVersion.version_string())

    # Create the linear solver with the GLOP backend.
    solver = pywraplp.Solver.CreateSolver("GLOP")
    if not solver:
        print("Could not create solver GLOP")
        return

    # Create the variables x and y.
    x_var = solver.NumVar(0, 1, "x")
    y_var = solver.NumVar(0, 2, "y")

    print("Number of variables =", solver.NumVariables())

    infinity = solver.infinity()
    # Create a linear constraint, x + y <= 2.
    constraint = solver.Constraint(-infinity, 2, "ct")
    constraint.SetCoefficient(x_var, 1)
    constraint.SetCoefficient(y_var, 1)

    print("Number of constraints =", solver.NumConstraints())

    # Create the objective function, 3 * x + y.
    objective = solver.Objective()
    objective.SetCoefficient(x_var, 3)
    objective.SetCoefficient(y_var, 1)
    objective.SetMaximization()

    print(f"Solving with {solver.SolverVersion()}")
    result_status = solver.Solve()

    print(f"Status: {result_status}")
    if result_status != pywraplp.Solver.OPTIMAL:
        print("The problem does not have an optimal solution!")
        if result_status == pywraplp.Solver.FEASIBLE:
            print("A potentially suboptimal solution was found")
        else:
            print("The solver could not solve the problem.")
            return

    print("Solution:")
    print("Objective value =", objective.Value())
    print("x =", x_var.solution_value())
    print("y =", y_var.solution_value())

    print("Advanced usage:")
    print(f"Problem solved in {solver.wall_time():d} milliseconds")
    print(f"Problem solved in {solver.iterations():d} iterations")


if __name__ == "__main__":
    init.CppBridge.init_logging("basic_example.py")
    cpp_flags = init.CppFlags()
    cpp_flags.stderrthreshold = True
    cpp_flags.log_prefix = False
    init.CppBridge.set_flags(cpp_flags)
    main()

・実行結果

Google OR-Tools version: 9.12.4544
Number of variables = 2
Number of constraints = 1
Solving with Glop solver v9.12.4544
Status: 0
Solution:
Objective value = 4.0
x = 1.0
y = 1.0
Advanced usage:
Problem solved in 0 milliseconds
Problem solved in 0 iterations

プログラムの理解にあたってはpywraplp.Solver.CreateSolver("GLOP")を用いてsolverオブジェクトを作成し、下記のように変数(Variables)や制約(Constraint)、目的関数(Objective)などを定義している点に着目しておくと良いです。

・変数(Variables)
→ x_var = solver.NumVar(0, 1, "x")のようにsolver.NumVarを用いて変数を定義

・制約(Constraint)
→ constraint = solver.Constraint(-infinity, 2, "ct")のようにconstraintオブジェクトを作成した後にconstraint.SetCoefficient(x_var, 1)のようにそれぞれの変数の制約を追加

・目的関数(Objective)
→ objective = solver.Objective()のようにobjectiveオブジェクトを作成した後にobjective.SetCoefficient(x_var, 3)のように変数に対し係数を定義

上記のように変数や制約、目的関数を定義した後にresult_status = solver.Solve()を元に最適化を実行します。最適化の結果は下記などを元に確認することができます。

・objective.Value() → 最適解における目的関数の値
・x_var.solution_value() → 最適解における変数の値

ネットワークフロー

最大フロー問題(Max Flow Problem)

最大フロー問題(Max Flow Problem)はネットワークのそれぞれのエッジのキャパシティが与えられた際にある点(start)からある点(goal)に流れる最大の流量を導出する問題です。OR-Toolsで最大フロー問題を解くにあたってはortools.graph.python.max_flow.SimpleMaxFlowクラスを用います。

"""From Taha 'Introduction to Operations Research', example 6.4-2."""
import numpy as np

from ortools.graph.python import max_flow


def main():
    """MaxFlow simple interface example."""
    # Instantiate a SimpleMaxFlow solver.
    smf = max_flow.SimpleMaxFlow()

    # Define three parallel arrays: start_nodes, end_nodes, and the capacities
    # between each pair. For instance, the arc from node 0 to node 1 has a
    # capacity of 20.
    start_nodes = np.array([0, 0, 0, 1, 1, 2, 2, 3, 3])
    end_nodes = np.array([1, 2, 3, 2, 4, 3, 4, 2, 4])
    capacities = np.array([20, 30, 10, 40, 30, 10, 20, 5, 20])

    # Add arcs in bulk.
    #   note: we could have used add_arc_with_capacity(start, end, capacity)
    all_arcs = smf.add_arcs_with_capacity(start_nodes, end_nodes, capacities)

    # Find the maximum flow between node 0 and node 4.
    status = smf.solve(0, 4)

    if status != smf.OPTIMAL:
        print("There was an issue with the max flow input.")
        print(f"Status: {status}")
        exit(1)
    print("Max flow:", smf.optimal_flow())
    print("")
    print(" Arc    Flow / Capacity")
    solution_flows = smf.flows(all_arcs)
    for arc, flow, capacity in zip(all_arcs, solution_flows, capacities):
        print(f"{smf.tail(arc)} / {smf.head(arc)}   {flow:3}  / {capacity:3}")
    print("Source side min-cut:", smf.get_source_side_min_cut())
    print("Sink side min-cut:", smf.get_sink_side_min_cut())


if __name__ == "__main__":
    main()

プログラムの把握にあたっては「start_nodes→end_nodesのエッジの容量」をcapacitiesで与えている点に着目すると良いです。

・実行結果

Max flow: 60

 Arc    Flow / Capacity
0 / 1    20  /  20
0 / 2    30  /  30
0 / 3    10  /  10
1 / 2     0  /  40
1 / 4    20  /  30
2 / 3    10  /  10
2 / 4    20  /  20
3 / 2     0  /   5
3 / 4    20  /  20
Source side min-cut: [0]
Sink side min-cut: [4, 1]

最適化の結果の出力にあたってはsmf.optimal_flow()で最大流量、solution_flows = smf.flows(all_arcs)によってそれぞれのエッジのFlowが得られる点に着目しておくと良いと思います。

最小費用フロー問題(Minimum-cost Flow Problem)

最小費用フロー問題(Minimum-cost Flow Problem)はネットワークのそれぞれのエッジの容量と輸送費用が与えられた際にコストが最小となる輸送フローを見つける問題です。OR-Toolsで最小費用フロー問題を解くにあたてはortools.graph.python.min_cost_flow.SimpleMinCostFlowクラスを用います。

"""From Bradley, Hax and Maganti, 'Applied Mathematical Programming', figure 8.1."""
import numpy as np

from ortools.graph.python import min_cost_flow


def main():
    """MinCostFlow simple interface example."""
    # Instantiate a SimpleMinCostFlow solver.
    smcf = min_cost_flow.SimpleMinCostFlow()

    # Define four parallel arrays: sources, destinations, capacities,
    # and unit costs between each pair. For instance, the arc from node 0
    # to node 1 has a capacity of 15.
    start_nodes = np.array([0, 0, 1, 1, 1, 2, 2, 3, 4])
    end_nodes = np.array([1, 2, 2, 3, 4, 3, 4, 4, 2])
    capacities = np.array([15, 8, 20, 4, 10, 15, 4, 20, 5])
    unit_costs = np.array([4, 4, 2, 2, 6, 1, 3, 2, 3])

    # Define an array of supplies at each node.
    supplies = [20, 0, 0, -5, -15]

    # Add arcs, capacities and costs in bulk using numpy.
    all_arcs = smcf.add_arcs_with_capacity_and_unit_cost(
        start_nodes, end_nodes, capacities, unit_costs
    )

    # Add supply for each nodes.
    smcf.set_nodes_supplies(np.arange(0, len(supplies)), supplies)

    # Find the min cost flow.
    status = smcf.solve()

    if status != smcf.OPTIMAL:
        print("There was an issue with the min cost flow input.")
        print(f"Status: {status}")
        exit(1)
    print(f"Minimum cost: {smcf.optimal_cost()}")
    print("")
    print(" Arc    Flow / Capacity Cost")
    solution_flows = smcf.flows(all_arcs)
    costs = solution_flows * unit_costs
    for arc, flow, cost in zip(all_arcs, solution_flows, costs):
        print(
            f"{smcf.tail(arc):1} -> {smcf.head(arc)}  {flow:3}  / {smcf.capacity(arc):3}       {cost}"
        )


if __name__ == "__main__":
    main()

プログラムの把握にあたっては前項で取り扱った最大フロー問題と同様に「start_nodes→end_nodes」のエッジの容量をcapacities、単位量あたりの輸送コストをunit_costsで与えている点に着目しておくと良いです。

Minimum cost: 150

 Arc    Flow / Capacity Cost
0 -> 1   12  /  15       48
0 -> 2    8  /   8       32
1 -> 2    8  /  20       16
1 -> 3    4  /   4       8
1 -> 4    0  /  10       0
2 -> 3   12  /  15       12
2 -> 4    4  /   4       12
3 -> 4   11  /  20       22
4 -> 2    0  /   5       0

また、最適化計算の出力にあたってはsmcf.optimal_cost()で最小費用、solution_flows = smcf.flows(all_arcs)で流量が得られることを確認しておくと良いと思います。

割当問題(Assignment Problem)の最適化

基本的な割当問題(Assignment Problem)

上記のようにワーカーがタスクに対応する際のコストが与えられる際のタスクの割当の最適化(コストの最小化)について以下確認します。OR-Toolsではこのような問題の最適化にあたってはortools.linear_solver.pywraplp.Solver.CreateSolverクラスを用います。

from ortools.linear_solver import pywraplp


def main():
    # Data
    costs = [
        [90, 80, 75, 70],
        [35, 85, 55, 65],
        [125, 95, 90, 95],
        [45, 110, 95, 115],
        [50, 100, 90, 100],
    ]
    num_workers = len(costs)
    num_tasks = len(costs[0])

    # Solver
    # Create the mip solver with the SCIP backend.
    solver = pywraplp.Solver.CreateSolver("SCIP")

    if not solver:
        return

    # Variables
    # x[i, j] is an array of 0-1 variables, which will be 1
    # if worker i is assigned to task j.
    x = {}
    for i in range(num_workers):
        for j in range(num_tasks):
            x[i, j] = solver.IntVar(0, 1, "")

    # Constraints
    # Each worker is assigned to at most 1 task.
    for i in range(num_workers):
        solver.Add(solver.Sum([x[i, j] for j in range(num_tasks)]) <= 1)

    # Each task is assigned to exactly one worker.
    for j in range(num_tasks):
        solver.Add(solver.Sum([x[i, j] for i in range(num_workers)]) == 1)

    # Objective
    objective_terms = []
    for i in range(num_workers):
        for j in range(num_tasks):
            objective_terms.append(costs[i][j] * x[i, j])
    solver.Minimize(solver.Sum(objective_terms))

    # Solve
    print(f"Solving with {solver.SolverVersion()}")
    status = solver.Solve()

    # Print solution.
    if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
        print(f"Total cost = {solver.Objective().Value()}\n")
        for i in range(num_workers):
            for j in range(num_tasks):
                # Test if x[i,j] is 1 (with tolerance for floating point arithmetic).
                if x[i, j].solution_value() > 0.5:
                    print(f"Worker {i} assigned to task {j}." + f" Cost: {costs[i][j]}")
    else:
        print("No solution found.")


if __name__ == "__main__":
    main()

プログラムの理解にあたってはcostsでタスクとワーカーのコストの表を定義している点にまず着目すると良いです。

・実行結果

Total cost = 265.0

Worker 0 assigned to task 3. Cost: 70
Worker 1 assigned to task 2. Cost: 55
Worker 2 assigned to task 1. Cost: 95
Worker 3 assigned to task 0. Cost: 45

実行結果からワーカーとタスクの対応とそれぞれのコストや総コストが確認できます。

包装(Packing)の最適化

ナップサック問題(Knapsack Problem)

OR-Toolsでナップサック問題(Knapsack Problem)を解くにあたってはortools.algorithms.python.knapsack_solver.KnapsackSolverクラスに基づくsolverを用います。

from ortools.algorithms.python import knapsack_solver


def main():
    # Create the solver.
    solver = knapsack_solver.KnapsackSolver(
        knapsack_solver.SolverType.KNAPSACK_MULTIDIMENSION_BRANCH_AND_BOUND_SOLVER,
        "KnapsackExample",
    )

    values = [
        # fmt:off
      360, 83, 59, 130, 431, 67, 230, 52, 93, 125, 670, 892, 600, 38, 48, 147,
      78, 256, 63, 17, 120, 164, 432, 35, 92, 110, 22, 42, 50, 323, 514, 28,
      87, 73, 78, 15, 26, 78, 210, 36, 85, 189, 274, 43, 33, 10, 19, 389, 276,
      312
        # fmt:on
    ]
    weights = [
        # fmt: off
      [7, 0, 30, 22, 80, 94, 11, 81, 70, 64, 59, 18, 0, 36, 3, 8, 15, 42, 9, 0,
       42, 47, 52, 32, 26, 48, 55, 6, 29, 84, 2, 4, 18, 56, 7, 29, 93, 44, 71,
       3, 86, 66, 31, 65, 0, 79, 20, 65, 52, 13],
        # fmt: on
    ]
    capacities = [850]

    solver.init(values, weights, capacities)
    computed_value = solver.solve()

    packed_items = []
    packed_weights = []
    total_weight = 0
    print("Total value =", computed_value)
    for i in range(len(values)):
        if solver.best_solution_contains(i):
            packed_items.append(i)
            packed_weights.append(weights[0][i])
            total_weight += weights[0][i]
    print("Total weight:", total_weight)
    print("Packed items:", packed_items)
    print("Packed_weights:", packed_weights)


if __name__ == "__main__":
    main()

・実行結果

Total value = 7534
Total weight: 850
Packed items: [0, 1, 3, 4, 6, 10, 11, 12, 14, 15, 16, 17, 18, 19, 21, 22, 24, 27, 28, 29, 30, 31, 32, 34, 38, 39, 41, 42, 44, 47, 48, 49]
Packed_weights: [7, 0, 22, 80, 11, 59, 18, 0, 3, 8, 15, 42, 9, 0, 47, 52, 26, 6, 29, 84, 2, 4, 18, 7, 71, 3, 66, 31, 0, 65, 52, 13]

プログラムでは荷物の価値に対応するvalueと荷物の重量に対応するweightsを定義した上でナップサックの最大の重量(容量)に対応するcapacitiesを定義し、solver.init(values, weights, capacities)やcomputed_value = solver.solve()を実行することで最適化を行います。

また、結果の出力にあたってはsolver.best_solution_contains(i)で計算した最適解(基本的には近似解)にi番目の荷物が含まれるかを調べられるので、プログラムではsolver.best_solution_containsを活用して得られた解の出力を行っている点に着目しておくと良いです。

ビンパッキング問題(Bin-Packing Problem)

OR-Toolsでビンパッキング問題(Bin-Packing Problem)を解くにあたってはortools.linear_solver.pywraplp.Solver.CreateSolverクラスに基づくsolverを用います。

from ortools.linear_solver import pywraplp


def create_data_model():
    """Create the data for the example."""
    data = {}
    weights = [48, 30, 19, 36, 36, 27, 42, 42, 36, 24, 30]
    data["weights"] = weights
    data["items"] = list(range(len(weights)))
    data["bins"] = data["items"]
    data["bin_capacity"] = 100
    return data


def main():
    data = create_data_model()

    # Create the mip solver with the SCIP(Solving Constraint Integer. Programs) backend.
    solver = pywraplp.Solver.CreateSolver("SCIP")

    if not solver:
        return

    # Variables
    # x[i, j] = 1 if item i is packed in bin j.
    x = {}
    for i in data["items"]:
        for j in data["bins"]:
            x[(i, j)] = solver.IntVar(0, 1, "x_%i_%i" % (i, j))

    # y[j] = 1 if bin j is used.
    y = {}
    for j in data["bins"]:
        y[j] = solver.IntVar(0, 1, "y[%i]" % j)

    # Constraints
    # Each item must be in exactly one bin.
    for i in data["items"]:
        solver.Add(sum(x[i, j] for j in data["bins"]) == 1)

    # The amount packed in each bin cannot exceed its capacity.
    for j in data["bins"]:
        solver.Add(
            sum(x[(i, j)] * data["weights"][i] for i in data["items"])
            <= y[j] * data["bin_capacity"]
        )

    # Objective: minimize the number of bins used.
    solver.Minimize(solver.Sum([y[j] for j in data["bins"]]))

    print(f"Solving with {solver.SolverVersion()}")
    status = solver.Solve()

    if status == pywraplp.Solver.OPTIMAL:
        num_bins = 0
        for j in data["bins"]:
            if y[j].solution_value() == 1:
                bin_items = []
                bin_weight = 0
                for i in data["items"]:
                    if x[i, j].solution_value() > 0:
                        bin_items.append(i)
                        bin_weight += data["weights"][i]
                if bin_items:
                    num_bins += 1
                    print("Bin number", j)
                    print("  Items packed:", bin_items)
                    print("  Total weight:", bin_weight)
                    print()
        print()
        print("Number of bins used:", num_bins)
        print("Time = ", solver.WallTime(), " milliseconds")
    else:
        print("The problem does not have an optimal solution.")


if __name__ == "__main__":
    main()

プログラムの理解にあたってはcreate_data_model関数のdata["weights"]で荷物の重量、data["bin_capacity"]でビンの容量(capacity)を指定していることに着目しておくと良いです。

・実行結果

Bin number 0
  Items packed: [0, 1, 2]
  Total weight: 97

Bin number 1
  Items packed: [3, 4, 5]
  Total weight: 99

Bin number 2
  Items packed: [6, 7]
  Total weight: 84

Bin number 3
  Items packed: [8, 9, 10]
  Total weight: 90


Number of bins used: 4
Time =  7  milliseconds

実行結果の理解にあたっては、「それぞれのビンのTotal weight」と「Items packedのインデックスに対応するdata["weights"]の和」の対応を確認しておくと良いと思います。

ルート最適化

セールスマン巡回問題(TSP; Travelling Salesperson Problem)

以下、当項ではortoolsのドキュメントのセールスマン巡回問題(TSP; Travelling Salesperson Problem)のプログラムを確認します。

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"] = [
        [0, 2451, 713, 1018, 1631, 1374, 2408, 213, 2571, 875, 1420, 2145, 1972],
        [2451, 0, 1745, 1524, 831, 1240, 959, 2596, 403, 1589, 1374, 357, 579],
        [713, 1745, 0, 355, 920, 803, 1737, 851, 1858, 262, 940, 1453, 1260],
        [1018, 1524, 355, 0, 700, 862, 1395, 1123, 1584, 466, 1056, 1280, 987],
        [1631, 831, 920, 700, 0, 663, 1021, 1769, 949, 796, 879, 586, 371],
        [1374, 1240, 803, 862, 663, 0, 1681, 1551, 1765, 547, 225, 887, 999],
        [2408, 959, 1737, 1395, 1021, 1681, 0, 2493, 678, 1724, 1891, 1114, 701],
        [213, 2596, 851, 1123, 1769, 1551, 2493, 0, 2699, 1038, 1605, 2300, 2099],
        [2571, 403, 1858, 1584, 949, 1765, 678, 2699, 0, 1744, 1645, 653, 600],
        [875, 1589, 262, 466, 796, 547, 1724, 1038, 1744, 0, 679, 1272, 1162],
        [1420, 1374, 940, 1056, 879, 225, 1891, 1605, 1645, 679, 0, 1017, 1200],
        [2145, 357, 1453, 1280, 586, 887, 1114, 2300, 653, 1272, 1017, 0, 504],
        [1972, 579, 1260, 987, 371, 999, 701, 2099, 600, 1162, 1200, 504, 0],
    ]
    data["num_vehicles"] = 1
    data["depot"] = 0
    return data


def print_solution(manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()} miles")
    index = routing.Start(0)
    plan_output = "Route for vehicle 0:\n"
    route_distance = 0
    while not routing.IsEnd(index):
        plan_output += f" {manager.IndexToNode(index)} ->"
        previous_index = index
        index = solution.Value(routing.NextVar(index))
        route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
    plan_output += f" {manager.IndexToNode(index)}\n"
    print(plan_output)
    plan_output += f"Route distance: {route_distance}miles\n"


def main():
    """Entry point of the program."""
    # 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)


    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)

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    )

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

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


if __name__ == "__main__":
    main()

・実行結果

Objective: 7293 miles
Route for vehicle 0:
 0 -> 7 -> 2 -> 3 -> 4 -> 12 -> 6 -> 8 -> 1 -> 11 -> 10 -> 5 -> 9 -> 0

上記のプログラムの理解にあたっては下記に着目すると良いです。

・Create Data
→ create_data_model関数に基づいてdata["distance_matrix"]、data["num_vehicles"]、data["depot"]を定義します。data["depot"]はルートのスタートとゴールに対応します。

・Routing Model
→ manager = pywrapcp.RoutingIndexManager...とrouting = pywrapcp.RoutingModel(manager) でRouting Modelの作成を行っています。

・distance_callback
→ data["distance_matrix"]の値などを出力します。

・Search Parameters
→ 検索パラメータはsearch_parameters = pywrapcp.DefaultRoutingSearchParameters()やsearch_parameters.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)を用いて定義します。``search_parameters.first_solution_strategyでは初期値を探索する戦略を指定しており、PATH_CHEAPEST_ARC`を用いる場合は重みの少ないエッジを繰り返して初期ルートの作成を行います。

・Solver
→ solution = routing.SolveWithParameters(search_parameters)を実行することによってソルバーを用いて最適化を行います。

配送計画問題(VRP; Vehicles Routing Problem)

当項ではortoolsのドキュメントの配送計画問題(VRP; Vehicles Routing Problem)のプログラムについて確認します。

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"] = [
        # fmt: off
      [0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354, 468, 776, 662],
      [548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674, 1016, 868, 1210],
      [776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164, 1130, 788, 1552, 754],
      [696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822, 1164, 560, 1358],
      [582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708, 1050, 674, 1244],
      [274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628, 514, 1050, 708],
      [502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856, 514, 1278, 480],
      [194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320, 662, 742, 856],
      [308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662, 320, 1084, 514],
      [194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388, 274, 810, 468],
      [536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764, 730, 388, 1152, 354],
      [502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114, 308, 650, 274, 844],
      [388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194, 536, 388, 730],
      [354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0, 342, 422, 536],
      [468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536, 342, 0, 764, 194],
      [776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274, 388, 422, 764, 0, 798],
      [662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730, 536, 194, 798, 0],
        # fmt: on
    ]
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")
    max_route_distance = 0
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        route_distance = 0
        while not routing.IsEnd(index):
            plan_output += f" {manager.IndexToNode(index)} -> "
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id
            )
        plan_output += f"{manager.IndexToNode(index)}\n"
        plan_output += f"Distance of the route: {route_distance}m\n"
        print(plan_output)
        max_route_distance = max(route_distance, max_route_distance)
    print(f"Maximum of the route distances: {max_route_distance}m")



def main():
    """Entry point of the program."""
    # 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 Distance constraint.
    dimension_name = "Distance"
    routing.AddDimension(
        transit_callback_index,
        0,  # no slack
        3000,  # vehicle maximum travel distance
        True,  # start cumul to zero
        dimension_name,
    )
    distance_dimension = routing.GetDimensionOrDie(dimension_name)
    distance_dimension.SetGlobalSpanCostCoefficient(100)

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    )

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

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)
    else:
        print("No solution found !")


if __name__ == "__main__":
    main()

・実行結果

Objective: 177500
Route for vehicle 0:
 0 ->  9 ->  10 ->  2 ->  6 ->  5 -> 0
Distance of the route: 1712m

Route for vehicle 1:
 0 ->  16 ->  14 ->  8 -> 0
Distance of the route: 1484m

Route for vehicle 2:
 0 ->  7 ->  1 ->  4 ->  3 -> 0
Distance of the route: 1552m

Route for vehicle 3:
 0 ->  13 ->  15 ->  11 ->  12 -> 0
Distance of the route: 1552m

Maximum of the route distances: 1712m

実行結果は制約条件のvehicle maximum travel distanceを2000に変えることで改善することも合わせて抑えておくと良いと思います。上記のVRPのプログラムは基本的には前項のTSPのプログラムに類似しているので以下ではTSPからの差分を中心に確認します。

・create_data_model
→ create_data_model関数でdata["num_vehicles"] = 4としている点に着目しておくと良いです。

・AddDimension
→ 下記のように車両毎の最大値を設定している点に注意しておくと良いです。

dimension_name = "Distance"
routing.AddDimension(
    transit_callback_index,
    0,  # no slack
    3000,  # vehicle maximum travel distance
    True,  # start cumul to zero
    dimension_name,
)
distance_dimension = routing.GetDimensionOrDie(dimension_name)
distance_dimension.SetGlobalSpanCostCoefficient(100)

容量制約付き車両ルート決定問題(CVRP; Capacited Vehicles Routing Problem)

当項では容量制約付き車両ルート決定問題(CVRP)のプログラムについて確認します。TSP・VRPと類似するかつ下記に記載のあるプログラムなのでプログラム全体については省略します。

・実行結果

Objective: 6208
Route for vehicle 0:
 0 Load(0) ->  7 Load(8) ->  3 Load(10) ->  4 Load(14) ->  1 Load(15) ->  0 Load(15)
Distance of the route: 1552m
Load of the route: 15

Route for vehicle 1:
 0 Load(0) ->  14 Load(4) ->  16 Load(12) ->  10 Load(14) ->  9 Load(15) ->  0 Load(15)
Distance of the route: 1552m
Load of the route: 15

Route for vehicle 2:
 0 Load(0) ->  12 Load(2) ->  11 Load(3) ->  15 Load(11) ->  13 Load(15) ->  0 Load(15)
Distance of the route: 1552m
Load of the route: 15

Route for vehicle 3:
 0 Load(0) ->  8 Load(8) ->  2 Load(9) ->  6 Load(13) ->  5 Load(15) ->  0 Load(15)
Distance of the route: 1552m
Load of the route: 15

Total distance of all routes: 6208m
Total load of all routes: 60

この問題ではそれぞれの車両に容量の制約(15)、それぞれの拠点に需要(0〜8)を設定しているので出力の形式もそれに合わせて用意されています。TSPやVSRとの差分については下記に着目しておくと良いです。

def create_data_model():
    ...
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]
    
    ...
    
def main():
    ...
    # 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",
    )

distanceとdemandについてcallbackが用意されていますが、callbackが「関数1に引数で渡して関数1の中で実行される関数」であることと合わせて確認しておくと良いと思います。

集荷と配達

当項では上図のように集荷と配達の組が定義されている場合(同じ車両で集荷と配達の組を通る必要あり)の最適化のプログラムについて確認します。

・実行結果

Objective: 226116
Route for vehicle 0:
 0 ->  13 ->  15 ->  11 ->  12 -> 0
Distance of the route: 1552m

Route for vehicle 1:
 0 ->  5 ->  2 ->  10 ->  16 ->  14 ->  9 -> 0
Distance of the route: 2192m

Route for vehicle 2:
 0 ->  4 ->  3 -> 0
Distance of the route: 1392m

Route for vehicle 3:
 0 ->  7 ->  1 ->  6 ->  8 -> 0
Distance of the route: 1780m

Total Distance of all routes: 6916m

プログラムの全体は上記に記載があるので以下ではVRPとの差分について確認します。まず、create_data_model関数に下記のような属性が追加される点に着目しておくと良いです。

def create_data_model():
    ...
    data["pickups_deliveries"] = [
        [1, 6],
        [2, 10],
        [4, 3],
        [5, 9],
        [7, 8],
        [15, 11],
        [13, 12],
        [16, 14],
    ]

また、main関数に下記のように「受け取りと宅配のリクエストを定義」が追加される点も抑えておく必要があります。

def main():
    ...
    for request in data["pickups_deliveries"]:
        pickup_index = manager.NodeToIndex(request[0])
        delivery_index = manager.NodeToIndex(request[1])
        routing.AddPickupAndDelivery(pickup_index, delivery_index)
        routing.solver().Add(
            routing.VehicleVar(pickup_index) == routing.VehicleVar(delivery_index)
        )
        routing.solver().Add(
            distance_dimension.CumulVar(pickup_index)
            <= distance_dimension.CumulVar(delivery_index)
        )

時間枠付き車両ルート選択問題(TRPTW; Vehicle Routing Problem with Time Windows)

当項では上図のように定義される時間枠に制約のある車両ルート選択問題(TRPTW; Vehicle Routing Problem with Time Windows)のプログラムについて確認します。

・実行結果

Objective: 71
Route for vehicle 0:
0 Time(0,0) -> 9 Time(2,3) -> 14 Time(7,8) -> 16 Time(11,11) -> 0 Time(18,18)
Time of the route: 18min

Route for vehicle 1:
0 Time(0,0) -> 7 Time(2,4) -> 1 Time(7,11) -> 4 Time(10,13) -> 3 Time(16,16) -> 0 Time(24,24)
Time of the route: 24min

Route for vehicle 2:
0 Time(0,0) -> 12 Time(4,4) -> 13 Time(6,6) -> 15 Time(11,11) -> 11 Time(14,14) -> 0 Time(20,20)
Time of the route: 20min

Route for vehicle 3:
0 Time(0,0) -> 5 Time(3,3) -> 8 Time(5,5) -> 6 Time(7,7) -> 2 Time(10,10) -> 10 Time(14,14) -> 0 Time(20,20)
Time of the route: 20min

Total time of all routes: 82min

プログラムの全体は上記に記載があるので以下ではVRPとの差分について確認します。まず、create_data_model関数に下記のような属性が追加される点に着目しておくと良いです。

def create_data_model():
    ...
    data["time_windows"] = [
        (0, 5),  # depot
        (7, 12),  # 1
        (10, 15),  # 2
        (16, 18),  # 3
        (10, 13),  # 4
        (0, 5),  # 5
        (5, 10),  # 6
        (0, 4),  # 7
        (5, 10),  # 8
        (0, 3),  # 9
        (10, 16),  # 10
        (10, 15),  # 11
        (0, 5),  # 12
        (5, 10),  # 13
        (7, 8),  # 14
        (10, 15),  # 15
        (11, 15),  # 16
    ]

上記はそれぞれのノードの対応可能な時間枠(Time Window)を定義しています。たとえば拠点1では[7, 12]、拠点2では[10, 15]の時間枠で訪れる必要があります。次に時刻に関するcallbackについて確認します。

def time_callback(from_index, to_index):
    """Returns the travel time between the two nodes."""
    # Convert from routing variable Index to time matrix NodeIndex.
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return data["time_matrix"][from_node][to_node]

transit_callback_index = routing.RegisterTransitCallback(time_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

時刻のコールバック関数は上記のようにdata["time_matrix"]に基づくtime_callback関数で定義されています。基本的にはdata["distance_matrix"]に基づくdistance_callbackと同様に理解すると良いです。また、時間枠の制約は下記のように追加されます。

time = "Time"
routing.AddDimension(
    transit_callback_index,
    30,  # allow waiting time
    30,  # maximum time per vehicle
    False,  # Don't force start cumul to zero.
    time,
)
time_dimension = routing.GetDimensionOrDie(time)

# Add time window constraints for each location except depot.
for location_idx, time_window in enumerate(data["time_windows"]):
    if location_idx == data["depot"]:
        continue
    index = manager.NodeToIndex(location_idx)
    time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])

# Add time window constraints for each vehicle start node.
depot_idx = data["depot"]
for vehicle_id in range(data["num_vehicles"]):
    index = routing.Start(vehicle_id)
    time_dimension.CumulVar(index).SetRange(
        data["time_windows"][depot_idx][0], data["time_windows"][depot_idx][1]
    )
for i in range(data["num_vehicles"]):
    routing.AddVariableMinimizedByFinalizer(
        time_dimension.CumulVar(routing.Start(i))
    )
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i)))

各拠点の時間枠の都合上それぞれの拠点で待つことを考慮したallow waiting timeが付与されていることに着目しておくと良いと思います。

リソースの制約

当項ではリソースの制約のあるVRPTW(Vehicle Routing Problem with Time Windows)のプログラムについて確認します。

・実行結果

Objective: 71
Route for vehicle 0:
0 Time(5, 5) -> 8 Time(8, 8) -> 14 Time(11, 11) -> 16 Time(13, 13) -> 0 Time(20,20)
Time of the route: 20min

Route for vehicle 1:
0 Time(0, 0) -> 12 Time(4, 4) -> 13 Time(6, 6) -> 15 Time(11, 11) -> 11 Time(14, 14) -> 0 Time(20,20)
Time of the route: 20min

Route for vehicle 2:
0 Time(5, 5) -> 7 Time(7, 7) -> 1 Time(11, 11) -> 4 Time(13, 13) -> 3 Time(14, 14) -> 0 Time(25,25)
Time of the route: 25min

Route for vehicle 3:
0 Time(0, 0) -> 9 Time(2, 3) -> 5 Time(4, 5) -> 6 Time(6, 9) -> 2 Time(10, 12) -> 10 Time(14, 16) -> 0 Time(25,25)
Time of the route: 25min

Total time of all routes: 90min

プログラムの全体は上記に記載があるので以下ではVRPとの差分について確認します。まず、create_data_model関数に下記のような属性が追加される点に着目しておくと良いです。

def create_data_model():
    ...
    data["time_windows"] = [
        (0, 5),  # depot
        (7, 12),  # 1
        (10, 15),  # 2
        (5, 14),  # 3
        (5, 13),  # 4
        (0, 5),  # 5
        (5, 10),  # 6
        (0, 10),  # 7
        (5, 10),  # 8
        (0, 5),  # 9
        (10, 16),  # 10
        (10, 15),  # 11
        (0, 5),  # 12
        (5, 10),  # 13
        (7, 12),  # 14
        (10, 15),  # 15
        (5, 15),  # 16
    ]
    ...
    data["vehicle_load_time"] = 5
    data["vehicle_unload_time"] = 5
    data["depot_capacity"] = 2