Quantum Optimization with Lagrangian Decomposition for Multiple-process Scheduling in Steel Manufacturing

Abstract

Steel manufacturing involves multiple processes, and optimizing the production schedule is essential for improving the quality, cost, and delivery of products. However, the corresponding optimization problems are usually NP-hard and generally intractable even for modern high-performance computers. Recent advances in quantum computers, such as those developed by D-Wave Systems, have opened new possibilities for handling this type of optimization problem. Nevertheless, a major obstacle to the application of currently available quantum computers is the relatively small number of quantum bits, which limits the feasible problem size. To overcome this obstacle, we propose an algorithm for solving optimization problems based on Lagrangian decomposition and coordination. D-Wave quantum computers can explore diverse solutions simultaneously, and we leverage this unique characteristic to obtain the upper bound for Lagrangian decomposition and coordination. We apply the proposed algorithm to a simplified problem of two-process production scheduling. By decomposing the problem into two stages for the processes, the number of steel products that can be handled increases from five to eight (60% increase). The optimal solutions are reached even for the extended case of the eight products. The proposed algorithm is a promising technique for enabling the application of quantum computers to real-world problems by thoroughly exploiting quantum bits.

1. Introduction

Steel is a fundamental material supporting the development of our society. In steel manufacturing, raw materials, mainly iron ore and coal, result in various types of steel products through a series of processes, including smelting, refining, continuous casting, rolling, annealing, and surface treatments.

For steel companies that produce various types of products, the optimization of production planning and scheduling for multiple processes is challenging. Over the multiple production processes of steel, the unit volume per batch gradually decreases throughout processing. For example, the typical volume upstream of steelmaking is approximately 300 tons, whereas the volume downstream when reaching surface treatment is approximately 10 tons. Therefore, a production unit in an upstream process includes different products or semi-products across downstream processes. Between all the process layers, the grouping of different sequence orders or semi-products should be jointly optimized. Such optimization problems must be solved under several constraints. For instance, two products or semi-products requiring very different processes or different specifications, such as chemical components, widths, and thicknesses, cannot be grouped. On the other hand, grouping products obtained from similar processing conditions is favorable for quality, productivity, and cost. Nevertheless, excessive grouping may lead to stock accumulation and delay because grouping implies waiting for other products that undergo similar processes before proceeding.

The tradeoff between grouping and due times should be optimized, but even a single process is difficult to optimize. When considering multiple processes, both the scale of the problem and complexity increase, as the grouping compatibility between any two products differs depending on the process. Such complexity hinders the development of an efficient algorithm for real production scheduling. A solution has been proposed in the previous work,¹⁾ where scheduling is solved using a two-step algorithm. First, the grouping size (lot size) is dynamically decided using Lagrangian decomposition and coordination (LDC). Then, a feasible schedule is generated heuristically considering the applicable constraints.

Recently, quantum computers have emerged to solve challenging combinatorial optimization problems. Their quantum superposition capability allows to search in parallel over numerous possible solution paths, likely allowing to obtain a better solution in a much shorter time than classical computers.

In particular, the quantum computers developed by D-Wave Systems are promising for quantum optimization because they offer thousands of quantum bits (qubits). A D-Wave quantum computer can solve quadratic unconstrained binary optimization (QUBO) problems, and several combinatorial optimization problems have already been addressed using these machines.^2,3,4,5,6)

Despite the promise of quantum computers, their application in industry is limited by the current inability to handle large-scale problems owing to the scarcity of qubits, lacking connections between qubits, and other factors. The connectivity between qubits is problematic, especially when a problem describes high correlations between a qubit and several other bits. To represent this type of correlation, many redundant qubits should be assigned to encode a single bit that represents the numerous connections. Scale limitations can be alleviated by adopting a hybrid approach that combines quantum annealers with classical computers for processing.^7,8) Therefore, we apply LDC to a production scheduling problem comprising two processes for a steel factory. Specifically, the optimization problem over the two processes is decomposed into two optimization problems for the processes, and the constraints between the problems are considered using LDC. We verify that LDC can determine the optimal solutions in all the experiments. For certain problems, LDC can reduce the required number qubits compared with conventional algorithms.

The remainder of this paper is organized as follows. In Sec. II, we define the problem of job scheduling for two processes and present the QUBO formulation required to use the D-Wave quantum computers. Section III details the proposed algorithm with focus on the novel LDC method. Experimental results are reported in Sec. IV. Finally, we discuss the study and findings in Sec. V and draw conclusions in Sec. VI.

2. Optimization Problem: Job Scheduling for Two Processes

We deal with a simplified job shop scheduling problem for two processes that represent the features of a variety of processes in steel manufacturing. The unit of a job corresponds to a steel product (e.g., coil, plate). All the products are sequentially processed in the first and second processes. A processing group and due time are assigned to each job. Optimization attempts to find the schedule that provides a finishing time of products close to their due times while minimizing the number of changes of processing groups during a production series to reduce costs.

2.1. Problem Definitions

We consider the following definitions in the optimization problem:

1) N products are finished through sequential manufacturing processes 1 and 2.

2) The time required to process one product in a process is constant for all the products and both processes. The time is thus normalized with this required time. Then, the problem is considered as distributing N products into several timeslots dedicated to processes 1 and 2, as shown in Fig. 1.

3) The due time from process 1 to process 2 is neglected.

