A Three-Stage Annealing Method Solving Slot-Placement Problems Using an Ising Machine

Ising machines are promising alternatives to solve combinatorial optimization problems, which search for their quasi-optimal solutions with high speed and high accuracy. However, the obtained solution much depends on the initial spin states, since the computation time is finite. Moreover, due to their probabilistic nature, they cannot always satisfy the constraints given to combinatorial optimization problems. In this paper, we propose a three-stage annealing method, targeting a slot-placement problem as a typical but difficult example of combinatorial optimization problems. The proposed method is composed of an initial process, an annealing process, and a correction process. The initial process and the correction process are executed by a classical computer while the annealing process is executed by an Ising machine. In the initial process, we give initial spin values that lead to a relatively good solution to the combinatorial optimization problem, which satisfies the given constraints. Then, the annealing process is executed by an Ising machine, and the solution obtained by the annealing process is further corrected to satisfy the constraints. The experimental results demonstrate that the proposed method reduces a minimum total weighted wiring length by 0.0898%–2.45% on average depending on the initial process methods used, compared to the existing method. The mean total weighted wiring length is reduced by 2.79%–7.08% on average depending on the initial process methods used.


I. INTRODUCTION A. ISING MACHINES
A combinatorial optimization problem is a problem to find a combination of variables that maximizes or minimizes an objective function while satisfying its given constraints. In combinatorial optimization problems, there exist NP-hard problems, in which it is very difficult to search for an optimal solution by using a conventional von Neumann architecture machine.
In recent years, various Ising machines, non-von Neuman architecture machine, using the Ising model [1], [2] have been studied to efficiently solve combinatorial optimization problems [3]- [12]. Ising machines search for their quasi-optimal solution with high speed and high accuracy, where a quasi-optimal solution refers to a solution that does not always maximize or minimize the objective function.
The associate editor coordinating the review of this manuscript and approving it for publication was Abdullah Iliyasu .
When we solve a combinatorial optimization problem using an Ising machine, we map it onto a model in statistical mechanics called an Ising model, or its equivalent quadratic unconstrained binary optimization (QUBO) model. Nowadays, various combinatorial optimization problems are mapped to Ising models or QUBO models and solved by an Ising machine such as in [13]- [15].

B. SLOT-PLACEMENT PROBLEM
The slot-placement problem [16], [17] is known as one of those difficult combinatorial optimization problems. Given lattice slots and several items connected by a certain number of wires, the goal of the slot-placement problem is to find an optimal item assignment to lattice slots.
In the slot-placement problem, every item must be assigned to a single slot and no more than one item cannot be assigned to a single slot. This is called a slot-placement constraint. Several items are connected by a certain number of wires. The wiring length is given by the Manhattan distance between the VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ slots where the items are assigned. Then we can define the total weighted wiring length over all the item pairs with wires. Minimizing the total weighted wiring length is the objective of the slot-placement problem (See Section III in detail). The slot-placement problem can be applied to many practical optimization problems such as labor shift scheduling and FPGA logic-block placement [18], [19]. It can be considered to be a variation of quadratic assignment problems, which are known as an NP-hard problem [20], and has been solved using simulated annealing (SA) [21] and branch-and-bound methods [22]. A method for mapping the slot-placement problem to the QUBO model and solving it using an Ising machine has also been proposed [16].
In this paper, we pick up the slot-placement problem as a typical but difficult example of combinatorial optimization problems for Ising machines.

C. HOW TO UTILIZE AN ISING MACHINE SOLVING COMBINATORIAL OPTIMIZATION PROBLEMS
Now we consider solving a combinatorial optimization problem by using an Ising machine. In practice, the obtained solution depends on the initial spin values and the optimal solution is not always obtained, since the computational time is finite. Moreover, the obtained solution can be a quasi-optimal solution and hence it cannot always satisfy the constraints given to the original combinatorial optimization problem. How to tackle those problems is one of the key issues to utilize Ising machines efficiently.
In [16], after the quasi-optimal solutions are obtained by using an Ising machine targeting the slot-placement problem, they are improved as post-processing using a classical computer, so that they satisfy the slot-placement constraint. However, no process is applied before the annealing process of the Ising machine, and all initial spin values are input as zero. As far as we know, there exist no previous researches where any initial process is applied before an Ising-machinebased annealing process.

