Homogeneous Formulation of Convex Quadratic Programs for Infeasibility Detection

Optimization algorithms are widely used in the control for constrained dynamic systems. In particular, the solution of convex Quadratic Programs (QPs) is a key ingredient of optimization algorithms. Convex QPs arise for example in the Model Predictive Control of linear systems, switched linear systems and nonlinear systems [1, 2]. The past two decades have witnessed the development of a number of algorithms for the solution of convex QPs arising in the context of optimization-based control. Initial work on QP algorithms were based on solving a system of equations representing the optimality conditions such as Interior Point Method (IPM) for QP [3], active-set methods [4, 5] and more recently, IPM for Second Order Cone Programs (SOCPs) [6]. In the last decade there has been interest in the development of first-order approaches such as gradient projection [7], dual gradient projection [8, 9], and splitting methods [10, 11, 12, 13] and iterative second-order approaches [14]. Recently, there has been work on employing these methods within Branch & Bound (B&B) algorithms for the solution of Mixed Integer Quadratic Programs (MIQPs) that arise in MPC for switched systems [15, 16, 17, 18].

A critical feature required of QP solvers for real-time applications is the ability to detect infeasibility and provide a graceful handling of the same. Infeasibility handling increases in importance when solving MIQPs in a B&B setting as it is expected that QPs resulting from the fixing of a subset of binary variables will likely be infeasible. Detecting such infeasible QPs and pruning the search tree in B&B is essential for computational efficiency of the MIQP solver.

Infeasibility of a QP is determined by the identification of a ray that is feasible for the dual and makes the dual objective unbounded. First-order primal-dual algorithms are capable of detecting infeasibility. Recent results on detecting infeasibility using first-order algorithms include [19, 20, 21]. A dual QP formulation, if it can be constructed explicitly, will allow the determination of a ray of infeasibility. However, such an explicit dual construction requires strong convexity of the QP. Active-set methods typically employ a feasibility phase to first determine a feasible point and then proceed to an optimization phase which produces the optimal solution. Failure of the feasibility phase allows to pronounce a QP as being infeasible [22]. The work performed in this phase is identical to that required for the optimization phase. When employing inexact linear algebra the first-phase is typically not employed and as a consequence the ability of such approaches to certify infeasibility is not clear. IPM algorithms, on the other hand, are not all capable of producing a certificate of infeasibility [23, 24]. IPMs based on Homogeneous Self-Dual (HSD) embedding [25, 26] are known to produce a certificate of infeasibility for Linear Programs (LPs). This technique has been implemented in MOSEK [27] and has also been extended to handle SOCPs. To detect infeasible QPs with IPMs, the QPs must first be formulated as a SOCP to which the HSD embedding is then applied. The ECOS solver [6] does exactly this for the solution of QPs arising in MPC and can detect infeasibility.

We present a novel embedding of the QP into a space that is one dimension higher. The embedding is produced by introducing an additional variable that is nonnegative and multiplies the affine terms in the objective and constraints. The objective contribution from the introduced variable is determined so the resulting Homogeneous Quadratic Program (HQP) satisfies the following properties: (i) HQP is always feasible, (ii) an optimal solution to the QP can always be recovered from the optimal solution to the HQP if the introduced variable is positive; and (iii) QP can be certified to be infeasible if the introduced variable is zero. More importantly, the HQP formulation can be presented to any QP solver to determine the optimal solution or infer infeasibility in a robust manner.

In this paper, we focus on applying an Infeasible IPM (IIPM) to solve the HQP. We are able to derive polynomial complexity of a standard IIPM when applied to HQP by adapting the analysis of IIPM for linear programs [23, Chapter 6]. Thus, we have an IIPM with polynomial iteration complexity for inferring optimality or infeasibility of QP without having to solve SOCPs.

The paper is organized as follows. §2 presents the QP and the assumptions. §3 presents the embedding of QP into a HQP and shows the equivalence of the formulations. An IIPM for solving HQP and the polynomial iteration complexity are described in §4. We present results in §5 when using the HQP formulation with IPOPT. Conclusions and directions for future work are provided in §6.

Notation. The set of reals is denoted by $\mathbb{R}$ and the set of vectors of dimension $n$ by $\mathbb{R}^{n}$. The set of $n\times n$ symmetric matrices is denoted by $\mathbb{S}^{n}$. For a matrix $A\in\mathbb{S}^{n}$, the notation $A\succcurlyeq(\succ)0$ denotes that A is positive semidefinite (definite). Given a vector $u\in\mathbb{R}^{n}$, $\textnormal{Diag}(u)$ denotes a diagonal matrix with the elements of the vector on the diagonal. The notation $I_{n}$ denotes the identity matrix of size $n$.

Consider the QP in the form

where $C\in\mathbb{S}^{n}$, $c\in\mathbb{R}^{n}$, $E\in\mathbb{R}^{m\times n}$, and $f\in\mathbb{R}^{m}$. Note that any QP formulation with finite lower or upper bound on each of the variables can be put into the form in (1) through a linear transformation of variables.

We make the following assumptions on QP (1).

Assumption 1 is not restrictive and can be easily satisfied by removing dependent rows if necessary. Assumption 2 requires that $C$ be positive definite when projected onto $Z$. This assumption readily holds for MPC formulations [19] and also for certain spectral relaxations of nonconvex MIQPs [28, 29]. Assumption 2 implies that the optimal solution to QP (1) is unique whenever QP (1) is feasible. Note that the assumptions do not preclude the infeasibility of QP (1). The basis $Z$ is typically computed using a QR factorization of $E^{T}$ [30].

