[PYTHON] I tried to find the optimal path of the dreamland by (quantum) annealing

Introduction

Dreamland ... Yes, it's Tokyo Disneyland, known as the Kingdom of Dreams and Magic.

I'm sure you've visited it once, but it's quite crowded on holidays, isn't it? In such a case, waiting time and travel distance are the burdens in the dreamland.

It would be ideal if there was a course that would shorten the waiting time and travel distance. However, when it comes to thinking for yourself, it is necessary to accurately grasp the location and waiting time of a huge number of attractions.

It's annoying, so it seems to break my dreams, but I made a web app that optimizes the route of the dreamland. There are various optimization methods, but I just tried quantum annealing, which I am dealing with in my university thesis.

If you want to try it first, please click here! https://tdr-planner.web.app/ The optimization process will take about 30 seconds, so please wait for a while m (_ _) m

Glossary

A lot of technical terms will appear, so I will explain them roughly in advance.

Quantum computer

Quantum computers are roughly classified into the following two types.

method	Feature
Gate method	General-purpose and can realize all kinds of quantum calculations (general-purpose machine)
Ising model method	Specializing in combinatorial optimization problems (dedicated machine)

At first glance, the general-purpose gate method seems to be superior, but the calculations that are guaranteed to be fast are still limited.

This time, we will consider Disney's course optimization problem as a combinatorial optimization problem and solve it using the Ising model method.

Combinatorial optimization problem

The ** combinatorial optimization problem ** that can be handled by the Ising model quantum computer is as follows.

The combinatorial optimization problem is to find the value (combination) of the variable that gives the best index (value) from among many options ** under various constraints.

Quoted from Annealing Cloud Web

There is a concrete example in the above citation source, so it will be easier to understand if you refer to it.

In the problem we are dealing with this time, we will find the combination that is most satisfying from the huge number of combinations of how to go around the attractions.

Ising model

My specialty suddenly increased, but I will explain it because it is unavoidable.

The Ising model quantum computer is a dedicated machine specializing in combinatorial optimization problems, so its input is limited to a specific form. The input format is a model used in statistical mechanics called the ** Ising model **.

When you input it into a quantum computer (Ising model method) in the form of an Ising model, it searches for a combination of variables that minimizes the energy of that model.

The concrete energy formula is expressed as follows.

H = \sum_{i \neq j}J_{ij}\sigma_{i}\sigma_{j} + \sum_{i}h_{i}\sigma_{i}\quad(\sigma = \pm 1)

Suddenly a complicated formula comes out, what? I'm feeling. I am not confident in the theoretical explanation, so I will omit it here. Just be aware that ** $ H $ is represented by a polynomial of up to quadratic order of $ \ sigma $ **.

When solving, enter the coefficient of $ \ sigma $ in matrix format, and it will search for the combination of $ \ sigma $ that minimizes $ H $. For the time being, the first goal of solving with a quantum computer is to express it with a polynomial of maximum quadratic order of $ \ sigma $.

QUBO

QUBO is another model similar to the Ising model, and the energy is expressed by the following formula.

H = \sum_{i \neq j}J_{ij}x_{i}x_{j} + \sum_{i}h_{i}x_{i}\quad(x \in 0,1)

The only difference from the Ising model is that the variable replaces $ \ sigma $, which takes $ \ pm 1 $, with $ x $, which takes $ 0,1 $.

QUBO will be used for this Disney course optimization.

Rough flow

Formulate Disney's course optimization as a combinatorial optimization problem with QUBO
Actually solve with a (quantum) computer
Try to make it a web application

Problem formulation

In formulating the problem, we will consider using the following figure.

The circles in the figure above correspond to the binary variable $ x_ {i, j} \ (1 \ leq i, j \ leq 4) $ in QUBO. "Which attraction to use" corresponds to the row, and "Which attraction to use" corresponds to the column. $ i $ Represents whether to use the attraction $ j $ as a matrix.

In the case of the above figure, the course is expressed in the order of Haunted Mansion → Pooh's Honey Hunt → Space Mountain → Splash Mountain.