D. OUR PROPOSAL
In this paper, we propose a three-stage annealing method solving the slot-placement problem. The proposed method is composed of three processes: an initial process, an annealing process, and a correction process. The initial process and the correction process are executed by a classical computer while the annealing process is executed by an Ising machine. In the initial process, we give initial values that give a relatively good solution to the slot-placement problem, 1 which satisfies the slot-placement constraints. Then, the annealing process is just executed by an Ising machine as usual. Finally, the correction process corrects the quasi-optimal solution 1 We can input initial values to spins in an Ising machine as follows: In semiconductor-based Ising machines [3], [8], [23], initial spin values are directly given to them using their API environment; and, in quantum annealers like D-Wave 2000Q [24], initial values can be given to spins using external magnetic fields. obtained by the annealing process so that it can satisfy the slot-placement constraints.
As the initial process, we propose three initial process methods: a pair-wise exchange method, a random exchange method, and a cluster-growth method. Given a random feasible slot-placement solution, the pair-wise exchange method repeatedly exchanges the paired items or moves the item to an empty slot which most reduces the objective function, until no further improvement is seen. The random exchange method randomly exchanges the paired items or moves the item to an empty slot for the fixed number of times, if improved. In the cluster-growth method, we first place the item which is most connected to other items to the center of the lattice slots. Then we place an item connected most to the items already assigned to a neighboring slot in a one-byone manner. The pair-wise exchange method gives a good initial solution but requires much time. The random exchange method runs within a limited time and gives an initial solution depending on its execution time. The cluster-growth method runs very fast but it may not give a good solution. Which method is better is depends on the situation where our proposed three-stage annealing method is utilized.
In the correction process, we correct the obtained quasi-optimal solution by an Ising machine so that it can satisfy the slot-placement constraint. We traverse the QUBO matrix horizontally and vertically and resolve any item assignment violations and/or slot assignment violations.
In our experiment, we apply our three-stage annealing method to various slot-placement problems and evaluate the solutions using an Ising machine hardware [23], [25], and confirm its effectiveness.

E. CONTRIBUTIONS OF THIS PAPER
The contributions of this paper are summarized as follows: 1) We propose a three-stage annealing method solving the slot-placement problem, which efficiently utilizes an Ising machine and obtains a feasible quasi-optimal solution to the slot-placement problem. 2) As the initial process, we propose a pair-wise exchange method, a random method, and a cluster-growth method. We also describe a correction process, which corrects an obtained quasi-optimal solution so that it satisfies the slot-placement constraint. 3) We applied the proposed three-stage annealing method to many types of slot-placement problems. Compared to the two-stage annealing method [16] not including any initial process, the proposed method reduces the minimum total weighted wiring length by approximately 2.42% on average in the case of using the pair-wise exchange method, approximately 2.45% on average in the case of using the random exchange method, and approximately 0.0898% on average in the case of using the cluster-growth method. The mean total weighted wiring length is reduced by approximately 6.92% on average in the case of using the pair-wise exchange method, approximately 7.08% on average in the case of using the random exchange method, and approximately 2.79% on average in the case of using the cluster-growth method.

F. ORGANIZATION OF THIS PAPER
This paper is organized as follows: Section II summarizes the related works; Section III defines the slot-placement problem and its constraints and objective function; Section IV introduces the Ising model and QUBO model and explains the slot-placement problem mapping to the QUBO model; Section V proposes a three-stage annealing method, where we particularly propose three initial processes, i.e., a pairwise exchange method, a random exchange method, and a cluster-growth method, and also a correction process; Section VI demonstrates the effectiveness of the proposed method using the Ising machine hardware; and Section VII summarizes this paper and gives several concluding remarks.

II. RELATED WORKS
In this section, we firstly summarize typical Ising-machineapplied combinatorial optimization problems. After that, we introduce several approaches to efficiently utilize an Ising machine to solve these combinatorial problems.

