【数理最適化】Or-Toolsで配送計画を最適化してみる（容量・時間制約付き配送経路最適化）

# 目的
- 数理最適化を使用して最適な配送計画を作成してみる
- 記事確認いただいた初学者の方が以下の内容を達成できる
    - 数理最適化の流れを理解する
    - サンプルコードを通して実装方法を理解する

# 課題設定
私は運送会社の配車組担当。弊社では主に遊園地等の行楽地にメーカー請けの食品を輸送しており、卸先は全て近隣のため、当日出庫・当日納品を行っている。これまでは私の経験によって配送計画表を作成していたが、私が高齢になってきたことや、後進もおらず、個人的な都合をつけ辛いことから、これを自動化する仕組みがあれば活用を検討したい。
以下の情報から、移動時間と作業時間の総和が最小になるようにルート組みできないか。

＜出庫地＞
- 千葉の倉庫から荷物を積み、千葉・東京の各行楽地（配送先）へトラックで輸送する
- 出庫場所ではトラックへの積作業（60min）が発生

＜卸地＞
- 配送先には卸可能時間（X:XX~X:XX）が存在する
- 配送先では卸作業がそれぞれ発生する

＜手配可能車両＞
- トラックは4t×4台（積載限界量2.4t×2、2.8t×2）
- 積載限界量を超えては積めない
- いずれもAM7:30~ 8時間を超えての稼働は不可

