Hybrid Quantum Benders’ Decomposition For Mixed-integer Linear Programming

MILP has been widely adopted optimization problems that include but are not limited to communication and networks. In the field of communication, for instance, it has been employed to the resource allocation in wavelength division multiple access (WDMA) [14], multi-access edge computing (MEC) [15], offloading for fog computing [16], etc. Furthermore, it plays an essential role in the field of network research, such as minimizing the network delay [17], scheduling virtual network re-configurations [18], and improving the supply chain network [19]. Besides civilian applications, it is no doubt that MILP has been adopted for military purposes such as assisting military pilot training [20]. These demonstrate that MILP is a powerful tool for either classical problems or future popular topics such as communication networks and resource allocation. Consider a MILP model as follows:

Here $x$ is a vector with binary variables and $A$ is its corresponding coefficient matrix in constraints. $y$ is a vector with non-negative continuous variables and $G$ is its corresponding coefficient matrix in constraints. $c^{\intercal}$ is the coefficient vector of variable $x$ in the objective function and $h^{\intercal}$ is the coefficient vector of variable $y$ in the objective function.

We first provide a brief introduction of the Benders’ Decomposition. In the MILP problem (1), for each possible choice of $\bar{x}\in X$, we find the best choice for $y$ by solving a linear program. So we regard $y$ as a function of $x$. Then we replace the contribution of $y$ to the objective with a scalar variable representing the value of the best choice for a given $\bar{x}$. We start with a crude approximation to the contribution of $y$ and then generate a sequence of dual solutions to tighten up the approximation. Hence, the original MILP problem in (1) can be written equivalently to a master problem as follows.

We denote $K$ and $J$ as the extreme points $u^{k}$ and extreme rays $r^{k}$ of the dual polyhedron $Q=\left\{u\in\mathbb{R}^{m}_{+}\mid G^{\intercal}u\geq h\right\}$ generated by the linear programming’s duality, which is called the subproblem. The subproblem is as follows.

In the subproblem, if the inner product between $(b-Ax)$ and any dual ray $r^{j^{\prime}}$ is negative then $z_{LP}(x)=-\infty$. Equivalently, in this situation, the dual problem of Problem (3) is infeasible. Then, $x$ does not allow a feasible solution to the original mixed-integer problem (1). Thus, we have a new feasibility cut, i.e.,

If $x$ satisfies (4), then we yield an extreme point $u^{k^{\prime}}$ and the value of $z_{LP}(x)$ is given by

Thus, problem (1) is written equivalently as problem (2) by denoting $z_{LP}(x)$ as a continuous real number variable $t$. Equation (5) posts a new optimality cut in master problem. In addition to the extreme points and extreme rays, we use $\hat{K}$ and $\hat{J}$ to denote the current known extreme points and extreme rays of $Q$, respectively, where $\hat{K}\subseteq K$ and $\hat{J}\subseteq J$.

Quantum annealers are able to solve the optimization problem in a QUBO formulation. To leverage state-of-art quantum annealers provided by D-Wave, the ILP problem has to be converted to the corresponding QUBO formulation. The definition of QUBO are (6) and (7). Let $f$:$\left\{0,1\right\}^{n}\rightarrow\mathbb{R}$ be a quadratic polynomial over binary variables $x$ of length $n$ and $q_{ij}$ is corresponding cost coefficient of $x_{i}x_{j}$. We get equation (6)

The QUBO problem consists of finding a binary vector $\mathbf{x}^{*}$ that is minimal with respect to $f$ among all other binary vectors, namely,

Here $\mathbf{x}$ is a vector of binary variables where the length of $n$, and $\mathbf{Q}$ is either an upper-diagonal matrix or a symmetric matrix. Since (6) is an unconstrained optimization model, we need to reformulate our constrained ILP as unconstrained QUBO by using penalties. Next, we get the optimal solution by finding the best penalty coefficients of the constraints. The principles of transforming classical constraints to their equivalent penalties are displayed in the TABLE I. Here $x_{1}$, $x_{2}$ and $x_{3}$ are binary variables. $s_{i}$ is a binary slack variable. $a_{l}$ is the coefficient for the corresponding slack variable. $b$ is a constant. $P$ is a user-defined penalty coefficient.

Now consider problem (2), the initial binary decision variables make up the vector $\mathbf{x}\in X$ with length of $n$. $\bar{t}\in\mathbb{R}$. In order to reformulate the master problem into the QUBO formulation, we needs to represent the continuous variable using binary bits. We use a binary vector $\mathbf{w}$ with length of $M$ bits to replace continuous variable $\bar{t}$. In general, $\bar{t}$ requires the binary numeral system assigning $M$ bits to replace continuous variable $\bar{t}$. Then we can recover the $\bar{t}$ by

In (8), $\bar{m}_{+}$ is the number of bits that assigned to represent positive integer part, $\underline{m}$ is the number of bits assigned to represent the positive decimal part, and $\bar{m}_{-}$ is the number of bits that represent the negative integer part.

Let’s reconsider problem (2) and replace $\bar{t}$ in (8). The new problem is only depend on binary decision variables $\mathbf{x}$ and $\mathbf{w}$. Then, the new master problem is

Then we create a new binary decision variable collection $\mathbf{x}^{\prime}=\left\{\mathbf{w},\mathbf{x}\right\}$ to set up our QUBO matrix, where the length of the binary vector $\mathbf{x}^{\prime}$ is $n+M$. According to the rule of setting up QUBO, we will transfer our objective and constraints one by one to be a symmetric matrix.

We deploy the D-Wave solver in our proposed heuristic algorithm to solve the model. In addition, we need to carefully tune the penalties for a decent QUBO model. An extremely large penalty may cause quantum annealer malfunctioning since it will explode the coefficients. Similarly, a soft penalty may make quantum annealer ignore the corresponding constraints. As long as the penalty is well-tuned, the quantum solver will give the right answer with a relatively high probability. Therefore, we summarize our proposed hybrid quantum Benders’ decomposition algorithm as Algorithm 1.

In our simulation, we consider a simple problem to test our proposed quantum algorithm, where $\mathbf{x}\in X$, $\mathbf{x}\in\left\{0,1\right\}^{2}$, $\mathbf{y}\in Y$, $\mathbf{y}\geq 0$

Fig. 2 shows how our proposed hybrid quantum Benders’ decomposition algorithm adds optimality and feasibility cuts to get the solution and corresponding $\bar{t}$. In the example, the algorithm takes 4 rounds to let $t$ converge. Four cuts are added iteratively in the space and tighten the feasible region. The algorithm can find the optimal solution in each round. Therefore, Fig. 2 demonstrates that our algorithm is reliable and efficient. In Fig. LABEL:fig:tutl, the hidden dashed line denotes that the lower bound in the corresponding round is negative infinity. As we can see in the graph, the upper bound and lower bound converge. Our algorithm only takes two rounds to find the non-negative infinity lower bound. This result proves that our proposed algorithm is mathematically consistent with the classic Benders’ decomposition algorithm. In other words, as long as the classic Benders’ decomposition algorithm can solve it, our algorithm can at least achieve the same.