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

The D-Wave quantum annealer has emerged as a novel computational architecture that is attracting significant interest, but there have been only a few practical algorithms exploiting the power of quantum annealers. Here we present a model predictive control (MPC) algorithm using a quantum annealer for a system allowing a finite number of input values. Such an MPC problem is classified as a non-deterministic polynomial-time-hard combinatorial problem, and thus real-time sequential optimization is difficult to obtain with conventional computational systems. We circumvent this difficulty by converting the original MPC problem into a quadratic unconstrained binary optimization problem, which is then solved by the D-Wave quantum annealer. Two practical applications, namely stabilization of a spring-mass-damper system and dynamic audio quantization, are demonstrated. For both, the D-Wave method exhibits better performance than the classical simulated annealing method. Our results suggest new applications of quantum annealers in the direction of dynamic control problems.


Model predictive control for finite input Systems using the D-Wave Quantum Annealer
Daisuke inoue * & Hiroaki Yoshida the D-Wave quantum annealer has emerged as a novel computational architecture that is attracting significant interest, but there have been only a few practical algorithms exploiting the power of quantum annealers. Here we present a model predictive control (MPC) algorithm using a quantum annealer for a system allowing a finite number of input values. Such an MPC problem is classified as a non-deterministic polynomial-time-hard combinatorial problem, and thus real-time sequential optimization is difficult to obtain with conventional computational systems. We circumvent this difficulty by converting the original MPC problem into a quadratic unconstrained binary optimization problem, which is then solved by the D-Wave quantum annealer. Two practical applications, namely stabilization of a spring-mass-damper system and dynamic audio quantization, are demonstrated. For both, the D-Wave method exhibits better performance than the classical simulated annealing method. Our results suggest new applications of quantum annealers in the direction of dynamic control problems.
Since quantum annealer 2000Q was released from D-Wave Systems Inc., research on quantum computing has rapidly progressed [1][2][3][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 various 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 conventional 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 predictive control (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 evaluation 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 evaluation function 13 . The particular challenge here is that sequential optimization is computationally intensive, and thus this technique has been conventionally applied only to systems with large time constants, such as chemical plants.
Although recent computational developments have gradually extended applicability to systems with small time constants, such as mechanical systems and electrical systems 12,14,15 , certain systems are still difficult 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 several well-known 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-Wave's quantum annealer 2000Q. So far, 2000Q has been used in several engineering fields [22][23][24][25][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. We 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 Toyota Central R&D Labs., Inc., Bunkyo-ku, Tokyo, 112-0004, Japan. *email: daisuke-inoue@mosk.tytlabs.co.jp open Results Model predictive control for finite input systems. Let us consider the following discrete-time linear system as a control target: ( ) denotes given signals;  ∈  y t ( ) denotes outputs of the system and , , n n n m n n 1 2 , and  ∈ ×  D m are known constant matrices that represent the dynamics of the system. We assume that system 1 is fully controllable and fully observable. In addition,  is a finite set composed of M vectors, and this set is defined as i i m  System 1 differs from the usual control system in that input u(t) at each time step only takes elements of the set  .
To control system 1, we design the following evaluation function: m m are positive definite symmetric matrices called weight matrices, which are the design parameters.
is also a design parameter representing the length of the evaluation section, which is called the prediction horizon. The problem of finding input sequence u(t) := [u(t) ⋯ u(t + N − 1)] ⊤ that minimizes the evaluation function 3, while observing states x(t) at time t, is called an optimal control problem. 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. We exploit the discrete nature of the input by transforming the above control problem into a quadratic unconstrained binary optimization (QUBO) problem. The Methods section provides details of the transformation. The evaluation function of the resulting QUBO problem takes the form 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 t ( ) are precisely defined in the Methods section. After arriving at this equivalent QUBO 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 two scenarios to evaluate 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.
• To validate the performance of the proposed method, we compare the following 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 converted 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 converted QUBO problem is solved using the quantum annealer 2000Q, which is physically located in British Columbia, Canada. We use HTTPS communication to command the execution of optimization from our local environment.
• The Solver API library provided by D-Wave was used for the quantum annealing method, and the anneal library included in the D-Wave Ocean software package was used for the simulated annealing method. Each method was implemented using programming language Python, and executed using a Windows machine with 3.40 GHz clock frequency and 16.0 GB memory capacity. In each method, the default parameters of the provided program were used.
www.nature.com/scientificreports www.nature.com/scientificreports/ Stabilization of spring-mass-damper system. Consider a continuous-time linear system as the control target of Eq. 1: where the values are in SI units. System 6 is a spring-mass-damper system with a negative 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: Considering the correspondence between system 6 and system 1, the following relation is obtained: where I denotes identity matrix. We design the weight matrices in evaluation function 3 as 5 For the set of inputs given in Eq. 2, we use and adopt N = 6 for the prediction horizon.
We show the generated input sequence in Fig. 1(a) and the time response of the states in Fig. 1(b,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,c), two states for the quantum annealing converge to the origin more rapidly than the simulated annealing. Comparison of the simulated annealing and the quantum annealing shows 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 values of the evaluation function 3, that is, The value of Eq. 15 in the quantum annealing is smaller than the corresponding value in the simulated annealing. The average value in each method is 4.22 × 10 8 and 5.62 × 10 8 , respectively, 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 , respectively, 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 value of the quantum annealing is 1.01 times that of the exact solution. Table 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 environment 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, respectively. To control the physical system given 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 above, 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 achieves input generation within 0.1 s. However, 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 achieved for control targets with small time constants, such as mechanical systems and electrical systems.
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 www.nature.com/scientificreports www.nature.com/scientificreports/ as digital signals 28 . We 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 , ( ) be the signal after quantization; that is, at each time 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 perceive sound with different sensitivities at different frequencies 29 , and its auditory characteristics are expressed as perception filters 30 . Among the available filters, we consider a linear time-invariant filter 31 : , and  ∈ × C n 1 are filter parameters. The quantization error after the signal passes through the perception filter is defined as and the evaluation function is defined as  www.nature.com/scientificreports www.nature.com/scientificreports/ Given perception filter P and audio signals {a(t)}, the problem of finding an input that minimizes the evaluation function is represented as an optimal control problem. Using Z-transform, Eq. 17 is expressed as the following state-space model: The correspondence between system 19 and system 1 is as follows: 1 2 and the correspondence between evaluation functions 18 and 3 is as follows: 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 following filter proposed in ref. 31  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: We adopt the following input sets: Comparison of the values obtained by evaluating Eq. 18 over the duration with three methods, that is, t 0 99 2 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 respectively, and the value of quantum annealing is 0.81 times that of the simulated annealing. The quantum annealing solutions have less variance; the www.nature.com/scientificreports www.nature.com/scientificreports/ variances of the simulated and quantum annealing are 1.94 and 1.78 × 10 −2 , respectively. When comparing the exact solution and the quantum annealing, the average value from the quantum annealing is as 2.1 times that of the exact solution, which is a conservative result when compared with the previous example. Table 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 even 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 evaluation 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 difference from the classical simulated annealing method. The latter employs an iterative 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 above, the library default value 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. To 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 evaluation function decreases as the iteration number is increased, and about one million steps are required to reach the level of the quantum annealing result. The average 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.
We 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 structure, in which closely connected 8-bit units are arranged vertically and horizontally 32 . For this reason, the variables of a given QUBO problem cannot be directly assigned to physical qubits. The method of converting the given graph structure to the chimera structure is called minor embedding, which is realized by expressing one logical variable with strongly coupled multiple physical qubits. The embedding problem itself is actually an NP-hard problem, and we used the API that D-Wave provides to perform embedding with heuristic algorithms 33 . Because multiple physical bits are sacrificed in minor embedding, the effective number of variables 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 variables 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   www.nature.com/scientificreports www.nature.com/scientificreports/ the original problem has a dense structure 34 . To check the sparseness of our formulated problem, we examine the value of all components of the QUBO matrix [J ij in Eq. 45] and plot them as histograms in Fig. 4. We 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 value of the QUBO matrix and normalize them such that . Many components have values close to zero, and qualitatively, 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 variables and the accuracy of the solution. D-Wave Systems Inc. recognizes this and has stated that they will adopt a new graph structure that increases the couplers per qubit in the next generation of machines 35 . Also, they plan to double the number of qubits in the next 2 years. Hence, the restriction on the problem size will soon be relaxed and the accuracy of solution will be improved in the future. We expect D-Wave'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 effective computational time of the quantum annealing compared with that of the simulated annealing, as shown in Tables 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 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
transformation of Mpc to QUBo. We first represent evaluation function 3 as a quadratic form for the input vector u(t). First, the dynamics of the system 1 is represented as with the following definitions:  Also, evaluation function 3 is expressed as by defining By substituting Eq. 28 into Eq. 34, we obtain the following quadratic form: where we assumed that Q is symmetrical and that c(t) is a constant independent of u(t).
The need to search the set  N makes it difficult to find the minima u(t) of the right side of Eq. 37. The set  N consists of M N elements (where M is the number of elements of the set  and N is the control horizon), and as N and M increase, the search becomes impractical. For this reason, when M and N are large, it is necessary to give up an exact solution and instead adopt an approximate solution. This study explores the quantum annealer as one of the approximation methods. However, because this method only accepts QUBO problems, we next need to convert Eq. 37 to a QUBO problem. For the set  defined by Eq. 2, we add the following assumptions: let M (the number of elements of sets  ) satisfy M = 2 L for some natural number 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 m ( 1, , ) j is a scaling parameter. For example, the set   in Eq. 14 and  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, respectively. The assumption that the input set  is expressed as Eq. 38 is made to express the evaluation function 3 in a QUBO form. We note that the expression in terms of Eq. 38 restricts the input set to the discrete values with a regular interval. Nevertheless, this is valid for various 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  , and its expression as a matrix is Eq. 37 can be written as a QUBO as follows: where ′ c t ( ) is a constant independent of b(t), and J(t) and h(t) are defined as 2 1 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.