# 参考
写経しつつ理解がお薦めです！
- [OR-Toolsを使って巡回セールスマン問題を解き、
効率的に沼津の聖地巡礼をする](https://jpn.pioneer/ja/piomatixlbs/useful/or-tools/)
- [第6回：OR-Toolsで巡回セールスマン問題を解く ～京都弾丸観光ツアーの作成を事例に～【ブレインパッドの数理最適化ブログ】](https://www.brainpad.co.jp/doors/contents/01_tech_2021-06-18-110000/)
- [Google公式ドキュメント](https://developers.google.com/optimization/routing/cvrp?hl=ja)

# 定式化
あまり数学に明るい方ではないですが（勉強中）、コスト関数と上記の条件を定式化してみます。

### コスト関数（目的関数）
```math
\displaylines{

t_{ij}:\quad Travel Time \quad and \quad Service Time\\
（i→j地点に移動した場合の移動時間＋i地点での作業時間）\\

x_{ij}:\quad 1 \quad if \quad moves,\quad 0 \quad otherwise：\\
（車両vがi→j地点に移動した場合に1、そうでない場合0を示す決定変数）\\

Min \sum_{v=1}^{V} \sum_{i=0}^{n} \sum_{j=0, j≠i}^{n}t_{ij}x_{ijv}\\
（移動時間+作業時間の総和を最小化）

}
```
### 制約
#### 車両の容量制約
```math
\displaylines{

q_{i}:\quad Order Quantity\\
（i地点に配送する荷量）\\

Q_{v}:\quad Vehicle Limit Quantity\\
（車両vの積載限界量）\\

\sum_{v=1}^{V} \sum_{i=0}^{n}q_{i}x_{ijv} ≤ Qv \\
（車両vに積載する量は車両vの積載限界量以下）

}
```

#### 配送先の到着時間枠制約
```math
\displaylines{

TwStart_{i} \quad TwEnd_{i}:\quad TimeWindow\\
（i地点の配送可能開始時間、終了時間）\\

ArriveAt_{iv}:\quad Vehicle Arrival Time\\
（車両vがi地点に到着する時間）\\

TwStart_{i} ≤ArriveAt_{iv}≤TwEnd_{i} \\
（車両vがi地点に到着するのは配送可能時間以上、終了時刻以下）

}
```

#### 車両の稼働時間制約
```math
\displaylines{

W_{v}:\quad WorkngTime\\
（車両vの最大稼働可能時間）\\

t_{ij}:\quad Travel Time \quad and \quad Service Time\\
（i→j地点に移動した場合の移動時間＋i地点での作業時間）\\

x_{ij}:\quad 1 \quad if \quad moves,\quad 0 \quad otherwise：\\
（車両vがi→j地点に移動した場合に1、そうでない場合0を示す決定変数）\\

\sum_{i=0}^{n} \sum_{j=0, j≠i}^{n}t_{ij}x_{ijv}≤W_{v}\\
（移動時間+作業時間の総和が車両vの最大稼働可能時間以下※一律値）

}
```

# 使用するライブラリ
```python
# 最適化ツール
!pip install ortools
# 位置情報から経路距離を求める
!pip install OSMnx
!pip install networkx

# インポート
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
import osmnx as ox
import networkx as nx
import numpy as np
import pandas as pd
import datetime
import warnings
warnings.simplefilter('ignore')
```

# データ作成
## 配送地点情報
「プロロジスパーク千葉」が積み場所、他は配送先になります。
配送先はなんとなく馴染みのある場所にしてみました。
```python
locations_demand = pd.DataFrame({'name':['プロロジスパーク千葉',
                                         '東京ドイツ村',
                                         'TDL',
                                         'TDS',
                                         '東京ソラマチ',
                                         '浅草花やしき',
                                         'ふなばしアンデルセン公園',
                                         '東京ドームシティ'],
                                'address':['千葉県千葉市稲毛区六方町２１０−２７',
                                           '千葉県袖ケ浦市永吉４１９',
                                           '千葉県浦安市舞浜１−１',
                                           '千葉県浦安市舞浜１−１３',
                                           '東京都墨田区押上１丁目１−２',
                                           '東京都台東区浅草２丁目２８−１',
                                           '千葉県船橋市金堀町５２５',
                                           '東京都文京区後楽１丁目３−６１'],
                                # 緯度経度
                                'lat':[35.65752312, 35.40588803, 35.63307077, 35.62710749, 35.71039879, 35.71560929, 35.76070174, 35.70553601],
                                'lon':[140.1386453, 140.0540462, 139.8803836, 139.885271 , 139.8115532, 139.7948823, 140.0616562, 139.7520476],
                                # 配送可能時間
                                'open':['6:00', '10:00', '8:00', '8:00', '7:30', '8:00', '8:00', '8:00'],
                                'close':['22:00', '10:00', '12:00', '12:00', '13:00', '11:00', '15:00', '15:00'],
                                # 地点での作業時間
                                'service_sec':[3600, 1800, 3600, 3600, 1800, 1800, 1800, 1800],
                                # 地点への配送荷量（g）
                                'demand_g':[0, 1000000, 2300000, 2400000, 1200000, 100000, 150000, 2800000]})

# openとcloseをsecondsに変換
locations_demand['open_time_sec'] = 0
locations_demand['close_time_sec'] = 0

for i in range(len(locations_demand)):
    opening = locations_demand['open'].loc[i]
    open_time_sec =3600*int(opening[ : opening.find(':')]) + 60*int(opening[ opening.find(':') + 1: ])
    
    closing = locations_demand['close'].loc[i]
    close_time_sec =3600*int(closing[ : closing.find(':')]) + 60*int(closing[ closing.find(':') + 1: ])
    
    locations_demand['open_time_sec'].loc[i] = open_time_sec
    locations_demand['close_time_sec'].loc[i] = close_time_sec

locations_demand
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2642161/b94e2be0-a435-5c29-b92c-ca79d211f5ad.png)

<img width="500" alt="image" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2642161/2e8eef90-8210-606f-ebb6-5fd8e6dc40e4.png">

## 各地点間の距離と移動時間の算出
```python
def culc_dist_durations(locations_demand):
    
    # 距離計算
    # https://osmnx.readthedocs.io/en/stable/user-reference.html#module-osmnx.graph
    G = ox.graph.graph_from_address(address = '千葉県千葉市稲毛区', 
                                   dist = 50000, 
                                   dist_type = 'network', 
                                   network_type = 'drive', 
                                   simplify = True,
                                   retain_all = False, 
                                   truncate_by_edge = False, 
                                   custom_filter = None)

    durations = []
    for i in locations_demand.index:
        start_point_lat = locations_demand['lat'].loc[i]
        start_point_lon = locations_demand['lon'].loc[i]
        start_node = ox.distance.nearest_nodes(G, start_point_lon,start_point_lat)

        durations_tmp = []
        for j in locations_demand.index:
            if i != j:

                end_point_lat = locations_demand['lat'].loc[j]
                end_point_lon = locations_demand['lon'].loc[j]
                end_node = ox.distance.nearest_nodes(G, end_point_lon,end_point_lat)

                dist_km = (nx.shortest_path_length(G, start_node, end_node, weight='length')/1000)

                # 経路距離を時速25kmで換算し移動時間とする
                dur_sec = int((3600*(dist_km/25)).round(0))
            else:
                dur_sec = 0

            durations_tmp.append(dur_sec)

        durations.append(durations_tmp)

    # グラフを初期化
    G=''

    # 移動時間を返す
    return durations
```
## インプットデータに整形
```python
# 車両の稼働開始時刻
start_time_at = '7:30'
start_time_at_sec = 3600*int(start_time_at[ : start_time_at.find(':')]) + 60*int(start_time_at[start_time_at.find(':') + 1 : ])

# 全車両に対し適用される限界稼働時間sec
time_limit_sec = 28800 # 8h=28800sで全箇所回り切れるか

# 許容する待機時間sec
slack_time_sec = 7200
list_name = locations_demand['name'].to_list() 
service_time = locations_demand['service_sec'].to_list() 


# 稼働開始がデポの開始時刻より遅い時
depot_open_time_at_sec = locations_demand['open_time_sec'].loc[0]
if start_time_at_sec >= depot_open_time_at_sec:
    
    open_close_time_tmp = [tuple(x) for x in locations_demand[['open_time_sec', 'close_time_sec']].values - start_time_at_sec]
    open_close_time = []
    for x,y in open_close_time_tmp:
        if x < 0:
            x = 0
        if y < 0:
            y = 0
        add = (x, y)
        open_close_time.append(add)

# 稼働開始がデポの開始時刻より早い時
else:
    open_close_time_tmp = [tuple(x) for x in locations_demand[['open_time_sec', 'close_time_sec']].values - depot_open_time_at_sec]
    open_close_time = []
    for x,y in open_close_time_tmp:
        if x < 0:
            x = 0
        if y < 0:
            y = 0
        add = (x, y)
        open_close_time.append(add)
    
data = {}
data['name'] = list_name

# data['time_matrix'] = culc_dist_durations(locations_demand)
# culc_dist_durations関数は情報取得に時間かかるのでコメントアウト
# 以下、上記のculc_dist_durationsの実行結果
data['time_matrix'] = [
  [0, 4856, 4071, 3985, 4601, 4901, 2361, 5365],
  [4857, 0, 7892, 7806, 8571, 8871, 6843, 9335],
  [4088, 7913, 0, 123, 2095, 2395, 3667, 2647],
  [4001, 7827, 122, 0, 2047, 2340, 3589, 2561],
  [4590, 8562, 2065, 2050, 0, 302, 3868, 865],
  [4855, 8827, 2327, 2311, 300, 0, 4107, 716],
  [2360, 6844, 3660, 3574, 3883, 4109, 0, 4667],
  [5335, 9307, 2588, 2573, 868, 703, 4614, 0]]

data['time_windows'] = open_close_time
data['service_time'] = service_time

data['depot'] = 0 # depotのインデックス
data['depot_open_time'] = locations_demand['open'].loc[0] # depotの開始時刻
data['depot_open_time_sec'] = depot_open_time_at_sec # depotの開始時刻 seconds表示
data['slack_time'] = slack_time_sec # 許容する待機時間
data['time_limit'] = time_limit_sec # 各車両の稼働限界時間
data['start_time'] = start_time_at # 車両の稼働開始時刻
data['start_time_sec'] = start_time_at_sec # 車両の稼働開始時刻 seconds表示

# 各地点への配送荷量
data['demands'] = list(locations_demand['demand_g'])
    
# 車両
data['vehicle_capacities'] = [2400000, 2400000, 2800000, 2800000] # 各車両の限界積載重量（g）
data['vehicles'] = 4 # 4台の車両を使用

data
```
data↓
```
{'name': ['プロロジスパーク千葉',
  '東京ドイツ村',
  'TDL',
  'TDS',
  '東京ソラマチ',
  '浅草花やしき',
  'ふなばしアンデルセン公園',
  '東京ドームシティ'],
 'time_matrix': [[0, 4856, 4071, 3985, 4601, 4901, 2361, 5365],
  [4857, 0, 7892, 7806, 8571, 8871, 6843, 9335],
  [4088, 7913, 0, 123, 2095, 2395, 3667, 2647],
  [4001, 7827, 122, 0, 2047, 2340, 3589, 2561],
  [4590, 8562, 2065, 2050, 0, 302, 3868, 865],
  [4855, 8827, 2327, 2311, 300, 0, 4107, 716],
  [2360, 6844, 3660, 3574, 3883, 4109, 0, 4667],
  [5335, 9307, 2588, 2573, 868, 703, 4614, 0]],
 'time_windows': [(0, 52200),
  (9000, 9000),
  (1800, 16200),
  (1800, 16200),
  (0, 19800),
  (1800, 12600),
  (1800, 27000),
  (1800, 27000)],
 'service_time': [3600, 1800, 3600, 3600, 1800, 1800, 1800, 1800],
 'depot': 0,
 'depot_open_time': '6:00',
 'depot_open_time_sec': 21600,
 'slack_time': 7200,
 'time_limit': 28800,
 'star（t_time': '7:30',
 'start_time_sec': 27000,
 'demands': [0, 1000000, 2300000, 2400000, 1200000, 100000, 150000, 2800000],
 'vehicle_capacities': [2400000, 2400000, 2800000, 2800000],
 'vehicles': 4}
 ```

# 最適化モデルを作成
### Routing Managerの作成
```python
manager = pywrapcp.RoutingIndexManager(
    len(data['time_matrix']),  # 訪問したい地点の数
    data['vehicles'],          # 使う車両の数
    data['depot']              # depotのindex
    )
```
### Routing Modelを作成（以下「モデル」）
```python
routing = pywrapcp.RoutingModel(manager) 
```
### 移動時間＋作業時間を返す関数をモデルに登録
```python
def time_callback(from_index, to_index):
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    travel_time = data['time_matrix'][from_node][to_node]
    service_time = data['service_time'][from_node]
    # 移動時間と地点の作業時間を足し合わせる
    return travel_time + service_time

transit_callback_index = routing.RegisterTransitCallback(time_callback)
```
### コストの設定
```python
# ArcCostの設定（本件では移動時間＋作業時間）
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index) 
```
### 配送先への荷量（＝需要量の設定）を返す関数をモデルに登録
```python
def demand_callback(from_index):
    from_node = manager.IndexToNode(from_index)
    return data['demands'][from_node]

demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
```
### 容量ディメンションを追加
```python
capacity = 'Capacity'
routing.AddDimensionWithVehicleCapacity(
    demand_callback_index,
    0,
    data['vehicle_capacities'],
    True,
    capacity)
```
### 時間枠ディメンションを追加
```python
time = 'Time'
routing.AddDimension(
    transit_callback_index,
    data['slack_time'],
    data['time_limit'], 
    False,
    time)
time_dimension = routing.GetDimensionOrDie(time)
```
### 各地点について時間枠の制約を追加する
```python
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(int(time_window[0]), int(time_window[1]))
```
### 稼働開始時間の時間枠の制約を各車両に追加する（depotの時間内に開始）
```python
depot_idx = data['depot']
for vehicle_id in range(data['vehicles']):
    index = routing.Start(vehicle_id)
    time_dimension.CumulVar(index).SetRange(
        int(data['time_windows'][depot_idx][0]), int(data['time_windows'][depot_idx][1]))

for i in range(data['vehicles']):
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.Start(i)))
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i))) 
```
### 地点に行かなかった（制約上、行けなかった）場合に加算されるコスト（ペナルティ）
```python
penalty = 10000000
for node in range(1, len(data['time_matrix'])):
    routing.AddDisjunction([manager.NodeToIndex(node)], penalty)
