Traffic signal optimization on a square lattice with quantum annealing

The spread of intelligent transportation systems in urban cities has caused heavy computational loads, requiring a novel architecture for managing large-scale traffic. In this study, we develop a method for globally controlling traffic signals arranged on a square lattice by means of a quantum annealing machine, namely the D-Wave quantum annealer. We first formulate a signal optimization problem that minimizes the imbalance of traffic flows in two orthogonal directions. Then we reformulate this problem as an Ising Hamiltonian, which is compatible with quantum annealers. The new control method is compared with a conventional local control method for a large 50-by-50 city, and the results exhibit the superiority of our global control method in suppressing traffic imbalance over wide parameter ranges. Furthermore, the solutions to the global control method obtained with the quantum annealing machine are better than those obtained with conventional simulated annealing. In addition, we prove analytically that the local and the global control methods converge at the limit where cars have equal probabilities for turning and going straight. These results are verified with numerical experiments.

signals installed at each intersection. To analytically handle this network, we consider a simplified situation in which two states are assumed for each signal: traffic is allowed in either the north-south direction or the east-west direction. The cars moving on the lattice are assumed to choose whether to make a turn or to go straight at an intersection with a given probability. We then formulate the signal operation problem as a combinatorial optimization problem. The objective function of the formulated problem is shown to be formally consistent with the Hamiltonian of the Ising model. The Ising model is a statistical ferromagnetism physics model that represents the behavior of a spin system, and it captures the relation between the microscopic state of spins and the macroscopic phenomena of magnetic phase transitions [35][36][37][38] . Importantly, the problem reformulated by means of the Ising model, with the aid of a graph embedding technique, is compatible with the class of problems that the 2000Q accepts; hence, one can apply the quantum annealing to solve the signal optimization problem.
By reformulating the problem using Ising minimization, this study makes three contributions to signal optimization. First, by performing numerical experiments, we confirm the engineering effectiveness of the proposed method using quantum annealing. Results of experiments using a large city consisting of 50 × 50 intersections show that the proposed method achieves high quality signal operation, compared with the results of a conventional local control method 39 . The reformulated optimization problem is also solved using a classical simulated annealing method, but the quantum annealing method is found to give a better solution in a specific parameter domain. Second, a theoretical correspondence between the local and global control methods is found. We analytically show that the conventional local control is consistent with the solution of the global signal optimization problem at the limit where the probability of cars going straight is equal to the probability of them turning. This result provides a theoretical basis for the numerical prediction of the previous study 39 , where the local control is found to cause phase transitions similar to those of the Ising model. The last contribution is the knowledge gained for the cooperative operation of traffic signals. Our numerical experiments show a strong correlation between a signal and its neighboring signals. In addition, a strong temporal correlation of signals emerges, that is, the signal display at a certain time is correlated with the displays in the previous several steps. This spatiotemporal correlation becomes stronger as the straight driving probability of cars increases. Our results suggest the necessity of signal cooperation for smooth traffic flow, with variation of cooperation strength depending on the rate at which vehicles drive straight.