4) If a product undergoes process 1 at timeslot s, it undergoes process 2 at timeslot s or later.

5) Considering manufacturing efficiency, N products are distributed into N consecutive timeslots for both processes 1 and 2.

6) Without loss of generality, process 1 handles products from the first timeslot, whereas process 2 may start later. Therefore, we prepare (N + e) timeslots for the process 2, where e depends on the due time, as shown in Fig. 1.

7) For every product, the due time and production groups of processes 1 and 2 are defined. The cost is counted when two adjacent products belong to different production groups. In real production, this cost is related to changes in the process conditions.

8) The formulated problem is aimed to find the optimal scheduling considering the tradeoff between the due time and group mismatch.

Fig. 1.

Simplified job scheduling for two processes. The valid timeslots are highlighted in gray.

2.2. Problem Formulation

Multiple-process scheduling can be defined as an integer programming problem (see the Appendix). Unfortunately, the existing quantum annealers only handle binary variables. In addition, we cannot take large time slots since the number of the available variables on the current D-Wave machines are limited. Thus, we formulate the scheduling problem to efficiently use the currently available quantum annealers. We detail the corresponding model below.

2.2.1. Parameters and Notation

The notation and parameters adopted in this study are detailed as follows:

p: Index of process, 1 ≤ p ≤ 2

i,j: Index of products, 1 ≤ i,j ≤ N

s: Index of timeslots, 1 ≤ s ≤ N

t_p,i: Time when product i undergoes process p

Δ: Gap between start times of processes 1 and 2

d_i: Due time of product i

L_p,i,j: Cost of group mismatch in process p; L_p,i,j = 0 if products i and j belong to the same group, and L_p,i,j = 1 if the products belong to different groups

w_g, w_e, w_l: Coefficients of costs

ρ, ρ′: Coefficients of penalty terms

2.2.2. Decision Variables

  

$$x_{p,i,s}=\{\begin{array}{ll}1,\hfill & \mathrm{i}\mathrm{f}\mathrm{\hspace{0.17em} } \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{\hspace{0.17em} } i\mathrm{\hspace{0.17em} } \mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{o}\mathrm{e}\mathrm{s}\mathrm{\hspace{0.17em} } \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{\hspace{0.17em} } p\mathrm{\hspace{0.17em} } \mathrm{a}\mathrm{t}\mathrm{\hspace{0.17em} } \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{t}\mathrm{\hspace{0.17em} } s\hfill \\ 0,\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}$$

2.2.3. Constraints

We define the following two constraints:   

$$\sum _{s=1}^N{x}_{p,i,s}=1{\hspace{1em}}^{\forall }i,{\mathrm{\hspace{0.17em} }}^{\forall }p,$$(1)

  

$$\sum _{i=1}^N{x}_{p,i,s}=1{\hspace{1em}}^{\forall }s,{\mathrm{\hspace{0.17em} }}^{\forall }p.$$(2)

Constraint (1) ensures that all products undergo all the processes once. Constraint (2) guarantees that only one product can be assigned to a timeslot. In addition, we introduce the following constraint for time t_p,i:   

$$t_{1,i}\le t_{2,i}{\hspace{1em}}^{\forall }i.$$(3)

Constraint (3) makes process 2 start only after a product has undergone process 1. Let us establish the definitions of t_p,i. When all the binary variables satisfy constraint (2), t_1,i and t_2,i can be written as   

$$t_{1,i}=\sum _{s=1}^{N}s{x}_{1,i,s}{\hspace{1em}}^{\forall }i,$$(4a)

  

$$t_{2,i}=\sum _{s=1}^N\left(s+\mathrm{\Delta }\right)x_{2,i,s}{\hspace{1em}}^{\forall }i.$$(4b)

Figure 2 illustrates our formulation and an example of the solution for N = 3 and Δ = 1. The gray areas correspond to x_p,i,s = 1. In the original definition of the multiple scheduling problem, process 2 has N + e timeslots. In our formulation, the timeslots in p = 2 range from 1 to N. Thus, our formulation reduces the number of variables. However, gap Δ between the start times of processes 1 and 2 should be considered in the scheduling problem, as explained in Sec. IV.

Fig. 2.

Diagram of proposed formulation and example of solution for N = 3 and Δ = 1.

2.2.4. Costs

We define two types of costs: due time and group mismatch. The costs of missing due times are defined as follows:   

$$Cos{t}_e\left({\mathit{x}}_2\right)=\sum _{i=1}^N\sum _{s=1}^{d_i-\mathrm{\Delta }-1}\left(d_i-\left(s+\mathrm{\Delta }\right)\right)x_{2,i,s},$$(5)

  

$$Cos{t}_l\left({\mathit{x}}_2\right)=\sum _{i=1}^N\sum _{s=d_i-\mathrm{\Delta }+1}^N\left(\left(s+\mathrm{\Delta }\right)-d_i\right)x_{2,i,s},$$(6)

where x_p denotes the binary variables for the p-th process, and Cost_e(x) (Cost_l(x)) corresponds to the cost when the products are processed earlier (later) than their due times. The cost of the group mismatch in the p-th process is defined as   

$$Cos{t}_g\left(\mathit{x}\right)=\sum _{p=1}^{2}Cos{t}_{g,p}\left({\mathit{x}}_p\right),$$(7)

where   

