An Active-Set Method for Quadratic Programming Based On Sequential Hot-Starts

Travis C. Johnson, Christian Kirches, Andreas Wächter
2015 SIAM Journal on Optimization  
A new method for solving sequences of quadratic programs (QPs) is presented. For each new QP in the sequence, the method utilizes hot-starts that employ information computed by an active-set QP solver during the solution of the first QP. This avoids the computation and factorization of the full constraint and Hessian matrices for all but the first problem in the sequence. The proposed algorithm can be seen as an extension of the iterative refinement procedure for linear systems to QP problems,
more » ... oupled with the application of an accelerated linear solver method that employs hot-started QP solves as a preconditioner. Local convergence results are presented. The practical performance of the proposed method is demonstrated on a sequence of QPs arising in nonlinear model predictive control and during the solution of a set of randomly generated nonlinear optimization problems using sequential quadratic programming. In these experiments, the method proves to be fairly reliable, despite the lack of global convergence guarantees. The results also show a significant reduction in the computation time for large problems with dense constraint matrices, as well as in the number of matrix-vector products. 1.1. Motivation. This research is motivated by a number of applications where a successive solution of QPs with similar data is required. (In this paper, we use the term "similar" loosely to express that the vectors and matrices, as well as the corresponding solutions of the QPs, are close to each other.) One example is nonlinear model-predictive control (NMPC), a numerical approach for optimally controlling a dynamic process (such as a chemical plant or a vehicle) in real-time. Here, at a given point τ in time, an optimal control action is computed as the solution of a QP that is obtained from a local quadratic model of the control problem and involves some linearization of the differential equations describing the process. (See section 6.1.3 for details.) The initial conditions and exogenous system parameters are chosen according to the actual or estimated state of the system at time τ . After a small time interval Δτ , a new control action is computed, now using the initial conditions and system parameter values corresponding to τ + Δτ . If Δτ is small and the state of the system has not changed very much, the QPs solved at τ and τ + Δτ as well as their solutions are similar. 969 The sequential quadratic programming (SQP) method for solving nonlinear programs (NLPs) represents another example. Consider an NLP of the form where the objective function f : R n −→ R and the constraint function c : R n −→ R m with m ≤ n are continuously differentiable. In a generic SQP method, the step d k at an iterate x k is obtained as the optimal solution of the QP (1.1) with = −x k , where g = ∇f (x k ) is the gradient of the objective function, c = c(x k ) is the residual of the constraints, A = ∇c(x k ) T is the Jacobian of the constraints, and W is (an approximation of) the Hessian of the Lagrangian function at x k for given multiplier estimates λ k for the equality constraints. Here, all vectors and matrices depend on the iterate x k , and, consequently, are different at each iterate of the method. However, if the iterates are close to each other, these quantities can be expected not to change very much. This is, for example, the case, when the SQP algorithm is close to convergence. Furthermore, suppose that we are interested in solving a sequence of NLPs (1.4) that differ only slightly in f (x) or c(x). If the optimal solution of the new NLP is close to the optimal solution of the previous one, the SQP method might require only a small number of iterations. The corresponding QPs are often similar not only to each other, but also across the different nonlinear problems. In this setting, it may be beneficial to solve the QPs arising during the SQP algorithm with the algorithm proposed in this paper, where the QP from the last SQP iteration of the first NLP is taken as the initial QP (1.2). This can be seen as a procedure for solving a sequence of similar NLPs using hot-starts. The solution of a sequence of closely related NLPs or QPs is also required during the execution of a branch-and-bound search for a mixed-integer nonlinear program (MINLP). Here, each node of the enumeration tree requires the solution of an NLP or QP relaxation, with different bound constraints. Moreover, during diving heuristics (see, e.g., [24, 23] and references therein) or strong-branching (see, e.g., [1, 23]), a succession of similar NLPs or QPs has to be solved. Structure of the article. This paper is organized as follows. Because our method crucially depends on the active-set method for QPs, we give a short summary in section 2 that introduces the notation and concepts necessary for the remained of this paper. In section 3, we briefly review the iterative refinement procedure for solving a linear system of equations. Reinterpreting this technique in the context of equality-constrained optimization in section 4.1, we make the connection to the use of hot-starts of QP (1.3) for the solution of the QPs (1.1). This is then generalized in section 4.2 to handle inequality constraints. An accelerated version is presented in section 5, where the solution of QP (1.3) is used as a preconditioner within an iterative linear solver method. In section 6.1, we explore the performance of the new QP solver in the context of an optimal control application. Section 6.2 examines the performance of the new method within an SQP framework applied to sequences of randomly generated NLPs with perturbed data. We conclude with some final remarks in section 7. Notation. Given a vector (or vector-valued function) x ∈ R n , we denote by x (i) the ith component of this vector. Given a set S ⊆ {1, . . . , n}, we denote by x S the vector composed from elements x (i) with indices i ∈ S, and S C denotes the complement Downloaded 11/09/15 to Redistribution subject to SIAM license or copyright; see T. C. JOHNSON, C. KIRCHES, AND A. WÄCHTER matrices in (4.13) and (5.6) are still nonsingular, and the sensitivity results required in the proof of Lemma 4.2 still apply. Numerical results. To examine the practical performance of the proposed approach, a prototype implementation of Algorithm 2, which we will refer to as iQP (inexact QP solver), was created in MATLAB R2012b. SQMR [17] was chosen as the iterative linear solver for the inner loop, because it exploits the symmetry of the matrix and allows indefinite preconditioners. For completeness, the detailed description of Algorithm 2 using SQMR is provided in Appendix A. We note that SQMR does not have theoretical convergence guarantees; in our implementation, SQMR is simply restarted if it breaks down, but this fall-back was triggered in our experiments very rarely. The QPs (4.7) and (5.4) were solved using the open-source active-set parametric QP solver qpOASES [11, 12] . All experiments were performed on an 8-core Intel-i7 3.4GHz 64bit Linux server with 32GB RAM. MATLAB was set to use only a single computational thread. We present two sets of numerical experiments. In section 6.1, Algorithm 2 is used to solve a sequence of QPs that arise in certain nonlinear model predictive control (NMPC) applications. In section 6.2, a sequence of randomly perturbed quadraticallyconstrained quadratic programs (QCQPs) is solved. The goal of these experiments is twofold. First, we explore the reliability of the new QP method in practice, given that convergence is not guaranteed. Second, we compare the performance of the iQP method (the hot-start approach) with that of a standard active-set solver, qpOASES. We refer to the latter as the warm-start approach to indicate that the solution of a new QP can be started from the optimal active set of the previous QP but that the KKT matrix in (2.2) has to be factorized from scratch because the entries changed. Because qpOASES uses dense linear algebra in its current implementation, our experiments are carried out for problems with dense matrices. Whether the new method requires overall less computation time for the solution of a new QP in a sequence depends on a number of factors. The warm-start approach requires the factorization of the KKT matrix in (2.2); this costs roughly O((n F ) 3 ) floating-point operations for dense matrices, where n F is the number of free variables. In addition, for each iteration of the active-set QP solver, in which one variable leaves or enters the active set, the linear system (2.2) is solved twice, and the factorization of the KKT matrix is updated for a new active set; this requires roughly O((n F ) 2 ) operations. On the other hand, the hot-start approach does not require a factorization with work O((n F ) 3 ), but the number of solves of the linear system (2.2) increases because one or two hot-started QPs have to be solved per iQP iteration, each of which might require several active-set changes, particularly in the first iQP iterations. Therefore, whether the new approach requires less computational effort depends on the number of iQP iterations and the relative cost of factorizing the KKT matrix in (2.2) versus the backsolve given the factorization. In some applications, the computation of the matrices A and W in (1.1) dominates the computational time. This is, for example, the case when the constraints arise from the integration of differential equations, as in the NMPC context (see section 6.1.3). Then, single matrix elements cannot be accessed individually, and the effort for computing the entire A matrix is equivalent to evaluating A rowwise (or columnwise) by computing n products of A (or m products of A T ) with unit vectors. Each product involves the computation of a directional derivative of the solution of the differential equation. Therefore, we also report the number of matrix-vector products involving both A · v and A T · v in our statistics. If this count is significantly smaller than n or Downloaded 11/09/15 to Redistribution subject to SIAM license or copyright; see
doi:10.1137/130940384 fatcat:6gtcmz2hhnarve4ab4yeoznd6a