Results
Traffic signal optimization problem. Consider L × L (L ∈ N) roads arranged in east-west and northsouth directions with a periodic boundary condition. Each road consists of two lanes, one in each direction. Traffic signals are located at each intersection to control the flow of vehicles traveling on the roads. The signal at each node i has one of two states: σ i = +1 , which allows vehicle flow only in the north-south direction, and σ i = −1 , which allows vehicle flow only in the east-west direction. Each car goes straight through each intersection at fixed probability a ∈ [0, 1] and otherwise turns to the left or right with equal probabilities, that is, (1 − a)/2 for each direction. Figure 1 illustrates this situation.
Reference 39 shows that the number of vehicles q ij ∈ R + in the traffic lane from intersection j to i evolves according to the following difference equation: where α := 2a − 1 , and s ij ∈ {±1} is the direction of the lane from node j to i; here, s ij = +1 denotes north-south and s ij = −1 denotes east-west. We note here that q ij is normalized by the numbers of cars passing per unit of time. Precisely, in terms of the mean flux of moving cars Q av and the dimensional time unit t , t = t * /�t and In the case of σ = +1 , the vehicles coming from the horizontal direction stop, and the vehicles coming from the vertical direction go straight at the rate of a, turn right at the rate of (1 − a)/2 , and turn left at the rate of (1 − a)/2 . The rate 1 − a shown for the horizontal direction is the sum of the vehicles from the two vertical directions. In the case of σ = −1 , the roles of the vertical and horizontal directions are reversed. This problem setting is basically following Ref. 39  www.nature.com/scientificreports/ q ij = q * ij /(Q av �t) , where t * is the dimensional time and q * ij is the number of vehicles in a lane. See "Methods" for the detailed derivation of Eq. (1). We define a quantity that represents the deviation of the north-south flow and the east-west flow at each intersection i as where N (i) represents the index of the four intersections adjacent to intersection i. Equation (2) transforms Eq. (1) into a time evolution equation for the flow bias x(t) as follows: where the flow bias vector is defined as x := [x 1 , . . . , x L 2 ] ⊤ and the signal state vector is defined as σ σ σ := [σ 1 , . . . , σ L 2 ] ⊤ . The matrix A ∈ R L 2 ×L 2 is the adjacent matrix of the periodic lattice graph.
Next, we define the following objective function to evaluate traffic conditions at each time step: where the first term on the right-hand side suppresses the flow bias during the next time step at each intersection, the second term prevents the traffic signal state at each intersection from switching too frequently, and η ∈ R + is a weight parameter for determining the ratio of the two terms. The traffic signal state σ i (t) at each time step is determined so that the objective function (4) is minimized; that is, we want to find the value of σ σ σ (t) that satisfies Ising formulation and optimization. Substituting Eq. (3) into Eq. (4) gives the following representation: where c(t) is a constant term that does not include σ σ σ (t) . By defining the variables we can represent the objective function (6) as follows: Equation (10) is a quadratic form with variables {±1} , which matches the Hamiltonian form of the Ising model 35 .
Hence, solving the signal optimization problem of the objective function (4) is regarded as equivalent to the problem of finding the spin direction σ i ∈ {±1} that minimizes the Ising Hamiltonian of Eq. (10). Because the Ising Hamiltonian is compatible with the class of problems that the 2000Q accepts, the quantum annealing can be applied to solve the signal optimization problem. We use a city consisting of 50 × 50 intersections to consider the signal operation problem, and we compare the results of numerical experiments on the following three methods for traffic control: • Local control, which determines the signal display at each time step with the following local rules: Equation (11) switches the display of the traffic signals to reduce the flow bias when the magnitude of the bias becomes larger than the threshold value θ ∈ R + at each intersection. To compare the local control with the optimal control, the value of the switching parameter θ is determined such that the common objective function (4) is minimized. For details, refer to "Methods" section.
• Optimal control with simulated annealing, which reduces Eq. (10) at each time step using the simulated annealing. The simulated annealing is an algorithm for finding a solution by examining the vicinity of the current solution at each step and probabilistically determining whether it should stay in the current state or  40 for details of the simulated annealing. We used the neal library provided by D-Wave for executing this algorithm. • Optimal control with quantum annealing, which reduces Eq. (10) by using the quantum annealing with the D-Wave 2000Q. Because the problem size exceeds the size of problems that 2000Q can solve, it is subdivided by the graph partitioning technique. We used the ocean library provided by D-Wave for executing this algorithm. See "Methods" for the detailed procedure. Figure 2 shows snapshots of the signal display at time t = 100 for α = 0.8 and η = 1.0 , where α is the parameter related to vehicle's straight driving probability and η is the weight parameter in the objective function (4). The flow bias distribution at the initial time x(0) is generated as random numbers following a uniform distribution in [−5.0, 5.0] , and the signal states at the initial time σ σ σ (0) are generated as random numbers following a binomial distribution of {±1} . In Fig. 2, blue dots mean that the cars are allowed to pass in the east-west direction, and red dots mean that the cars are allowed to pass in the north-south direction. We observe the synchronization of proximity signals under optimal control (see Fig. 2b,c), while the two direction states are distributed rather uniformly under local control (see Fig. 2a). The correlation of proximity signal states is quantitatively analyzed in "Discussion". Figure 3a plots the time evolution of the Hamiltonian of Eq. (10) for each method in the case of α = 0.8 and η = 1.0 . In all three methods, the signals change rapidly over time to reduce the Hamiltonian. The value of the Hamiltonian in the steady state is the smallest in the quantum annealing method, followed by the simulated annealing method, and it is the largest under local control. That is, the optimal control using the quantum annealing exhibits the best performance among the these methods. An attempt to compare with the exact solution has also been made for the same simulation using Gurobi package. Although the full exact solution for the entire time interval was not feasible in a reasonable time because of the large number of variables (2500 variables), the Hamiltonian averaged over first three steps showed the same order of accuracy. In Fig. 3, the response of   www.nature.com/scientificreports/ the quantum annealing and the simulated annealing is more oscillatory than that of the local control. This is because the objective function of Eq. (4) only contains states up to one step ahead. An optimal value at one time is not necessarily consistent with the optimal values for long time behavior, resulting in more oscillatory response. If we use an objective function including more than two steps ahead, the oscillatory phenomenon should be suppressed, although the latter makes the formulation more complex, hindering the direct use of the quantum annealing.
We examine the effect of changing the parameter α , the vehicle's straight driving probability, on the Hamiltonian of Eq. (10). The time average of the Hamiltonian of Eq. (10), denoted as H , is plotted in Fig. 3b. As α approaches zero, the values of the Hamiltonian for the local and optimal control methods converge to a common value. This suggests that the local control gives the solution to the signal optimization problem at the limit of α → 0 . The validity of this conjecture is explored in "Discussion". In the interval of α ∈ [0.2, 0.8] , the Hamiltonian under optimal control is smaller than that under local control, showing that the optimum control method exhibits performance better than that of the local control in this range. However, in the simulated annealing method at α > 0.8 , the value of the Hamiltonian is larger than that under the local control method, suggesting that the simulated annealing does not reach the global optimal solution. Conversely, under the quantum annealing method, the value of the Hamiltonian is smaller than those under the other two methods, which means that the solution is closer to the global optimum. Here, we briefly discuss the slightly better value of H for the simulated annealing in a parameter domain of α ∈ [0.2, 0.8] , than that for the quantum annealing. In the range of large values of α , obtaining an exact solution of Eq. (4) is hard because of the high impact of the quadratic term. Actually in this parameter range, the quantum annealing gives better optimization results than the simulated annealing. On the other hand, regardless of problem to be solved, the quantum annealing generally contains stochastic fluctuations in the solutions 28,32 . When the parameter α is in the intermediate range where the difficulty inherent in the optimization problem is moderate, both the simulated annealing and the quantum annealing give high quality solutions, but the simulated annealing gives slightly better solutions than the quantum annealing because the relative strength of stochastic fluctuations is large.