$$Cos{t}_{g,p}\left({\mathit{x}}_p\right)=\sum _{i,j=1}^N\sum _{s=1}^{N-1}{L}_{p,i,j}{x}_{p,i,s}{x}_{p,j,s+1}.$$

Thus, we define the costs for processes 1 and 2 as   

$$Cos{t}_1\left({\mathit{x}}_1\right)=w_{g}Cos{t}_{g,1}\left({\mathit{x}}_1\right),$$

  

$$Cos{t}_2\left({\mathit{x}}_2\right)=w_{g}Cos{t}_{g,2}\left({\mathit{x}}_2\right)+w_{e}Cos{t}_e\left({\mathit{x}}_2\right)+w_{l}Cos{t}_l\left({\mathit{x}}_2\right).$$

The total cost is given by   

$$Cost\left(\mathit{x}\right)=Cos{t}_1\left({\mathit{x}}_1\right)+Cos{t}_2\left({\mathit{x}}_2\right).$$(8)

2.2.5. QUBO Formulation

To use the current D-Wave machines, we express constraints (1), (2), and (3) in a quadratic form over x. Equality constraints (1) and (2) can be represented using penalty terms as⁹⁾   

$$\sum _{s=1}^N{x}_{p,i,s}=1\mathrm{\hspace{0.17em} }{ }^{\forall }i,{\mathrm{\hspace{0.17em} }}^{\forall }p\leftrightarrow \sum _{p=1}^2\sum _{i=1}^N{\left(\sum _{s=1}^N{x}_{p,i,s}-1\right)}^2,$$(9)

  

$$\sum _{i=1}^N{x}_{p,i,s}=1\mathrm{\hspace{0.17em} }{ }^{\forall }s,{\mathrm{\hspace{0.17em} }}^{\forall }p\leftrightarrow {\sum _{p=1}^2\sum _{s=1}^N\left(\sum _{i=1}^N{x}_{p,i,s}-1\right)}^2.$$(10)

We express the penalty terms for the p-th process as follows:   

$$Penalt{y}_p\left({\mathit{x}}_p\right)=\left\{{\sum _{t=1}^N\left(\sum _{i=1}^N{x}_{p,i,t}-1\right)}^2+{\sum _{i=1}^N\left(\sum _{t=1}^N{x}_{p,i,t}-1\right)}^2\right\}.$$(11)

To deal with inequality constraints, various techniques are required. A simple approach is to introduce the following alternative expression to constraint (3):   

$$\sum _{i=1}^N\sum _{s,s^{\prime }=1}^N\mathrm{\Theta }\left(s-\left(s^{\prime }+\mathrm{\Delta }\right)\right)x_{1,i,s}{x}_{2,i,s^{\prime }},$$(12)

where Θ(z) is the Heaviside step function, with Θ(z) = 1 if z > 1 and Θ(z) = 0 otherwise. From the definition of Eqs. (4a) and (4b), [s−(s′+Δ)] corresponds to (t_1,i−t_2,i). Thus, constraint (11) can be interpreted as the penalty when process 2 starts earlier than process 1.

Substituting Eqs. (9), (10), and (11) into Eq. (8), we obtain the following QUBO model:   

$$\begin{array}{l}E\left(\mathit{x}\right)=Cost\left(\mathit{x}\right)+\rho \sum _{p=1}^{2}Penalt{y}_p\left({\mathit{x}}_p\right)\\ \hspace{1em}\hspace{1em}+{\rho }^{\prime }\sum _{i=1}^N\sum _{s,s^{\prime }=1}^N\mathrm{\Theta }\left(s-\left(s^{\prime }+\mathrm{\Delta }\right)\right)x_{1,i,s}{x}_{2,i,s^{\prime }}.\end{array}$$(13)

3. Methods

3.1. Quantum Annealing

Quantum annealing (QA) is a metaheuristic approach to solve optimization problems.¹⁰⁾ A QA system consists of target and driver parts. The target part corresponds to the QUBO model to be optimized, while the driver part introduces quantum fluctuations into the system and creates a superposition state that contains all the patterns of solutions. To start QA, we set the driver part as dominant and explore diverse solutions. Then, we gradually reduce the influence of the driver part, and the system evolves to a solution. If QA is operated slowly enough and the system suffers no external interferences, we can reach the optimal solution. In fact, QA has outperformed classical simulated annealing under ideal conditions.^11,12) Unfortunately, the QA implementation in an actual quantum device is affected by random noise, and it is difficult to conduct ideal QA. Therefore, current quantum annealers establish samplers that stochastically generate diverse approximated solutions.¹³⁾

The D-Wave 2000Q quantum annealer has approximately 2000 physical qubits implemented in a chimera graph. D-Wave 2000Q can solve a QUBO model defined as   

$$E_{\text{QUBO}}={\mathit{x}}^{T}Q\mathit{x},$$

where x is a vector of binary variables and Q is a real-valued matrix. However, the connectivity between qubits on the chimera graph is limited. Consequently, D-Wave 2000Q cannot represent the QUBO problem directly, and the embedding and unembedding become necessary.^14,15) The embedding technique maps the QUBO problem onto the D-Wave hardware. After QA, we use the unembedding to retransform the annealing results into the original QUBO representation. Thus, the embedding and unembedding allow to solve various QUBO problems on D-Wave 2000Q. As embedding requires many physical bits, the problem size that can be computed reduces dramatically. For example, when the QUBO problem is given by a fully connected graph, D-Wave 2000Q can only compute up to 64 variables. Furthermore, a previous study reported that the success probability of the embedding algorithm reduces dramatically as the problem size increases.¹⁵⁾

