{'fraction_non_alphanumeric': 0.04194254130922909, 'fraction_numerical': 0.014134377338937187, 'mean_word_length': 4.3626119172584135, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Hybrid Quantum and Classical Computing Algorithms for Mixed-Integer Programming

Quantum computing is emerging as a new computing resource that could be superior to conventional computing for certain classes of optimization problems. However, in principle, most existing approaches to quantum optimization are intended to solve unconstrained binary programming problems, while mixed-integer linear programming is of most interest in practice. We attempt to bridge the gap between the capability of quantum computing and real-world applications by developing a new approach for mixed-integer programming. The approach applies Benders decomposition to decompose the mixed-integer programming into binary programming and linear programming sub-problems, which are solved by a noisy intermediate-scale quantum processor and conventional processor, respectively. The algorithm is provably able to reach the optimal solution of the original mixed-integer programming problem. The algorithm is tested on a D-Wave 2000Q quantum processing unit and is shown to be effective for small-scaled test cases. We also test the algorithm on a mixed-integer programming inspired by power system applications. Many insights are drawn from the numerical results for both the capabilities and limitations of the proposed algorithm.


I. INTRODUCTION
Recent years have have witnessed exciting progress in the field of quantum computing (QC). Many QC algorithms have been developed to solve various types of problems, optimization being among them. For integer factorization, Shor's algorithm is exponentially faster than the best-known classical algorithms [1]. Grover's algorithm [2] provides a quadratic speedup in database search problems. The Harrow-Hassidim-Lloyd (HHL) algorithm [3] provides an exponential speedup over classical computing (CC) in solving linear equations. And quantum machine learning ( [4] and references therein ), built largely on the aforementioned QC algorithms, is also emerging as a popular research direction. However, many of these quantum algorithms with provable speedups are contingent upon the as-yet unrealized construction of fault-tolerant quantum processors [5]. The potential of QC is immense, but the classes of optimization problems that nearterm, noisy intermediate-scale quantum (NISQ) processors are able to address remain restricted [6]. This paper is intended to increase the impact of NISQ optimization by leveraging its capability with classical co-processing to solve mixed-integer programming (MIP) problems, which are of interest in a wide range of real-world applications.
Many practical control and optimization problems can be formulated as MIPs because the controls (or decision variables in the context of optimization) are a mix of discrete and continuous actions. MIPs occur in various applications such as the classical traveling salesman or routing problems [7], [8], supply chain planning [9], [10], energy market operation [11], and power system expansion planning [12]. Finding an optimal solution to a MIP in principle requires complete combinatorial enumeration, leading to exponential complexity. In addition to state-of-the-art branch-and-cut (or branch-and-bound) methods, various specialized algorithms have been developed to efficiently solve some special cases of MIP. This reflects the importance and challenges involved in MIP.
Quantum processors are able to directly represent exponentially many solutions on linearly scaling hardware resources through qubit state superposition and leverage a uniquely quantum phenomenon known as amplitude amplification in order to soften the hard exponential scaling of combinatorial enumeration. Some variants of the amplitude amplification mechanism are thought to be able to soften the timecomplexity scaling exponent for combinatorial optimization so effectively that they outclass the best possible classical heuristics, opening the door to a discussion of quantum advantage. However, most near-term quantum optimization heuristics rely on solving Quadratic Unconstrained Binary Optimization (QUBO) problem formulations. Any optimization problem for which there is an analytic expression can be reformulated as a QUBO, and many formal mathematical and computational problems have been recast in QUBO form, including the minimal dominating set formulation of optimal phasor measurement unit placement on the power grid [13], the Markov decision process [14], and many NP difficult problems [15]. With well established CC algorithms for convex optimization, decomposing a MIP problem into a series of one or more QUBO and convex optimization sub-problems, which are respectively solved by QC and CC, is a promising way to leverage QC capabilities. There are only a few works that engage in this hybrid QC-CC concept. In [16], an alternating direction method of multipliers (ADMM)-based hybrid QC-CC algorithm was developed to solve a class of MIP with constraints represented as equalities coupling binary and continuous variables. However, the ADMM-based methods can not guarantee the convergence due to the non-smooth nature of MIP. In [17], hybrid QC-CC algorithms were proposed for job-shop scheduling problems, cell formation manufacturing problems, and vehicle routing problems. The latter two problems are non-convex even without considering the discrete binary variables. Therefore, some heuristic methods instead of algorithms with exact ones were developed for them in [17].
In this paper, we propose applying Benders decomposition [18] for a hybrid QC-CC algorithm that is provably convergence for a class of MIP. The equivalence between the original MIP and a binary programming problem is derived. The binary programming sub-problem can involve a large number of constraints. In order to avoid this problem, we follow the Benders decomposition iterative procedure by only adding one constraint at a time to gradually tighten the optimization search space. With a sufficient number of back-and-forth iterations of solving the QUBO and linear programming (LP) (used to identify and add constraints to QUBO) sub-problems, the algorithm in principle can find an optimal solution. However, there are many challenges for the practical implementations. For small-scale MIPs, the algorithm can find the optimal solution with 2000Q D-Wave System, while for a more engaging formulation of MIP inspired by power systems applications, the algorithm can only reach a near global optimal solution. Limitations of reformulating inequality constraints to QUBOs and current quantum annealing hardware post some challenges for the practicability of the algorithm, which will be detailed in the numerical studies.
The rest of the paper is organized as follows. In section II, we begin with a simplified form of MIP in order to review the core ideas of Benders decomposition. Section III extends the results in section II for general MIP problems. Section IV presents our numerical results.
Notation: Let R and R + respectively denote the set of real and non-negative real numbers; N denotes the set of natural numbers. Given a set S, |S| is the cardinality of it. For a scalar a ∈ R, a is the largest integer less than a.

II. HYBRID ALGORITHM FOR SIMPLIFIED MIXED-INTEGER PROGRAMMING
We consider a special case of MIP in this section, shown in the following min y∈{0,1} ny ,z∈R nz where n y and n z are respectively the number of binary and continuous variables, f : {0, 1} ny → R is a quadratic function (can be convex or non-convex), b ∈ R m b , d ∈ R m d , m b and m d are the number of equality and inequality constraints, and A y , A z , C are in proper dimensions. We first decompose (1) into two subproblems such that one is a pure binary programming problem that can be solved by QC, and the other can be efficiently solved by conventional computers. Let us rewrite (1) in the following equivalent form where q(y) = min We can view (3) as a feasibility problem. In the context of optimization theory [19], primal optimization (3) has the associated dual problem formulation. Let α ∈ R m b , β ∈ R m d + , λ ∈ R nz + be the dual variables associated with constraints (1b), (1c) and z ∈ R nz + respectively. The dual problem of (3) is given as which is equivalent to the following Because Slater's condition holds for the linear programming (3), the duality gap between the primal problem, (3), and dual problem, (4), is zero. This implies that the optimal value of (4) should be zero if (3) is feasible. On the other hand, if (3) is infeasible, then the optimal value of (4) is infinity. The optimal value of (4) should be non-negative either way. Optimization (4) has the optimal value non-negative only if α A z + β C − λ = 0. Otherwise, there exists a z such that (4) has the optimal value −∞. We can then rewrite (4) in the following equivalent form The optimal value of (5) is either zero or infinity, which respectively indicates (3) is feasible or infeasible for the given y. The y such that (3) is feasible should make the optimum of (5) zero. Optimization (5) with y given has the optimum being zero if and only if where Γ is the set of extreme rays of the set To this end, we can conclude that (3) is feasible if and only if (6) holds. Replacing q(y) in (2) by (6) gives an equivalent formulation of (1), shown as follow min y∈{0,1} ny f (y), s.t. (6) holds.
Optimization (7) is a quadratic binary programming problem and in principle can be solved by QC. The potential issue is that the number of extreme rays, or |Γ|, may incur too many inequality constraints, which in turn would require a prohibitive number of qubits.
In the context of Benders decomposition, optimization (7) without constraint (6) is known as the master problem, and (6) collects the set of "cuts" (or constraints) that may gradually add to the master problem. If the optimal solution y 0 of the unconstrained (7) has (5) unbounded, it only indicates that y 0 is not a feasible solution, but it does not provide information on which cut to add. For convenience of determining which cut to add to the master problem, we define (8), which is the same as (5) except that bounds on α and β are added. max |α|≤α, |β|≤β,β≥0 where α ∈ R m b + and β ∈ R m d + . Both α and β are finite. The optimal solution of (8) is clearly bounded. In addition, the optimal solution of (8) should be one of the extreme rays in Γ. Repeatedly solving (8) leads to a sequence of cuts to add to unconstrained (7). Algorithm 1 summarizes the steps described above.
Find (α t , β t ) by solving (8) for the given y t 5: Conclude y t is the optimal solution 10: Solve (3) with given y t Find the optimal z 11: Set ξ = 1 12: end if 13: end while Note that in the worst case, one would need to run |Γ| iterations before finding the optimal solution. One alternative may be setting an upper bound on the number of iterations. If Algorithm 1 is unable to find the optimal solution within the maximal number of iterations, one can assume that the incumbent solution is close to the global optimum and use that to assist the conventional computer to solve (1).

III. GENERAL MIXED-INTEGER PROGRAMMING
We now generalize (1) by adding terms linear in the continuous variables to the objective function. The equalities in (1b) are replaced by inequalities.
Following the steps in the last section, we can rewrite optimization (9) in the following min y∈{0,1} ny ,q∈R where Γ r and Γ p respectively are the set of extreme rays and extreme points of the set (7), (10) is not a binary programming because it has one continuous scalar variable q. This is a trade-off that we need to make when we consider a more general class of MIP. We follow a similar Benders decomposition process to solve (10). Consider the following dual optimization (similar logic to (8) The optimal solution of (11) should be in the set Γ r or Γ p . Similar to the previous section, repeatedly solving (11) leads to a sequence of cuts to add to the unconstrained (10). However, any added cut depends on whether the solution of (11), (α , β ), is in Γ r or Γ p . We use the following simple logic to determine which set the solution is in where α (k) is the k th entry of α and similar notation applies for other vector variables. If an entry of α (or β ) reaches the upper bound, then it is in Γ r if α and β are chosen large enough, because the optimal points in Γ r are unbounded if we do not add the bounds α and β. On the contrary, the extreme points are in the bounded set and constitute a strict subset of {(α, β) | |α| < α, |β| < β}. The logic shown in (12) is accurate as long as we have α and β large enough.
Algorithm 2 summarizes the steps for finding the optimal point of (9). The main difference between Algorithm 1 and 2 lies in the form of the cuts. The cuts added in Algorithm 1 are known as "feasibility cuts" because those cuts are added for the purpose of enforcing a feasible solution for the original optimization (1). The intuition behind the feasibility cut is that the continuous variables, z, only affect the feasibility of the optimization because they only appear in the constraint of (1). Optimization (9) that Algorithm 2 targets to solve has z simultaneously in the constraint set and objective function. As a result, in the course of iterative steps of Algorithm 2, we may add feasibility cuts (similar to Algorithm 1) or "optimality cuts" to reflect appearance of z in the constraint set and the objective function. Feasibility and optimality cuts are very similar, the only difference between them being the dependency on the variable q as shown in line 10 and 13 of Algorithm 2.

IV. REFORMULATION OF BINARY PROGRAMMING
Both Algorithms 1 and 2 pose a binary programming problem to be solved via quantum computing. The binary programming problem becomes constrained when cuts are added during each of the while loop iterations. However, because most existing quantum optimization heuristics are designed only to accept QUBO problems, we need to recast hard constraints as soft, penalty functions added to the Algorithm 2 Solve the optimization (y t , p t ) = argmin y∈{0,1} ny ,q∈R f (y) + γq, Find (α t , β t ) by solving (11) for the given y t

5:
Extract the optimal value of the primal problem: Extract the optimal value of (11) given as if s pt < s dt then 9: if (α t , β t ) ∈ Γ r then 10:  Conclude (y t , q t ) is the optimal solution 19: Find the optimal z with the given y t 20: Set ξ = 1 21: end if 22: end while objective function. We will focus on a penalization method for Algorithm 1 for simplicity. The logic of the penalization method for Algorithm 2 is very similar.
In the following we write down in explicit form the optimization in line 3 of Algorithm 1 whereÂ t ∈ R t×ny andb t ∈ R t . The number of rows for bothÂ t andb t increases by one whenever t increases, which can be equivalently interpreted as an added cut. We consider penalizing each cut by a quadratic objective function shown in the following where a k ∈ R ny is the k th row of the matrixÂ t , w t ∈ R + and w x k ∈ R + are penalty weights, and x k (j) is the j th element of x k . The ideal penalty weights w t and w x k have the following properties: (i) for any y such that the inequalities in (13) hold, there exists feasible x k , k = 1, · · · , t, such that (14b) is zero; (ii) all solutions violating the constraint (13) are not the optimal y for (14). If (i) and (ii) hold true, then (14) and (13) are equivalent in the sense that they have the same optimal y. We first analyze condition (i). If a k andb t (k) have integer values, then we can always choose w x k = 1 for all k = 1, · · · , t to ensure that for any y feasible for (13), there exists a x k such that a k y −b t (k)+w x k nx k j=1 2 (j−1) x k (j) = 0 as long as n x k is sufficiently large. When the values of a k andb t (k) are not integer, we choose a small value for w x k so that there exists a x k such that a k y −b t (k)+w x k nx k j=1 2 (j−1) x k (j) = 0 for any feasible y. In this case, we have a trade-off between the accuracy and the number of qubits. If w x k is chosen very small for accuracy purposes, a large n x k is necessary, which implies a large number of qubits. Currently we choose w x k in heuristic ways, but a rigorous way of choosing w x k may be attainable. For condition (ii), there always exists a large enough w t so that (ii) holds. However, in the context of QC, if w t is overwhelmingly large compared to the actual objective function f (y), then noise could impair the QC from reaching the ground state, i.e., optimal solution. As a result, QC may reach a solution that satisfies the penalty constraints while f (y) is not minimized. Again we compromise by choosing w t in heuristic ways.
V. NUMERICAL STUDIES There are many potential applications that can be formulated as (1) or (9) so that the proposed algorithm can become useful. For example, planning of transmission systems, cooptimization of grid and transportation. In this section, we consider two different problem formulations for the algorithm validation purpose. One is a small-scale MIP and the other is a more engaging and realistic MIP for a unit commitment problem for power systems. For the small-scale example, we implement the proposed algorithm on an actual D-Wave system while for the unit commitment problem, we use a simulated annealing sampler as an alternative to the D-Wave system because of current size and connectivity limitations.

A. Small-Scale Test Case for D-Wave System
We validate the proposed algorithm by running small-scale MIP optimizations on a 2000Q D-Wave quantum processing unit (QPU). The 2000Q D-Wave QPU is constructed from a 2,048 qubit (16, 4)-chimera hardware graph topology consisting of a 16 × 16 array of 8−qubit chimera unit cells, each of which is organized as a K 4,4 complete bipartite graph. More details can be found in [20]. We accessed the D-Wave system by their Leap 2 project.
The objective function for the binary variables takes the linear form f (y) = g y y.
Both C and d are zero in this test case. Using our proposed decomposition approach, the original MIP is decomposed into a binary programming and LP, which are solved in the alternating sequence described in Algorithm 1. The D-Wave system solves binary programming and classical computing solves LP. We set the penalty strength (w t ) at 50, granularity constant (w x k ) at 0.01, and n x k = 10 for every t and k. In addition to those parameters, one needs to define the chain strength for the D-Wave system, which gauges the degree to which multiple physical qubits act as a single logical qubit. On the D-Wave system, the optimization problem is embedded onto its physical Chimera hardware graph. The embedding process may group multiple qubits together to represent one logical variable in the QUBO. Chain strength controls the coupling strength between the grouped qubits. In our simulation setup, we set the chain strength at 75 for all t to ensure it outweighs the penalty strength. Under this numerical setup, Algorithm 1 finds the optimal solution in two iterations. We next simulate a more general MIP with an additional term on the objective function shown in (9) with g = [20, 2, 1, 4] . The remaining constraints are the same as (15). We use 30 qubits to approximate the continuous scalar q for QC (Line 3 in Algorithm 2), shown in the following where w q = 0.01 and x(j) ∈ {0, 1} for all j = 1, · · · , 30. We set penalty strength and chain strength respectively at 2 and 3 in this case. Algorithm 2 reaches the optimal solution in four iterations. Note that the penalty and chain strengths are much smaller in our second test case compared to the first one. The reason is subtle and is related to the upper bounds on α in (8) and (12) 1 . For Algorithm 1, we set every entry of α at 1 and can always decide whether to add a cut or not by checking whether the objective value of (8) is strictly positive. However, for Algorithm 2, we need to differentiate feasibility and optimality cuts, in this case, entries of α are set at some arbitrary large values. The large α are translated to large coefficients for the feasibility cut defined in Line 10 of Algorithm 2. To prevent the penalty functions for the feasibility cuts from outweighing the objective function, we normalize all the cuts in Algorithm 2 by dividing their coefficients by the square root of the optimal value of (12). Even with the presence of the normalization, the coefficients of these cuts are still relatively large compared to our first test case, so we choose much smaller penalty and chain strengths.

B. Unit Commitment Test Case
In this section, we test Algorithm 1 on a unit commitment problem for power systems applications. The formulation is similar to standard optimal power flow (OPF) with the main difference being in the additional option for every generator i ∈ N G to switch "on/off" based on the total power demand, where N G is the set of generators. Those switching decisions 1 No β is defined for the numerical examples because we do not impose any linear inequality constraint (1c) or (9c) in the numerical study. naturally introduce binary variables. The other main difference from the standard OPF is that we consider planning for a certain time horizon, T hours, with a sampling time of one hour (τ = 1, · · · , T ). The switching between consecutive timestamps is constrained by the minimal locked time of the generators. That is, every generator has to stay in "on" (or "off") state before it is allowed to change state. Without loss of generality, we assume that the minimal "on" ("off") time, T on ∈ N (T of f ∈ N), is strictly larger than one. Otherwise, the minimal "on" ("off") time is not binding for the unit commitment problem. Capturing the locked time constraints introduces additional binary variables. Such a unit commitment problem can be either formulated in the form of (1) or (9).
1) Formulation for a 6-Bus System: Given a 6-bus system shown in Figure 1, our goal is to minimize the generation cost subject to the OPF constraints. In the following we go through all the OPF constraints and reformulate them to the ones for the unit commitment problem (with additional binary decisions). First of all, every generator has the upper and lower bounds on the generation limits given as for all τ = 1, · · · , T , where p i is the active power generation of generator i, P i ∈ R + and P i ∈ R + respectively are the lower and upper bounds of the generation. Define y G,i (·) ∈ {0, 1} for every i ∈ N G such that y G,i (·) = 0 corresponds to generator i being in the "on" state (y G,i (·) = 1 for the "off" state). The constraint above is reformulated as for all τ = 1, · · · , T , where M is an arbitrary scalar larger than max{P 1 , · · · , P | N G | }. Equation (18a) and (18b) applies well-known Big M method so that (17) is active only when y G,i (·) = 0. Equation (18c) ensures that p i (τ ) = 0 if y G,i (τ ) = 1. We next apply similar methods to reformulate the ramp constraints where P r,i is the maximal ramp rate. The reformulated (19) is given as for all τ = 2, · · · , T , where y on,i (·) ∈ {0, 1} with y on,i (τ − 1) = 1 indicating that generator i should stay in the "on" state at time τ , y of f,i (·) ∈ {0, 1} with y of f,i (τ − 1) = 1 indicating that generator i should stay in the "off" state at time τ . Equation (20) has the ramp constraint (19) active if y on,i (τ −1) = y of f,i (τ −1) = 0. We will next post constraints such that y G,i (τ ) = y G,i (τ − 1) = 0 if y on,i (τ − 1) = 1 and y G,i (τ ) = y G,i (τ − 1) = 1 if y of f,i (τ − 1) = 1, so no need to have the ramp constraints when y on,i (τ − 1) = 1 or y of f,i (τ − 1) = 1. For every i ∈ N G , y on,i (·), y of f,i (·), y G,i (·) are coupled by the following constraints: at the "on" state if y on,i (τ − 1) = 1. Similar logic applies to (21b). Equations (21c) and (21d) respectively have y on,i (τ ) = 1 and y of f,i (τ ) = 1 when generator i switches at time τ . Equations (21e) and (21f) respectively keep y on,i (τ ) and y of f,i (τ ) stay at 1 before the associated generator is allowed to switch. All the constraints above are about nodal power injections. The correlations between the buses (nodes) are defined by the graph of the network and the power flow equations. Define N and E respectively as the set of buses (nodes) and transmission lines (edges). Constraints associated with the line power flows are given as for all τ = 1, · · · , T , where p ik (·) is the active line power flow for line {i, k} ∈ E, P ik (·) is the upper bound of it, A is the incidence matrix associated with the power network graph (with an arbitrary direction for every line), B is a diagonal matrix with each of the diagonal element being associated with the impedance of a line, P (·) ∈ R | N | collects the nodal active power injections and P L (·) ∈ R | E | collects p ik for all {i, k} ∈ E, Θ(·) ∈ R | N | collects the phase angle of all the buses. Equation (22c) corresponds to DC power flow formulation which assumes purely inductive lines so the power flow equations can be approximated in the linear form. More details about the DC power flow can be found in [21]. Because our formulation requires all the continuous variables to be nonnegative, we change all the constraints in (22) to ensure the positiveness of the variables: where P L is the collection of P ik for all {i, k} ∈ E. All the constraints for the unit commitment problems have been included in (18), (20), (21) and (23). Our goal is to minimize the generation cost given as τ =1,··· ,T i∈N G c i p i (τ ), where c i ∈ R + for all i ∈ N G . To this end, the complete formulation of the unit commitment problem is given as s.t. (18), (20), (21) and (23), where y and z are the collections of all the binary and continuous variables, respectively. Unit commitment problem (24) fits directly into the formulation of (9). In fact, one can also choose to reformulate (24) to (1) by compromising on the size of the problem. Specifically, we can introduce slack continuous variables to transform all the inequality constraints to equality constraints. The objective function of the continuous variables is approximated by binary variables with the method similar to (16). In this moderately sized test case, the additional slack variables do not pose major a challenge, so we choose to reformulate (24) to (1) for both simpler coding and demonstration purposes.     2) Results with Binary Programming Solver: We first test Algorithm 1 against the unit commitment problem formulated in the previous subsection with the numerical setup detailed in Tables I-IV. Both the binary programming and convex optimization of Algorithm 1 are solved by PuLP toolkit for python. The upper bound of every entry of α, α, is set at 1000. We compare two scenarios with different planning time horizons (2 and 3 hours). As shown in Figure 2, Algorithm 1 finds the optimal solutions in both scenarios. The number of iterations it took to find the optimal solution grows with the problem size. The scenario with a 3 hour planning horizon has roughly 1.5 times more of the constraints and decision variables compared to the one with 2 hour, which results in roughly 100 more iterations (cuts) needed to reach the optimal solution. The simulation results give a rough picture on how the number of cuts grow with respect to the problem size. Figure 3 shows that the optimal value of (8) decays gradually to zero as opposed to the abrupt jump to zero in Figure 2. The optimal value of (8) serves as a good indicator of how close Algorithm 1 is to the optimum.
3) Simulated Annealing Sampler: Simulated annealing algorithm is based on the technique of cooling metal from a high temperature to improve its structure (annealing) [22]. We choose a simulated annealing sampler for the unit commitment problem largely due to hardware constraints of current D-Wave systems. Results in the last subsection show that Algorithm 1 needs to add 104 constraints (cuts) before reaching the optimal solution. Each additional inequality constraint requires introducing concomitant slack binary variables as discussed in Section IV. The number of variables added depends on the desired floating point precision. Even if we choose a relatively large floating point precision at 1 (granularity constant w x k = 1 in (14)), around 15 additional binary variables are needed for each new constraint. The binary program in Algorithm 1 will therefore grow to a scale of around 2000 binary variables in the last few iterations. Such binary programs are not solvable by the D-Wave 2000Q system that we have access to because embedding the binary program in the D-Wave hardware graph requires physical qubit overhead that results in far more than the 2000 logical qubits needed for the 2000 binary variable problem. The simulated annealing sampler follows similar heuristics to the D-Wave system although it is run on classical computers and relies on "thermal" hopping rather than quantum tunneling to traverse a problem's energy landscape. We therefore view simulated annealing as a good proxy for an actual D-Wave processor for algorithm validation purposes. The view aligns with the literature well [23], [24].
We compare the performance of Algorithm 1 with two different implementations: one uses PuLP to solve the binary programs involved in Algorithm 1 as shown in the last subsection; the other uses simulated the annealing sampler. Both implementations are tested on the unit commitment problem with T = 2. For simulated annealing, we set the penalty strength at 2 (one need not set the chain strength because it is only required for an actual D-Wave system). In addition to the penalty strength, other parameters, including number of reads (N read ) and number of sweeps (N sweep ), have a significant impact on simulated annealing results. N read is the number of runs of the simulated annealing algorithm (analogous to the number of runs of quantum annealing). N sweep is the number of fixed temperature sampling sweeps for each run of the simulated annealing algorithm. As the size of the optimization problem grows, one requires larger values of those parameters for similar precision of the solutions. Therefore, we take the approach of gradually increasing N read and N sweep as the number of iterations increases in Algorithm 1. Specifically, we set N read = 80 and N sweep = 5000 at t = 0 and allow them to grow exponentially with the rate of 0.02, given as where t is the iteration number of the algorithm.
As shown in Figure 4, PuLP obtains the optimal solution in 104 steps while simulated annealing fails to find the optimum with 120 steps. The optimal solution that PuLP found has G1 and G2 in the "on" state and G3 in the "off" state for all time. The solution that simulated annealing reached has all the generators in the "on" states. Both solutions are feasible with a marginal difference between the optimal values, which are respectively 5986 and 6021. Simulated annealing version still found a reasonably good solution with less than 1% off the optimum. A similar conclusion can be drawn by observing the evolution of the optimal values of (8) in Figures 5 and 6, which show that the simulated annealing approaches the optimal Fig. 4: Convergence to the optimal point solution as the number of iterations increase. The main reason that simulated annealing fails to find the optimal solution is the violation of the constraints (cuts). Recall that in simulated annealing, we penalize the inequality constraints by quadratic functions given by (14b). We quantify the violations by the following The violations as a function of t shown in Figure 7 are within a small margin but fail to drop further to the exact optimal solution as the PuLP case. This reflects one of the challenges as the number of constraints (cuts) grows. Adding a constraint in the form of an extra quadratic penalty function effectively diminishes the weight share of the previous penalty functions and the original objective function. This is not an issue from the analytical standpoint because as long as the penalty weights are non-zero, the QUBO always has the same optima as the original constrained binary programming as discussed in Section IV. However, numerically, especially for simulated annealing or QPU that rely heavily on the probabilistic distribution of samples, it is challenging to find penalty weights that guarantee distinction between the global optima from near optimal solutions. We have many trials of various penalty weights, but they all generate similar results in the sense that reaching near global optima in the last few iterations. Using the metric in equation (8) (illustrated in Figure 6) for early termination on a near optimum is a reasonable way to implement Algorithms 1 and 2.
We last want to comment on the runtime for the simulated annealing. Each emulation of a result from quantum annealer needs N read · N sweep number of samples, which translates to the order of 10 7 (or more) samples from iteration 100 onward. Therefore, the runtime for generating Figure 7 is around 37 hours (did not time exactly) on a personal laptop with 3.5GHz CPU and 8GB RAM. This is a lot longer than many existing solvers. However, simulated annealing sampler only serves the purpose of generating similar results of quantum annealer. We expect the actual runtime on quantum annealer to be a lot faster, and advantages over MIP standard solvers could be observed in large-scale MIP in future endeavors.

VI. CONCLUSION
In this paper, we proposed applying Benders decomposition for a hybrid QC-CC algorithm. The algorithms presented take an iterative approach that gradually closes the gap the optimal solution. Small-scale simulation results demonstrate that the algorithms can find the global optimum with proper selection of D-Wave quantum annealer parameters. For a more realistic MIP that is inspired by unit commitment applications, we identify some limitations of the effectiveness of the algorithm mainly due to the random nature of the annealing mechanism and the requirement that hard constraints be softened in such approaches. In addition to the challenges of the random nature and methods of softening the constraints, another important improvement need of the proposed algorithm is characterizing the bound of the number of iterations required for finding a reasonably good solution. Our near future work will focus more on some tangible tasks, including rigorous methods for identifying penalty weights for robust numerical performance of the hybrid algorithm. We will also look towards additional real-world MIP applications and refine the algorithms for those applications.