Model Predictive Control for Finite Input Systems using the D-Wave Quantum Annealer

Model Predictive Contr ol f or Finite Input Systems using the D-W a ve Quantum Annealer Daisuke Inoue 1,* and Hiroaki Y oshida 1 1 T oy ota Central R&D Labs., Inc., Bunky o-ku, T okyo , 112-0004, Japan * daisuke-inoue@mosk.tytlabs .co .jp ABSTRA CT The D-W av e quantum annealer has emerged as a no vel computational architecture that is attracting significant interest, b ut there hav e been only a f ew practical algorithms exploiting the pow er of quantum annealers. Here we present a model predictiv e control (MPC) algorithm using a quantum annealer f or a system allo wing a finite number of input v alues. Such an MPC prob lem is classified as a non-deter ministic polynomial-time-hard combinatorial problem, and thus real-time sequential optimization is difficult to obtain with con ventional computational systems. W e circumvent this difficulty b y con ver ting the original MPC prob lem into a quadratic unconstrained binar y optimization problem, which is then solved by the D-Wa ve quantum annealer . T wo practical applications , namely stabilization of a spring-mass-damper system and dynamic audio quantization, are demonstrated. F or both, the D-Wa ve method e xhibits better performance than the classical simulated annealing method. Our results suggest new applications of quantum annealers in the direction of dynamic control prob lems. Introduction Since quantum annealer 2000Q was released from D-W ave Systems Inc., research on quantum computing has rapidly progressed 1 – 4 . The quantum annealer is a computational system specialized for solving combinatorial optimization problems, and it is expected to solve these problems with high accuracy and speed compared to classical computer systems 5 , 6 . Examples of using quantum annealers are expected to greatly expand in the future because the combinatorial optimization problem appears in v arious engineering fields, such as logistics 7 , finance 8 , transportation 9 , and social infrastructure 10 . Still, the class of problems that can be handled by quantum annealers is restricted compared with con ventional computers, and the current limitation on problem size further hinders the expansion of practical applications. In this context, we shed light on the potential for application of quantum annealers to a model pr edictive contr ol (MPC) , which is one of the modern control algorithms. In the MPC, the future trajectory of the controlled states of a target system is predicted using a dynamic model of the target, and the input is sequentially generated so as to minimize the ev aluation function for the future state trajectory 11 , 12 . Compared with other classical control methods, this method has flexibly to reflect the designer’ s control requirements in the ev aluation function 13 . The particular challenge here is that sequential optimization is computationally intensi ve, and thus this technique has been con v entionally applied only to systems with large time constants, such as chemical plants. Although recent computational dev elopments hav e gradually extended applicability to systems with small time constants, such as mechanical systems and electrical systems 12 , 14 , 15 , certain systems are still dif ficult to control by means of the MPC. A system that allows only a predetermined finite number of input values (here refer to as a finite input system ) is one example of such systems. The simplest systems consisting of on/off 16 input, and systems with finite digital bit input 17 are also included in this class of systems. Finding the optimal input for such systems becomes drastically difficult as the size of the problem increases, because this problem is classified as a non-deterministic polynomial-time (NP)-hard combinatorial optimization problem. Although there are se veral well-kno wn approximation methods, such as simulated annealing 18 , tabu search 19 , and genetic algorithms 20 , the computational costs for all of these methods are still too large 21 . In the present study , we propose a method to solve the MPC problem by using the quantum annealing with D-W a ve’ s quantum annealer 2000Q. So far , 2000Q has been used in se veral engineering fields 22 – 26 , but to the best of our knowledge, this is first time that the 2000Q is applied to dynamic control problems in which an optimization must be obtained sequentially at high speed. W e first give a method to transform the original MPC problem into a quadratic unconstrained binary optimization (QUBO) problem, which is the only class of problem that the 2000Q can solv e. W e then apply the proposed method to two scenarios. W e first consider the stabilization of a spring-mass-damper system , which is one of the most basic control problems to stabilize a physical system. The second scenario is dynamic quantization of an audio signal , in which the audio signal is dynamically quantized considering human auditory characteristics. In both scenarios, we verify that the performance is better than the case with the simulated annealing, which is one of the classical approximation algorithms. Results Model Predictive Control f or Finite Input Systems Let us consider the following discrete-time linear system as a control tar get: x ( t + 1 ) = Ax ( t ) + B 1 u ( t ) + B 2 a ( t ) , y ( t ) = Cx ( t ) + Du ( t ) + a ( t ) , (1) where t ∈ N denotes discrete time steps; x ( t ) ∈ R n denotes states; u ( t ) ∈ U ⊂ R m denotes inputs; a ( t ) ∈ R ` denotes given signals; y ( t ) ∈ R ` denotes outputs of the system and A ∈ R n × n , B 1 ∈ R n × m , B 2 ∈ R n × ` , C ∈ R ` × n , and D ∈ R ` × m are known constant matrices that represent the dynamics of the system. W e assume that system ( 1 ) is fully controllable and fully observable. In addition, U is a finite set composed of M vectors, and this set is defined as U = { u i | u i ∈ R m , i = 1 , . . . , M } . (2) System ( 1 ) differs from the usual control system in that input u ( t ) at each time step only takes elements of the set U . T o control system ( 1 ), we design the following e v aluation function: H = t + N − 1 ∑ k = t  y ( k ) > Qy ( k ) + u ( k ) > Ru ( k )  , (3) where Q ∈ R ` × ` and R ∈ R m × m are positi ve definite symmetric matrices called weight matrices , which are the design parameters. N ∈ N is also a design parameter representing the length of the ev aluation section, which is called the pr ediction horizon . The problem of finding input sequence u ( t ) : = [ u ( t ) · · · u ( t + N − 1 )] > that minimizes the ev aluation function ( 3 ) , while observing states x ( t ) at time t , u ∗ ( t ) : = arg min u ( t ) ∈ U N H , (4) is called an optimal contr ol pr oblem . In the MPC, Eq. ( 4 ) is solved at each time step and the given input is applied to the control target at each time step. This is done in the present study by using the quantum annealer 2000Q. W e exploit the discrete nature of the input by transforming the abov e control problem into a quadratic unconstr ained binary optimization (QUBO) problem. The Methods section provides details of the transformation. The ev aluation function of the resulting QUBO problem takes the form H = b ( t ) > J ( t ) b ( t ) + h ( t ) > b ( t ) + c 0 ( t ) , (5) with the binary design variable b ( t ) ∈ { 0 , 1 } mLN ( L is a natural number satisfying M = 2 L ). The matrix J ( t ) , vector h ( t ) , and function c 0 ( t ) are precisely defined in the Methods section. After arriving at this equi valent Q UBO expression, which is compatible with quantum annealers, we are able to solve the optimal control problem using the quantum annealer 2000Q. In the following discussion, we consider tw o scenarios to ev aluate the performance of the proposed method: • Stabilization of spring-mass-damper system: For an unstable physical system, optimal stabilization control is formulated as an MPC problem. • Dynamic quantization of audio signal: Audio quantization that considers human auditory characteristics is formulated as an MPC problem. T o validate the performance of the proposed method, we compare the follo wing three solutions: • Exact solution: A brute-force search is used to find the minimum value of Eq. ( 4 ) by enumerating all possible combinations of variables. • Approximate solution using the simulated annealing: Eq. ( 4 ) is con verted into a QUBO problem and solved using the simulated annealing, which is one of the classical optimization methods for combinatorial optimization. See Ref. 18 for the detailed explanation of the simulated annealing. • Approximate solution using the quantum annealing: The con verted Q UBO problem is solved using the quantum annealer 2000Q, which is physically located in British Columbia, Canada. W e use HTTPS communication to command the ex ecution of optimization from our local en vironment. The Solver API library provided by D-W ave was used for the quantum annealing method, and the anneal library included in the D-W ave Ocean software packag e was used for the simulated annealing method. Each method was implemented using programming language Python, and executed using a W indows machine with 3.40 GHz clock frequency and 16.0 GB memory capacity . In each method, the default parameters of the pro vided program were used. 2/ 11 Stabilization of Spring-Mass-Damper System Consider a continuous-time linear system as the control target of Eq. ( 1 ): ˙ x ( τ ) =  0 1 − 10 1  x ( τ ) +  0 1  u ( τ ) , x ( 0 ) =  10 − 10  , (6) y ( τ ) = x ( τ ) , (7) where the values are in SI units. System ( 6 ) is a spring-mass-damper system with a negati ve damper coefficient, which is a controllable but unstable system 27 . By discretizing system ( 6 ) with zero-order hold of the sampling period 0 . 1 s , the following discrete-time linear system is obtained: x ( t + 1 ) = ˆ Ax ( t ) + ˆ Bu ( t ) , (8) y ( t ) = x ( t ) , (9) ˆ A : =  0 . 94872313 0 . 1034271 − 1 . 03427103 1 . 05215023  , (10) ˆ B : =  0 . 00512769 0 . 1034271  . (11) Considering the correspondence between system ( 6 ) and system ( 1 ), the following relation is obtained: A = ˆ A , B 1 = ˆ B , B 2 = 0 , C = I , D = 0 , (12) where I denotes identity matrix. W e design the weight matrices in ev aluation function ( 3 ) as Q = 10 5  1 0 0 1  , R = 1 . (13) For the set of inputs gi ven in Eq. ( 2 ), we use ˆ U = {− 60 , − 45 , − 30 , − 15 , 0 , 15 , 30 , 45 } (14) and adopt N = 6 for the prediction horizon. W e sho w the generated input sequence in Fig. 1 -(a) and the time response of the states in Fig. 1 -(b) and (c). First, for the input in Fig. 1 -(a), the proposed method using the quantum annealing produces an input that is closer to the exact solution than that of the simulated annealing. Also, in Fig. 1 -(b) and (c), two states for the quantum annealing con ver ge to the origin more rapidly than the simulated annealing. Comparison of the simulated annealing and the quantum annealing sho ws that the latter method is more likely to output a solution closer to the global optimum solution. In Fig. 1 -(d), we show the total v alues of the e valuation function ( 3 ), that is, ˆ H : = 99 ∑ t = 0  y ( t ) > Qy ( t ) + u ( t ) > Ru ( t )  . (15) The v alue of Eq. ( 15 ) in the quantum annealing is smaller than the corresponding v alue in the simulated annealing. The a verage value in each method is 4 . 22 × 10 8 and 5 . 62 × 10 8 , respectiv ely , and the value in the quantum annealing is suppressed to 0 . 75 times that in the simulated annealing. The standard deviations are 1 . 37 × 10 8 and 5 . 19 × 10 8 , respecti vely , which means the variance of the quantum annealing solution is smaller than that of the simulated annealing. When the quantum annealing method is compared with the exact solution method, the average v alue of the quantum annealing is 1 . 01 times that of the exact solution. T able 1 shows the elapsed time for 100 time steps in each method to find the input. Quantum annealing takes 96 . 8 s to perform the optimization. Most of this time is occupied by communication time between the local en vironment and Canada, where the 2000Q is located. The calculation time on the 2000Q is 7 . 95 × 10 − 1 s , which is 0 . 8% of the total time, and the time required for the annealing operation is as short as 2 . 0 × 10 − 3 s . In contrast, the calculation time for the simulated annealing and the exact solution are 5 . 51 × 10 − 1 s and 2 . 46 × 10 2 s , respectiv ely . T o control the physical system giv en in Eq. ( 6 ) in real time, it is desirable to generate the input at each step at a time shorter than the sampling period of 0 . 1 s . In the numerical example abov e, the time required for one step of input generation in each method is 2 . 46 s in the brute-force search, 5 . 51 × 10 − 3 s in the simulated annealing, and 0 . 97 s in the quantum annealing. Among them, only the simulated annealing achie ves input generation within 0 . 1 s . Howe ver , the proposed method takes 7 . 95 × 10 − 3 s to generate inputs when communication time is not counted. This suggests that if the quantum annealing can be performed in the local environment, real-time control is achie ved for control targets with small time constants, such as mechanical systems and electrical systems. 3/ 11 0 20 40 60 80 100 t − 60 − 40 − 20 0 20 40 u ( t ) Exact Sim ulated Annealing Quan tum Annealing 0 20 40 60 80 100 t − 5 0 5 10 x 1 ( t ) Exact Sim ulated Annealing Quan tum Annealing (a) (b) 0 20 40 60 80 100 t − 20 − 10 0 x 2 ( t ) Exact Sim ulated Annealing Quan tum Annealing Exact SA QA 0 2 4 6 ˆ H × 10 8 (d) (c) Figure 1. Result of the stabilization of a spring-mass-damper system. (a) Input u ( t ) generated with each method. (b) Time response of the first component of the state x ( t ) . (c) T ime response of the second component of the state x ( t ) . In (a)–(c), the input generated with the quantum annealing is compared to the input from the simulated annealing and the brute-force search. (d) V alue of ev aluation function gi ven in Eq. ( 15 ). The box-and-whisker plots are created by performing 10 numerical simulations. Dynamic Quantization of Audio Signal Here, we consider the problem of audio signal quantization, in which real-time input signal generation is not as important as in physical systems. Audio quantization is not only important from the viewpoint of data compression, but it is also indispensable when analog signals are expressed as digital signals 28 . W e first show that the quantization to improve the naturalness of human auditory signals is attributed to the optimal control problem for the finite input systems. Consider the problem of quantizing the audio signal { a ( t ) } , t ∈ N , a ( t ) ∈ R . Let { u ( t ) } , t ∈ N , u ( t ) ∈ R be the signal after quantization; that is, at each time t ∈ N , u ( t ) ∈ U is satisfied for a pre-fix ed set U = { u i | u i ∈ R , i = 1 , . . . , M } consisting of M elements. The goal of quantizing audio signal is to reduce the data size while retaining the naturalness of the human voice content. Human ears perceiv e sound with different sensiti vities at dif ferent frequencies 29 , and its auditory characteristics are expressed as per ception filters 30 . Among the av ailable filters, we consider a linear time-in variant filter 31 : P ( z ) = 1 + ¯ C ( zI − ¯ A ) − 1 ¯ B , (16) where ¯ A ∈ R n × n , ¯ B ∈ R n × 1 , and ¯ C ∈ R 1 × n are filter parameters. The quantization error after the signal passes through the perception filter is defined as y ( t ) : = P ( z )( a ( t ) − u ( t )) , (17) and the ev aluation function is defined as H = t + N − 1 ∑ k = t ( y ( k )) 2 . (18) 4/ 11 T able 1. Elapsed time for each method. Method Elapsed T ime [s] Exact solution 2 . 46 × 10 2 Simulated Annealing 5 . 51 × 10 − 1 Quantum Annealing (T otal) 9 . 68 × 10 1 QA (Computation T ime) 7 . 95 × 10 − 1 QA (Annealing T ime) 2 . 0 × 10 − 3 (a) (b) 10 2 10 3 10 4 0 20 Magnitude (dB) 10 2 10 3 10 4 F requency (Hz) − 90 − 45 0 Phase (deg) exact SA QA 0 5 10 15 20 ¯ H Figure 2. Result of dynamic audio quantization. (a) Bode plot of perception filter giv en in Eq. ( 22 ). (b) V alue of ev aluation function giv en in Eq. ( 27 ) . The solution from the quantum annealing (green) is compared to that from the simulated annealing (orange) and exact solution (blue). The box-and-whisker plot is created by performing 10 numerical simulations. Gi ven perception filter P and audio signals { a ( t ) } , the problem of finding an input that minimizes the ev aluation function is represented as an optimal control problem. Using Z-transform, Eq. ( 17 ) is expressed as the following state-space model: x ( t + 1 ) = ¯ Ax ( t ) + ¯ B ( a ( t ) − u ( t )) , y ( t ) = ¯ Cx ( t ) + a ( t ) − u ( t ) . (19) The correspondence between system ( 19 ) and system ( 1 ) is as follows: A = ¯ A , B 1 = − ¯ B , B 2 = ¯ B , C = ¯ C , D = − 1 , (20) and the correspondence between ev aluation function ( 18 ) and ( 3 ) is as follo ws: Q = 1 , R = 0 . (21) In this way , the problem of quantizing the signal to minimize evaluation function ( 18 ) is expressed as an optimal control problem. For the numerical simulation, we use the follo wing filter proposed in Ref. 31 (with a sampling frequency of 44 . 1 kHz): P ( z ) = 1 + z − 1 2 . 245 − 0 . 664 z − 1 1 − 1 . 335 z − 1 + 0 . 644 z − 2 . (22) The Bode plot for this filter is depicted in Fig. 2 -(a). Eq. ( 22 ) is essentially a low-pass filter , which also enhances frequencies in the band around 4 kHz. The parameters when expressing Eq. ( 22 ) in the state space are as follows: ¯ A =  1 . 335 − 0 . 644 1 . 0 0  , (23) ¯ B =  1 . 0 0  , (24) ¯ C =  2 . 245 − 0 . 664  . (25) 5/ 11 10 0 10 2 10 4 10 6 Iteration n um b er 5 × 10 8 6 × 10 8 7 × 10 8 8 × 10 8 ˆ H Sim ulated Annealing Quan tum Annealing Figure 3. Iteration number for the simulated annealing versus the v alue of the ev aluation function gi ven in Eq. ( 15 ) . At each value of iteration number , 10 computations are performed (blue dot) and the average v alue is plotted (blue solid line). The corresponding value of the e valuation function obtained with the quantum annealing method is also shown by the green dashed line. W e adopt the following input sets: ¯ U = {− 1 . 0 , − 0 . 5 , 0 , 0 . 5 } . (26) W e use N = 10 for the prediction horizon. For the pre-quantized audio signal { a ( t ) } , we use a random number sequence of 100 steps following a uniform distrib ution on [ − 1 , 1 ) . Comparison of the values obtained by e v aluating Eq. ( 18 ) over the duration with three methods, that is, ¯ H : = 99 ∑ t = 0 ( y ( t )) 2 , (27) is shown in Fig. 2 -(b). Comparison of simulated and quantum annealing results shows that the value of Eq. ( 27 ) is smaller for the quantum annealing. The average values are 13 . 9 and 11 . 2 respectiv ely , and the value of quantum annealing is 0 . 81 times that of the simulated annealing. The quantum annealing solutions ha ve less variance; the v ariances of the simulated and quantum annealing are 1 . 94 and 1 . 78 × 10 − 2 , respectiv ely . When comparing the exact solution and the quantum annealing, the av erage v alue from the quantum annealing is as 2 . 1 times that of the e xact solution, which is a conserv ati ve result when compared with the previous e xample. T able 2. Elapsed time for each method. Method Elapsed T ime [s] Exact solution 8 . 61 × 10 2 Simulated Annealing 3 . 74 × 10 − 1 Quantum Annealing (T otal) 8 . 6 × 10 1 QA (Computation T ime) 7 . 91 × 10 − 1 QA (Annealing T ime) 2 . 0 × 10 − 3 T able 2 shows the elapsed time required for the optimization after 100 steps in each method. The quantum annealing takes 8 . 6 × 10 1 s to perform the optimization. As mentioned in the previous example, HTTPS communication accounts for most of the time, and the 2000Q calculation time is 7 . 91 × 10 − 1 s , only 0 . 9 % of the total time. In addition, the time taken for the annealing itself is e ven shorter, 2 . 0 × 10 − 3 s . In contrast, the calculation time for the simulated annealing and the exact solution is 3 . 74 × 10 − 1 s and 8 . 61 × 10 2 s , respectively . The proposed method using the quantum annealing enables higher-speed quantization than the brute-force search and higher performance than the simulated annealing. Discussion In the 2000Q, the ev aluation function is minimized using a natural phenomenon called a quantum effect, and the solution is not updated by the algorithmic iteration, and indeed this is the essential dif ference from the classical simulated annealing method. 6/ 11 (a) (b) 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 J ij (absolute v alue) 0 20 40 60 80 F requency 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 J ij (absolute v alue) 0 25 50 75 100 125 F requency Figure 4. Histogram of the absolute value of the QUBO matrix. (a) Stabilization of spring-mass-damper system. (b) Dynamic audio quantization. The latter employs an iterativ e algorithm, where the neighborhoods of the current solution are examined at each time step, with probabilistically judging whether to update the solution or to remain in the current state. The number of iterations is a design parameter (in the numerical experiment abo ve, the library default v alue of 1000 is used as the number of iterations). Clearly the larger the number of iterations is, the more accurate the obtained solution is, but the longer the required computational time becomes. T o examine this trade-off, we show in Fig. 3 the relationship between the iteration number of the simulated annealing method and the quality of the solution in the spring-mass-damper system. The value of the ev aluation function decreases as the iteration number is increased, and about one million steps are required to reach the lev el of the quantum annealing result. The av erage computational time of 1 . 0 × 10 6 iterations was 3 . 9 seconds, which is 4 . 9 times longer than 7 . 95 × 10 − 1 second, the time required for the quantum annealing. W e also mention the restriction of the current quantum annealing using the 2000Q, which is equipped no more than 2,048 qubits. Each qubit is not coupled with all the other qubits, but instead the assemblage has a chimera structur e , in which closely connected 8-bit units are arranged vertically and horizontally 32 . For this reason, the variables of a gi ven QUBO problem cannot be directly assigned to physical qubits. The method of con verting the gi ven graph structure to the chimera structure is called minor embedding , which is realized by e xpressing one logical v ariable with strongly coupled multiple physical qubits. The embedding problem itself is actually an NP-hard problem, and we used the API that D-W av e provides to perform embedding with heuristic algorithms 33 . Because multiple physical bits are sacrificed in minor embedding, the effecti ve number of v ariables decreases depending on the density of the original problem, that is, the number of non-zero elements of a QUBO matrix. In the chimera structure, N 2 / 4 physical qubits are needed for representing a fully connected N -variable problem, which means that the maximum number of v ariables that the 2000Q can use is as small as 64 when the original problem has a fully connected structure. This implies that mLN ≤ 64 must be satisfied for the dimension b ( t ) in Eq. ( 5 ) . In addition to the limited number of variables, it is known that the quality of the solution degrades when the original problem has a dense structure 34 . T o check the sparseness of our formulated problem, we e xamine the value of all components of the QUBO matrix [ J i j in Eq. ( 45 ) ] and plot them as histograms in Fig. 4 . W e show the result of stabilization of the spring-mass-damper system in Fig. 4 -(a) and of dynamic audio quantization in Fig. 4 -(b). In both cases, we plot the absolute v alue of the QUBO matrix and normalize them such that max i , j J i j = 1 . Many components hav e values close to zero, and qualitati vely , 51 % of components are smaller than 0 . 1 for both cases, and the distribution is nearly exponential. This implies that the problem considered in this study is relatively sparse and matches the structure of the problems that the 2000Q can solve. Still, the sparse chimera structure limits the number of v ariables and the accuracy of the solution. D-W ave Systems Inc. recognizes this and has stated that they will adopt a ne w graph structure that increases the couplers per qubit in the next generation of machines 35 . Also, the y plan to double the number of qubits in the ne xt 2 years. Hence, the restriction on the problem size will soon be relaxed and the accuracy of solution will be improv ed in the future. W e expect D-W ave’ s next-generation machine to be capable of solving MPC problems with higher performance than the current generation machines. One of the important results obtained in the present study is short ef fectiv e computational time of the quantum annealing compared with that of the simulated annealing, as shown in T ables 1 and 2 . In the MPC, optimization must be performed sequentially at each time step for input calculation, and the performance deteriorates significantly when the input generation is not synchronized with the control cycle 36 . If the problem at each time step is NP-hard, as in this study , this performance 7/ 11 degradation becomes a critical problem. Our results of using the quantum annealing to obtain a solution within a control period suggest the possibility of expanding the classes of systems that can be controlled in real time. Methods T ransformation of MPC to QUBO W e first represent ev aluation function ( 3 ) as a quadratic form for the input vector u ( t ) . First, the dynamics of the system ( 1 ) is represented as y ( t ) = A x ( t ) + B 1 u ( t ) + B 2 a ( t ) , (28) with the following definitions: y ( t ) : =    y ( t ) . . . y ( t + N − 1 )    ∈ R ` N , (29) a ( t ) : =    a ( t ) . . . a ( t + N − 1 )    ∈ R ` N , (30) A : =      C CA . . . CA N − 1      ∈ R ` N × n , (31) B 1 : =       D 0 · · · 0 CB 1 D . . . . . . . . . . . . . . . 0 CA N − 2 B 1 · · · CB 1 D       ∈ R ` N × mN , (32) B 2 : =       1 0 · · · 0 CB 2 1 . . . . . . . . . . . . . . . 0 CA N − 2 B 2 · · · CB 2 1       ∈ R N × N . (33) Also, ev aluation function ( 3 ) is expressed as H = x ( t ) > Qx ( t ) + u ( t ) > Ru ( t ) (34) by defining Q : =       Q 0 · · · 0 0 Q . . . . . . . . . . . . . . . 0 0 0 · · · Q       ∈ R ` N × ` N , (35) R : =       R 0 · · · 0 0 R . . . . . . . . . . . . . . . 0 0 0 · · · R       ∈ R mN × mN . (36) By substituting Eq. ( 28 ) into Eq. ( 34 ), we obtain the following quadratic form: H = u ( t ) > ( B > 1 QB 1 + R ) u ( t ) + 2 ( A x ( t ) + B 2 a ( t )) > QB 1 u ( t ) + c ( t ) , (37) 8/ 11 where we assumed that Q is symmetrical and that c ( t ) is a constant independent of u ( t ) . The need to search the set U N makes it dif ficult to find the minima u ( t ) of the right side of Eq. ( 37 ) . U N consists of M N elements (where M is the number of elements of the set U and N is the control horizon), and as N and M increases, the search becomes impractical. F or this reason, when M and N are large, it is necessary to gi ve up an exact solution and instead adopt an approximate solution. This study explore the quantum annealer as one of the approximation methods. Howe ver , because this method only accepts QUBO problems, we next need to conv ert Eq. ( 37 ) to a QUBO problem. F or the set U defined by Eq. ( 2 ) , we add the follo wing assumptions: let M (the number of elements of sets U ) satisfy M = 2 L for some natural number L ∈ N , and for all i , let u i ∈ U satisfy u i =    K 1  − 2 L − 1 b 1 , L + ∑ L − 1 ` = 1 2 ` − 1 b 1 ,`  . . . K m  − 2 L − 1 b m , L + ∑ L − 1 ` = 1 2 ` − 1 b m ,`     , (38) where b j ,` ∈ { 0 , 1 } ( j = 1 , . . . , m , ` = 1 , . . . , L ) is a binary variable that takes a value of either 0 or 1 and K j ∈ R ( j = 1 , . . . , m ) is a scaling parameter . For e xample, the set ˆ U in Eq. ( 14 ) and ¯ U in Eq. ( 26 ) are expressed by Eq. ( 38 ) using m = 1 , L = 3 , and K 1 = 15 and m = 1 , L = 2 , and K 1 = 0 . 5 , respecti vely . The assumption that the input set U is expressed as Eq. ( 38 ) is made to express the ev aluation function ( 3 ) in a QUBO form. W e note that the expression in terms of Eq. ( 38 ) restricts the input set to the discrete v alues with a regular interv al. Nevertheless, this is v alid for v arious practical systems, and wider classes of problems should be dealt with by means of appropriate variable transformations, which is indeed included in our future studies. Eq. ( 38 ) can be viewed as a linear mapping from { 0 , 1 } mL to U , and its expression as a matrix is u i = E b , (39) E : =       E 1 0 · · · 0 0 E 2 . . . . . . . . . . . . . . . 0 0 0 · · · E m       ∈ R m × mL , (40) E i : =  K i · · · K i 2 L − 1 − K i 2 L  ∈ R 1 × L , (41) b : =  b 1 , 1 · · · b 1 , L · · · b m , 1 · · · b m , L  > ∈ { 0 , 1 } mL . In addition, by defining following block matrices, E : =       E 0 · · · 0 0 E . . . . . . . . . . . . . . . 0 0 0 · · · E       ∈ R mN × mLN , (42) b ( t ) : =    b ( t ) . . . b ( t + N − 1 )    ∈ { 0 , 1 } mLN , (43) Eq. ( 37 ) can be written as a QUBO as follo ws: H = b ( t ) > J ( t ) b ( t ) + h ( t ) > b ( t ) + c 0 ( t ) , (44) where c 0 ( t ) is a constant independent of b ( t ) and J ( t ) and h ( t ) are defined as J ( t ) : = E > ( B > 1 QB 1 + R ) E , (45) h ( t ) : = 2 ( A x ( t ) + B 2 a ( t )) > QB 1 E . (46) In this way , the problem of determining the right-hand side of Eq. ( 3 ) is expressed as a QUBO that minimizes Eq. ( 44 ) on b ( t ) . Once the QUBO is solved by the 2000Q, the original input u ( t ) is recovered by substituting solution b ( t ) into Eq. ( 39 ). 9/ 11 References 1. Kado waki, T . & Nishimori, H. Quantum Annealing in the T ransverse Ising Model. Phys. Rev. E 58 , 5355 (1998). 2. Johnson, M. W . et al. Quantum Annealing with Manufactured Spins. Natur e 473 , 194 (2011). 3. O’Gorman, B., Babb ush, R., Perdomo-Ortiz, A., Aspuru-Guzik, A. & Smelyanskiy , V . Bayesian network structure learning using quantum annealing. The Eur. Phys. J. Special T op. 224 , 163–188 (2015). 4. V enturelli, D., Marchand, D. J. & Rojo, G. Quantum annealing implementation of job-shop scheduling. arXiv pr eprint arXiv:1506.08479 (2015). 5. McGeoch, C. C. & W ang, C. Experimental Evaluation of an Adiabiatic Quantum System for Combinatorial Optimization. In Pr oceedings of the A CM International Confer ence on Computing F r ontiers , 23 (A CM, 2013). 6. Denche v , V . S. et al. What Is the Computational V alue of Finite-Range Tunneling? Phys. Rev. X 6 , 031015 (2016). 7. Sbihi, A. & Eglese, R. W . Combinatorial optimization and green logistics. Annals Oper. Res. 175 , 159–175 (2010). 8. Chinchuluun, A. & P ardalos, P . M. A survey of recent de velopments in multiobjecti ve optimization. Annals Oper. Res. 154 , 29–50 (2007). 9. T oth, P . & V igo, D. The V ehicle Routing Pr oblem (SIAM, 2002). 10. Frangopol, D. M. & Liu, M. Maintenance and management of civil infrastructure based on condition, safety , optimization, and life-cycle cost ∗ . Struct. infr astructur e engineering 3 , 29–41 (2007). 11. Rawlings, J. B. T utorial ov erview of model predicti ve control. IEEE Contr ol. Syst. 20 , 38–52 (2000). 12. Morari, M. & Lee, J. H. Model Predicti ve Control: Past, Present and Future. Comput. & Chem. Eng. 23 , 667–682 (1999). 13. K ouro, S., Cortés, P ., V argas, R., Ammann, U. & Rodríguez, J. Model predicti ve control—A simple and po werful method to control power con verters. IEEE T ransactions on industrial electr onics 56 , 1826–1838 (2009). 14. Qin, S. J. & Badgwell, T . A. A Surve y of Industrial Model Predictive Control T echnology. Contr ol. Eng. Pract. 11 , 733–764 (2003). 15. Lin, S., De Schutter , B., Xi, Y . & Hellendoorn, H. Fast model predicti ve control for urban road networks via MILP. IEEE T ransactions on Intell. T ransp. Syst. 12 , 846–856 (2011). 16. W onham, W . M. & Johnson, C. Optimal bang-bang control with quadratic performance index. J. Basic Eng. 86 , 107–115 (1964). 17. Bemporad, A. & Morari, M. Control of Systems Integrating Logic, Dynamics, and Constraints. Automatica 35 , 407–427 (1999). 18. Suman, B. & Kumar , P . A Survey of Simulated Annealing as a T ool for Single and Multiobjective Optimization. J. operational r esear ch society 57 , 1143–1160 (2006). 19. Glov er , F . & Laguna, M. T abu Search. In Handbook of Combinatorial Optimization , 2093–2229 (Springer, 1998). 20. Mitchell, M. Handbook of Genetic Algorithms. Artif. Intell. 100 , 325–330 (1998). 21. V azirani, V . V . Appr oximation Algorithms (Springer Science & Business Media, 2013). 22. O’Malley , D., V esselinov , V . V ., Alexandro v , B. S. & Alexandro v , L. B. Nonnegativ e/binary matrix factorization with a D-W av e quantum annealer . arXiv preprint arXiv:1704.01605 (2017). 23. Ohzeki, M., Okada, S., T erabe, M. & T aguchi, S. Optimization of neural networks via finite-value quantum fluctuations. Sci. r eports 8 , 9950 (2018). 24. T ran, T . T . et al. A hybrid quantum-classic al approach to solving scheduling problems. In Ninth Annual Symposium on Combinatorial Sear ch (2016). 25. Nev en, H., Denchev , V . S., Rose, G. & Macready , W . G. Training a Binary Classifier with the Quantum Adiabatic Algorithm. arXiv:0811.0416 [quant-ph] (2008). 0811.0416 . 26. Nev en, H., Denche v , V . S., Rose, G. & Macready , W . G. Training a Large Scale Classifier with the Quantum Adiabatic Algorithm. arXiv:0912.0779 [quant-ph] (2009). 0912.0779 . 27. Cairano, S. D., Bemporad, A., K olmanovsk y , I. V . & Hrov at, D. Model predictive control of magnetically actuated mass spring dampers for automotiv e applications. Int. J. Contr ol. 80 , 1701–1716 (2007). 28. Mitchell, J. L. Introduction to digital audio coding and standards. J. Electr on. Imaging 13 , 399 (2004). 10/ 11 29. Zwicker , E. & Fastl, H. Psychoacoustics: Facts and models. Springer Sci. & Bus. Media (2013). 30. Pohlmann, K. C. Principles of Digital Audio (McGraw-Hill Professional, 2000), 4th edn. 31. W annamaker , R. Psychoacoustically Optimal Noise Shaping. J. Audio Eng. Soc. 40 , 611–620 (1992). 32. Boothby , T ., King, A. D. & Roy , A. Fast clique minor generation in chimera qubit connectivity graphs. Quantum Inf. Pr ocess. 15 , 495–508, DOI: 10.1007/s11128- 015- 1150- 6 (2016). 33. Cai, J., Macready , W . G. & Roy , A. A practical heuristic for finding graph minors. arXiv:1406.2741 [quant-ph] (2014). 1406.2741 . 34. Hamerly , R. et al. Experimental in vestigation of performance dif ferences between coherent Ising machines and a quantum annealer . Sci. Adv. 5 , DOI: 10.1126/sciadv .aau0823 (2019). https://adv ances.sciencemag.org/content/5/5/eaau0823.full.pdf . 35. W . Johnson, M. Future Hardware Directions of Quantum Annealing. In Qubits Eur ope 2018 D-W ave Users Confer ence (Munich, 2018). 36. Zhang, W ., Branicky , M. S. & Phillips, S. M. Stability of networked control systems. IEEE contr ol systems magazine 21 , 84–99 (2001). Ackno wledgements The authors would like to thank Dr . Akihisa Okada, Dr . T adayoshi Matsumori, Dr . Y uji Ito and Dr . Kiyosumi Kidono of T oyota Central R&D Labs. Inc., and Dr . Y oshiki Matsuda of Fixstars Corporation for the useful discussions. 11/ 11