num_people = 30
num_cars = 10
car_capacity = 6

import numpy as np
from geopy.distance import geodesic
    
    
def generate_problem(
    num_people: int,
    num_cars: int,
    latitude_range: tuple[float, float],
    longitude_range: tuple[float, float],
    seed=999,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    rng = np.random.default_rng(seed)
    
    # 利用者と車の座標を生成
    people_coords = rng.uniform(
        low=np.array([latitude_range[0], longitude_range[0]]),
        high=np.array([latitude_range[1], longitude_range[1]]),
        size=(num_people, 2),
    )
    car_coords = rng.uniform(
        low=np.array([latitude_range[0], longitude_range[0]]),
        high=np.array([latitude_range[1], longitude_range[1]]),
        size=(num_cars, 2),
    )
    
    # 距離行列 d を作成 (d[i, k] が利用者 i と車 k の距離)
    distances = np.zeros((num_people, num_cars))
    for i in range(num_people):
        for k in range(num_cars):
            distances[i, k] = geodesic(
                people_coords[i], car_coords[k]
            ).m  # 単位はメートル
    
    return people_coords, car_coords, distances

# 船橋駅周辺
latitude_range = (35.675500, 35.76)
longitude_range = (139.9, 140.08)
    
people_coords, car_coords, distances = generate_problem(
    num_people, num_cars, latitude_range, longitude_range
)

import folium
    
    
def plot_people_and_cars(
    people_coords: np.ndarray,
    car_coords: np.ndarray,
    latitude_range: tuple[float, float],  # 描画範囲 (緯度)
    longitude_range: tuple[float, float],  # 描画範囲 (経度)
):
    m = folium.Map(
        [sum(latitude_range) / 2, sum(longitude_range) / 2],
        tiles="OpenStreetMap",
        zoom_start=12,
    )
    
    for latitude, longitude in people_coords:
        folium.Marker(
            location=(latitude, longitude),
            icon=folium.Icon(icon="user", prefix="fa", color="orange"),
        ).add_to(m)
    
    for latitude, longitude in car_coords:
        folium.Marker(
            location=(latitude, longitude),
            icon=folium.Icon(icon="car", prefix="fa", color="green"),
        ).add_to(m)
    
    return m

plot_people_and_cars(people_coords, car_coords, latitude_range, longitude_range)

average = distances.mean()
std = distances.std()
distances: np.ndarray = (distances - average) / std

from amplify import VariableGenerator
    
gen = VariableGenerator()
q = gen.array("Binary", num_people, num_cars)

from amplify import sum
    
alpha = 1.0
    
# 第一項: 利用者の移動距離を短くする
distance_objective = (distances * q).sum()
    
# 第二項: 乗車率の二乗和を大きくする
occupancies = q.sum(axis=0) / car_capacity  # 乗車率
occupancy_objective = (occupancies**2).sum()
    
objective = distance_objective - alpha * occupancy_objective

from amplify import one_hot
    
one_person_one_car_constraints = one_hot(q, axis=1)

from amplify import less_equal
    
car_capacity_constraints = less_equal(q, car_capacity, axis=0)

one_person_one_car_constraints *= np.max(distances) + alpha * 1
car_capacity_constraints *= np.max(distances) + alpha * 2 / car_capacity

from amplify import Model
    
model = objective + one_person_one_car_constraints + car_capacity_constraints

from amplify import FixstarsClient, solve
from datetime import timedelta
    
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=2000)  # タイムアウトは 2000 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境等で使用する場合は、Fixstars Amplify AE のアクセストークンを入力してください。
    
result = solve(model, client)
if len(result) == 0:
    # 実行可能解が見つかっていなければ例外を投げる
    raise RuntimeError("No feasible solution was found.")

q_values = q.evaluate(result.best.values)

# 車インデックスのリスト (利用者インデックス順)
car_indices = (q_values @ np.arange(num_cars)).astype(int)
# (利用者インデックス, 車インデックス) のタプルのリスト
car_assignment = list(enumerate(car_indices))
    
for person_index, car_index in car_assignment:
    print(f"利用者 {person_index}: 車 {car_index}")