Discussion
Performance analysis of quantum annealing. The performance of the D-Wave 2000Q is known to vary depending on the structure of the problem. In particular, when the matrix J in Eq. (8) has a sparse structure, the accuracy of the solution is improved 21 . To check the sparseness of our formulated problem, we examine the value of all components of J in Eq. (8). First, expanding J yields the following expression: where the number of non-zero elements in each column of A is 4, because it is equal to the number of degrees of each node in the lattice graph (see the green nodes in Fig. 4a). Also, the number of non-zero elements in each column of A ⊤ A is 9 because it coincides with the number of nodes which are connected with the reference node via two edges in the lattice graph (see the orange nodes in Fig. 4a). Thus, the number of all non-zero elements in J is expressed as 13L 2 . From this, we calculate S J (L) , the sparseness of matrix J, defined as the ratio of the number of 0-valued elements and the number of all elements in the matrix: where we confirm that S J (L) → 1 as L → ∞ . In Fig. 4b, we plot S J (L) given in Eq. (13), to show that the sparseness of matrix J increases as increasing city size. This allows us to expect that the performance of the D-Wave Local and optimal control correspondence. As shown in Fig. 5, when the parameter α of Eq. (1) is sufficiently small, the local control of Eq. (11) approaches the optimal control that is the solution of Eq. (5). When α ≈ 0 is valid, the term associated with α in Eq. (10) can be ignored, yielding Because J in Eq. (14) is a diagonal matrix, the first term σ σ σ (t) ⊤ Jσ σ σ (t) on the right-hand side of Eq. (10) is a constant that does not depend on σ σ σ . Therefore, the minimizer of H(σ σ σ (t)) is determined depending only on the sign of h in Eq. (15), that is, for all i = 1, . . . , L 2 . By transforming Eq. (16), we obtain for all i = 1, . . . , L 2 . The control method of Eq. (17) is equivalent to the local control (11) in Ref. 39 . Because α = 0 ⇔ a = 0.5 holds, this optimality means that an appropriate vehicle turning rate autonomously eases the flow bias in the local control laws. In addition, the occurrence of this magnetic transition for the signal display, stated in Ref. 39 , is consistent with the fact that the local control in Eq. (11) actually minimizes the Ising Hamiltonian in Eq. (10). However, note that the optimality of the local control is valid only when α ≈ 0 , but not when α → 1 , where the phase transition occurs.
Signal synchronization analysis. To analyze the signal correlation observed in Fig. 2, we calculate the magnetization, which is regarded as an important quantity in the Ising model: In the Ising model, this value represents the spin bias of the entire system, and it is an indicator of ferromagnetic transitions in the system. Figure 5a shows the time variation of magnetization m(t) . The value of magnetization remains small under local control, whereas it becomes significantly larger under both optimal control methods (simulated annealing and quantum annealing). For each method, at α = 0.8 , the response of the magnetization oscillates or fluctuates around zero. To confirm this observation, the time average of the magnetization of Eq. (18), denoted as m , is plotted in Fig. 5b. Here, the ferromagnetic transition at α → 1 , that is, the finite value of m , is observed for the magnetization under local control, which was originally reported in Ref. 39 . Also, under optimal control, the time average of the magnetization m takes a large value when α → 1 , which shows that a ferromagnetic transition similar to that under local control occurs under optimal control.
In addition to the ferromagnetic transition, the large amplitudes observed under optimal control are indeed a quantification of the synchronization of proximity signals observed in Fig. 2. For further analysis of this (14) J ≈ (1 + η)I, www.nature.com/scientificreports/ synchronization, we also evaluate two types of autocorrelation functions. Figure 6a shows the autocorrelation function obtained from the time-series data of the signal state σ i (t) for t ∈ [0, 200] . Here, the autocorrelation function is computed at all intersections, and the average value is displayed in Fig. 6a. Under local control, there is a negative correlation peak around t = 3 , which means that the signals switch approximately every 3 time steps. In contrast, under optimal control, the negative correlation peak is in the interval of t = [10,15] steps, and the same state is maintained for a time longer than that under local control. In general, excessive signal switching is undesirable from a traffic engineering standpoint, and the optimization-based method overcomes this issue.
Next, Fig. 6b shows the correlation between the display of signals at one intersection and another intersection, with the distance between the intersections as a parameter. Here, the correlation function is calculated for all the intersections for fixed time t = 100 , and the average value thereof is plotted. In Fig. 6b, the distance is normalized to make the distance of adjacent intersections equal to 1. There is almost no correlation between adjacent signals under local control, while there is a positive correlation of up to 4-6 adjacent intersections under optimal control. Then, we extract quantities from these correlation functions to investigate the effect of α . First, considering that both the temporal and spatial autocorrelations in Fig. 6 decay while oscillating, both functions are fitted with the following equation: where represents the damping rate coefficient, ω represents the vibration frequency coefficient, and z ∈ R + represents different variables, i.e., the time t for the time autocorrelation function and the distance between intersections for the spatial autocorrelation function. Figure 7a plots ω values obtained by fitting Eq. (19) to the time autocorrelation, as a function of α . Under local control, the vibration frequency is ω ≈ 1 regardless of the value of α , while ω decreases as increasing α under optimal control. This suggests that the frequency of signal switching reduces as the vehicle straight driving rate increases in order to guarantee optimality. In view of the large difference in ω between the local control and the optimization-based controls for large values of α , we expect that optimization-based signal controls are particularly effective in preventing excessive switching for high vehicle straight driving rates. Next, we show in Fig. 7b the value of obtained by fitting Eq. (19) to the spatial   Effect of parameter η. Here we examine the effect of parameter η , which controls the priority of the smoothness of the entire traffic flow to the signal switching frequency in the objective function of Eq. (4). The objective function is designed such that the priority is to smooth the flow of the car for small value of η , and inversely for a large value of η , preventing excessive signal switching is prior. The time average of the objective function H is obtained for various values of η ; η ∈ {0.125, 0.25, 0.5, 1, 2, 4, 8} . We show the results in Fig. 8 for the simulated annealing ( H SA ) and the quantum annealing ( H QA ), where the ratios H QA /H SA is plotted; H QA /H SA < 1 means that the quantum annealing is better than the simulated annealing, and vice versa. The quantum annealing method shows better performances for η larger than 0.5, and the simulated annealing is better for η smaller than 0.5. This suggests that the quantum annealing works better when the priority is on preventing excessive signal switching.

