Here, we explain the details of the ARP. In the event of a disaster, such as an earthquake, an on-site witnessing will be conducted by adjusters dispatched by the insurance company to damaged buildings to confirm the damage. Particularly, in the case of large-scale disasters, a single insurance company has to assess tens of thousands of cases, requiring the dispatch of hundreds of personnel from the office each day. Therefore, an optimization of the route (i.e., minimizing the time required for witnessing, including travel time) is an extremely important problem from both a time and cost perspective. This problem can be regarded as the vehicle routing problem (VRP) [17] because multiple persons visit multiple locations and optimize their routes. In the setting where multiple vehicles visit multiple customers, VRP is the problem of finding a route that minimizes the sum of costs among the routes that satisfy the following constraints:

• Each customer is visited exactly once by exactly one vehicle.

Here, the interpretation of the cost as the time required for witnessing, the customer as the building to be witnessed, the vehicle as the adjuster, and depot as the dispatcher’s office, respectively, VRP is basically consistent with the problem considered. However, there are some additional constraints in the ARP. First, the adjuster who can witness the building is determined according to the type of the building. The level of difficulty of witnessing a building depends on the type of building materials. Basically, three types of building materials (wood, steel, and reinforced concrete) are considered. Therefore, there are three types of buildings with three levels of difficulty: “hard,” “normal,” and “easy.” Adjusters cannot witness all of these types of buildings, and the buildings that can be witnessed are determined according to their capabilities. We refer to this capability as the adjuster’s spec. Three types of specs are defined, “high,” “middle,” and “low,” from the highest to lowest. Buildings with a hard difficulty level can be witnessed by an adjuster with a high spec, and buildings with a normal difficulty level can be witnessed by an adjuster with a high or middle spec. In optimizing the route, we should consider the matching of the adjuster’s spec and the difficulty level of witnessing the building. In addition, it is necessary to be able to start witnessing at the desired time submitted in advance by each building. This desired time is specified by time zone, for example, AM, PM1, and PM2, and each adjuster must be able to witness to satisfy the desired time. We also consider the time required for witnessing, which is determined according to the total area of the building (i.e., it increases in proportion to the area), and the total operating time includes the total travel time and the time required for witnessing at each building. Furthermore, even when all the aforementioned constraints are satisfied, if the number of buildings witnessed by one adjuster is large, the witnessing cannot be completed in time. Therefore, it is necessary to limit the number of buildings one adjuster can witness.

In summary, the ARP is the problem of finding a route that minimizes the total operating time among routes that satisfy the following constraints:

1. Every adjuster leaves from certain office and must return to the same office at the end of the last witnessing.

As mentioned previously, the usual VRP does not satisfy constraints 3, 4, and 5. Therefore, the ARP is an extension of the VRP. Other VRP extensions include the vehicle routing problem with time windows (VRPTW) satisfying constraint 3 [18,19] and the capacitated vehicle routing problem (CVRP) satisfying constraint 5 [20].

As the VRP is in the class NP, we consider solving the ARP using a quantum computer (quantum annealing).

In this study, we propose a method for solving the ARP using quantum annealing. The following decision variables are needed for route optimization: x n , i k = 1 when adjuster k witness building n for the i - th time 0 otherwise . (1)

Let K be the number of adjusters, N the number of buildings, and I the number of buildings witnessed by one adjuster. The size of the optimization problem is K × N × I. Quantum annealing can only solve problems of small size when three parameters (k, n, i) are introduced because the size of the problem is limited by hardware constraints, such as the number of qubits and connectivity [21]. Further information about these constraints is presented in the following paragraph. We divided the ARP into two subproblems to reduce the number of decision variables: clustering and routing phases. The clustering phase determines which building each adjuster will witness, and the routing phase determines the order each adjuster will witness. A method for dividing the VRP into clustering and routing phases has already been proposed [22,23], and we adopt a similar approach for the ARP. To apply quantum annealing, we formulate each phase in QUBO.