3.2. Lagrangian Decomposition and Coordination

LDC is a relaxation and decomposition method for solving constrained optimization problems.^16,17) Specifically, LDC generates an approximated relaxation problem by incorporating some of constraints into Cost(x). In addition, the accuracy of LDC is guaranteed by the duality gap between its upper and lower bounds. LDC has been applied in various production scheduling problems.^1,18,19,20) We use LDC for the scheduling problem addressed in this study.

We apply LDC to inequality constraint (3). With Lagrange multiplier λ, the QUBO formulation for the relaxation problem is given by   

$$E_{\text{LDC}}\left(\mathit{x}\right)=Cost\left(\mathit{x}\right)+\rho \sum _{p=1}^{2}Penalt{y}_p\left({\mathit{x}}_p\right)+\sum _{i=1}^N{\lambda }_i\left(t_{1,i}-t_{2,i}\right).$$(14)

Considering the definitions in Eqs. (4a), (4b), and (8), we can separate E_LDC(x) into the following two QUBO models:   

$$E_1\left({\mathit{x}}_1\right)=Cos{t}_1\left({\mathit{x}}_1\right)+\rho Penalt{y}_1\left({\mathit{x}}_1\right)+\sum _{i=1}^N{\lambda }_i{t}_{1,i},$$(15a)

  

$$E_2\left({\mathit{x}}_2\right)=Cos{t}_2\left({\mathit{x}}_2\right)+\rho Penalt{y}_2\left({\mathit{x}}_2\right)-\sum _{i=1}^N{\lambda }_i{t}_{2,i}.$$(15b)

As t_p,i is a linear over x_p, E₁(x₁) and E₂(x₂) are given as QUBO formulations. Thus, LDC allows to solve E₁(x₁) and E₂(x₂) separately as follows:   

$${\mathit{x}}_p^{\text{*}}=\underset{{\mathit{x}}_p}{\text{argmin}}{E}_p\left({\mathit{x}}_p\right).$$(16)

To obtain a feasible solution, the appropriate Lagrangian multipliers are searched. To this end, we use the subgradient method as follows:   

$$\begin{array}{c}{\lambda }_i^{new}={\lambda }_i+\alpha max\left(t_{1,i}^{*}-t_{2,i}^{*},0\right),\end{array}$$(17)

where parameter α > 0 determines the step size. In addition, we can obtain $t_{1,i}^{*}$ ($t_{2,i}^{*}$) by substituting ${\mathit{x}}_1^{\text{*}}$ (${\mathit{x}}_2^{\text{*}}$) into Eq. (4a) [(4b)]. As the convergence criterion, we use the following duality gap:   

$$duality\mathrm{\hspace{0.17em} } gap=\frac{\left|Cost\left({\mathit{x}}^{\text{*}}\right)-Cost\left({\mathit{x}}^{\text{feas}}\right)\right|}{\left|Cost\left({\mathit{x}}^{\text{*}}\right)\right|},$$

where Cost(x^*) is the lower bound, x^* represents {${\mathit{x}}_1^{\text{*}}$, ${\mathit{x}}_2^{\text{*}}$} obtained from Eq. (16), and upper bound Cost(x^feas) can be defined by a feasible solution x^feas, which is often computed using a heuristic algorithm.

LDC provides various advantages. In the original QUBO model, E(x), the number of variables is 2N², while LDC requires only N² variables because the total QUBO model is divided into E₁(x₁) and E₂(x₂). Let us consider a more general case with N_p processes. Even in this case, the number of variables in the LDC method is N², whereas the original QUBO model would require N_pN² variables. Thus, LDC can considerably reduce the number of variables for solving optimization problems. In addition, LDC can be run in parallel. By solving E₁(x₁) and E₂(x₂) in parallel, we can efficiently search for a feasible solution.

3.3. Proposed Algorithm

We propose a hybrid LDC algorithm that combines LDC and QA. The main process in our algorithm uses D-Wave machines to determine the lower and upper bounds. As mentioned above, D-Wave machines are samplers that generate diverse solutions stochastically. Therefore, we input QUBO model E_p(x_p) into the D-Wave machines to then obtain samples {${\mathit{x}}_p^{\nu }$}, where ν is the index of each sample. From the samples obtained by QA, we define the lowest-cost solution for the p-th process as   

$${\mathit{x}}_p^{\text{cost}}=\underset{x_p\in \left\{x_p^{\nu }\right\}}{\text{argmin}}{E}_p\left({\mathit{x}}_p\right).$$(18)

The lowest cost, Cost₁(${\mathit{x}}_1^{\text{cost}}$) + Cost₂(${\mathit{x}}_2^{\text{cost}}$), corresponds to the lower bound. In LDC, the combination of ${\mathit{x}}_1^{\text{cost}}$ and ${\mathit{x}}_2^{\text{cost}}$ may not be feasible initially. However, the combination of samples in {${\mathit{x}}_1^{\nu }$} and {${\mathit{x}}_2^{\nu }$} can be feasible. We define feasible solution x^feas as   

