Ride sharing¶
このチュートリアルで取り扱う問題は、集合型ライドシェアと呼ばれる問題です。
集合型ライドシェアとは、複数の利用者がいくつかの大型駐車場に集合し、
同じ車に乗って同じ目的地を目指す形式のライドシェアを指します。
(他の形式として巡回型ライドシェアというのもありますが、今回は名前のみの紹介とします。)
今回は、同じ目的地を持つ複数の人物と利用可能な車(駐車場)の候補が与えられている場合に、各人の車までの移動距離と使用する車の台数をなるべく小さくするような人と車の割り当てを考えます。
イジングマシンで実行可能なモデルとして定式化を行い、最小化問題として割り当てを求めていきます。
一つの駐車場に複数台の車がある場合は、乗車可能人数がその駐車場にある車の乗車可能人数の合計と等しい一台の車だけが駐車してあるとみなし、 駐車場と車が一対一に対応していることを仮定します。
定式化¶
まず、定式化に必要な定数・変数を定義します。
定数¶
- $N$:ライドシェア利用者数
- $M$:使用可能な車の台数
- $C$:車一台あたりの乗車可能人数
- $D$:$ik$成分$(d_{ik})$が利用者$i$と車$k$の間の距離となるような行列
変数¶
- $q_{ik}\in\{0,1\}\quad(i\in\{1,\dots,N\}, k\in\{1,\dots,M\})$
人$i$が車$k$に乗るかどうかを表現するバイナリ変数($q_{ik}=1\Leftrightarrow$人$i$が車$k$に乗る) - $y_{lk}\in\{0,1\}\quad(l\in\{0,\dots,C\},k\in\{1,\dots,M\})$
$\sum_ly_{lk}=\sum_iq_{ik}$を満たす変数(乗車可能人数に関する制約を表現するために使用)
続いて、変数が満たすべき制約を考えます。
制約条件¶
- 一人必ず一台の車に乗る
$\sum_{k=1}^Mq_{ik}=1(\forall i\in\{1,\dots,N\})$ - 実際の乗車人数が乗車可能人数を上回らない
$\sum_{i=1}^Nq_{ik}\leq C(\forall k\in\{1,\dots,M\})$
最後に、今回の目的である、
- 利用者はなるべく近い場所にある車を利用する
- なるべく少ない車で移動する
を満足できるような目的関数を考えます。
目的関数¶
- 利用者がなるべく無駄な移動をしない
$\text{minimize}\quad\sum_{i,k}d_{ik}q_{ik}$ - 車の使用台数をなるべく少なく抑えたい$\space\Rightarrow\space$一台あたりの乗車率を最大化する
$\text{maximize}\quad\sum_{k}\left(\sum_i\frac{q_{ik}}{C}\right)^2$
これら二項目を考慮すると、以下の目的関数が考えられます。
$$\sum_{i,k}d_{ik}q_{ik}-\alpha\sum_{k}\left(\sum_i\frac{q_{ik}}{C}\right)^2$$
Note¶
$\alpha$は車の使用台数をどのくらい重視するかを決定するパラメーターで、$\alpha>0$を満たすとします。
$\alpha$が0に近いほど移動距離を小さくするように最適化が行われ、$\alpha$の値が大きいほど車の使用台数を小さくするように最適化が行われます。
$\alpha$が大きい場合には移動距離に関する項が軽視されることにもなるため、可視化結果は$\alpha$が小さい方が綺麗なものになります。
まとめ¶
以上から、集合型ライドシェア問題は以下のようなQUBOとして定式化することができます。
$$ \begin{align} H&=H_{\rm cost}+H_{\rm costraint}\\ H_{\rm cost}&= \sum_{i,k}d_{ik}q_{ik}-\alpha\sum_{k}\left(\sum_i\frac{q_{ik}}{C}\right)^2\\ H_{\rm constraint} &= k_1\sum_{i=1}^N\left(\sum_{k=1}^Mq_{ik}-1\right)^2+k_2\sum_{k=1}^M\left(\sum_{i=1}^Nq_{ik}-\sum_{l=0}^Cy_{lk}\right)^2 \end{align} $$
ここで、$k_1, k_2$は制約の強さを決定する定数です。
解の実行可能性を担保するために、制約を破ることによって目的関数が改善しないように定数の大きさを設定する必要があります。
今回は少なくとも以下の不等式が成り立てば良いです。導出の詳細については省略します。
$$ \begin{align} k_1&>{\rm max}\left(− {\rm min\space}d_{ik}+ \frac{2C − 1}{C^2}\alpha,\space {\rm max\space}d_{ik}−\frac{2C − 1}{C^2}\alpha \right)\\ k_2&>\frac{2C − 1}{C^2}\alpha \end{align} $$
問題の実装¶
車と利用者の位置が入力データとして必要なため、それぞれの位置(緯度、経度)をランダムに生成する関数を作ります。
import numpy as np
from geopy.distance import geodesic
def generate_problem(
lon_range,
lat_range,
parking,
ncars=None,
npeople=None,
C=None,
lb=1,
ub=160,
seed=1,
):
"""
車の数、人の数、車の定員をランダムに決定した後、
車の数+人の数の点の座標を生成し、座標を元に距離行列を生成する関数
Params
------
lon_range : list
車と人の座標をサンプルする際の経度の範囲
lat_range : list
車と人の座標をサンプルする際の緯度の範囲
parking : list
使用可能な駐車場の座標が入ったリスト
ncars : int (default=None)
実際に使用する駐車場の数を指定する整数
ncars <= len(parking)を満たすように設定しなければならない
npeople : int (default=None)
ライドシェアの利用者数
C : int (default=None)
車の乗車可能人数の上限
lb : int (default=1)
ランダムな整数値を生成する時の下限
seed : int (default=None)
乱数のシード
Retuens
-------
ncars : int
実際に使用する駐車場の数を指定する整数
ncars <= len(parking)を満たすように設定しなければならない
npeople : int
ライドシェアの利用者数
D : np.ndarray
利用者、車の間の距離を成分に持つ行列
C : int
車の乗車可能人数の上限
ind2coord : dict
車(人)のインデックスと座標を紐づける辞書
key : int
車(人)のインデックス
value : list
座標を表すリスト
"""
np.random.seed(seed)
if ncars is None or (isinstance(ncars, int) and ncars > len(parking)):
if isinstance(ncars, int) and ncars > len(parking):
print(
f"Maximum value of ncars is {len(parking)}.\n ncars : {ncars} -> {len(parking)}."
)
ncars = len(parking)
if npeople is None:
npeople = np.random.randint(lb, ub)
if C is None:
C = np.random.randint(npeople // ncars + 1, npeople + 2)
if ncars * C < npeople:
print("Fail to create valid problem.\nPlease retry after changing random seed.")
return None, None, None, None, None
n = ncars + npeople
ind2coord = dict()
tmp = [
parking[i][::-1] for i in np.random.choice(len(parking), ncars, replace=False)
]
for i in range(ncars):
ind2coord[i] = (tmp[i][0], tmp[i][1])
for i in range(ncars, n):
lon = np.random.uniform(lon_range[0], lon_range[1])
lat = np.random.uniform(lat_range[0], lat_range[1])
tmp.append((lon, lat))
ind2coord[i] = (lon, lat)
D = np.zeros((n, n))
for i in range(n):
for j in range(n):
D[i, j] = geodesic(tmp[i][::-1], tmp[j][::-1]).m
return ncars, npeople, D, C, ind2coord
また可視化のため、車と利用者の座標を入力すると、それらを地図上にプロットする関数を作ります。
import folium
_colors = [
"green",
"orange",
"blue",
"pink",
"red",
"purple",
"darkblue",
"cadetblue",
"darkred",
"lightred",
"darkgreen",
"lightgreen",
"lightblue",
"gray",
"darkpurple",
]
def simple_plot(coord, ncars):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
tmp = list(coord.items())
for j, x in enumerate(tmp):
if j < ncars:
folium.Marker(
location=x[1][::-1],
icon=folium.Icon(icon="car", prefix="fa", color=_colors[0]),
).add_to(m)
else:
folium.Marker(
location=x[1][::-1],
popup="person",
icon=folium.Icon(icon="user", prefix="fa", color=_colors[1]),
).add_to(m)
return m
車の位置の候補を以下のように定め、先ほど定義したgenerate_problem
関数で利用者数、利用人数、車の乗車可能人数、利用者と車の位置を生成します。simple_plot
関数を用いてそれらを可視化します。
# 船橋駅周辺
lon = (139.9, 140.08)
lat = (35.675500, 35.76)
# 9箇所
parking = [
(35.67699938102926, 140.0434199237448),
(35.68494726920934, 139.99303731029542),
(35.68604762650153, 140.01831984588475),
(35.69720660219214, 139.98034538800417),
(35.6981824540223, 140.00360550271415),
(35.698774929464875, 139.9982410856558),
(35.700029569368, 139.98558105961536),
(35.70599837320516, 139.93269833544272),
(35.71199204224218, 140.0415316476293),
]
ncars, npeople, D, C, index2coordinate = generate_problem(lon, lat, parking, seed=0)
simple_plot(index2coordinate, ncars)
print(ncars, npeople, C)
二次多項式モデルの構築¶
次に必要なQUBO変数を定義します。$N$人の利用者に対して$M$台の車を用意するので、$N\times M$の二次元配列としてQUBO変数を以下のように定義します。
from amplify import BinarySymbolGenerator
gen = BinarySymbolGenerator()
q = gen.array(npeople, ncars)
続いて、目的関数、制約条件を定義していきます。
まず、目的関数の距離に関係する項と車の台数に関係する項のオーダーを揃えるため、以下の関数を用いて距離行列の要素の平均を0、分散を1に直します。
この操作を行うことで、$\alpha$の値を問題に依存せずに決定することができるようになります。
from math import sqrt
def regularizeDistance(D):
average = D.mean(axis=0, keepdims=True)
std = D.std(axis=0, keepdims=True, ddof=0)
return (D - average) / std
D = regularizeDistance(D)
次に目的関数を定義します。
sum_poly
関数を用いてQUBO変数を含む多項式を表現します。
目的関数は
$$\sum_{i,k}d_{ik}q_{ik}-\alpha\sum_{k}\left(\sum_i\frac{q_{ik}}{C}\right)^2$$
です。前半が移動距離に関係する項、後半が乗車率に関係する項です。
from amplify import sum_poly, einsum
def setObjective(q, ncars, npeople, D, C, alpha=1):
"""目的関数"""
# 各利用者の移動距離に関係する項
distance_cost = sum_poly(D[ncars:, :ncars] * q)
# 各車の乗車率に関係する項
ride_rate_cost = ((q.sum(axis=0) / C) ** 2).sum()
cost = distance_cost - alpha * ride_rate_cost
return cost
また、制約式は以下のように表現できます。
one_hot
関数を用いて制約(1)$\sum_{k=1}^Mq_{ik}=1(\forall i\in\{1,\dots,N\})$を、
less_equal
関数を用いて制約(2)$\sum_{i=1}^Nq_{ik}\leq C(\forall k\in\{1,\dots,M\})$を表現しています。
各制約項の強さの目安として、
$$
\begin{align}
k_1&>{\rm max}\left(− {\rm min\space}d_{ik}+
\frac{2C − 1}{C^2}\alpha,\space
{\rm max\space}d_{ik}−\frac{2C − 1}{C^2}\alpha
\right)\\
k_2&>\frac{2C − 1}{C^2}\alpha
\end{align}
$$
と冒頭で記載しました。
今回はこれを満たす係数として、
$$
\begin{align}
k_1&=2+\frac{2\alpha}{C}+1\\
k_2&=\frac{2}{C}\alpha+1
\end{align}
$$
を選択しました。
これらが冒頭の説明で与えた条件を満たしていることは、以下のようにして確認できます。
まず、距離行列は正規化されているため、$\max|d_ik|\leq 1$が成り立ちます。従って、
$$
\max\left(− \min d_{ik}+
\frac{2C − 1}{C^2}\alpha,\space
\max d_{ik}−\frac{2C − 1}{C^2}\alpha
\right) < 2 + \frac{2C − 1}{C^2}\alpha
$$
となります。
また、$C\geq 1$を仮定しているため、$\frac{2C-1}{C^2} < \frac{2C}{C^2} = \frac{2}{C}$が成り立ちます。
以上から、上のように定めた$k_1,k_2$が条件を満たしていることがわかりました。
from amplify.constraint import one_hot, less_equal
def setConstraints(q, ncars, npeople, C, k1=None, k2=None, alpha=1):
"""小規模問題の制約式を設定する関数"""
if k2 is None:
k2 = 2 * alpha / C + 1
if k1 is None:
k1 = (2 + 2 * alpha / C) + 1
# 一人一台の車に乗車するという制約(1)
allocate_constraints = [one_hot(q[i]) for i in range(npeople)]
# 一台の車に乗れるのはC人以下だという制約(2)
capacity_constraints = [less_equal(sum_poly(q[:, j]), C) for j in range(ncars)]
constraints = k1 * sum(allocate_constraints) + k2 * sum(capacity_constraints)
return constraints
上記目的関数と制約条件を足し合わせることによって最終的なQUBOを得ます。
cost = setObjective(q, ncars, npeople, D, C)
constraints = setConstraints(q, ncars, npeople, C)
model1 = cost + constraints
イジングマシンの実行¶
イジングマシンのクライアントをFixstarsClient
に設定、ソルバーを作成して以下のように問題を解きます。
from amplify import Solver
from amplify.client import FixstarsClient
client = FixstarsClient()
client.parameters.timeout = 2000 # 制限時間
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # ローカル環境で使用する場合は、Fixstars Amplify AEのアクセストークンを入力してください
solver = Solver(client)
result = solver.solve(model1)
続いて、得られた解を確認していきます。
q.decode()
を用いてもとの変数に代入することができます。
if len(result) == 0:
# 実行可能解が見つかっていなければ例外を投げる
raise RuntimeError("No feasible solution was found.")
q_values = q.decode(result[0].values)
最後に、以下の関数を用いて得られた割り当てを可視化します。
def plot_result(coord, q_values):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
npeople = len(q_values)
ncars = len(q_values[0])
columns = ["latitude", "longitude", "size", "name"]
data = {label: list() for label in columns}
answer = dict()
for i in range(npeople):
car = np.where(np.array(q_values[i]) == 1)[0][-1]
if car not in answer:
answer[car] = []
answer[car].append(i + ncars)
for k in range(ncars):
status = "active"
car_loc = coord[k]
if k in answer:
tmp = answer[k]
x = [coord[p][0] for p in tmp] + [car_loc[0]]
y = [coord[p][1] for p in tmp] + [car_loc[1]]
else:
x = car_loc[:1]
y = car_loc[1:]
status = "empty"
folium.Marker(
location=[y[-1], x[-1]],
popup=f"cluster{k}",
icon=folium.Icon(icon="car", prefix="fa", color=_colors[k % len(_colors)]),
).add_to(m)
for a, b in zip(y[:-1], x[:-1]):
folium.Marker(
location=[a, b],
popup=f"person{k}",
icon=folium.Icon(
icon="user",
prefix="fa",
color="white",
icon_color=_colors[k % len(_colors)],
),
).add_to(m)
return m
plot_result(index2coordinate, q_values)
発展的な話題(問題の分割)¶
現在、アニーリングマシンで解ける問題のサイズは限られているため、2クラスのクラスタリングによって問題を分割することを考えます。
問題のビット数が設定した値を下回るまでクラスタリングを繰り返し、得られた各クラスタについて冒頭の問題を解くことで、
計算時間の削減やビット数の問題を解決することを狙っています。
やりたいこと¶
以下のフローチャートに沿った最適化が最終的な目的です。
定式化¶
定式化は以下の通りです。
定数¶
- $N$:ライドシェア利用者数
- $M$:使用可能な車の台数
- $D$:$ik$成分$(d_{ik})$が利用者(車)$i$と利用者(車)$k$の間の距離となるような行列
変数¶
$q_{k}\in\{0,1\}\quad(k\in\{1,\dots,M,\dots,M+N\})$
人(or 車)$k$がどのクラスタに属するかを表すバイナリ変数
($q_{k}=1\Leftrightarrow$人(or 車)$k$はクラスタ1に属する)
制約¶
- なるべく均等に分けたい
$\sum_{k=1}^Mq_k=\frac{M}{2}$
$\sum_{k=M+1}^{M+N}q_k=\frac{N}{2}$
目的関数¶
-
互いに近くにいる(ある)人・車は同じクラスタに属する
互いに遠くにいる(ある)人・車は異なるクラスタに属する$\text{minimize}\quad\sum_{i,j}d_{ij}(2q_i-1)(2q_j-1)$
以上から、2クラスのクラスタリングは以下のようなQUBOとして定式化することができます。
$$ \begin{align} H&=H_{\rm cost}+H_{\rm constraint}\\ H_{\rm cost}&=\sum_{i,j}d_{ij}(2q_i-1)(2q_j-1)\\ H_{\rm constraint}&=k_1\left(\sum_{k=1}^Mq_k-\frac{M}{2}\right)^2+k_1\left(\sum_{k=M+1}^{M+N}q_k-\frac{N}{2}\right)^2 \end{align} $$ 解の実行可能性を保証するために、定数$k_1$は不等式$k_1>2{\space\rm max}\sum_jd_{ij}$を満たす必要があります。
実装¶
以上の定式化を元に、amplifyを用いて実装していきます。
まずは変数を定義します。
n = ncars + npeople
q = BinarySymbolGenerator().array(n)
print(D.shape)
print(q.shape)
続いて目的関数を用意します。
from amplify import einsum
cost = einsum("ij,i,j->", D, (2 * q - 1), (2 * q - 1))
制約条件は以下の通りです。
from amplify.constraint import equal_to
# 分割後の台数が元の台数の半分になるための制約
car_constraints = equal_to(sum_poly(q[:ncars]), ncars // 2)
# 分割後の人数が元の台数の半分になるための制約
people_constraints = equal_to(sum_poly(q[ncars:n]), npeople // 2)
# 制約の強さを設定
k1 = 2 * int(D.sum(axis=1).max()) + 3
constraints = car_constraints + people_constraints
constraints *= k1
最終的なモデルは以下のようになります。
model_split = cost + constraints
イジングマシンの実行¶
これまでと同様に、イジングマシンを実行して問題を解いていきます。
result = solver.solve(model_split)
if len(result) == 0:
# 実行可能解が見つからなければ例外を投げる
raise RuntimeError("No feasible solution was found.")
else:
# 実行可能解が見つかればそれらの目的関数値を順番に表示する
for solution in result:
energy = solution.energy
values = solution.values
print(f"energy = {energy}")
続いて、解を元にして距離行列と座標を分割するための関数を定義します。
def devide(q, D, coord, result):
"""クラスタリングの結果を分割する関数"""
energy, values = result[0].energy, result[0].values
q_values = q.decode(values)
cluster1 = np.where(np.array(q_values) == 1)[0]
cluster2 = np.where(np.array(q_values) == 0)[0]
nc1 = len(cluster1)
nc2 = len(cluster2)
D1 = np.zeros((nc1, nc1))
D2 = np.zeros((nc2, nc2))
C1 = dict()
C2 = dict()
for i in range(nc1):
C1[i] = coord[cluster1[i]]
for j in range(nc1):
D1[i][j] = D[cluster1[i]][cluster1[j]]
for i in range(nc2):
C2[i] = coord[cluster2[i]]
for j in range(nc2):
D2[i][j] = D[cluster2[i]][cluster2[j]]
return D1, D2, C1, C2
次に分割結果を地図上に描写する関数を定義します
def plot_splited_problem(coord: list, ncars: list):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
for i in range(len(ncars)):
tmp = list(coord[i].items())
for j, x in enumerate(tmp):
if j < ncars[i]:
folium.Marker(
location=x[1][::-1],
popup=f"cluster{i}",
icon=folium.Icon(icon="car", prefix="fa", color=_colors[i]),
).add_to(m)
else:
folium.Marker(
location=x[1][::-1],
popup=f"person{i}",
icon=folium.Icon(
icon="user", prefix="fa", color="white", icon_color=_colors[i]
),
).add_to(m)
return m
plot_splited_problem
関数を用いて、分割後の座標をプロットします。
D1, D2, cluster1, cluster2 = devide(q, D, index2coordinate, result)
plot_splited_problem([cluster1, cluster2], [ncars // 2, ncars - ncars // 2])
まとめ¶
最後に、以下のフローチャート(再掲)に沿って問題を分割$\Rightarrow$分割後の問題を解くという一連の流れを実装していきます。
まずは問題を分割するためのモデルを作成する関数を定義します。
def splitProblem(q, ncars, npeople, D, k1=None):
"""問題を分割するためのモデルを作成する関数"""
n = ncars + npeople
if k1 is None: # 係数が大きすぎると逆に問題がありそうなので、なるべく小さく設定
k1 = 2 * int(max([sum(D[i]) for i in range(n)])) + 3
half_cars = ncars // 2
half_emp = npeople // 2
cost = einsum("ij,i,j->", D, (2 * q - 1), (2 * q - 1))
constraints = equal_to(sum_poly(q[:ncars]), half_cars) + equal_to(
sum_poly(q[ncars:n]), half_emp
)
model = cost + k1 * constraints
return model, half_cars, half_emp
続いて、生成された各子問題をモデル化するための関数を定義します。
def construct(ncars, npeople, D, C, k1=None, k2=None, alpha=1):
"""分割後の小規模な問題のためのモデルを作成する関数"""
D = regularizeDistance(D)
q = BinarySymbolGenerator().array(npeople, ncars)
cost = setObjective(q, ncars, npeople, D, C, alpha=alpha)
constraints = setConstraints(q, ncars, npeople, C, k1=k1, k2=k2, alpha=alpha)
model = cost + constraints
return model, q
次に、子問題の最適化結果を統合するための関数を定義します。
def make_data(coord, q_values, C, last=0, data=None, nframe=1):
if data is None:
columns = ["latitude", "longitude", "size", "name", "time", "color"]
data = {label: list() for label in columns}
npeople = len(q_values)
ncars = len(q_values[0])
answer = dict()
for i in range(npeople):
car = np.where(np.array(q_values[i]) == 1)[0][-1]
if car not in answer:
answer[car] = []
answer[car].append(i + ncars)
loc = [[], []]
for k in range(ncars):
status = "active"
car_loc = coord[k]
if k in answer:
tmp = answer[k]
x = [coord[p][0] for p in tmp] + [car_loc[0]]
y = [coord[p][1] for p in tmp] + [car_loc[1]]
else:
x = car_loc[:1]
y = car_loc[1:]
status = "empty"
loc[0] += y
loc[1] += x
for i in range(nframe):
data["latitude"] += list(
map(lambda a: ((nframe - i) * a + y[-1] * i) / nframe, y)
)
data["longitude"] += list(
map(lambda a: ((nframe - i) * a + x[-1] * i) / nframe, x)
)
data["size"] += [0.5] * (len(x) - 1) + [3]
data["name"] += [f"group{k+last}({status})"] * len(x)
data["time"] += [i] * len(x)
data["color"] += [_colors[k + last]] * len(x)
return data
def plot_final_result(data):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
df = pd.DataFrame(data)
for _name in data["name"]:
tmp = df[df["name"] == _name]
x = list(tmp["longitude"])
y = list(tmp["latitude"])
folium.Marker(
location=[y[-1], x[-1]],
popup=f"cluster_{_name}",
icon=folium.Icon(icon="car", prefix="fa", color=list(tmp["color"])[-1]),
).add_to(m)
for a, b, c in zip(y[:-1], x[:-1], list(tmp["color"])[:-1]):
folium.Marker(
location=[a, b],
popup=f"person_{_name}",
icon=folium.Icon(icon="user", prefix="fa", color="white", icon_color=c),
).add_to(m)
return m
続いて、分割後の問題の大きさを制御するパラメータ、$\alpha$(車の台数に関する項の強さを設定)を管理するnamedtuple
を作成します。
今回は、問題を50ビットまで分割することにします。
from collections import namedtuple
Parameter = namedtuple("Config", ("bit_size", "alpha"))
param = Parameter(bit_size=50, alpha=10)
以下のsolve
関数を使って問題を解きます。
まず問題サイズがある程度小さくなるまで問題の分割を行い、次に分割後の問題を解いていく流れです。
小規模問題に必要なビット数は、(車の台数) * (利用者数 + 乗車可能人数)のため、
この値が事前に設定した値(今回は50)を上回っている限り2クラスのクラスタリングを行います。
from collections import deque
from IPython.display import display
import pandas as pd
def solve(ncars, npeople, D, C, coord, param, debug=False):
print("問題設定")
display(simple_plot(coord, ncars))
client = FixstarsClient()
client.parameters.timeout = 200 # 制限時間
solver = Solver(client)
print("問題を分割しています...", end="")
queue = deque([(ncars, npeople, D, coord)])
while (queue[0][1] + C) * queue[0][0] > param.bit_size:
(ncars, npeople, D, coord) = queue.popleft()
q = BinarySymbolGenerator().array(ncars + npeople)
model, ncars1, npeople1 = splitProblem(q, ncars, npeople, D)
result = solver.solve(model)
if len(result) == 0:
raise RuntimeError("No feasible solution was found.")
D1, D2, C1, C2 = devide(q, D, coord, result)
queue.append((ncars1, npeople1, D1, C1))
queue.append((ncars - ncars1, npeople - npeople1, D2, C2))
print("完了")
print("結果を描写しています...", end="")
m = plot_splited_problem([x[-1] for x in queue], [x[0] for x in queue])
display(m)
print("完了")
print("分割後の問題を解いています...")
client = FixstarsClient()
client.parameters.timeout = 1000 # 制限時間
solver = Solver(client)
index = 0
last = 0
data = None
while queue:
index += 1
(ncars, npeople, D, coord) = queue.pop()
model, q = construct(ncars, npeople, D, C, alpha=param.alpha)
result = solver.solve(model)
if len(result) == 0:
raise RuntimeError("No feasible solution was found.")
print(f"子問題{index}が解けました")
q_values = q.decode(result[0].values)
data = make_data(coord, q_values, C, data=data, last=last)
last += ncars
print("結果を描写しています...")
m = plot_final_result(data)
display(m)
実際に解いてみましょう。
problem = generate_problem(lon, lat, parking, C=12, seed=0)
solve(*problem, param)