Constraints 4 and 5 are considered in the clustering phase because the building to be witnessed by the adjuster is determined in this phase. Specifically, for constraint 4, we interpret the adjuster and building matching as bipartite graph matching and introduce constraints into the clustering phase. For constraint 5, we introduce an element count constraint into the clustering. Constraint 3 is considered in the routing phase because it affects the order of witnessing by the adjuster. Although the routing phase optimizes the route of a single adjuster to witness multiple buildings, and thus can be regarded the same as the traveling salesman problem (TSP) as an optimization problem, we introduce the concept of time into the TSP and introduce constraints to realize routing that satisfies the desired time. When the problem is divided into two phases, the probability that no feasible solution exists in the routing phase increases if the desired time zones of the buildings are concentrated in the same cluster in the clustering phase or if only buildings with large areas are clustered together. Therefore, we introduce constraints in the clustering phase such that they are properly distributed to ensure the diversity of cluster elements. Furthermore, we conducted numerical experiments using D-Wave to confirm the correctness of our formulation.

Our contributions can be summarized as follows:

• Applications of quantum annealing to the real-world problem: Our ARP is based on the VRP and extended to a realistic setting. We formulated the problem in a form amenable to quantum annealing and confirmed correctness through numerical experiments.

We review the theoretical backgrounds of quantum annealing and D-Wave machines that implement quantum annealing. Several solvers can perform QUBO, such as Digital Annealer [25] and CMOS annealer [26]. In particular, quantum annealing machines have been increasingly applied in recent years, and there are high expectations that both the speed of solution seeking and the quality of the solution will exceed those of classical computers.

Quantum annealing is a meta-heuristic algorithm for optimization problems introducing quantum effects, proposed by Kadowaki and Nishimori [2]. It can solve optimization problems in QUBO or Ising form. Because QUBO belongs to the NP-hard, quantum annealing can be applied to an extremely wide range of optimization problems.

In quantum annealing, the objective function of the optimization problem under consideration is represented by a quantum system Hamiltonian H ̂ obj . Then, a quantum fluctuation H ̂ q is introduced, and the ground state of H ̂ obj is searched. Thus, quantum annealing deals with the following Hamiltonian H ̂ qa : H ̂ qa t = A t H ̂ obj + B t H ̂ q , (2) where A(t)⩾0, B(t)⩾0 and the control shall be performed as in A(0) = 0, B(τ) = 0 with the annealing start time of t = 0 and end time of t = τ. Let ℏ be the reduced Planck constant, ψ ( t ) the spin wavefunction ( k ( t ) means the kth excited state), Δ _k (t) the energy gap between the state k and ground state, and s = t/τ. By time-evolving H ̂ qa with the Schrödinger equation, i ℏ d d t ψ t = H ̂ qa t ψ t , (3) and with the asymptotic expansion of adiabatic condition ² , τ ≫ k s | H ̂ ̇ qa s | 0 s Δ k s 2 , (4) H ̂ qa follows the ground state at each moment, and finally, the ground state of H ̂ obj is obtained. As the ground state of H ̂ obj is the solution to the original optimization problem, it can be solved using quantum annealing. This is an outline of the algorithm.

Physical devices that enable quantum annealing have been developed by D-Wave Systems [3], which recently launched a quantum cloud computing service called D-Wave Leap. In the D-Wave machine, the target quantum system is a spin system, and the Ising model and transverse magnetic field are used as H ̂ obj and H ̂ q , respectively. Quantum annealing machines are more stable than quantum gate-based machines in controlling quantum devices such as noise, and the machine at the first time of the implementation announcement had only 128 qubits (D-Wave One), increasing to 2,048 qubits in 2017 (D-Wave 2000Q) and 5,436 qubits in September 2020 (Advantage).

As the computational principle of quantum annealing is the transverse-field Ising model, if the combinatorial optimization problem we wish to solve can be formulated in the Ising model, the solution can be obtained by providing information on its coefficients as input. Conversely, when we design a combinatorial optimization problem for a real-world problem, decision variables such as x i ∈ 0,1 are often used, and the input forms of the objective function that D-Wave can optimize correspond to both of these, called the Ising and QUBO forms, respectively, and they are equivalent by variable transformation.

The spin–spin interaction in the transverse-field Ising model assumes a fully connected state. Therefore, realizing the fully connected state between spins as a physical device must construct a complete graph with all spins as nodes on the hardware. Because it is currently difficult to construct such a graph, D-Wave introduces a special graph for the connection of spins, called a chimera graph in D-Wave 2000Q and a Pegasus graph in Advantage [21]. Whatever optimization problem we design, embedding operations are essentially required when performing optimization computations with D-Wave, and the number of installed qubits and the size of the problem to be solved do not generally match.