import branca.colormap as cm
    
    
def plot_result(
    car_assignment: list[tuple[int, int]],
    people_coords: np.ndarray,
    car_coords: np.ndarray,
    latitude_range: tuple[float, float],  # 描画範囲 (緯度)
    longitude_range: tuple[float, float],  # 描画範囲 (経度)
):
    # マップを用意
    m = folium.Map(
        [sum(latitude_range) / 2, sum(longitude_range) / 2],
        tiles="OpenStreetMap",
        zoom_start=12,
    )
    
    # 色の用意
    colormap = cm.linear.Set1_09.scale(0, len(car_coords)).to_step(len(car_coords))  # type: ignore
    
    # 車のプロット (k 番目の車を色 k で塗る)
    for car_index, (latitude, longitude) in enumerate(car_coords):
        folium.Marker(
            location=(latitude, longitude),
            popup=f"car {car_index}",
            icon=folium.Icon(
                icon="car", prefix="fa", color="white", icon_color=colormap(car_index)
            ),
        ).add_to(m)
    
    # 利用者のプロット (車 k に乗る利用者を色 k で塗る)
    for person_index, car_index in car_assignment:
        latitude, longitude = people_coords[person_index]
        folium.Marker(
            location=(latitude, longitude),
            popup=f"car {car_index}",
            icon=folium.Icon(
                icon="user",
                prefix="fa",
                color="white",
                icon_color=colormap(car_index),
            ),
        ).add_to(m)
    
    return m

plot_result(car_assignment, people_coords, car_coords, latitude_range, longitude_range)

num_people = 60
num_cars = 18
    
people_coords, car_coords, distances = generate_problem(
    num_people, num_cars, latitude_range, longitude_range
)
    
plot_people_and_cars(people_coords, car_coords, latitude_range, longitude_range)

num_clusters = 3

gen = VariableGenerator()
x = gen.array("Binary", num_people, num_clusters)
y = gen.array("Binary", num_cars, num_clusters)

from amplify import einsum, Poly
    
objective: Poly = einsum("ij,ik,jk->", distances, x, y)  # type: ignore

from amplify import one_hot, equal_to
    
# x の各行に対する one-hot 制約
people_row_constraints = one_hot(x, axis=1)
    
# y の各行に対する one-hot 制約
car_row_constraints = one_hot(y, axis=1)
    
