import numpy as np
from geopy.distance import geodesic
    
# 車両数
nvehicle = 2
    
# 都市の数
ncity = 15
avg_cities_per_vehicle = ncity // nvehicle
    
# 乱数シードの固定
seed = 12345
rng = np.random.default_rng(seed)
    
# 各都市における配送需要（重量）をランダムに決定
demand = rng.integers(1, 100, size=ncity)
demand_max = np.max(demand)
demand_mean = demand.mean()
    
# 全体的な需要に合わせ、車両の積載可能量 Q を設定する。
capacity = int(demand_max) + int(demand_mean) * avg_cities_per_vehicle

# 座標の取り得る範囲を設定
lat_range = [35.7014, 35.968]
lon_range = [139.34, 140.04]
    
# デポと各都市の座標をランダムに決定
ind2coord = [
    (
        rng.uniform(lon_range[0], lon_range[1]),
        rng.uniform(lat_range[0], lat_range[1]),
    )
    for i in range(ncity + 1)
]
    
# 2都市間の座標距離行列 D
distance_matrix = np.array(
    [
        [geodesic(coord_i[::-1], coord_j[::-1]).m for coord_j in ind2coord]
        for coord_i in ind2coord
    ]
)

_colors = [
    "green",
    "orange",
    "blue",
    "red",
    "purple",
    "pink",
    "darkblue",
    "cadetblue",
    "darkred",
    "lightred",
    "darkgreen",
    "lightgreen",
    "lightblue",
    "darkpurple",
]

import folium
    
    
def plot_solution(coord: list[tuple], title: str, best_tour: dict = dict()):
    l = len(coord)
    center = [
        sum(lat for _, lat in coord) / l,
        sum(lon for lon, _ in coord) / l,
    ]
    m = folium.Map(center, tiles="OpenStreetMap", zoom_start=10)
    folium.Marker(
        location=coord[0][::-1],
        popup=f"depot",
        icon=folium.Icon(icon="car", prefix="fa"),
    ).add_to(m)
    
    _color = _colors[1]
    if best_tour:
        for k, tour in best_tour.items():
            _color = _colors[k % len(_colors)]
            for city in tour:
                if city == 0:
                    continue
    
                folium.Marker(
                    location=coord[city][::-1],
                    popup=f"person{k}",
                    icon=folium.Icon(
                        icon="user", prefix="fa", color="white", icon_color=_color
                    ),
                ).add_to(m)
            folium.vector_layers.PolyLine(  # type: ignore
                locations=[coord[city][::-1] for city in tour], color=_color, weight=3
            ).add_to(m)
    else:
        for k, node in enumerate(coord):
            if k == 0:
                continue
            folium.Marker(
                location=node[::-1],
                popup=f"customer{k}",
                icon=folium.Icon(
                    icon="user", prefix="fa", color="white", icon_color=_color
                ),
            ).add_to(m)
    
    title = f"<h4>{title}</h4>"
    m.get_root().html.add_child(folium.Element(title))  # type: ignore
    
    return m

# 埼玉・千葉・東京近辺
title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}"
plot_solution(ind2coord, title)

def upperbound_of_tour(capacity: int, demand: np.ndarray) -> int:
    max_tourable_cities = 0
    for w in sorted(demand):
        capacity -= w
        if capacity >= 0:
            max_tourable_cities += 1
        else:
            return max_tourable_cities
    return max_tourable_cities

from amplify import VariableGenerator
    
gen = VariableGenerator()
    
# 積載可能量から1台の車両が訪問できる都市の最大数
max_tourable_cities = upperbound_of_tour(capacity, demand)
    
x = gen.array("Binary", shape=(max_tourable_cities + 2, ncity + 1, nvehicle))

x[0, 1:, :] = 0
x[-1, 1:, :] = 0
    
x[0, 0, :] = 1
x[-1, 0, :] = 1

from amplify import one_hot
    
one_trip_constraints = one_hot(x[1:-1, :, :], axis=1)

one_visit_constraints = one_hot(x[1:-1, 1:, :], axis=(0, 2))

from amplify import einsum
    