We conclude this section with a statement of the conditions for optimality and infeasibility of QP (1).

A point $y^{\star}$ minimizes QP (1) if there exist multipliers $\nu^{\star}\in\mathbb{R}^{m}$, $\xi^{\star}\geq 0\in\mathbb{R}^{n}$, satisfying the first-order optimality conditions

The conditions (2) are necessary and sufficient for a minimizer of QP (1) under Assumption 2.

The QP (1) is infeasible if there exist $\nu^{\circ}\in\mathbb{R}^{m}$, $\xi^{\circ}\in\mathbb{R}^{n}$ satisfying

The conditions in (3) are obtained from applying a Theorem of the Alternative [31] to the constraints in (1b)-(1c).

Consider embedding QP (1) into a Homogeneous Quadratic Program (HQP) as

where $\tau\in\mathbb{R}$ is an additional nonnegative variable and $\theta>0$ is a parameter. Observe that the objective (4a) is obtained by multiplying the linear term $c^{T}y$ in (1a) with $\tau$ and appending with the term $(\theta/2)(\tau^{2}-\tau)$. The equality constraints (4b) is obtained by multiplying the right-hand side of (1b) with $\tau$. §3.1 provides conditions on $\theta$ that ensures HQP is convex. §3.2 shows how the parameter $\theta$ can be computed efficiently. Finally, §3.3 shows that solving HQP allows to recover a solution to QP or declare infeasibility.

For the remainder of the paper, we assume that the HQP has the form

where $x=\left(\begin{smallmatrix}y\\ \tau\end{smallmatrix}\right)$, $Q=\left(\begin{smallmatrix}C&c^{T}\\ c&\theta\end{smallmatrix}\right)$, $q=\left(\begin{smallmatrix}0\\ -\theta\end{smallmatrix}\right)$, $A=\left(\begin{smallmatrix}E&-f\end{smallmatrix}\right)$. We assume the following for the HQP (17).

The above HQP can be identified with original QP (1) when $f=0$ and the Assumptions 1-2 imply satisfaction of Assumptions 3-4. If $f\neq 0$ §3 shows how to cast (1) satisfying Assumptions 1-2 into a HQP of the form in (17) satisfying Assumptions 3-4.

In the rest of the section, we show how a standard IIPM for LPs can be extended to HQP (17). The IIPM also enjoys polynomial iteration complexity.

We present preliminary numerical experiments with the HQP formulation using the IPM solver IPOPT [32]. IPOPT  implements the standard IPM for QPs i.e. the algorithm cannot determine a certificate of infeasibility. In our experiments, we generated random instances of QP (1) with a single equality constraint with $C=I_{n}$, $c=e$, $E\in\mathbb{R}^{1\times n}$, $f=-1$. The coefficients in the single equality are restricted to be nonnegative and generated at random. It is easy to deduce that there exists no $y\geq 0$ such that $Ey<0$ for $E\geq 0$. Hence, the instances are all infeasible. We presented the QP formulation directly to IPOPT. The HQP formulation  (4) is derived for each such random instance and is also presented for solution to IPOPT. We generated 10 random instances of different sizes $n\in\{10,25,50\}$.

In the first set of experiments, we set the IPOPT  option mehrotra_algorithm=no. In this case, the predictor-corrector algorithm of Mehrotra is disabled. When solving the QP formulation, IPOPT  always terminates after 8-10 iterations with the indication that the problem might be infeasible. When solving the HQP formulation, IPOPT  always terminates with an optimal solution of $0$ as specified in Theorem 2 in about 16-20 iterations.

In the second set of experiments, we enabled the predictor-corrector algorithm in IPOPT  by setting mehrotra_algorithm=yes. In this setting IPOPT  hits the iteration limit on every single instance when solving the QP formulation. The dual infeasibility blows up to infinity in every instance. When solving the HQP formulation IPOPT  always terminates with an optimal solution of $0$ in about 8-10 iterations.

In the third set of experiments, we set $f=1$. In this case the QP instances are always feasible. We verified that both the QP and HQP formulations solved the problem to optimality. Further, the solutions to QP can be recovered from the HQP could be recovered according to Theorem 1. Thre is no significant difference in the number of iterations taken for convergence using either formulation. However, the computational time per iteration for HQP is higher than that of QP formulation. This is likely due to the matrix in (21) being dense. This can be rectified using a tailored Schur-complement-based implementation for the step computation. We will investigate this in a future work.

We presented a homogeneous formulation of QP that allows for robust detection of infeasibility in QPs. We also presented an infeasible IPM for the solution of the HQP and showed polynomial iteration complexity. In the context of IPMs, the paper leaves open a number of future avenues for exploration. Firstly, the sparsity of the linear systems in the step computation of IIPM are now affected by the sparsity of $c,f$ which can be dense. To effectively handle this a tailored linear algebra solution is required. Secondly, the linear system $E$ may have structure such as in the case of MPC. Exploiting this through a decomposition approach such as in [3] is critical to improving the computational efficiency. Thirdly, the user of predictor-corrector steps will be crucial for improving the practial performance. We will explore these in a future work. The homogeneous formulation can be readily presented to any QP algorithm. We will also investigate the performance of other QP algorithms on the proposed formulation.