```
### 計算のパラメータ
```python
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.local_search_metaheuristic = (
    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH) # 指定時間で厳密解を探索

search_parameters.time_limit.seconds = 10 # 計算時間を10sで指定
```
### 結果出力
```python
 # 結果の表示
def print_solution(data, manager, routing, solution):
    time_dimension = routing.GetDimensionOrDie('Time')
    total_time = 0
    total_load = 0
    for vehicle_id in range(data['vehicles']):
        print('-------------------------------------------------------------------------------------')
        index = routing.Start(vehicle_id)
        _vehicle_capacity_kg = int(data['vehicle_capacities'][vehicle_id]/1000)
        print('Route for vehicle {0}  Limit Capacity {1}kg:'.format(vehicle_id, _vehicle_capacity_kg))
        total_time_route = 0
        total_load_route = 0
        if data['start_time_sec'] >= data['depot_open_time_sec']:
            
            while not routing.IsEnd(index):
                time_var = time_dimension.CumulVar(index)
                h = int(data['start_time'][ : data['start_time'].find(':')])
                m = int(data['start_time'][data['start_time'].find(':') + 1 : ])
                print( '(到着{1}  出発{2})  {0}'.format(
                    data['name'][manager.IndexToNode(index)],
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var) + data['service_time'][manager.IndexToNode(index)])
                )
                )
                total_load_route += data['demands'][manager.IndexToNode(index)]

                index = solution.Value(routing.NextVar(index))

            time_var = time_dimension.CumulVar(index)
            h = int(data['start_time'][ : data['start_time'].find(':')])
            m = int(data['start_time'][data['start_time'].find(':') + 1 : ])
            print('(到着{1}  出発{2})  {0}'.format(
                data['name'][manager.IndexToNode(index)],
                datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var)) # 戻り時は作業時間を加算しない
            )
            )
        
        else:
            while not routing.IsEnd(index):
                time_var = time_dimension.CumulVar(index)
                h = int(data['depot_open_time'][ : data['depot_open_time'].find(':')])
                m = int(data['depot_open_time'][data['depot_open_time'].find(':') + 1 : ])
                print( '(到着{1}  出発{2})  {0}'.format(
                    data['name'][manager.IndexToNode(index)],
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var) + data['service_time'][manager.IndexToNode(index)])
                )
                )
                total_load_route += data['demands'][manager.IndexToNode(index)]

                index = solution.Value(routing.NextVar(index))

            time_var = time_dimension.CumulVar(index)
            h = int(data['depot_open_time'][ : data['depot_open_time'].find(':')])
            m = int(data['depot_open_time'][data['depot_open_time'].find(':') + 1 : ])
            print('(到着{1}  出発{2})  {0}'.format(
                data['name'][manager.IndexToNode(index)],
                datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var)) # 戻り時は作業時間を加算しない
            )
            )

        total_time_route = solution.Min(time_var)
        total_time += total_time_route
        print(total_time_route)

        total_load += total_load_route
        print('Total time of route: {}'.format(datetime.timedelta(seconds=total_time_route)))
        print('Total load of route: {}kg'.format((total_load_route)/1000)) 

    print('-------------------------------------------------------------------------------------')
    print(f'Objective: {solution.ObjectiveValue()}')
    print('Total time of all routes: {}'.format(datetime.timedelta(seconds=total_time))) 