Future improvements.
Here we discuss three possible improvements of the results obtained in this study.
First, we expect that the solution is improved by using the most recently released D-Wave's machine. The D-Wave machine used in this study has 2048 qubits that are connected with the chimera structure, in which closely connected 8-qubit units are arranged 41 . Since the chimera structure is sparser than the fully-connected structure, the representation of arbitrary Ising problems requires a process called minor-embedding to map logical variables to physical qubits. This process however significantly reduces the number of available qubits, and also deteriorates the computational accuracy. Very recently, D-Wave has updated the hardware with a new graph structure called the pegasus structure. The number of qubits has increased from 2048 to 5024, and the maximum number of connections in the graph structure has increased from 8 to 15 42 . These improvements allow us to deal with much larger problems, and to realize efficient embedding with sacrificing less physical qubits, than the previous D-Wave machine. For the proposed method, this enhancement will significantly reduce the number of divisions of Hamiltonian (see "Methods" for details), which contributes to fast and high-accurate computations.
Next, adjusting the hyper-parameters of the solver could improve the performance of our method. The D-Wave machine contains a few hyper-parameters, such as, number of samplings, chain strength, and post processing. We left most of the parameters at their default values because we focus on examining ability of the quantum annealing to solve the traffic signal optimization problem. However, as these hyperparameters affect the optimized result, more careful tuning of these parameters may achieve faster and more accurate calculations. An error mitigation scheme proposed in Ref. 43 could enhance the performance. We remark here that the problem formulated in this study has the form of Ising model that is solvable by using several dedicated computers other than the D-Wave machine [20][21][22][23] . Since the development of these dedicated machines is expected to further accelerate, the proposed framework for traffic flow control will be more generally available in the future.
We finally discuss a further improvement toward application to a real city. The parameter α in our model is expected to be relatively large in a city with many rational players, because in a grid-like city, each vehicle can reach its destination from any starting point with only one right or left turn. In our experiments, the size of the grid L is empirically determined as L = 50 so that it would be comparable to the size of typical grid cities in the world (Kyoto, Japan; Barcelona, Spain; La Plata, Argentina, etc.). It is however desirable to identify these parameters in advance using real-world data. Since the constant probability of each vehicle driving straight ahead and the grid topology of the city are both idealistic assumptions, our traffic signal control method has to be further improved by relaxing these assumptions. (1). Let q * ij (t) be the numbers of cars that exist between the intersections i and j and t be the minimum time interval at which a signal is switched. We denote by Q av the average flow rate of cars passing during t . Then the change in the numbers of cars from time t * to the next time t * + t is represented as By normalizing Eq. (20) with t := t * /�t and q := q * /(Q av �t) , we obtain the following equation: which is essentially identical to Eq. (1). In this paper, we consider the result of solving Eq. (21). The dimensional time and the actual numbers of cars are recovered with inverse transformation of the above normalization.
We remark on the numbers of cars, speed, and minimum signal switching interval on the model. First, a solution of the model in Eq. (1) is valid for an arbitrary numbers of cars. For example, the solution for q * (0) = γ q * (0) with some γ ∈ R + , is obtained by setting Q av = γ Q av because the average flow rate is defined by "vehicle density" × "average speed". Second, let us consider the case of the average speed multiplied by γ . The average flow rate Q av should then be Q av = γ Q av , while q * (0) remains the same. Therefore, while Eq. (21) normalized by Q av does not apparently change, the initial value should be appropriately adjusted q ij (0) = q * 0 ij /(Q av �t) = q ij (0)/γ . Similarly, for the case of �t := γ �t , the normalized Eq. (21) does not apparently change, but the initial value should be q ij (0) = q * 0 ij /(Q av �t) = q ij (0)/γ.
Parameter identification for objective function. As stated in "Discussion", a direct correspondence between the optimal control and local control is established for small values of α , with the apparent relation θ = η between the local control switching constant θ in Eq. (11) and the optimal control weight parameter η in Eq. (4). To make a systematic comparison for an arbitrary value of α , however, we still need to construct a protocol to determine the values of θ and η . The strategy is described as follows. Given a value of η , we select a value of θ , denoted by θ , from a candidate set via the following auxiliary numerical analysis: 1. For one value of θ in the set , numerical simulation using local control (11) is performed to obtain time series data x(t) and σ σ σ (t) . The value of the objective function (4) using the given η is calculated from the obtained time series data. This time average is denoted as H (θ).