$${\mathit{x}}^{\text{feas}}=\underset{x_1\in \left\{x_1^{\nu }\right\},\mathrm{\hspace{0.17em} } x_2\in \left\{x_2^{\nu }\right\}}{\text{argmin}}{E}_1\left({\mathit{x}}_1\right)+E_2\left({\mathit{x}}_2\right)s.t.\hspace{1em}{t}_{1,i}\le t_{2,i}{\hspace{1em}}^{\forall }i.$$(19)

Thus, we can obtain the lower and upper bounds from sampling on the D-Wave machines. The proposed hybrid LDC algorithm is detailed as follows:

1) Initialize Lagrangian multipliers as λ_i = 0 ^∀i

2) Perform embedding

3) Initialize counter t = 1

4) Compute QUBO matrices

5) Input the QUBO matrices and an embedding result into D-Wave 2000Q to obtain samples {${\mathit{x}}_1^{\nu }$} and {${\mathit{x}}_2^{\nu }$}

6) Compute the lowest-cost solutions as ${\mathit{x}}_1^{\text{cost}}$ and ${\mathit{x}}_2^{\text{cost}}$

7) Compute feasible solution x^feas and the upper bound

8) Check convergence: when one of the following criteria is satisfied, terminate the iterative process:

 i) t > 20

 ii) duality gap < 0.01

 iii) upper bound does not improve in 10 consecutive steps.

9) Update the multipliers as λ_i = λ_i + αmax($t_{1,i}^{*}$ − $t_{2,i}^{*}$,0) ^∀i

10) Update the counter as t = t + 1

11) Iterate steps 4–10 until convergence

To obtain accurate results, lower-cost solutions should be sampled. We use two techniques:¹⁴⁾ 1) energy minimization method for unembedding, which proceeds while computing a local QUBO model; 2) greedy steepest descent (GSD) to modify the samples obtained from QA on the D-Wave machines while minimizing a local QUBO model.

4. Experiments and Results

4.1. Experimental Setup

We tested the performance of the proposed hybrid algorithm on artificial setups. Table 1 lists the setup of the five- and eight-product problems, whose product groups and due times were generated randomly. For comparison, we evaluated the performance of solving the original QUBO model, E(x), using D-Wave 2000Q. We also compared the results obtained by solving the integer linear programming formulation (see the Appendix) with the CPLEX software.²¹⁾ We set the annealing time to 20 μs and generated 1000 samples on D-Wave 2000Q. The coefficients of the penalty terms were given by ρ, ρ′ = 5 × max(w_g, w_e, w_l), and those of the costs were set as w_e = 3, w_d = 1, and w_g = {4, 10, 100}. Step size parameter α in the subgradient method was set to 0.01.

Table 1. Experimental setup for five- and eight-product problem.

Product	1	2	3	4	5
Due time	5	1	4	2	3
Group 1	2	3	4	2	4
Group 2	4	6	6	8	1

Product	1	2	3	4	5	6	7	8
Due time	7	1	5	8	6	2	3	4
Group 1	3	2	4	3	4	1	3	3
Group 2	7	2	7	8	3	2	3	7

4.2. Dependence on Gap Time Δ between Processes

We use parameter Δ to reduce the number of variables by representing the gap between the start times of processes 1 and 2, which must be determined before scheduling. Table 2 details the Δ-dependence of each cost in the five-product problem in terms of Cost^*, $Cos{t}_e^{*}$, $Cos{t}_l^{*}$, and $Cos{t}_g^{*}$, which indicate the costs in the best feasible solution obtained by the hybrid LDC algorithm. The values highlighted in boldface in rows Cost^* are consistent with the optimal solutions obtained by CPLEX. Table 2 shows that Cost^* is minimized at Δ = 1 for every w_g. This is because the due times start from t = 1 in our problem. When Δ shifts from 1, the values of Cost_e + Cost_l increase. Thus, the hybrid LDC algorithm can achieve the optimal solution for the appropriate value of Δ.

Table 2. Δ-dependence of costs.

	Δ = 0	Δ = 1	Δ =2	Δ = 3	Δ = 4	Δ = 5
wg = 4
Cost*	29	24	35	50	65	80
$Cos{t}_e^{*}$	5	0	0	0	0	0
$Cos{t}_l^{*}$	0	0	5	10	15	20
$Cos{t}_g^{*}$	6	6	5	5	5	5
wg = 10
Cost*	59	58	65	80	95	110
$Cos{t}_e^{*}$	6	2	0	0	0	0
$Cos{t}_l^{*}$	1	2	5	10	15	20
$Cos{t}_g^{*}$	5	5	5	5	5	5
wg = 100
Cost*	509	508	515	530	545	566
$Cos{t}_e^{*}$	6	2	0	0	0	0
$Cos{t}_l^{*}$	1	2	5	5	15	20
$Cos{t}_g^{*}$	5	5	5	10	5	5

4.3. Accuracy and Computation Time

To investigate the typical performance of the hybrid LDC, we performed five independent trials. To evaluate the algorithm performance, we used the average relative error and computation times. The relative error is defined as   

$$\text{error}=\frac{\left|Cost\left({\mathit{x}}^{\text{opt}}\right)-Cost\left({\mathit{x}}^{\text{feas}}\right)\right|}{\left|Cost\left({\mathit{x}}^{\text{opt}}\right)\right|},$$