# x の各列の和が N/K に等しい
people_col_constraints = equal_to(x, num_people // num_clusters, axis=0)
    
# y の各列の和が M/K に等しい
car_col_constraints = equal_to(y, num_cars // num_clusters, axis=0)
    
constraints = (
    np.max(distances)
    * (num_people // num_clusters)
    * (
        people_row_constraints
        + people_col_constraints
        + car_row_constraints
        + car_col_constraints
    )
)

clustering_model = Model(objective, constraints)

client = FixstarsClient()
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
client.parameters.timeout = timedelta(milliseconds=1000)  # タイムアウト 1 秒
    
result = solve(clustering_model, client)
    
if len(result) == 0:
    # 実行可能解が見つからなければ例外を投げる
    raise RuntimeError("No feasible solution was found.")

x_values = x.evaluate(result.best.values)
y_values = y.evaluate(result.best.values)
    
people_cluster_indices = (x_values @ np.arange(num_clusters)).astype(int)
car_cluster_indices = (y_values @ np.arange(num_clusters)).astype(int)
    
for person_index, cluster_index in enumerate(people_cluster_indices[:5]):
    # 長いので最初の 5 個だけ表示
    print(f"利用者{person_index}: クラスター{cluster_index}")
    
for person_index, cluster_index in enumerate(people_cluster_indices[:5]):
    # 長いので最初の 5 個だけ表示
    print(f"車{person_index}: クラスター{cluster_index}")

def plot_clusters(
    people_coords,
    car_coords,
    people_cluster_indices,
    car_cluster_indices,
    latitude_range,
    longitude_range,
):
    m = folium.Map(
        [sum(latitude_range) / 2, sum(longitude_range) / 2],
        tiles="OpenStreetMap",
        zoom_start=12,
    )
    
    colors = ["red", "blue", "lightgray"]
    assert len(colors) == num_clusters
    
    for person_index, cluster_index in enumerate(people_cluster_indices):
        latitude, longitude = people_coords[person_index]
        folium.Marker(
            location=(latitude, longitude),
            popup=f"cluster{cluster_index}",
            icon=folium.Icon(icon="user", prefix="fa", color=colors[cluster_index]),
        ).add_to(m)
    
    for car_index, cluster_index in enumerate(car_cluster_indices):
        latitude, longitude = car_coords[car_index]
        folium.Marker(
            location=(latitude, longitude),
            popup=f"cluster{cluster_index}",
            icon=folium.Icon(icon="car", prefix="fa", color=colors[cluster_index]),
        ).add_to(m)
    
    return m
    
    
plot_clusters(
    people_coords,
    car_coords,
    people_cluster_indices,
    car_cluster_indices,
    latitude_range,
    longitude_range,
)

# 各クラスターに所属する利用者のインデックスのリストを作成
clustered_people_indices = [[] for k in range(num_clusters)]
for person_index, cluster_index in enumerate(people_cluster_indices):
    clustered_people_indices[cluster_index].append(person_index)
    
# 各クラスターに所属する車のインデックスのリストを作成
clustered_car_indices = [[] for k in range(num_clusters)]
for car_index, cluster_index in enumerate(car_cluster_indices):
    clustered_car_indices[cluster_index].append(car_index)
    
# 距離行列を分割
clustered_distances = [
    distances[clustered_people_indices[k]][:, clustered_car_indices[k]]
    for k in range(num_clusters)
]

def assign_cars(distances):
    num_people = distances.shape[0]
    num_cars = distances.shape[1]
    
    # 問題の正規化
    average = distances.mean()
    std = distances.std()
    distances = (distances - average) / std
    
    #  変数配列の作成
    gen = VariableGenerator()
    q = gen.array("Binary", num_people, num_cars)
    
    # 目的関数の構築
    alpha = 1.0  # 目的関数の重みパラメータ
    distance_objective = (distances * q).sum()  # 目的関数第一項
    occupancies = q.sum(axis=0) / car_capacity  # 乗車率
    occupancy_objective = (occupancies**2).sum()  # 目的関数第二項
    objective = distance_objective - alpha * occupancy_objective
    
    # 制約条件の構築
    one_person_one_car_constraints = one_hot(q, axis=1)
    car_capacity_constraints = less_equal(q, car_capacity, axis=0)
    # 制約重みを設定
    one_person_one_car_constraints *= np.max(distances) + alpha * 1
    car_capacity_constraints *= np.max(distances) + alpha * 2 / car_capacity
    
    # モデルの構築
    model = objective + one_person_one_car_constraints + car_capacity_constraints
    
    # ソルバークライアントの設定
    client = FixstarsClient()
    client.parameters.timeout = timedelta(milliseconds=2000)  # タイムアウトは 2000 ms
    # client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境等で使用する場合は、Fixstars Amplify AE のアクセストークンを入力してください。
    
    # 求解
    result = solve(model, client)
    q_values = q.evaluate(result.best.values)
    car_indices = (q_values @ np.arange(num_cars)).astype(int)
    
    # (利用者のインデックス, 車のインデックス) からなるリストを返却
    return list(enumerate(car_indices))

index_pairs_list = []
    
for d in clustered_distances:
    index_pairs_list.append(assign_cars(d))
    
for k, index_pairs in enumerate(index_pairs_list):
    print(f"クラスター {k}: ")
    for person_index, car_index in index_pairs:
        print(f"  利用者{person_index}: 車{car_index}")

car_assignment = [
    (clustered_people_indices[k][person_index], clustered_car_indices[k][car_index])
    for k, index_pairs in enumerate(index_pairs_list)
    for person_index, car_index in index_pairs
]
    
for person_index, car_index in car_assignment:
    print(f"利用者{person_index}: 車{car_index}")

plot_result(car_assignment, people_coords, car_coords, latitude_range, longitude_range)

	車 0	車 2
利用者 0	1	0
利用者 1	1	0
利用者 2	0	1
利用者 3	0	1

ライドシェア¶

問題設定¶

データ¶

定式化¶

変数¶

制約条件¶

目的関数¶

Note¶

問題の作成¶

組合せ最適化の実行¶

問題の正規化¶

変数¶

目的関数¶

制約条件¶

1 台の車を割り当てる制約の重み¶

車の定員を超えない制約の重み¶

組合せ最適化モデルの構築¶

組合せ最適化ソルバーの実行¶

発展的な話題（問題の分割）¶

問題設定¶

定式化¶

変数¶

制約¶

目的関数¶

実装¶

イジングマシンの実行¶