Optimal annealing time τ should be long enough to give a sufficiently high probability of finding the ground state (Eq. 4). At the same time, it should be as short as possible for practical purposes. Therefore, τ is a trade-off value. Since the optimal τ depends on the structure of the problem we wish to solve [27], a parameter study should essentially be performed. However, as discussed in Section 3, we used D-Wave’s hybrid solver system (HSS) for the solver, where τ cannot be treated as a variable parameter, so we experimented numerically using τ as a fixed value of 20 µs.

If the clustering target is a building and each cluster is applied to each adjuster, clustering is synonymous with determining the group of buildings to be witnessed by each adjuster. Thus, defining x n k ∈ { 0,1 } is a decision variable, 1 if a building n belongs to a certain cluster k and 0 if it does not. The simplest Hamiltonian that classifies buildings close to each other into the same cluster can be described as follows: H vanilla x = ∑ k = 0 K − 1 ∑ u v ∈ E D u v x u k x v k + ∑ n = 0 N − 1 ∑ k = 0 K − 1 x n k − 1 2 , (5) where N is the total number of buildings, K is the total number of clusters, and D _uv is the weight (distance) between u and v connected by an edge E. The first term is the distance cost term, and the second term is the constraint term for any building to always belong to one of the clusters.

In ARP clustering, it is necessary to add a cost term that accounts for variance (H _div), a difficulty-matching constraint term between adjuster and building (H _matching), and a constraint term on the number of elements in the cluster (H _num). Thus, the Hamiltonian of our clustering takes the following form: H clustering = H vanilla + H div + H matching + H num . (6) In the following sections, we show how these three conditions can be formulated.

After the clustering phase is completed, the goal is to find the optimal route in each cluster. Following Feld et al. [23], we solve the TSP for every generated cluster. However, we need to find the shortest time route that satisfies the desired time for each building. Therefore, we should introduce time into TSP. Several studies have introduced time to the TSP [30,31]. We adopt the time-scheduled traveling salesman problem (TS-TSP) proposed by Irie et al. [32] as a route optimization model. In addition, we introduce the constraint for the desired time in the TS-TSP. We begin with a brief review of the TS-TSP in the following section and describe the Hamiltonian formulation of our routing phase.

In the numerical experiments, multiple cases were set up and verified for the ARP. The details of the settings are described as follows. The target area was chosen to be near Kawajiri Station in Kumamoto Prefecture, Japan, as an example of major earthquake damage in the past (Figure 1). First, we defined the size of the problem in the clustering phase as one instance such that 100 visits are divided into 20 clusters. The instances were configured with information on the buildings to be visited (location, distance between buildings, difficulty level for assessment, required time, and desired time for witnessing) and the percentage of corresponding adjuster specs for each instance. To account for the effect of the different geographic distributions of instances, we prepared 10 instances, as shown in Figure 2. As a practical matter, the numbers of adjusters and their specs are limited, and considering the matching constraints between specs and difficulty, it becomes obvious that if there is a bias in the difficulty levels, a feasible solution in the routing will not be obtained. Therefore, the adjuster’s spec ratio and the assessment difficulty ratio were set to be completely common (fixed) in each instance. In contrast, the time required for witnessing and the desired time for witnessing were set to change from case to case. The common settings for case-independent instances are as in Table 1. Furthermore, we defined the base case as shown in Table 2 for the ratio of time desired for witnessing. For each of the settings in Table 2, we defined cases as a derivation of the base case with each of the parameters changed except the depot location. Numerical experiments were also conducted for these cases (details are given in Section 3.5).

