Nonlinear Programming Strategies on High-Performance Computers Jia Kang, Naiyuan Chiang, Carl D. Laird, and Victor M. Zavala Abstract—We discuss structured nonlinear programming problems arising in control applications, and we review software and hardware capabilities that enable the efficient exploitation of such structures. We focus on linear algebra parallelization strategies and discuss how these interact and influence high- level algorithmic design elements required to enforce global convergence and deal with negative curvature in a nonconvex setting. I. NLP FRAMEWORK AND NOTATION We consider the general nonlinear programming problem (NLP) of the form min f (x) (I.1a) s.t. c(x)=0 (λ) (I.1b) x ≥ 0 (ν ) (I.1c) Here, x ∈< n are the primal variables, λ ∈< m are the dual variables for the equality constraints, and ν ∈< n + are the dual variables for the bounds. The objective function f (x): < n →< and constraints c(x): < n →< m are assumed to be twice continuously differentiable and are allowed to be nonconvex (the linear algebra concepts presented also hold for convex problems). General inequality constraints can be handled by introducing slack variables. To solve the NLP (I.1a), we use a primal-dual interior- point framework. This framework is particularly attractive because it enables modular implementations of linear algebra kernels. The barrier subproblem solved in the IP framework is given by min ϕ μ (x) := f (x) - μ n X j=1 ln x (j) s.t. c(x)=0, (λ) (I.2) where μ> 0 is the barrier parameter and x (j) denotes the j th component of vector x. To solve each barrier problem we apply Newton’s method to the Karush-Kuhn-Tucker system: ∇ x ϕ μ (x)+ ∇ x c(x)λ =0 (I.3a) c(x)=0 (I.3b) N. Chiang and V. M. Zavala are with the Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439, USA. E-mail: {nychiang, vzavala}@mcs.anl.gov. J. Kang, is a Ph.D. student in the Department of Chemical Engineering, Texas A&M University, College Station, TX. Email: jkang@tamu.edu. C. D. Laird is an Associate Professor in the School of Chemical Engineering, Purdue University, West Lafayette, IN. Email: carllaird@purdue.edu. Preprint ANL/MCS-P5359-0615. with the implicit restriction x ≥ 0. This restriction is typically enforced using a fraction-to-the-boundary rule [39]. We denote the variables at iterate k as (x k ,λ k ). In a line- search setting [16, 47], the primal search direction d k and the dual updates λ + k are typically obtained by solving the augmented system:  W k (δ w ) J T k J k -δ c I m  d k λ + k  = -  g k c k  . (I.4) Here, c k := c(x k ), J k := ∇ x c(x k ) T , g k := ∇ x ϕ μ k , W k (δ w ) := H k +Σ k + δ w I n , H k := ∇ xx L(x k ,λ k ), L(x k ,λ k ) := ϕ μ (x k )+ λ T k c(x k ), and Σ k := X -1 k V k with V k := diag(ν k ), and X k := diag(x k ). The term I n denotes an identity matrix of dimension n. We also define the primal and dual regularization parameters δ w ,δ c ≥ 0 and, to enable compact notation, we define the augmented matrix M k (δ w ,δ c ) :=  W k (δ w ) J T k J k -δ c I m  . (I.5) In a trust-region setting [9], the primal step is typically decomposed as d k = n k + t k , where n k is the normal component satisfying J k n k = -c k and t k is the tangential component computed from the augmented system,  W k (0) J T k J k  t k λ + k  = -  g k + W k (0)n k 0  . (I.6) We note that regularization is not added directly in the augmented matrix as in the line-search setting. The reason is that regularization is treated implicitly in trust-region procedures. The structure of the augmented systems in both settings, however, is the same. Consequently, unless stated otherwise, we will refer to the more general augmented system (I.4). Many methods can be considered for the parallel solution of (I.4) including iterative linear algebra or tailored decomposition of problem-specific structure; however, the possibility of nonconvexity may impact the algorithms se- lected. A key difference between convex and nonconvex NLPs is the presence of negative curvature. The presence of negative curvature indicates that the Newton step d k might not correspond to a minimum of the associated quadratic programming problem. In a line-search setting this is an im- portant issue because the Newton step cannot be guaranteed to provide a descent direction for the objective function when the constraint violation is sufficiently small. In particular, we seek that d T k g k < 0 whenever c k ≈ 0. In a trust-region setting, on the other hand, we need to guarantee that the