Ride sharing¶
The problem we are dealing with in this tutorial is called collective ridesharing.
Collective ridesharing refers to a form of ridesharing in which multiple users gather at several large parking lots and ride in the same car to the same destination(there is another type of ridesharing called traveling ridesharing, but it will not be discussed here).
Here, given multiple people with the same destination and available cars, we will find the allocation of people and cars such that the travel distance to the parking lot for each person and the number of cars to be used are as small as possible.
We formulate the problem as a model that can be run on an Ising machine, and find the allocation as a minimization problem.
Formulation¶
First, we define the constants and variables necessary for the formulation.
Constant¶
- NN:Number of rideshare users
- MM:Number of available cars
- CC:Number of available seats per car
- DD:Matrix such that ikikcomponent(dik)(dik) is the distance between user ii and car kk
Variables¶
- qik∈{0,1}(i∈{1,…,N},k∈{1,…,M})qik∈{0,1}(i∈{1,…,N},k∈{1,…,M})
Binary variables representing whether or not person ii rides in car kk (qik=1⇔qik=1⇔ person ii rides in car kk) - ylk∈{0,1}(l∈{0,…,C},k∈{1,…,M})ylk∈{0,1}(l∈{0,…,C},k∈{1,…,M})
Binary variables that satisfy ∑lylk=∑iqik∑lylk=∑iqik (used to express constraints on the number of passengers)
Then, we consider the constraints where the variables must satisfy.
Constraints¶
-
Each person always rides in one car.
∑Mk=1qik=1(∀i∈{1,…,N})∑Mk=1qik=1(∀i∈{1,…,N}) -
The actual number of passengers does not exceed the number of available seats.
∑Ni=1qik≤C(∀k∈{1,…,M})∑Ni=1qik≤C(∀k∈{1,…,M})
Finally, we will consider an objective function that satisfies the followings:
- Users use a car that is as close as possible to their location.
- Users travel with as few cars as possible.
Objective function¶
- Users should avoid unnecessary travel as much as possible. minimize∑i,kdikqikminimize∑i,kdikqik
- We want to minimize the number of cars used as much as possible ⇒ ⇒ Maximize the number of passengers per car. maximize∑k(∑iqikC)2maximize∑k(∑iqikC)2
Considering these two items, the following objective function can be set.
∑i,kdikqik−α∑k(∑iqikC)2∑i,kdikqik−α∑k(∑iqikC)2
Note¶
Let α>0α>0 be the parameter that determines how much importance is placed on the number of
cars in
use.
The closer the value of αα is to 00, the more the optimization places emphasis on minimizing
the
distances traveled by users. The greater the value of αα is, the more the optimization places
emphasis on minimizing the number of cars used.
If αα is large, the term regarding the distance traveled will be less important. The
visualization
result thus will be cleaner if αα is small.
Summary¶
From the above, the collective ridesharing problem can be formulated as the following Ising model.
H=Hcost+HconstraintHcost=∑i,kdikqik−α∑k(∑iqikC)2Hconstraint=k1N∑i=1(M∑k=1qik−1)2+k2M∑k=1(N∑i=1qik−C∑l=0ylk)2H=Hcost+HconstraintHcost=∑i,kdikqik−α∑k(∑iqikC)2Hconstraint=k1N∑i=1(M∑k=1qik−1)2+k2M∑k=1(N∑i=1qik−C∑l=0ylk)2(1)(2)(3)
k1,k2k1,k2 are constants that determine the strength of the constraints.
In order to ensure the feasibility of the solution, the size of the constant must be set so that the
objective function is not improved by violating the constraint.
In the present case, at least the following inequality should hold. The details of the derivation are
omitted.
k1>max(−min dik+2c−1c2α, max dik−2c−1c2α)k2>2c−1c2αk1>max(−min dik+2c−1c2α, max dik−2c−1c2α)k2>2c−1c2α(4)(5)
Implementing the problem¶
Since we need the positions of the cars and the users as input data, we create a function to randomly generate their positions (latitude and longitude).
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,
):
"""
A function that randomly determines the number of cars, the number of people, and the capacity of cars,
then generates the coordinates of the points of the number of cars + the number of people, and generates a distance matrix based on the coordinates.
Params
------
lon_range : list
Range of longitude when sampling coordinates of cars and people
lat_range : list
Latitude range for sampling the coordinates of a car and a person
parking : list
List containing coordinates of available parking spaces
ncars : int (default=None)
Integer specifying the number of parking spaces actually used
Must be set so that ncars <= len(parking)
npeople : int (default=None)
Number of rideshare users
C : int (default=None)
Maximum number of people that can ride in a car
lb : int (default=1)
Lower limit for generating random integer value
seed : int (default=None)
Seed of random numbers
Retuens
-------
ncars : int
Integer specifying the actual number of parking spaces to be used.
Must be set to satisfy ncars <= len(parking)
npeople : int
Number of rideshare users
D : np.ndarray
Matrix whose components are distances between users and cars
C : int
Maximum number of people that can ride in a car
ind2coord : dict
dictionary that associates car (person) indices with coordinates
key : int
index of car (person)
value : list
list of coordinates
"""
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:
raise ValueError(
"Fail to create valid problem.\nPlease retry after changing random seed."
)
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
For visualization purposes, we also create a function that plots the coordinates of cars and users on a map when they are input.
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=9)
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
After we define the candidate locations of the cars as follows, we use generate_problem
function defined earlier to generate the number of users, the number of people using the cars, the
number
of available seats in the cars, and the locations of the users and the cars. The
simple_plot
function is used to visualize them.
# Near Funabashi Station
lon = (359, 360)
lat = (51.5, 52)
# 9 locations
parking = [
(51.67699938102926, 359.6434199237448),
(51.60494726920934, 359.19303731029542),
(51.55604762650153, 359.41831984588475),
(51.69720660219214, 359.98034538800417),
(51.4581824540223, 359.30360550271415),
(51.408774929464875, 359.7982410856558),
(51.800029569368, 359.18558105961536),
(51.75599837320516, 359.83269833544272),
(51.65199204224218, 359.2415316476293),
]
ncars, npeople, D, C, index2coordinate = generate_problem(lon, lat, parking, seed=0)
simple_plot(index2coordinate, ncars)
print(ncars, npeople, C)
Building a quadratic polynomial model¶
Next, we define the necessary QUBO variables. Since we will have MM cars for NN users, we will define the QUBO variables as a M×NM×N two-dimensional array as follows:
from amplify import VariableGenerator
gen = VariableGenerator()
q = gen.array("Binary", npeople, ncars)
Then, we define the objective function and constraints. First, in order to align the order of the terms related to the distance and the number of cars in the objective function, we use the following function to adjust the mean of the elements of the distance matrix to 0 and the variance to 1.
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)
The next step is to define the objective function.
The amplify.sum
function is used to represent polynomials containing QUBO variables. The
objective function is as follows:
∑i,kdikqik−α∑k(∑iqikC)2
The first half of the term is related to the distance traveled and the second half is related to the occupancy rate of each car.
from amplify import sum
def setObjective(q, ncars, npeople, D, C, alpha=1):
"""Objective function"""
# Term related to the distance traveled by each user
distance_cost = sum(D[ncars:, :ncars] * q)
# Term related to the occupancy rate of each vehicle
ride_rate_cost = ((q.sum(axis=0) / C) ** 2).sum()
cost = distance_cost - alpha * ride_rate_cost
return cost
The constraint equation can be expressed as follows.
The one_hot
function is used to express constraint(1) ∑Mk=1qik=1(∀i∈{1,…,N}) and the less_equal
function is used to express constraint(2)
∑Ni=1qik≤C(∀k∈{1,…,M}).
As a guide to the weight of each constraint term, we introduced
k1>max(−min dik+2C−1C2α, max dik−2C−1C2α)k2>2C−1C2α
in the beginning. Here, as weights satisfying the above equations,
k1=2+2αC+1k2=2Cα+1
was chosen.
We can verify that they satisfy the above conditions as follows. First, since the distance matrix is normalized, max|dik|≤1 holds. Therefore,
max(−mindik+2C−1C2α, maxdik−2C−1C2α)<2+2C−1C2α
Also, since C≥1 is assumed, 2C−1C2<2CC2=2C.
From the above, we see that K1,K2 determined above satisfy the conditions.
from amplify import one_hot, less_equal
def setConstraints(q, ncars, npeople, C, k1=None, k2=None, alpha=1):
"""Functions to set constraint equations for small-scale problems"""
if k2 is None:
k2 = 2 * alpha / C + 1
if k1 is None:
k1 = (2 + 2 * alpha / C) + 1
# Constraint(1) that each person rides in one car
allocate_constraints = [one_hot(q[i]) for i in range(npeople)]
# Constraint(2) that no more than C people can fit in one car
capacity_constraints = [less_equal(sum(q[:, j]), C) for j in range(ncars)]
constraints = k1 * sum(allocate_constraints) + k2 * sum(capacity_constraints)
return constraints
The final QUBO equation is obtained by adding up the above objective function and constraints.
cost = setObjective(q, ncars, npeople, D, C)
constraints = setConstraints(q, ncars, npeople, C)
model1 = cost + constraints
Running the Ising machine¶
Set the Ising machine client to FixstarsClient
, create a solver, and solve the problem
as
follows:
# Use in the solving part
from amplify import FixstarsClient, solve
from datetime import timedelta
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=2000) # timeout is 2000 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
# Solve the problem
result = solve(model1, client)
Then, we will check the obtained solutions. You can use q.evaluate
to substitute them to
the
original variables.
if len(result) == 0:
raise RuntimeError("No feasible solution was found.")
q_values = q.evaluate(result.best.values)
Finally, we visualize the resulting assignments using the following function:
def plot_result(coord, q_values):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=9)
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)
Developmental topics (Splitting the problem)¶
Currently, the size of the problem that can be solved by an annealing machine is limited, so we
consider
dividing the problem by two classes of clustering.
By repeating the clustering process until the number of bits in the problem is below a set value, and
solving the problem at the beginning for each cluster obtained,
we aim to reduce the computation time and solve the problem of the number of bits.
Purpose¶
The ultimate goal is to solve this problem with following clustering technique:
- Divide the original problem into clusters until each cluster becomes the specified problem size.
- Solve each divided problem.
- Combine the divided solutions.
Formulation¶
The formulation is as follows:
Constant¶
- N:Number of rideshare users
- M:Number of available cars
- D:matrix that the ik component (dik) is the distance between user (car) i and user (car) k
Variables¶
qk∈{0,1}(k∈{1,…,M,…,M+N})
Binary variable representing which cluster a person (or car) k belongs to
(qk=1⇔ people (or cars) k belong to cluster 1)
Constraints¶
-
Divide it as evenly as possible
∑Mk=1qk=M2
∑M+Nk=M+1qk=N2
Objective function¶
-
People/cars that are (are) near each other belong to the same cluster
-
People/cars that are (are) far away from each other belong to different clusters
minimize∑i,jdij(2qi−1)(2qj−1)
From the above, the clustering of the two classes can be formulated as the following Ising model.
H=Hcost+HconstraintHcost=∑i,jdij(2qi−1)(2qj−1)Hconstraint=k1(M∑k=1qk−M2)2+k1(M+N∑k=M+1qk−N2)2 To guarantee the feasibility of the solution, the constant k1 must satisfy the inequality k1>2 max∑jdij.
Implementation¶
Based on the above formulation, we implement it using AMPLIFY.
The first step is to define the variables.
n = ncars + npeople
q = VariableGenerator().array("Binary", n)
print(D.shape)
print(q.shape)
Then, prepare the objective function.
from amplify import einsum
cost = einsum("ij,i,j->", D, (2 * q - 1), (2 * q - 1))
The constraints are as follows:
from amplify import equal_to, sum
# Constraint for the number of cars after the split to be half of the original number
car_constraints = equal_to(sum(q[:ncars]), ncars // 2)
# Constraint for the number of people after the split to be half of the original number.
people_constraints = equal_to(sum(q[ncars:n]), npeople // 2)
# Set the strength of the constraint
k1 = 2 * int(D.sum(axis=1).max()) + 3
constraints = car_constraints + people_constraints
constraints *= k1
The final model will look like this.
model_split = cost + constraints
Running the Ising machine¶
As before, we run the Ising machine to solve the problem.
result = solve(model_split, client)
if len(result) == 0:
# Throw an exception if no viable solution is found
raise RuntimeError("No feasible solution was found.")
else:
# If feasible solutions are found, display the objective function value corresponding to the best solution
objective = result.best.objective
values = result.best.values
print(f"objective = {objective}")
Then we define a function to divide the distance matrix and coordinates based on the solution.
def divide(q, D, coord, result):
"""Function to split the result of clustering"""
objective, values = result.best.objective, result.best.values
q_values = q.evaluate(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
Next, we define a function that will draw the result of the segmentation on the map
def plot_split_problem(coord: list, ncars: list):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=9)
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
Use the plot_split_problem
function to plot the coordinates after splitting.
D1, D2, cluster1, cluster2 = divide(q, D, index2coordinate, result)
plot_split_problem([cluster1, cluster2], [ncars // 2, ncars - ncars // 2])
We have implemented functions for dividing the problem. Let us implement more functions below to solve the divided problem.
The first step is to define a function to create a model for the divided problems.
def splitProblem(q, ncars, npeople, D, k1=None):
"""Function to divide the problem and create a model"""
n = ncars + npeople
if (
k1 is None
): # Set the coefficient as small as possible, since large coefficients may cause a problem
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(q[:ncars]), half_cars) + equal_to(
sum(q[ncars:n]), half_emp
)
model = cost + k1 * constraints
return model, half_cars, half_emp
Then, define a function to model each small-scale problem generated.
def construct(ncars, npeople, D, C, k1=None, k2=None, alpha=1):
"""Function to create a model for a small-scale problem after partitioning"""
D = regularizeDistance(D)
q = VariableGenerator().array("Binary", 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
Next, we define a function to integrate the optimization results of the small-scale problem.
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
import pandas as pd
def plot_final_result(data):
m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=9)
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
Then, we create a namedtuple
that manages the parameter α (which sets the
strength of
the term regarding the number of cars) that controls the size of the problem after splitting.
Here, we divide the problem up to 50 bits.
from collections import namedtuple
Parameter = namedtuple("Config", ("bit_size", "alpha"))
param = Parameter(bit_size=50, alpha=10)
We solve the problem by using the solve
function defined below.
The first step is to divide the problem into parts until the size of the problem becomes small
enough,
and then solve the divided problem.
The number of bits required for the small-scale problem is (number of cars) * (number of users +
number of
possible passengers),
so as long as this value is greater than a predetermined value (50 in this case), two classes will be
clustered.
from collections import deque
from IPython.display import display
def solve_ride_share(ncars, npeople, D, C, coord, param, debug=False):
print("Problem setting")
display(simple_plot(coord, ncars))
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=200) # timeout is 200 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
print("Splitting the problem...", 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 = VariableGenerator().array("Binary", ncars + npeople)
model, ncars1, npeople1 = splitProblem(q, ncars, npeople, D)
result = solve(model, client)
if len(result) == 0:
raise RuntimeError("No feasible solution was found.")
D1, D2, C1, C2 = divide(q, D, coord, result)
queue.append((ncars1, npeople1, D1, C1))
queue.append((ncars - ncars1, npeople - npeople1, D2, C2))
print("Completed")
print("Describing the results...", end="")
m = plot_split_problem([x[-1] for x in queue], [x[0] for x in queue])
display(m)
print("Completed")
print("Solving the problem after the split...")
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=1000) # timeout is 1000 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
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 = solve(model, client)
if len(result) == 0:
raise RuntimeError("No feasible solution was found.")
print(f"The small-scale problem {index} was solved")
q_values = q.evaluate(result.best.values)
data = make_data(coord, q_values, C, data=data, last=last)
last += ncars
print("Describing the results...")
m = plot_final_result(data)
display(m)
Let's solve it.
problem = generate_problem(lon, lat, parking, C=12, seed=0)
solve_ride_share(*problem, param)