Furthermore, by setting up numerical experiments as described previously, H _div and H _matching in Eq. 13 can be written down specifically, respectively. Thus, the Hamiltonian H clustering sim ( x ) of the clustering phase in the numerical experiment is as follows: H clustering sim x       = A ∑ k = 0 K − 1 ∑ u v ∈ E D u v x u k x v k + B ∑ n = 0 N − 1 1 − ∑ k = 0 K − 1 x n k 2 + D ∑ k = 0 K − 1 M − ∑ n = 0 N − 1 x n k 2 + C ∑ k = 0 K − 1 ∑ n ∈ AM e n x n k 2 + ∑ n ∈ PM 1 e n x n k 2 + ∑ n ∈ PM 2 e n x n k 2 + E ∑ n ∈ Hard ∑ k ∈ High x n k − 1 2 + ∑ n ∈ Normal ∑ k ∈ High ∪ Middle x n k − 1 2 , (21) where N = 100, K = 20, and M = 5 correspond to the total number of visits, the total number of clusters, and the number of elements per cluster, respectively, and D _μν represent the distance matrix between visits μ and ν. We used A = 1, B = D = E = max(D _uv ) × K × 2, C = 60 for each coefficient.

Similarly, the desired time constraint in Eq. 20 can be specifically written down, and the Hamiltonian H route sim ( x ) for the routing phase in the numerical experiment becomes H route sim x = ∑ 1 ⩽ a ≠ b ⩽ N 1 ⩽ τ ⩽ T − 1 n a b − μ ρ × x τ + e a b , a x τ , b + λ ∑ 1 ⩽ δ τ ⩽ e a b − 1 x τ + δ τ , a x τ , b + λ ∑ 1 ⩽ a < b ⩽ N 1 ⩽ τ ⩽ T x τ , a x τ , b + ∑ 1 ⩽ a ⩽ N 1 ⩽ τ ≠ τ ′ ⩽ T x τ , a x τ ′ , a + λ ∑ 1 ⩽ a ⩽ N 1 < τ ⩽ e a , d e p o t − 1 x τ , a + λ ∑ 1 ⩽ a ⩽ N T − e depot , a ⩽ τ ⩽ T x τ , a + λ { ∑ a ∈ AM ∑ τ ∈ AM x τ , a − 1 2 + ∑ b ∈ PM 1 ∑ τ ∈ PM 1 x τ , b − 1 2 + ∑ c ∈ PM 2 ∑ τ ∈ PM 2 x τ , c − 1 2 } , (22) where N = 5 and T = 32 ⁴ correspond to the number of buildings per cluster and the number of time intervals, respectively, and ∑ _a∈AM means taking the sum for building a such that the desired time for witnessing is the AM zone (the same for PM1 and PM2). n _ab in the first term can be obtained using the distance matrix D _ab and Eq. 15. We took λ = 40 for the coefficient.

Here are some of the clustering results for each instance in the base case using the Hamiltonian in Eq. 21, accounting for the matchings of adjusters and buildings and dispersion in the required and desired time for witnessing.

Figure 3 shows ⁵ a portion of the clustering results for instance No. 1 of the base case. The red point cloud in Figure 2, an example of an urban area, color-coded by clusters to visualize the location relationships and the data for each of the same clusters, is shown in Table 3. The distances between buildings in the cluster are short because of the high density near the depot, but the buildings closest to each other are necessarily not in the same cluster, which means the difficulty and spec correspondence and dispersion of desired and required time for witnessing are considered. For clusters other than those shown in Figure 3 and Table 3, the matching conditions for adjusters (clusters) and buildings were satisfied, and the required time and desired time for witnessing were appropriately distributed. There was no solution that did not satisfy the constraints (we call such a solution “break solution”) for either the SA or HSS solutions. No break solution occurred in all other instances, including sparsely populated areas, confirming that a feasible solution was obtained during the entire clustering phase.

Here are some of the results of route optimization using Eq. 22 for the results of the clustering described in Section 3.2. In the routing phase, unlike the clustering phase, some cases occurred where feasible solutions were not obtained, so we also addressed such solutions.

Figure 4 visualizes the solution of route optimization for one cluster belonging to instance No. 1, corresponding to the top five clusters in Tables 3, 4 summarizes the travel time in the travel order for each visit. As shown in the chart, the route departs from the depot and returns to the depot within 8 h and is such that the desired time for witnessing each visit is satisfied. In addition, it was confirmed that all the route optimization solutions for the other clusters in instance No. 1 are also feasible.