A. ISING-MACHINE-APPLIED COMBINATORIAL PROBLEMS
Firstly, we summarize typical combinatorial optimization problems solved by an Ising machine. In [13], a rectangle packing problem is solved using an Ising machine. Given a set of rectangles, the problem arranges them into a small-sized area without overlapping. The problem is efficiently mapped onto the Ising model and solved by the Ising machine.
In [14], [26], a graph partitioning problem is solved using an Ising machine. Given a graph composed of vertices and edges, the problem is to minimize the number of edges connecting the two partitioned vertex sets, such that the number of elements in the two partitioned vertex sets must be equal. The problem is efficiently mapped onto the Ising model and solved by the Ising machine.
In [15], an induced subgraph isomorphism problem is solved using an Ising machine. Given two graphs, a substitution matrix X with binary variables as elements is introduced, which represents the mapping relation between the two graphs, and the problem is mapped to the QUBO model. The problem has two constraints: one is that there exists only one (+1) binary variable in the row direction of X , and the other is that there exists only one (+1) binary variable in the column direction of X . The objective of the problem is to find a substitution matrix X that satisfies XA 2 X T = A 1 , given the adjacency matrices A 1 and A 2 of the two graphs.
In [16], a slot-placement problems is solved using an Ising machine. This problem is mapped to the QUBO model using t-by-m matrix with binary variables as elements, where t shows the number of slots and m shows the number of items. The objective of this problem is to minimize the total weighted wiring length under the slot-placement constraint.
In [27], a number partitioning problem is solved using an Ising machine. Given a set of numbers, the problem is to partition the set so that the sum of the elements of the two partitioned sets is equal. The problem is efficiently mapped onto the Ising model and solved by the Ising machine.
In [28], a traveling salesman problem is solved using an Ising machine. The problem is mapped to the QUBO model using an n-by-n matrix with binary variables as elements, where n refers to the number of cities to be traversed. The problem has two constraints: one is that a person cannot visit multiple cities at the same time, and the other is that the same city cannot be visited multiple times. The objective of the problem is to minimize the total travel cost in a route that visits all cities while satisfying the above two constraints.
In [29], a nurse scheduling problem is solved using an Ising machine. This problem has three constraints: upper and lower limit of the number of breaks, the number of nurses in duty for each shift slot, and upper and lower limit of time interval between two days of duty. The objective of the problem is to find an optimal schedule while satisfying the three constraints.
All the above combinatorial optimization problems except for [16] are directly solved by Ising machines. No preprocesses nor post-processes are applied to them to further improve the obtained solution. Due to the probabilistic nature of Ising machines, the obtained solution cannot always become an optimal solution and hence it cannot always satisfy the constraints given to the original combinatorial optimization problem.

B. HOW TO EFFICIENTLY UTILIZE AN ISING MACHINE
In order to efficiently solve combinatorial optimization problems, several approaches have been proposed to utilize an Ising machine and a classical computer.
In [30], a black-box optimization approach is proposed. Black-box optimization is composed of a machine-learning process by a classical computer and annealing process by an Ising machine. Black-box optimization is effective for particular applications as in material design, but is not applicable for general combinatorial optimization problems since its target is machine-learning-based optimization and it requires a large number of iterations.
In [16], an efficient two-stage annealing method utilizing an Ising machine and a classical computer is proposed to solve combinatorial optimization problems. Firstly, an Ising machine is applied to solve a combinatorial optimization problem. After that, a classical computer improves the obtained solution so that it satisfies the constraints given to the original combinatorial optimization problem. This approach can be effectively applied to the slot-placement problem and considered to be one of the state-of-the-art methods proposed before, since no other methods are proposed combining Ising machines and classical computers to satisfy the constraints as post-processing, as far as we know. VOLUME 9, 2021 However, the obtained solutions by this method much depend on initial spin values in the Ising machine used, and the accuracy of the solutions may not be sufficient.
Based on these discussions, hereinafter, we pick up the slot-placement problem as an example in Sections III and IV and propose an effective annealing method combining an Ising machine and classical computer in Section V.

III. FORMULATION OF THE SLOT-PLACEMENT PROBLEM
In this section, we define and formulate the slot-placement problem as well as its constraints and objective function. Let Let w(c i , c j ) be the number of wires between two items is also zero when i = j. When an item c i is assigned to any slot, the slot is denoted by s(c i ). The weighted wiring length between items c i and c j is denoted by w(c i , c j ) × l(s(c i ), s(c j )). The total weighted wiring length (TWWL), L, is defined by Eqn. (2) is the objective function of the slot-placement problem.
In addition, every item must be assigned to a single slot, which is called item assignment constraint and no more than one item cannot be assigned to a single slot, which is called slot assignment constraint. The slot-placement constraint refers to both of the item assignment constraint and slot assignment constraint.
Hence, we define the slot-placement problem as follows: Definition 1: Given t = h × v slots and m items which are connected by wires, the slot-placement problem is to find the slot assignment of items that minimizes the total weighted wiring length, L, while satisfying the slot-placement constraint. Fig. 1 shows an example of the slot-placement problem of h = v = 3 and m = 5. As in Fig. 1(a), we have five items and 3 × 3 slots, and some items are connected to each other with wires. For example, the items 1 and 2 are connected by three wires, and the items 1 and 4 are connected by two wires. Fig. 1(b) shows the solution to the slot-placement problem. For example, the item 1 and item 4 are assigned to the slot 1 and the slot 2, respectively. The total weighted wiring length is calculated as: (3)

