【Python】ortoolsを用いた最適化

ortoolsの概要とサンプルコードの実行

ortoolsの概要

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

ortoolsのインストール

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

$ pip install ortools

ortoolsのサンプルコードの実行

\begin{align}
\mathrm{Objective} \quad &: 3x+y \\
\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() → 最適解における変数の値

ルート最適化

セールスマン巡回問題(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; Travelling Salesperson Problem)