2.
Step 1 is performed for all θ in to find θ that minimizes the time average H , that is, θ = argmin θ∈�H (θ).
We plot the result of the above procedure in Fig. 9. Figure 9a shows H against θ when η is fixed as η = 1.0 . When α = 0 , H is a convex function and indeed θ ≈ η is satisfied. However, for larger values of α , H becomes non-convex, and particularly for α = 0.995 , the relation θ = η no longer holds. Figure 9b shows the value of θ that minimizes H versus η for the interval η ∈ [0.0, 3.0] . When α = 0 , the linear relation θ = η approximately holds, but when α = 0 , this relation breaks down and some discontinuities appear. These discontinuities correspond to the changes in the local minima observed in Fig. 9a.

Implementation on D-wave 2000Q.
All experiments in this study are conducted on a Linux computer with 64 GB of memory and a clock speed of 3.70 GHz. All methods are implemented using the programming language Python (version 3.7). www.nature.com/scientificreports/ We use DW_2000Q_VFYC_5 as a machine solver with the aid of D-Wave's ocean library for the actual implementation. Here, the VFYC solver partially emulates some qubits that are temporarily unavailable because of hardware failures 44 . The number of samplings is specified through a parameter named num_reads, which we set 100 in all experiments. The validity of this parameter setting is confirmed by preliminary experiments using several candidate parameters.
For embedding the logical variables to the physical qubits on the D-Wave machine, we use a tool called minorminer (Apache license 2.0), which is a heuristic embedding method in ocean library 45 . We perform embedding operation every time when the problem is sent to 2000Q in order to average out the bias of the embedding quality.
For the simulated annealing, neal solver in ocean library is used. This solver also allows us to specify the number of samplings through a parameter called num_reads, which we set 100, i.e., the same value as the one in the quantum annealing.
For the chimera structure in the 2000Q, N 2 /4 physical qubits are necessary for embedding N-variable problem for the worst case, which means that the maximum number of variables that the 2000Q is capable of handling is as small as 64. This implies that L 2 ≤ 64 ⇔ L ≤ 8 must be satisfied for the number of roads L in our problem setting. A method exists for solving a problem that exceeds the size limitation: to divide the Hamiltonian variable of Eq. (10) into several groups and minimize the approximate Hamiltonian for each group. We define the traffic signal state vector of the jth group as σ σ σ j := [σ i 1 , σ i 2 , . . . , σ i m ] ⊤ , where i 1 , i 2 , . . . , i m are subscripts of variables included in the jth group. Then, we define the Hamiltonian of the group j as where J jj is a matrix extracting the (j, j) th components of matrix J in Eq. (10). Similarly, h j is a vector obtained by extracting the jth component of h. The index j represents the set of variables not belonging to group j. One naive approximation is to regard the variables outside group j as constant. This allows the annealing machine to deal with a Hamiltonian exceeding the limitation, but at the same time this approach degrades the control performance. To reduce such errors, the variables having a large interaction should be in the same group, and the variable interaction between different groups should be small. Such a problem is called a graph partitioning problem, which is known to be an NP-hard problem, but there are some approximation methods with adequate accuracy. For the actual implementation, we used the Metis software (Apache license 2.0), which is a widely used solver for graph partitioning problems, to break up the large-scale problem into several groups having fewer than 64 variables 46 . Figure 10 shows the result of the graph partitioning of the city of L = 50 into 42 groups using Metis, where we certainly see that the adjacent intersections, i.e., the strongly interacting variables, are included in the same group.
We evaluate the effect of this graph partitioning method. To this end, the simulated annealing optimization on a system with L = 8 is performed, and results are compared between the no-partitioned and quadruple-partitioned cases. Figure 11 shows the time average of the objective function for various α . The values with partitioning are larger than those without partitioning, where the difference between these values represents the error caused by partitioning. The error increases with larger α , indicating that the larger the straight driving rate of vehicles, the more the partitioning has a negative effect. In Fig. 3, the quantum annealing has advantage at large α , and thus a higher performance signal control should be achieved once the D-wave without partitioning is realized. (22) H j (σ σ σ j (t)) := σ σ σ j (t) ⊤ J jj σ σ σ j (t) + (h j + σ σ σ¯j(t) ⊤ J¯j j )σ σ σ j (t),