where x^opt is the optimal solution obtained by CPLEX. In addition, we measured three computation times: t_cloud, t_unemb, and t_LDC. Total cloud time t_cloud involves internet latency, QA, and other processes of D-Wave 2000Q. Time t_unemb is related to unembedding in the LDC iterative process. The total computation time in the hybrid LDC algorithm is denoted by t_LDC and involves t_cloud, t_unemb, and the computation times for GSD and other processes.

Tables 3 and 4 list the experimental results of the five- and eight-product problems for Δ = 1, respectively. For comparison, we show the average relative error, t_cloud, and t_unemb when solving E(x) using D-Wave 2000Q. Table 3 shows that the average errors obtained by the hybrid LDC algorithm are zero for every w_g. Hence, the hybrid LDC achieves optimal solutions in all the five trials. D-Wave 2000Q also finds the optimal solution in all experiments except for w_g = 4. In addition, t_cloud, and t_unemb of D-Wave 2000Q are shorter than those of the proposed hybrid LDC algorithm, which requires multiple sampling operations on D-Wave 2000Q to update the multipliers. Table 4 shows that the hybrid LDC algorithm achieves the optimal solution in all trials. Unfortunately, the direct implementation in D-Wave 2000Q cannot solve the eight-product problem because the embedding fails. In the eight-product problem, the number of variables is 128 because the QUBO model, E(x), requires 2N² variables. As reported in the previous work,¹⁵⁾ the success probability of embedding reduces as the problem size increases. Consequently, in the eight-product problem, the embedding fails to map E(x) onto the hardware graph of D-Wave 2000Q. In the hybrid LDC algorithm, the QUBO model is separated into E₁(x₁) and E₂(x₂). Therefore, the number of variables is 64 even for the eight-product problem because E₁(x₁) and E₂(x₂) require N² variables. Owing to this reduction, the hybrid LDC algorithm can solve larger problems than using the direct implementation in D-Wave 2000Q.

Table 3. Average relative error and computation times in five-product problem.

	wg = 4				wg = 10				wg = 100
	Error	tcloud	tunemb	tLDC	Error	tcloud	tunemb	tLDC	Error	tcloud	tunemb	tLDC
Hybrid LDC	0.0	17.05	4.38	25.45	0.0	18.32	4.17	26.45	0.0	16.89	4.17	25.07
D-Wave 2000Q	0.07	1.74	0.49	–	0.0	2.32	0.49	–	0.0	2.24	0.49	–

Table 4. Average relative error and computation times in eight-product problem.

	wg = 4				wg = 10				wg = 100
	Error	tcloud	tunemb	tLDC	Error	tcloud	tunemb	tLDC	Error	tcloud	tunemb	tLDC
Hybrid LDC	0.0	25.60	13.71	50.90	0.0	47.67	27.07	97.65	0.0	24.79	13.97	50.23
D-Wave 2000Q	–	–	–	–	–	–	–	–	–	–	–	–

5. Discussion

The proposed hybrid LDC algorithm can successfully optimize scheduling problems for up to eight products. However, much larger scheduling optimization should be achieved to solve practical problems. In 2020, the next-generation machine D-Wave Advantage quantum computer was released.²²⁾ It can handle up to 118 variables for the QUBO problem being expressed as a fully connected graph. Moreover, a quantum–classical hybrid solver, Leap Hybrid, was launched in 2020,²³⁾ and it can solve up to 20000-variable problems on a fully connected graph. Thus, with the development of improved hardware and software, the proposed hybrid LDC algorithm can solve more practical scheduling problems in the steel industry.

Here, we discuss the computation times in our algorithm. In our experiments, most of the hybrid LDC time, t_LDC, corresponds to the total cloud time t_cloud. Although QA is performed in 20 μs, t_cloud includes the processing time for generating 1000 samples, waiting for queued results, and communication delays. Moreover, the hybrid LDC algorithm requires multiple iterations of D-Wave sampling. As a result, the total cloud time takes a few tens of seconds. In addition, unembedding has a considerable impact on t_LDC as the problem size increases. This is because we use the energy minimization method for the unembedding to obtain high-quality solutions. The energy minimization method performs the unembedding while optimizing a local QUBO model, consequently increasing the computation time. However, because the current D-Wave machines are implemented on a sparse graph, the embedding and unembedding are essential. When the hardware on a denser graph is realized in the future, our approach will be greatly accelerated.

6. Conclusion

A novel hybrid LDC algorithm is proposed to increase the feasible size of a combinatorial optimization problem using quantum computers. The unique feature of D-Wave quantum computers to generate various solutions in a single run is exploited for the first time to calculate the upper bound of the total cost during each LDC update. We validated the effectiveness of the proposed algorithm for scheduling two processes typically found in steel manufacturing. Optimization was successful even when considering eight products, which establish a currently intractable problem without decomposition. The proposed algorithm is applicable to other types of combinatorial optimization problems, and it will likely contribute to scaling-up various problems to enable real-world applications.

References