IV. QUBO MODEL MAPPING OF THE SLOT-PLACEMENT PROBLEM
In this section, we firstly introduce the Ising model and QUBO model. After that, we describe QUBO model mapping to the slot-placement problem.

A. ISING MODEL AND QUBO MODEL
The Ising model [1], [2] is a fundamental model in statistical mechanics, which describes the behavior of an entire system of microscopic elements called spins, depending on the interaction between the spins and the external magnetic field acting on each spin. As in Fig. 2, the Ising model is defined on an undirected graph G = (V , E), where V is a set of vertices and E is a set of edges between the vertices. When two vertices i, j ∈ V are connected, (i, j) ∈ E denotes the connected edge between the vertices i and j. Let σ i be the spin placed on the vertex i. σ i takes either (+1) or (−1), where (+1) is the upward spin and (−1) is the downward spin. Let J i,j be the interaction coefficient between the two spins σ i and σ j and h i be the external magnetic field coefficient acting on The smaller the value of H is, the more stable the state of the Ising model is. The state giving the minimum value of H is called the ground state.
The quadratic unconstrained binary optimization (QUBO) model [31] is an equivalent to the Ising model, which uses a binary variable n i that takes the value of 0 or 1 instead of using a spin σ i . The Hamiltonian H of the QUBO model is expressed as follows: where J i,j is the interaction coefficient acting between the two binary variables n i and n j , h i is the external magnetic field coefficient acting on the binary variable n i , and const is a constant. Eqn. (4) can be converted to Eqn. (5) by using Eqn. (6) below: In the rest of this paper, we focus on QUBO models. Note that, QUBO models can be equivalently converted into Ising models as described above.

B. QUBO MODEL MAPPING OF THE SLOT-PLACEMENT PROBLEM
In this subsection, we describe slot-placement problem mapping to the QUBO model based on [15], [16]. Let m be the number of items, t = h × v be the number of slots, c i (1 ≤ i ≤ m) be the i-th item, and s a (1 ≤ a ≤ t) be a-th slot. We define a binary variable x a,i as follows: x a,i = 1 (if c i is assigned to the slot s a ) 0 (otherwise).

1) OBJECTIVE FUNCTION
The objective function in the slot-placement problem is given by Eqn. (2). Using the binary variables x a,i and x b,j defined above, Eqn. (2) can be re-written as follows: Let h a be the minimum value of H A . h a depends on the slot-placement problem instance.

2) ITEM ASSIGNMENT CONSTRAINT
The item assignment constraint shows that every item c i must be assigned to a single slot. Fig. 3(b) shows an example. In Fig. 3(b), the green check indicates that the constraint is satisfied and the red checks indicate that the constraint is violated. The constraint can be formulated as follows: Then we introduce the Hamiltonian H B below: H B takes the minimum value of zero if and only if all items satisfy Eqn. (9).

3) SLOT ASSIGNMENT CONSTRAINT
The slot assignment constraint shows that no more than one item cannot be assigned to every slot. Fig. 3(c) shows an example. In Fig. 3(c), the green checks indicate that the constraint is satisfied and the red check indicates that the constraint is violated. The constraint can be formulated as follows: Then we introduce the Hamiltonian H C below: H C takes the minimum value of t/4 if and only if all slots satisfy Eqn. (11).
The Hamiltonian H takes the minimum value of h a + αt/4, if and only if the ground state is found and the binary variables   giving the ground state shows an optimal solution. How to set up the hyper-parameter α will be discussed in VI-C.
Note that the number of binary variables required to express Eqn. (13) becomes (m × t), where m shows the number of items and t shows the number of slots.