I'm thinking of only 4 attractions, but there are 65536 combinations of 4 × 4 = 16 $ x $ in 2 to the 16th power. The quantum computer finds the combination of $ x $ that minimizes the energy $ H $ indicated by QUBO from these 65536 combinations and returns it in the matrix format as shown in the above figure. I will.

In other words, if the expression that minimizes $ H $ in the ideal course can be expressed using $ x $ (second order or less), it can be optimized by a quantum computer. In general, when we think of $ H $, we think of objectives and constraints as the terms that make up $ H $.

Item	role
Objective	The target term that you want to maximize or minimize. Set so that the better the result, the smaller the value.
Constraint	The minimum term if the constraint is met. If the constraint is violated, the value is set to be large as a penalty.

After that, we will consider $ H $ based on ↑. However, the number of attractions used (number of rows) and the number of attractions (number of columns) are generalized to $ n and m $, respectively. (For convenience, $ n = m = 4 $ in the figure)

Purpose

At the beginning, I mentioned that I wanted to shorten the waiting time and travel distance, so that seems to be the purpose, but when I actually go to Disneyland, I don't think there are people who say, "I want to wait within 3 hours in total!"

In reality, I think the purpose is "I want to ride a lot of attractions!" However, quality as well as quantity must be considered. I don't want to go to Minnie's House 10 times because the waiting time is short. Therefore, we set 10 levels of satisfaction for each attraction and aim to ** maximize total satisfaction **.

Assuming that the satisfaction level of the attraction $ j $ is $ s_j \ \ (1 \ leq j \ leq n, 1 \ leq s_j \ leq 10) $, the total satisfaction level $ S $ is as follows.

S = \sum_{j=1}^{m}\sum_{i=1}^{n}s_{j}x_{i,j}

It looks like this in the figure.

This time, the purpose is to maximize the total satisfaction **, but since the objective should be set so that the better the result, the smaller the value, the inverted sign is the objective $ H_ {aim} $. Set as.

H_{aim} = -\sum_{j=1}^{m}\sum_{i=1}^{n}s_{j}x_{i,j}

Constraint A

Now that the objectives have been decided, what happens if we just optimize them for the time being? There is a high probability that such a solution will be obtained.

That's right. Satisfaction will be maximized if you use all of them. However, this means that you will be using four attractions at the same time, and the course has not been established.

Therefore, set a restriction that ** the use of multiple attractions at the same time is prohibited **. This means that there is always only one blue circle on each line. This constraint is a one-hot constraint that appears frequently in QUBO and is expressed by the following formula.

H_{a} = \sum_{i=1}^{n}\left(1-\sum_{j=1}^{m}x_{i,j}\right)^2

It looks like this in the figure.

In the figure above, there are two blue circles on the second line and $ H_ {a} = 1 $ due to a constraint violation. If the constraint is met, $ H_ {a} $ has a minimum value of 0.

Constraint B

What if we optimize with the objective and constraint A? You will probably get a solution like this. (* Numbers in parentheses are satisfaction levels)

Although it is established as a course by constraint A, it is a bold course with 4 consecutive Splash Mountains. It's only natural that Splash Mountain has the highest satisfaction level. However, no matter how popular the attraction is, it is not a good course to ride four times.

Therefore, if you use the attraction multiple times, set a constraint that the satisfaction level will decrease by 3 each time you use it. For example, if you use Splash Mountain 3 times, your satisfaction will decrease, such as 10 for the first time, 7 for the second time, and 4 for the third time. This way, the third satisfaction is 4, so Space Mountain (6) and Haunted Mansion (5) will be chosen.

Although it is a constraint, the objective term $ H_ {aim} $ is modified because it is related to the calculation of total satisfaction.

The policy is to calculate the total satisfaction level of one attraction and sum it up for all attractions. The number of times $ k_ {a} $ used for an attraction $ a $ is

k_{a}=\sum_{i=1}^{n}x_{i,a}