weight_sums = einsum("j,ijk->ik", demand, x[1:-1, 1:, :])

from amplify import less_equal, ConstraintList
    
capacity_constraints: ConstraintList = less_equal(
    weight_sums,  # type: ignore
    capacity,
    axis=0,
    penalty_formulation="Relaxation",
)

from amplify import Poly, einsum
    
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
    
# 経路の総距離
objective: Poly = einsum("pq,ipk,iqk->", distance_matrix, x[:-1], x[1:])  # type: ignore

from amplify import Model
    
constraints = one_trip_constraints + one_visit_constraints + capacity_constraints
constraints *= np.max(distance_matrix)  # 重みの設定
    
model = Model(objective, constraints)

from amplify import FixstarsClient
from datetime import timedelta
    
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=2000)  # タイムアウト2秒
# client.token = "API トークンを入力してください"

from amplify import solve
    
    
result = solve(model, client)
if len(result) == 0:
    raise RuntimeError("Some of the constraints are not satisfied.")
    
x_values = result.best.values  # 目的関数の値が最も低い解の変数値の取り出し

# one-hot な変数テーブルを辞書に変換。key：車両インデックス, value：各車両が訪問した都市の順番が入ったリスト
def onehot2sequence(solution: np.ndarray) -> dict[int, list]:
    nvehicle = solution.shape[2]
    sequence = dict()
    for k in range(nvehicle):
        sequence[k] = np.where(solution[:, :, k])[1]
    return sequence

def process_sequence(sequence: dict[int, list]) -> dict[int, list]:
    new_seq = dict()
    for k, v in sequence.items():
        v = np.append(v, v[0])
        mask = np.concatenate(([True], np.diff(v) != 0))
        new_seq[k] = v[mask]
    return new_seq

solution = x.evaluate(x_values)  # 結果が入ったnumpy配列
sequence = onehot2sequence(
    solution
)  # one-hot な変数テーブルを辞書に変換。key：車両インデックス, value：各車両が訪問した都市の順番が入ったリスト
best_tour = process_sequence(sequence)  # 上の辞書からデポへの余計な訪問を取り除く
print(f"Cost: {result.solutions[0].objective}")  # 目的関数値を表示
print(*best_tour.items(), sep="\n")  # 求めた解を表示

title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}, cost={result.solutions[0].objective:.2f}"
plot_solution(ind2coord, title, best_tour)

$x_{i,j,1}$	depot	A	C	E	F
0(出発)	1	0	0	0	0
1	0	1	0	0	0
2	0	0	1	0	0
3	0	0	0	1	0
4	0	0	0	0	1
5	1	0	0	0	0
6	1	0	0	0	0
7	1	0	0	0	0
8	1	0	0	0	0
9(到着)	1	0	0	0	0

$x_{i,j,2}$	depot	B	G	I
0(出発)	1	0	0	0
1	1	0	0	0
2	1	0	0	0
3	1	0	0	0
4	0	1	0	0
5	0	0	1	0
6	0	0	0	1
7	1	0	0	0
8	1	0	0	0
9(到着)	1	0	0	0

$x_{i,j,3}$	depot	D	H
0(出発)	1	0	0
1	0	1	0
2	0	0	1
3	1	0	0
4	1	0	0
5	1	0	0
6	1	0	0
7	1	0	0
8	1	0	0
9(到着)	1	0	0

容量制約つき運搬経路問題 (Capacitated Vehicle Routing Problem, CVRP)¶

定式化¶

定数¶

変数¶

制約条件¶

目的関数¶

具体例を使った制約式の説明¶

実装¶

問題の作成・可視化¶

訪問都市の上限の計算¶

二次多項式のモデルの構築¶

制約条件1: 全ての車両はデポ（都市 0）から出発し、デポに到着する¶

制約条件2: 車両$k$が$i$番目に訪問する場所は一箇所だけ¶

制約条件3: 都市$j$はいずれかの車両によってちょうど一回だけ訪問される¶

制約条件4: 車両 $k$ が運ぶ荷物の総重量は $Q$ 以下¶

目的関数¶

クライアントの設定¶

結果取得・可視化¶