V. THREE-STAGE ANNEALING METHOD UTILIZING AN ISING MACHINE A. THREE-STAGE ANNEALING STRATEGY
In general, the annealing process by an Ising machine tries to minimize the Hamiltonian expressed by Eqn. (4) or Eqn. (5). However, the solution obtained by the Ising machine does not necessarily become a ground-state solution and hence does not necessarily minimize Eqn. (4) or Eqn. (5). Particularly, the constraint terms such as Eqn. (10) and Eqn. (12) in the slot-placement problem cannot be minimized and hence those constraints cannot necessarily be satisfied.
Based on this discussion, we firstly introduce a correction process as in [16] after the annealing process by an Ising machine is over. The correction process is executed by a classical computer and corrects slightly the quasi-optimal solution obtained by the annealing process so that it satisfies the constraints to the original combinatorial optimization problem.
The annealing process globally searches for the solution space represented by the Ising model or QUBO model and reaches the ground state regardless of its initial spin values theoretically if it takes the infinite amount of time [32], [33].
However, in practice, the obtained solution depends on the initial spin values and the optimal solution is not always obtained, since the computational time is finite. There is a trade-off between time and solution accuracy, and finding a quasi-optimal solution in a short time is as demanding as finding an exact optimal solution in practice. Therefore, by appropriately setting the initial spin values or binary variables, we can obtain a better quasi-optimal solution in a short time by using an Ising machine. Furthermore, due to the existence of constraints in combinatorial optimization problems, the solution space can be multimodal and it may take much time to obtain a global optimal solution [34]. By providing initial values close to the optimal solution, we can improve the probability to reach a better quasi-optimal solution since we can start the search which may be closer to one of them.
Based on this discussion, we secondly introduce an initial process before the annealing process by an Ising machine starts. The initial process is executed by a classical computer and gives initial values for spins or binary variables so that the annealing process can obtain a better solution.
In summary, Fig. 4 shows the proposed three-stage annealing method including the initial process and the correction process. (Stage 1) First, the initial process gives an initial solution that satisfies the constraints to the combinatorial optimization problem, which may be close to the optimal solution. (Stage 2) Then, we perform the annealing process on an Ising machine [3], [4], [8], [23], [24] as usual. (Stage 3) Finally, the correction process corrects the quasi-optimal solution obtained by the annealing process so that it satisfies the constraints to the combinatorial optimization problem. Note that, the time complexity of the proposed initial processes (Stage 1) and correction process (Stage 3) can be calculated as in the subsequent subsections, while the annealing time (Stage 2) using an Ising machine much depends on the Ising machine hardware used and its parameter settings (See Section VI-D for the discussion on execution time).
In the rest of this section, we propose an initial process and also describe a correction process targeting the slotplacement problem.

B. INITIAL PROCESS
As an initial process to the slot-placement problem, we propose a pair-wise exchange method (Section V-B1), a random exchange method (Section V-B2), and a cluster-growth method (Section V-B3).

1) PAIR-WISE EXCHANGE METHOD
We propose a pair-wise exchange method as an initial process to the slot-placement problem. Firstly, we generate a random feasible slot-placement solution, where a feasible solution shows the slot placement satisfying the slot-placement constraint. This is done by just assigning every item to a slot randomly without overlapping. Next, we select two slots and exchange the items assigned to them. If no item is assigned to one of the selected slots, we move the item to the empty slot. If no items are assigned to both the selected slots, no move is done. We try all those exchanges and moves and accept the one which most reduces the total weighted wiring length. This process is repeated until the total weighted wiring length no longer improves.
Algorithm 1 shows the pair-wise exchange method. By performing the pair-wise exchange method, a set of initial binary variables to the QUBO model is obtained.
The pair-wise exchange method starts with a random solution that satisfies the slot-placement constraint and greedily exchanges two items or moves one item to an empty slot to reduce the total weighted wiring length. Although it may lead to a locally optimal solution, it can provide a relatively good initial solution. Since the pair-wise exchange method repeats the exchange and move process until it well converges to a local minima, it requires several amount of time. In some cases, it may become longer than the annealing time (See the discussion in Section VI).
The time complexity of the pair-wise exchange method is calculated as follows: The single item exchange or item move requires O(m) time to re-evaluate the total weighted wiring length, where m shows the number of items. Since we try all item exchanges and moves, we require O(m × t 2 ) time in every iteration, where t shows the number of slots. Assume that we repeat this process N pwe times, i.e., after N pwe iterations, no further improvement of the total weighted wiring length is seen. N pwe depends on every slot-placement problem instance. Then the time complexity of the pair-wise exchange method becomes O(N pwe × mt 2 ).

Algorithm 1 Pair-Wise Exchange Method
Generate a random feasible slot-placement solution; while (TWWL is improved) do Try all item exchanges and moves; Accept the one which most reduces TWWL; end while return (Slot-placement solution at this time)