It will be. Since the satisfaction level of attraction $ a $ decreases by 3 each time it is used, it can be regarded as an arithmetic progression with the first term $ s_a $ tolerance-3, and the total $ S_ {a} $ is expressed by the following formula. I will. (Sum of arithmetic progression)

S_{a}=\frac{1}{2}k_{a}(2s_{a}+(k_{a}-1)\times(-3))= k_{a}(s_{a} - \frac{3}{2}k_{a} + \frac{3}{2})

It looks like this when using Splash Mountain three times.

The sum of all attractions is the total satisfaction, so

S = \sum_{j=1}^{m}S_{j} = \sum_{j=1}^{m} \left( \sum_{i=1}^{n}x_{i,j} \left( s_j-\frac{3}{2}\sum_{i=1}^{n}x_{i,j}+\frac{3}{2} \right) \right)

It will be. Set the inverted sign as the objective term in the same way as [First set objective term](# objective).

H_{aim}= -\sum_{j=1}^{m} \left( \sum_{i=1}^{n}x_{i,j} \left( s_j-\frac{3}{2}\sum_{i=1}^{n}x_{i,j}+\frac{3}{2} \right) \right)

Constraint C

What happens if you optimize using the newly set objectives? The following is an example of the optimal solution.

Now that the course has become quite decent, let's think about the distance traveled. The positional relationship of each attraction and the course in the above figure are as follows.

Since it is a schematic diagram, it differs from the actual positional relationship, but the magnitude relationship of the travel distance is correctly expressed. If the order of Space Mountain and Haunted Mansion is reversed, the travel distance will be shorter.

As mentioned above, it is ideal to shorten the travel distance, but it is not realistic to specifically want to keep the total travel distance within 〇〇m.

Therefore, set a constraint that the total required time is within XX hours **. If you calculate the travel time and waiting time and express this constraint, you will maximize the total satisfaction within the limited time, so the travel time and waiting time should be shortened naturally.

In expressing this constraint, the following three are defined.

--Scheduled stay time (minutes): $ t_ {max} $ --Cost of attraction $ j $ (= waiting time + required time) (minutes): $ c_ {j} \ (1 \ leq j \ leq n) $ --Travel time (minutes) between two attractions $ i, j $: $ d_ {i, j} \ (1 \ leq i, j \ leq n) $

Divide the total travel time into ** "total cost of attractions to use" ** + ** "total travel time" **.

First, the total cost of the attraction to be used, $ C $, can be expressed as follows in the same way as First set objective. I just replaced the satisfaction $ s $ with the cost $ c $.

C = \sum_{j=1}^{m}\sum_{i=1}^{n}c_{j}x_{i,j}

Next, the total travel time $ D $ is [Traveling salesman problem](http://blog.brainpad.co.jp/entry/2017/04/20/160000#%E4%BE%8B3-%E5%B7 % A1% E5% 9B% 9E% E3% 82% BB% E3% 83% BC% E3% 83% AB% E3% 82% B9% E3% 83% 9E% E3% 83% B3% E5% 95% 8F It can be expressed as follows in the same way as% E9% A1% 8C). (I gave up because it was difficult to express in the figure ... If you want to understand, please check the traveling salesman problem)

D = \sum_{t=1}^{n-1}\sum_{i=1}^{m}\sum_{j=1}^{m}d_{i,j}x_{t,i}x_{t+1,j}

Set the difference between the total travel time $ (C + D) $ and the estimated stay time $ t_ {max} $ as the constraint $ H_ {b} $.

H_{b}=\sum_{j=1}^{m}\sum_{i=1}^{n}c_{j}x_{i,j}+\sum_{t=1}^{m-1}\sum_{i=1}^{n}\sum_{j=1}^{n}d_{i,j}x_{t,i}x_{t+1,j} -t_{max}

$ H_ {b} $ serves as a constraint because the value increases as a penalty when the total required time exceeds the planned stay time. However, the shorter the total travel time, the smaller $ H_ {b} $ will be, so it will be easier to choose a course with a shorter total travel time. This will be adjusted with the hyperparameters described later.

Constraint D

Optimized by adding the constraint $ H_ {b} $, it looks like this.

You have come to go to the Haunted Mansion first, considering the distance properly. However, I'm worried about the continuation of Splash Mountain.

Computerly, if you ride continuously, the travel time between them will be 0, so if you use it multiple times, it will be continuous. Even if it's OK to ride twice, you want to avoid continuous riding if possible.

Therefore, set a constraint that ** the same attraction should not be used consecutively **. This corresponds to the fact that blue circles are not adjacent in the column direction, and is expressed by the following formula.

H_{c} = \sum_{j=1}^{m}\sum_{i=1}^{n-1}x_{i,j}x_{i+1,j}

It looks like this in the figure.

In the above figure, there is a part where the blue circle is adjacent in the first column and it is $ H_ {c} = 1 $ because of the constraint violation. If the constraint is met, $ H_ {c} $ will have a minimum value of 0.

Constraint E

Since I have considered four attractions so far, it suddenly becomes a course like Splash Mountain. But in a real course, you have to enter from the entrance and come back to the entrance.

Therefore, we set a constraint that ** use the entrance at the beginning and end of the course **. The entrance is also considered as an attraction and added as a line. Satisfaction is set to 0 because it is not necessary to use it in the middle of the course. Assuming the entrance is in the first column, this constraint is expressed by the following formula.

H_{d} = 2 - x_{1,1} - x_{n,1}

It looks like this in the figure. If the first $ (x_ {1,1}) $ and the last $ (x_ {5,1}) $ in the first column are blue circles, the constraint is satisfied.

If the constraint is met, $ H_ {d} $ will have a minimum value of 0.

Formulation summary

It's been a long time, but now I'm looking for a fairly reasonable course. Below is a summary of Disney's course optimization formulations.

Input data

variable	Contents
n	Number of attractions used in one course (how many you want to ride)
m	Number of attractions (Tokyo Disneyland: 34）
t_{max}	Estimated stay time(Minutes)
s_{j}	AttractionsjSatisfaction(1 \leq j \leq n, 1 \leq s_{j} \leq 10)
c_{j}	AttractionsjCost(Minutes)　(1 \leq j \leq n)
d_{i,j}	Attractionsi,jTravel time between(Minutes)　(1 \leq i,j \leq n)

Purpose

--Maximizing total satisfaction

Constraint

--The use of multiple attractions at the same time is prohibited ――If you use the attraction multiple times, your satisfaction will decrease by 3 each time you use it. --Total time required within the planned stay time --Continuous use of the same attraction is prohibited --Use the entrance at the beginning and end of the course

Energy function

\begin{align}
H_{aim}&=-\sum_{j=1}^{m}\left(\sum_{i=1}^{n}x_{i,j}\left(s_j-\frac{3}{2}\sum_{i=1}^{n}x_{i,j}+\frac{3}{2}\right)\right) \\
H_{a}&=\sum_{i=1}^{n}\left(1-\sum_{j=1}^{m}x_{i,j}\right)^2 \\
H_{b}&=\sum_{j=1}^{m}\sum_{i=1}^{n}c_{j}x_{i,j}+\sum_{t=1}^{n-1}\sum_{i=1}^{m}\sum_{j=1}^{m}d_{i,j}x_{t,i}x_{t+1,j} -t_{max} \\
H_{c}&=\sum_{j=1}^{m}\sum_{i=1}^{n-1}x_{i,j}x_{i+1,j} \\
H_{d}&=2- x_{1,1}-x_{n,1}
\end{align}

The final energy function $ H $ looks like this: $ \ Alpha, \ beta, \ gamma, \ delta $ are positive hyperparameters that represent constraint weights.

H = H_{aim} + \alpha H_{a} + \beta H_{b} + \gamma H_{c} + \delta H_{d}

Try to solve with (quantum) computer

Let's actually solve it using the formulated energy function. This time, we will focus on 10 types of attractions that are likely to be popular, selected by dogmatism and prejudice.

First, prepare the following input data.

For the waiting time, I referred to the here site. Uses the average daily wait time for each attraction on Saturday-Sunday levels of normal congestion. (Waiting time, but is time-series data that changes from moment to moment, be carried out at any time for calculation simplicity we assume a constant waiting time.) The time required will be the one posted on Disney's official website. The cost is waiting time + required time, and satisfaction is my favorite.

When I was having trouble finding data on the travel time between attractions, a friend who loves Disney cooperated and calculated the distance between attractions. When asked how it was calculated, he said that he set the branch point of each attraction and passage as a node on Google Maps and found the shortest distance by Dijkstra's algorithm. Moreover, since there are passages that only the cast can pass through with Google Maps alone, it is the first time to make a map of only the passages that guests can pass through based on Google Maps. It's scary to like Disney in science ... ↑ A map created by a friend to find the shortest distance between attractions in Dijkstra. Does anyone know that it's Disneyland by looking at this ...

Now that the input data is ready, convert it to the format to be input to the quantum computer. Originally, I managed to expand $ H $ to create a coefficient matrix (QUBO matrix) for the linear and quadratic terms of $ x $, but this time I generated a QUBO matrix with highly readable code. I used the python library pyqubo.

import numpy as np
import pandas as pd
from pyqubo import Array, Sum, Num, Constraint, Placeholder

cost_data = pd.read_csv('tdl_cost.csv')
distance_data = pd.read_csv('tdl_distance.csv', header=None)

s = cost_data['value'].to_numpy()
c = cost_data['cost'].to_numpy()
d = (distance_data / 80).to_numpy() #Travel time assuming a walking speed of 80 m / min(Minutes)Convert to
n = 12
m = len(c)
t_max = 60 * 12
x = Array.create('x', (n, m), 'BINARY')

#Objective
H_aim = -Sum(0, m, lambda j: Sum(0, n, lambda i: x[i,j]) * (s[j] - 1.5 * Sum(0, n, lambda i: x[i,j]) + 1.5))

#Use of multiple attractions at the same time is prohibited
H_a = Constraint(Sum(0, n, lambda i: (1 - Sum(0, m, lambda j: x[i,j])) ** 2), 'alpha')

#Total travel time within the planned stay time
C = Sum(0, m, lambda j: Sum(0, n, lambda i: (c[j] * x[i,j])))
D = Sum(0, n-1, lambda t: Sum(0, m, lambda i: Sum(0, m, lambda j: (d[i,j] * x[t,i] * x[t+1,j]))))
H_b = Constraint(C + D - t_max, 'beta')

#Continuous use of the same attraction is prohibited
H_c = Constraint(Sum(0, m, lambda j: Sum(0, n-1, lambda i: x[i,j] * x[i+1,j])), 'gamma')

#Use the entrance at the beginning and end of the course
H_d = Constraint(2 - x[0,0] - x[n-1,0], 'delta')

#Energy function(I also weighted the objectives)
H = 30 * H_aim + Placeholder('alpha') * H_a + Placeholder('beta') * H_b + Placeholder('gamma') * H_c + Placeholder('delta') * H_d

#Constraint weighted hyperparameters
params = {
    'alpha': 60,
    'beta': 2,
    'gamma': 60,
    'delta': 100
}

#Generate QUBO matrix
model = H.compile()
qubo, offset = model.to_qubo(feed_dict=params)

You can get the optimized result by throwing the qubo generated in the last line to the quantum computer. There is a machine called D-Wave as a quantum computer that can be generally used, but D-Wave has a limitation on the QUBO model that can be solved, so this problem is not suitable. Hmm.

This time, I used Fujitsu Digital Annealer used in university research. The digital annealer is composed of digital circuits inspired by quantum phenomena, and is not exactly a quantum computer. So, it says, "I tried to find the optimal path of the dream country by quantum annealing", but I'm sorry it's a title scam m (_ \ ) m Please recognize that it is a simulation of quantum annealing m ( _) m

The code that uses the digital annealer cannot be posted because it uses a private library, but the result is as follows.

The total travel time is over 10 minutes, but it is acceptable. Tuning the hyperparameters should improve it. I haven't used only Pooh's honey hunt, but it's a reasonable result because it's the fourth highest cost but the lowest satisfaction.

When I draw the course on the Disney map, it looks like the following. The arrow is a straight route, so it is not accurate, but it has become a lean course that goes around the park.

Web application

I had a Disney-loving friend who cooperated with me check the validity of the course, but each time I changed the input on the local PC and recalculated, it was troublesome to send the result, so I made it a web application I tried it.

https://tdr-planner.web.app/

There is also DisneySea data so you can try both Land and Sea!

Back end

I used Google Cloud Functions because I only need one simple API that performs optimization processing and returns the course.

It is basically the same as the python implementation of [here](#solve with a quantum computer), satisfaction $ (s) $, number of attractions used $ (m) $, estimated stay time $ (t_ { max}) $ is input from HTTP request.

It's the essential annealing process, but I had to think of an alternative because I couldn't use a research machine in my private application. Thankfully, SA (Simulated Annealing) is implemented in pyqubo, so I used this on Cloud Functions.

You can easily simulate the generated QUBO matrix with the following code.

Of course, this is not quantum annealing m (_ \ _) m

from pyqubo import solve_qubo

raw_solution = solve_qubo(qubo)
decoded_solution, broken, energy = model.decode_solution(raw_solution, vartype='BINARY', feed_dict=params)

I made the API that returns the following json by shaping the result obtained in ↑ in a nice way.

{
    "totalScore": 72.0,
    "totalAttractionTime": 699,
    "totalTravelTime": 26,
    "totalTime": 725
    "course": [
        {
            "attraction": "entrance",
            "attractionId": 0,
            "waitingTime": 0,
            "requiredTime": 0.0,
            "startTime": "09:00",
            "endTime": "09:00",
            "x": 655,
            "y": 830
        },
        {
            "attraction": "Monsters Inc",
            "attractionId": 31,
            "waitingTime": 85,
            "requiredTime": 4.0,
            "startTime": "09:02",
            "endTime": "10:31",
            "x": 815,
            "y": 740
        },
        {
            "attraction": "Space mountain",
            "attractionId": 33,
            "waitingTime": 87,
            "requiredTime": 3.0,
            "startTime": "10:35",
            "endTime": "12:05",
            "x": 1005,
            "y": 600
        },
        ...
        {
            "attraction": "entrance",
            "attractionId": 0,
            "waitingTime": 0,
            "requiredTime": 0.0,
            "startTime": "21:13",
            "endTime": "21:13",
            "x": 655,
            "y": 830
        }
    ]
}

front end

I only needed the time and satisfaction input form and the ability to hit the API to display the results, so I made it quickly with Nuxt + Vuetify, which I often see recently, and hosted it as a SPA on firebase. (The favicon is Nuxt and the feeling of Vuetify is totally embarrassing ...) Ideally, I wanted to draw a route on Disney's official map. When I contacted Tokyo Disney Resort for the time being, it was normally NG. However, it should be called Disney, and he also politely responded to the playful college students who asked me to use the official map image.

Please refrain from posting the map image of Tokyo Disneyland Sea on the website of individuals. We are sorry for the inconvenience that we will give you such an answer even though you have contacted us, but please understand.

in conclusion

Since I will do quantum annealing in my thesis, my original motivation was to solve various combinatorial optimization problems and get used to it. However, I set this problem because I thought it would be nice to have a realistic and interesting problem if I could solve it anyway.

At first, I intended to do it as lightly as a break from my thesis, but my friend cooperated with me astonishingly (the person who calculated the distance), and the result was more decent than I had imagined, so I'm out. I thought I'd keep it together.

Actual Disney involves more complicated elements, so the following issues are piled up as an application.

--The waiting time is not treated as time series data (the same waiting time is used no matter what time you go) --Not considering Fastpass ――We do not consider anything other than attractions such as lunch, dinner and shows

In short, it's still a fucking app, so I'm thinking of evolving it in my spare time.

We are looking for other problems to solve with quantum annealing, so if you have any interesting ideas, please write them in the comments!