• 1)   S.  Tomiyama,  S.  Hojo,  O.  Yamaguchi and  K.  Muramatsu: J. Jpn. Ind. Manag. Assoc., 70 (2020), 185, (in Japanese). https://doi.org/10.11221/jima.70.185
• 2)   E.  Boyda,  S.  Basu,  S.  Ganguly,  A.  Michaelis,  S.  Mukhopadhyay and  R. R.  Nemani: PLoS One, 12 (2017), e0172505. https://doi.org/10.1371/journal.pone.0172505
• 3)   F.  Neukart,  G.  Compostella,  C.  Seidel,  D.  von Dollen,  S.  Yarkoni and  B.  Parney: Front. ICT, 4 (2017), 29. https://doi.org/10.3389/fict.2017.00029
• 4)   D.  Venturelli,  D. J. J.  Marchand and  G.  Rojo: arxiv:1506.08479, (2016), https://arxiv.org/abs/1506.08479.
• 5)   K.  Ikeda,  Y.  Nakamura and  T. S.  Humble: Sci. Rep., 9 (2019), 12837. https://doi.org/10.1038/s41598-019-49172-3
• 6)   M.  Ohzeki,  A.  Miki,  M. J.  Miyama and  M.  Terabe: Front. Comput. Sci., 1 (2019), 9. https://doi.org/10.3389/fcomp.2019.00009
• 7)   M.  Ohzeki: Sci. Rep., 10 (2020), 3126. https://doi.org/10.1038/s41598-020-60022-5
• 8)   K.  Yonaga,  M. J.  Miyama and  M.  Ohzeki: arxiv:2012.06119, (2020), https://arxiv.org/abs/2012.06119.
• 9)   A.  Lucas: Front. Phys., 2 (2014), 5. https://doi.org/10.3389/fphy.2014.00005
• 10)   T.  Kadowaki and  H.  Nishimori: Phys. Rev. E, 58 (1998), 5355. https://link.aps.org/doi/10.1103/PhysRevE.58.5355
• 11)   R.  Martoňák,  G. E.  Santoro and  E.  Tosatti: Phys. Rev. E, 70 (2004), 057701. https://link.aps.org/doi/10.1103/PhysRevE.70.057701
• 12)   G. E.  Santoro,  R.  Martoňák,  E.  Tosatti and  R.  Car: Science, 295 (2002), 2427. https://doi.org/10.1126/science.1068774
• 13)   M.  Benedetti,  J.  Realpe-Gómez,  R.  Biswas and  A.  Perdomo-Ortiz: Phys. Rev. A, 94 (2016), 022308. https://link.aps.org/doi/10.1103/PhysRevA.94.022308
• 14)  D-Wave Systems, Inc.: D-Wave Ocean Software Documentation, https://docs.ocean.dwavesys.com/en/stable/index.html, (accessed 2021-11-22).
• 15)   J.  Cai,  W. G.  Macready and  A.  Roy: arxiv:1406.2741, (2014), https://arxiv.org/abs/1406.2741.
• 16)   A. M.  Geoffrion: Approaches to Integer Programming (Mathematical Programming Studies Vol. 2), ed. by M. L. Balinski, Springer, Berlin-Heidelberg, (1974), 82. https://doi.org/10.1007/BFb0120690
• 17)   M. L.  Fisher: Manage. Sci., 27 (1981), 1. https://doi.org/10.1287/mnsc.27.1.1
• 18)   P. B.  Luh,  D. J.  Hoitomt,  E.  Max and  K. R.  Pattipati: IEEE Trans. Robot. Autom., 6 (1990), 687. https://doi.org/10.1109/70.63271
• 19)   D. J.  Hoitomt,  P. B.  Luh and  K. R.  Pattipati: IEEE Trans. Robot. Autom., 9 (1993), 1. https://doi.org/10.1109/70.210791
• 20)   Y.  Zhang,  P. B.  Luh,  K.  Yoneda,  T.  Kano and  Y.  Kyoya: Proc. 1997 IEEE Int. Conf. on Control Applications, (USA), IEEE, Piscataway, NJ, (1997), 429. https://doi.org/10.1109/CCA.1997.627618
• 21)  IBM: IBM CPLEX Optimizer, https://www.ibm.com/analytics/cplex-optimizer, (accessed 2021-11-22).
• 22)   C.  McGeoch and  P.  Farré: The D-Wave Advantage System: An Overview, D-Wave Technical Report Series, D-Wave Systems, Inc., Burnaby, (2020), 14-1049A-A, https://www.dwavesys.com/media/s3qbjp3s/14-1049a-a_the_d-wave_advantage_system_an_overview.pdf, (accessed 2021-11-22).
• 23)   C.  McGeoch,  P.  Farré and  W.  Bernoudy: D-Wave Hybrid Solver Service + Advantage: Technology Update, D-Wave Technical Report Series, D-Wave Systems, Inc., Burnaby, (2020), 14-1048A-A, https://www.dwavesys.com/media/m2xbmlhs/14-1048a-a_d-wave_hybrid_solver_service_plus_advantage_technology_update.pdf, (accessed 2021-11-22).

Appendix

As a comparison with the results obtained from the D-Wave quantum computer, the same scheduling problems for the two processes were resolved using integer linear programming implemented in the commercial CPLEX software. This appendix presents the formulation of the corresponding problem. The parameters and notations are the same as those indicated in Sec. II unless otherwise specified.

In addition to the decision variables for real products x_p,i,k (p =1,2;1 ≤ i ≤ N), we prepare variables for dummy products x_p,i,k (p = 2;N + 1 ≤ i ≤ N + e) for process 2 (p = 2) to facilitate the formulation. The variable range of k is thus 1 ≤ k ≤ N for p = 1 and 1 ≤ k ≤ N + e for p = 2. Another variable y, which represents the consecutive processing of a product, is introduced to keep the formulation linear, unlike the QUBO model that is quadratic.   