Algorithm 2 Random Exchange Method
Generate a random feasible slot-placement solution; for (Iteration N ) do Try item exchange or item move randomly; if (old TWWL > new TWWL) then Accept the trial above; end if end for return (Slot-placement solution at this time)

2) RANDOM EXCHANGE METHOD
Next, we propose a random exchange method as an initial process to the slot-placement problem. Firstly, we generate a random feasible slot-placement solution. Then, secondly, we try to exchange the items assigned to the slots or move the item to an empty slot for several times. When the total weighted wiring length is reduced, we accept the exchange or move. Algorithm 2 shows the random exchange method. The number N of iterations is given beforehand. By performing the random exchange method, a set of initial binary variables to the QUBO model is obtained.
The random exchange method starts with a random solution that satisfies the slot-placement constraint and randomly exchanges two items or moves an item to an empty slot for a certain number of iterations to reduce the total weighted wiring length. Depending on the number of iterations, a relatively good solution can be obtained as the initial solution. Note that, since the number of iterations can be given in the proposed random exchange method, the CPU time required to the initial process can be adjusted.
The time complexity of the random exchange method is calculated as follows: In the same way as the discussion in the pair-wise exchange method, we require O(m × t 2 ) time in every iteration. Since the iteration is repeated N times, time complexity of the random exchange method becomes O (N × mt 2 ).

3) CLUSTER-GROWTH METHOD
Lastly, we propose a cluster-growth method as an initial process to the slot-placement problem. Firstly, we assign an item i which is connected most to other items onto the center of the lattice slots. Next, among the items that have not yet been assigned to slots, we pick up the item j which is most connected to the item i and assign j to the slot neighboring to i. Similarly, we pick up an unassigned item which has the maximum connections to the items already assigned and VOLUME 9, 2021

Algorithm 3 Cluster-Growth Method
Assign an item which has the most connections to other items to the center of the lattice slots; for (m − 1) do Select an unassigned item j which is most connected to the items already assigned; Assign the item j to the empty slot neighboring to the occupied slots; end for return (Slot-placement solutioin at this time) assign it to an empty slot neighboring to the occupied slots. This process is repeated until all the items are assigned to slots. Algorithm 3 shows the cluster-growth method. By performing the cluster-growth method, a set of initial binary variables to the QUBO model is obtained.
Unlike the pair-wise exchange method and random exchange method, the cluster-growth method does not calculate the total weighted wiring length explicitly during its process. Therefore, the cluster-growth method runs fast even when the scale of the problem becomes larger.
The time complexity of the cluster-growth method is calculated as follows: When we find out the item most connected to other ones, we require O(m 2 ) time, where m shows the number of items. Since we assign an item to a slot in a oneby-one manner, the time complexity of the cluster-growth method becomes O(m 3 ).

C. CORRECTION PROCESS
When an annealing process in Fig. 4 outputs a slot-placement solution which does not satisfy the slot-placement constraint, the correction process corrects the slot-placement solution so that it satisfies the slot-placement constraint. Now the QUBO matrix {x a,i } refers to the matrix arranging the binary variables given by Eqn. (7) as in Fig. 3. According to [15], [16], we simply apply the processes below as the correction process:

1) DELETE REDUNDANT ITEMS FROM QUBO MATRIX COLUMNS AND ROWS
When two or more (+1) binary variables exist in any column in the QUBO matrix as in the first column of Fig. 3(b), we select one of the (+1) binary variables which minimizes the total weighted wiring length and set the other binary variables in this column to zero. We traverse the QUBO matrix from the first column to the last column to resolve the slot direction overlap.
After that, when two or more (+1) binary variables exist in any row in the QUBO matrix as in the second row of Fig. 3(c), we select one of the (+1) binary variables which minimizes the total weighted wiring length and set the other binary variables in this row to zero. We traverse the QUBO matrix from the first row to the last row to resolve the item direction overlap.
The time complexity of this process is calculated as follows: When verifying that a given QUBO matrix satisfies the slot-placement constraint, we require O(mt) time, where m and t show the number of items and the number of slots, respectively. This is because we traverse the QUBO matrix just twice in the column direction and row direction. Assume that we find out the violation of item assignment constraint in the QUBO matrix column. We require O(m 2 + t) time to try the operation that one of the binary variables is fixed to (+1) and the others to zero in this column, to flip binary variables in this column and re-evaluate the total weighted wiring length. Since we have totally m items and we try the operation for these items at worst, we require O(m 3 + mt) time. The similar discussion holds for the violation of slot assignment constraint. Totally, the time complexity of this process becomes O(m 3 + mt).