Conversely, only in instance No. 9, the brown point cloud in Figure 2, break solutions occurred in 2 out of 20 clusters. Figure 5 visualizes one of the clusters in instance No. 9, which is one of the solutions for which the constraint was broken. Table 5 summarizes the travel time for each building in order of travel; the desired time zone for witnessing the fourth visit is AM, but as a solution, it arrives at PM1, which breaks the desired time constraint. Similarly, for the other break solution, the break of the desired time constraint occurred in a cluster in a sparsely populated area far from the depot, resulting in a solution with a relatively long total travel distance. For these solutions, optimization with only the constraint term after excluding the cost term, the first term in Eq. 22, did not break, suggesting that the trade-off between travel time minimization and desired time constraints is more pronounced in the case of long total travel distances.

In the base case, we measured processing times in the two phases of the numerical experiment, the clustering and routing phases, and compared the differences between the solvers. The processing time was defined as the QUBO computation and solving time for both clustering and routing, both for one instance (for clustering, the process of dividing 100 buildings into 20 clusters of 5 elements each; for routing, the routing process for 20 clusters). The results are shown in Table 6. The processing time variation in clustering is small for HSS but large for routing. HSS, including quantum annealing, is superior in clustering, whereas SA is superior in the routing phase. The reason for these differences is that HSS performs quantum annealing while dividing QUBO according to the problem size, which is effective for large problems. However, extra processing time is presumably required for too small problems. The advantages of HSS will become even more effective when the problem size per processing unit increases during the routing phase.

For the break solution in the routing described previously, no feasible solution was obtained by HSS instead of SA. The break solution occurred even when the entire routing phase was solved with HSS from the beginning, and the same was true when both clustering and routing phases were solved consistently with HSS. Focusing especially on the clustering phase, where the superiority of the HSS solution was confirmed in terms of processing time, the percentage of solutions with the smaller energy value is shown in Table 7. The largest energy difference occurred at instance No. 6, yellow dots in Figure 2, where the HSS solution was 4.6% smaller than the SA solution in terms of magnitude. However, this level of energy difference is considered to have a small impact on the routing phase solution because there are cases where the break solution occurs in the routing phase, regardless of whether the HSS or SA solution is used in the clustering phase.

In this paper, clustering is performed in the form of assigning a building to each adjuster, so even if a feasible solution is obtained in the clustering phase, the impact of the solution will strongly affect the next routing phase. To reduce this effect, for example, let us consider the case where the number of elements in one cluster increased. In this form of TS-VRP, a route optimization in which several people visit dozens of buildings, it may be possible to reduce the number of break solutions in the HSS solver and reduce the processing time. Mitigating the constraints that the solution of the clustering phase imposes on the next routing phase is a topic for future work.

A list of changes from the base case for each case is shown in Table 8. The application of the results of five numerical experiments, each with SA for all cases, is shown in Table 9. Here, the processing time is calculated by the total processing time of the clustering and routing phases according to Table 6. The percentage of break solutions means the percentage of the total 200 route optimization solutions. In all cases, the changes in each setting value affected the processing time and solution quality. However, especially for case #4, although the processing time was reduced, the quality dropped significantly due to many broken solutions in the final routing phase caused by the increase in unit time in route optimization. There are instances of urban and sparsely populated areas in the numerical experiments. Particularly, in the former, where buildings to visit are relatively close to each other, a too large unit time would have a negative impact on the routing optimization. Therefore, an appropriate value for the problem setting should be set like the unit time of route optimization. It was also found that when there is no significant bias in the desired time zone, as in #1 and #2, a feasible solution is sufficient in the area within a radius of 16 km. For cases such as #5–#7, where there is a significant bias in the desired time zone, the Q2Q approach requires measures to ease the constraints on the routing solution by adjusting the size of the solution passed between the two phases, as described in Section 3.4. As also discussed in Section 3.4, the HSS has advantages, at least in clustering. Therefore, when HSS is used for severe cases such as #2, #7, and #8, it may be possible to suppress break solutions.

Conversely, it is also important in real-world problems to consider multiple targets for optimization. In this paper, we aimed to optimize total travel time and total witnessing time, but various optimization targets depend on the needs, such as optimization to minimize the number of adjusters to be dispatched or optimization to complete all witnessing in a minimum number of operating days at the expense of operating time. These will be discussed in future prospects.