```
### 最適化モデル計算実行
```python
solution = routing.SolveWithParameters(search_parameters)

# 解が見つかれば結果出力
if solution:
    print_solution(data, manager, routing, solution)
# 解が見つからぬ場合
else:
    print('No solution found.') 
```
### 出力結果
制約条件を満たしていることがわかります。
- いずれの車両も8h内の稼働時間に収まっている
- いずれの配送先にも時間枠内に到着できている
- 各車両の限界積載量を超えていない

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2642161/07dc8be4-72ef-1948-5de3-2c4d0e405b3a.png)
同色のラインが同一車両のルートです。
<img width="500" alt="image" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2642161/c8d82402-e2ea-7800-b8f4-45cb415682b2.png">

ちなみにこれを5h内の稼働時間という制約で実行すると・・・
```python
time_limit_sec = 18000
```
出力結果
![スクリーンショット 2024-10-16 175741.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2642161/a1b95acc-0d28-a6bf-4a0c-563fed5d10bb.png)

- Objectiveが2千万弱と二か所の配送先に行けていないことでペナルティ（2箇所×10000000）が加算されている
- 行けなかった配送先は「東京ドームシティ」「TDL」
- vehicle=3は積載残2650kg、稼働時間残2.2h程度。TDL分は積載できないのか？
    - ふなばしアンデルセン公園からTDLまでの移動時間は3660s、卸作業が3600s、TDLから発地までの移動時間は4088s、合計3.15hは必要になる。稼働時間の残時間では配送できないため、積載できない


## 最適化モデルコードまとめ
```python
def main(data):
    

    # Routing Managerを作成
    manager = pywrapcp.RoutingIndexManager(
        len(data['time_matrix']),  # 訪問したい地点の数
        data['vehicles'],          # 使う車両の数
        data['depot']              # depotのindex
    )

    # Routing Modelを作成
    routing = pywrapcp.RoutingModel(manager) 


    def time_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        travel_time = data['time_matrix'][from_node][to_node]
        service_time = data['service_time'][from_node]
        # 移動時間と配送先での滞在時間を足し合わせる
        return travel_time + service_time


    transit_callback_index = routing.RegisterTransitCallback(time_callback)

    # Arc Costの設定
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index) 

    # demandの設定
    def demand_callback(from_index):
        # 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)

    # 容量ディメンションを追加
    capacity = 'Capacity'
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # null capacity slack
        data['vehicle_capacities'],
        True,  # start cumul to zero
        capacity)
    
    # 時間枠ディメンションを追加
    time = 'Time'
    routing.AddDimension(
        transit_callback_index,
        data['slack_time'], 
        data['time_limit'], 
        False,               
        time)
    time_dimension = routing.GetDimensionOrDie(time)

    # 各地点について時間の制約を追加する
    for location_idx, time_window in enumerate(data['time_windows']):
        index = manager.NodeToIndex(location_idx)
        time_dimension.CumulVar(index).SetRange(int(time_window[0]), int(time_window[1]))

    # 稼働開始時間の時間枠の制約を各車両に追加する
    depot_idx = data['depot']
    for vehicle_id in range(data['vehicles']):
        index = routing.Start(vehicle_id)
        time_dimension.CumulVar(index).SetRange(
            int(data['time_windows'][depot_idx][0]),
            int(data['time_windows'][depot_idx][1])) 

    for i in range(data['vehicles']):
        routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.Start(i)))
        routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i))) 
        
    # 地点に行かなかった（制約上、行けなかった）場合に加算されるコスト（ペナルティ）
    penalty = 10000000
    for node in range(1, len(data['time_matrix'])):
        routing.AddDisjunction([manager.NodeToIndex(node)], penalty)

    # 計算時のパラメータ
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)

    search_parameters.time_limit.seconds = 10

     # 結果の表示
    def print_solution(data, manager, routing, solution):
        time_dimension = routing.GetDimensionOrDie('Time')
        total_time = 0
        total_load = 0
        for vehicle_id in range(data['vehicles']):
            print('-------------------------------------------------------------------------------------')
            index = routing.Start(vehicle_id)
            _vehicle_capacity_kg = int(data['vehicle_capacities'][vehicle_id]/1000)
            print('Route for vehicle {0}  Limit Capacity {1}kg:'.format(vehicle_id, _vehicle_capacity_kg))
            total_time_route = 0
            total_load_route = 0
            if data['start_time_sec'] >= data['depot_open_time_sec']:

                while not routing.IsEnd(index):
                    time_var = time_dimension.CumulVar(index)
                    h = int(data['start_time'][ : data['start_time'].find(':')])
                    m = int(data['start_time'][data['start_time'].find(':') + 1 : ])
                    print( '(到着{1}  出発{2})  {0}'.format(
                        data['name'][manager.IndexToNode(index)],
                        datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                        datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var) + data['service_time'][manager.IndexToNode(index)])
                    )
                    )
                    total_load_route += data['demands'][manager.IndexToNode(index)]

                    index = solution.Value(routing.NextVar(index))

                time_var = time_dimension.CumulVar(index)
                h = int(data['start_time'][ : data['start_time'].find(':')])
                m = int(data['start_time'][data['start_time'].find(':') + 1 : ])
                print('(到着{1}  出発{2})  {0}'.format(
                    data['name'][manager.IndexToNode(index)],
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var)) # 戻り時は作業時間を加算しないで表示
                )
                )

            else:
                while not routing.IsEnd(index):
                    time_var = time_dimension.CumulVar(index)
                    h = int(data['depot_open_time'][ : data['depot_open_time'].find(':')])
                    m = int(data['depot_open_time'][data['depot_open_time'].find(':') + 1 : ])
                    print( '(到着{1}  出発{2})  {0}'.format(
                        data['name'][manager.IndexToNode(index)],
                        datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                        datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var) + data['service_time'][manager.IndexToNode(index)])
                    )
                    )
                    total_load_route += data['demands'][manager.IndexToNode(index)]

                    index = solution.Value(routing.NextVar(index))

                time_var = time_dimension.CumulVar(index)
                h = int(data['depot_open_time'][ : data['depot_open_time'].find(':')])
                m = int(data['depot_open_time'][data['depot_open_time'].find(':') + 1 : ])
                print('(到着{1}  出発{2})  {0}'.format(
                    data['name'][manager.IndexToNode(index)],
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Min(time_var)),
                    datetime.timedelta(hours=h, minutes=m, seconds=solution.Max(time_var)) # 戻り時は作業時間を加算しないで表示
                )
                )

            total_time_route = solution.Min(time_var)
            total_time += total_time_route
            print(total_time_route)

            total_load += total_load_route
            print('Total time of route: {}'.format(datetime.timedelta(seconds=total_time_route)))
            print('Total load of route: {}kg'.format((total_load_route)/1000)) 

        print('-------------------------------------------------------------------------------------')
        print(f'Objective: {solution.ObjectiveValue()}')
        print('Total time of all routes: {}'.format(datetime.timedelta(seconds=total_time))) 


    # 最適化計算実行
    solution = routing.SolveWithParameters(search_parameters)

    if solution:
        print_solution(data, manager, routing, solution)
    else:
        print('No solution found.') 


main(data)
```