2) ADD ITEMS TO QUBO MATRIX COLUMNS
After deleting several items as above, there may exist no (+1) binary variables in a column in the QUBO matrix as in the second column of Fig. 3(b). In this case, we flip one of the binary variables in the column which minimizes the total weighted wiring length. We traverse the QUBO matrix from the first column to the last column and add an item if necessary to satisfy the item assignment constraint. The time complexity of this process is calculated as follows: In the worst case, we try to flip almost every binary variable in the QUBO matrix and re-evaluate its total weighted wiring length. Hence, the time complexity of this process becomes O(mt × m) time.
By applying the above processes, the slot-placement solution can satisfy the slot-placement constraint and hence we can always obtain a feasible slot-placement solution.

VI. EXPERIMENTAL EVALUATIONS
We have evaluated the proposed three-stage annealing method applying to various slot-placement problems.

A. COMPARISON METHOD
We have compared the following two methods to demonstrate the effectiveness of the proposed method:

1) OUR PROPOSED METHOD
The first one is the proposed method. In the proposed method, we use the classical computer environment of CentOS Linux 7.7.1908 and Intel Xeon Gold 6148 CPU processor for the initial process and the correction process. We run the pairwise exchange (PWE) method, the random exchange (RE) method, and the cluster-growth (CG) method as an initial process. The number of iterations in the random exchange method is set to 10,000. ''Ours (PWE)'' refers to the proposed method with the pair-wise exchange method. ''Ours (RE)'' refers to the proposed method with the random exchange method. ''Ours (CG)'' refers to the proposed method with the cluster-growth method.
As an Ising machine, we use one of the semiconductorbased Ising machines, called Digital Annealer (DA) Unit 2 [23], [25], for the annealing process. 2 The number of iterations in the annealing process in DA is set to 1,000,000 (default). In DA, we use the parallel tempering mode (PT mode), where 128 quasi-optimal solutions are obtained simultaneously. We apply the correction process to these 128 solutions if necessary and obtain the 128 feasible slot-placement solutions. Fig. 5 summarizes the experimental flow of the proposed method.

2) BASELINE METHOD [16] (TWO-STAGE ANNEALING METHOD)
For comparison, we run the existing state-of-the-art method [16] which is composed of the annealing process and the correction process. Before starting the annealing process, we initialize all the binary variables to zero. In this method, we use the same classical computer environment as Ours. We also use a Digital Annealer Unit 2 for the annealing 2 The Ising machine that we use is implemented on ASIC and realized by hardware [23], [25]. It can deal with a maximum of 8,192 binary variables and hence we set up the problem instances as described in Table 1 where we can use at most 7,500 binary variables. process. The number of iterations in the annealing process is set to 1,000,000. In the same way, we use PT mode, where 128 quasi-optimal solutions are obtained simultaneously. We apply the correction process to these 128 solutions if necessary and obtain the 128 feasible slot-placement solutions.
Note that, if we do not apply the correction process, the obtained slot-placement solutions cannot necessarily satisfy the slot-placement constraint. In fact, if we do not apply the correction process in our experiments, a maximum of 98% solutions violate the slot-placement constraint (See Table 6 in detail). Hence, we compare our proposed method to the two-stage annealing method above.

B. PROBLEM SETTING
We prepare 4 × 4 slots to 10 × 10 slots. For every p × p slot instance, we prepare three item instances, m = (p × p) × (1/2), m = (p × p) × (3/4), and m = (p × p). Note that, we did not prepare an item instance m = 100 for 10 × 10 slot instance, because the number of variables exceeds the upper limit of the Ising machine that we use for the annealing process. The items are randomly connected by wires such that w(c i , c j ) ∈ [0, 10] for any two items c i and c j .
For every slot-placement problem instance, we randomly generate 10 different problems. For each of 10 different problems, we obtained 128 feasible slot-placement solutions by our methods or the baseline method [16]. As in Fig. 5, we picked up the best solution (Min. value) which minimizes the total weighted wiring length over these 128 solutions. We also calculated the mean value (Avg. value) of the total weighted wiring length over 128 solutions. Then, we averaged these Min. values and Avg. values over 10 different problems and summarized them.