$$y_{p,i_1,i_2,s}=\{\begin{array}{ll}1,\hfill & \mathrm{i}\mathrm{f}\mathrm{\hspace{0.17em} } \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{s}\mathrm{\hspace{0.17em} } i_1\mathrm{\hspace{0.17em} } \mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\hspace{0.17em} } i_2\mathrm{\hspace{0.17em} } \mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{o}\mathrm{\hspace{0.17em} } \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{\hspace{0.17em} } \hfill \\ \hfill & p\mathrm{\hspace{0.17em} } \mathrm{a}\mathrm{t}\mathrm{\hspace{0.17em} } \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{t}\mathrm{s}\mathrm{\hspace{0.17em} } s\mathrm{\hspace{0.17em} } \mathrm{a}\mathrm{n}\mathrm{d}\mathrm{\hspace{0.17em} } \left(s+1\right),\mathrm{\hspace{0.17em} } \mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{y}\hfill \\ 0,\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}$$(12)

The constraints for x are the same as those in Eqs. (1) and (2), respectively:   

$$\sum _{i=1}^N{x}_{1,i,s}=1\begin{array}{c}{\hspace{1em}}^{\forall }s\in \left\{1,2,\dots ,N\right\},\end{array}$$(13)

  

$$\sum _{s=1}^N{x}_{1,i,s}=1\begin{array}{c}{\hspace{1em}}^{\forall }i\in \left\{1,2,\dots ,N\right\},\end{array}$$(14)

  

$$\sum _{i=1}^{N+e}{x}_{2,i,s}=1{\hspace{1em}}^{\forall }s\in \left\{1,2,\dots ,N+e\right\},$$(15)

  

$$\sum _{s=1}^{N+e}{x}_{1,i,s}=1{\hspace{1em}}^{\forall }i\in \left\{1,2,\dots ,N+e\right\}.$$(16)

From the definition of y, the following constraints are derived:   

$$\sum _{i_2=1}^N{y}_{1,i_1,i_2,s}=x_{1,i_1,s}{\hspace{1em}}^{\forall }{i}_1\in \left\{1,2,\dots ,N\right\},{\mathrm{\hspace{0.17em} }}^{\forall }s\in \left\{1,2,\dots ,N-1\right\},$$(17)

  

$$\sum _{i_1=1}^N{y}_{1,i_1,i_2,s-1}=x_{1,i_2,s}{\hspace{1em}}^{\forall }{i}_2\in \left\{1,2,\dots ,N\right\},{\mathrm{\hspace{0.17em} }}^{\forall }s\in \left\{2,3,\dots ,N\right\},$$(18)

  

$$\begin{array}{c}\begin{array}{l}\sum _{i_2=1}^{N+e}{y}_{2,i_1,i_2,s}=x_{2,i_1,s}{\hspace{1em}}^{\forall }{i}_1\in \left\{1,2,\dots ,N+e\right\},\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}{}^{\forall }s\in \left\{1,2,\dots ,N+e-1\right\},\end{array}\end{array}$$(19)

  

$$\begin{array}{l}\sum _{i_1=1}^{N+e}{y}_{2,i_1,i_2,s-1}=x_{2,i_2,s}{\hspace{1em}}^{\forall }{i}_2\in \left\{1,2,\dots ,N+e\right\},\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}{}^{\forall }s\in \left\{2,3,\dots ,N+e\right\},\end{array}$$(20)

The time to process a product can be expressed as   

$$t_{p,i}=\sum _{s=1}^{N}s{x}_{p,i,s}{\hspace{1em}}^{\forall }i\in \left\{1,2,\dots ,N\right\}.$$(21)

The constraint for the processing order is the same as in Eq. (3):   

$$t_{1,i}\le t_{2,i}{\hspace{1em}}^{\forall }i\in \left\{1,2,\dots ,N\right\}.$$(22)

Finally, we set the constraint that a dummy product does not stay between two real products. This constraint can be represented by the number of transitions from a dummy product to a real product being at most one:   

$$\sum _{s=1}^{N+e-1}\sum _{i_1=N+1}^{N+e}\sum _{i_2=1}^N{y}_{2,i_1,i_2,s}+\sum _{i=1}^N{x}_{2,i,1}=1.$$(23)

The second term in the left-hand side of the equation represents the case in which a sequence starts with a real product. In this case, the sequence does not experience a transition from a dummy product to a real product.

The total cost is the sum of the three costs for late delivery, early delivery, and production group mismatch:   

$$\begin{array}{l}Cost\left(\mathit{x}\right)=w_d\sum _{i=1}^N\sum _{s=d_i+1}^{N+e}\left(s-d_i\right)x_{2,i,s}+w_e\sum _{i=1}^N\sum _{s=1}^{d_i-1}\left(d_i-s\right)x_{2,i,s}\\ +w_g\sum _{s=1}^{N-1}{\sum }_{(i_1,i_2)\in L{C}_1}{y}_{1,i_1,i_2,s}+w_g\sum _{s=1}^{N+e-1}{\sum }_{(i_1,i_2)\in L{C}_2}{y}_{2,i_1,i_2,s},\end{array}$$(24)

where LC_p (p = 1,2) are sets of pairs of two products i₁ and i₂ (i₁, i₂ ∈{1, 2, …, N}) that belong to different production groups.