C. SEARCH FOR THE OPTIMAL HYPER-PARAMETER
Before we evaluate our proposed method, we have searched for the optimal hyper-parameter α value described in VOLUME 9, 2021 Eqn. (13). In this search, we did not use any initial process and set zero to all the binary values at first. We used the correction process to obtain a feasible solution. We prepared 10 different problems for every slot-placement problem instance and searched for the hyper-parameter value in the range of α ∈ [10, 600] for 4 × 4, 5 × 5, and 6 × 6 slot instances and in the range of α ∈ [100, 700] for 7 × 7 slot instance. We also used PT mode as in Section VI-A and we obtained 128 feasible solutions at a time for every problem. We calculated and summarized averaged Avg. value and Min. value as described in Section VI-B for every slot-placement problem instance. Fig. 6 shows the results. In most of the cases, Avg. value and Min. value are minimized at the particular α value, which gives the optimal hyper-parameter. This is because of the following reason: if the α value is too small, the constraint terms of H B and H C of Eqn. (13) cannot be minimized and hence the non-feasible solutions are obtained. Then the correction process corrects those solutions and their total weighted wiring lengths are not globally minimized. If the α value is too large, the constraint terms of H B and H C of Eqn. (13) can be minimized and the obtained solutions can satisfy the slot-placement constraint. But the first term H A of Eqn. (13) is not fully minimized and hence their total weighted wiring lengths are not minimized.
According to the results of Fig. 6, Fig. 7 plots the optimal hyper-parameter value versus the number of variables. As in Fig. 7, the optimal hyper-parameter value is  proportional to the number of variables used and it is expressed by: (14) where N v shows the number of variables and α opt shows the optimal hyper-parameter value. Based on this discussion, we summarize the hyper-parameter values in Table 1 for every slot-placement problem instance. Table 2 summarizes the averaged minimum and mean values (Min. values and Avg. values) of the total weighted wiring length over 10 different problems for 4 × 4 to 10 × 10 slot instances by using our methods and the baseline method [16]. The numbers in parentheses indicate the improvement of our methods over [16]. From Table 2, the minimum values are   to [16]. The mean values are improved by approximately 6.92% on average in the case of Ours (PWE), 7.08% on average in the case of Ours (RE), and 2.79% on average in the case of Ours (CG). The results indicate that Ours is superior to [16] in almost all the cases. The total weighted wiring length obtained by Ours (PWE) and Ours (RE) is almost the same. In any cases except for small examples, the total weighted wiring length is reduced by introducing the initial process, which definitely indicates that the three-stage annealing process is effective to the slot-placement problem. Table 3 shows the total weighted wiring length after the initial process and the averaged minimum value (Min. value) after the correction process. The numbers in parentheses indicate the improvement of Min. value over ''after initial process.'' The averaged minimum total weighted wiring length obtained after the correction process are improved approximately up to 4.4% in the case of Ours (PWE), 2.9% in the case of Ours (RE), and 12% in the case of Ours (CG) compared to those after initial process. These results demonstrate that the annealing process successfully improves the total weighted wiring length, even after the initial solution is given by the initial process. Table 4 summarizes the averaged execution time for each stage and Table 5 summarizes the total averaged execution time for multiple stages. The initial process itself requires several amounts of time, but by introducing the initial process, the time required for the correction process is much reduced. Particularly, Ours (CG) requires less time than Ours (PWE) and Ours (RE), and the overall execution time of Ours (RE) and Ours (CG) are comparable to [16]. Table 6 shows the average number of the correction process runs. As described, we obtain 128 solutions simultaneously and hence the correction process can be applied to each of them, i.e., the maximum number of the correction process runs becomes 128. If all the solutions obtained by the annealing process are feasible ones, then the number of the correction process runs becomes zero. Table 6 clearly indicates that, by introducing the initial process, the annealing process tends to output the feasible solutions and hence the average number of the correction process runs is much decreased.

D. EVALUATION RESULTS
Overall, the three-stage annealing process is effective enough to the slot-placement problem in terms of reducing the objective function, satisfying the slot-placement constraint, and computational time.

VII. CONCLUSION
In this paper, we have proposed a three-stage annealing method solving the slot-placement problem. The experimental results demonstrate that the proposed method reduces the minimum total weighted wiring length by 0.0898%-2.45% on average depending on the initial process methods used, compared to the existing method. The mean total weighted wiring length is reduced by 2.79%-7.08% on average depending on the initial process methods used. Totally, the proposed method is effective to the slot-placement problem. In the future, we will apply the proposed approach to other combinatorial problems and further confirm its effectiveness.