Adjoint-based SQP Method with Block-wise quasi-Newton Jacobian Updates for Nonlinear Optimal Control

TL;DR

This paper introduces an adjoint-based SQP method with block-wise quasi-Newton Jacobian updates for nonlinear MPC, offering improved convergence and computational efficiency in solving complex non-convex optimization problems.

Contribution

It presents a novel local convergence analysis and an efficient implementation of an adjoint-based SQP algorithm using block-structured TR1 quasi-Newton updates for nonlinear optimal control.

Findings

Convergence analysis confirms local convergence properties.

The proposed method reduces computational load in implicit discretizations.

Simulation case studies demonstrate improved performance in nonlinear MPC.

Abstract

Nonlinear model predictive control~(NMPC) generally requires the solution of a non-convex optimization problem at each sampling instant under strict timing constraints, based on a set of differential equations that can often be stiff and/or that may include implicit algebraic equations. This paper provides a local convergence analysis for the recently proposed adjoint-based sequential quadratic programming~(SQP) algorithm that is based on a block-structured variant of the two-sided rank-one~(TR1) quasi-Newton update formula to efficiently compute Jacobian matrix approximations in a sparsity preserving fashion. A particularly efficient algorithm implementation is proposed in case an implicit integration scheme is used for discretization of the optimal control problem, in which matrix factorization and matrix-matrix operations can be avoided entirely. The convergence analysis results as well as the computational performance of the proposed optimization algorithm are illustrated for two simulation case studies of nonlinear MPC.

Figures

Equations

$$\underset{x(\cdot ),u(\cdot )}{\min }$$

$$x_0-{\hat{x}}_0=0,$$

$$0=f(\dot{x}(t),x(t),u(t)),$$

$$p(x(t),u(t))\le 0,$$

$$0=f(\dot{x}(\tau ),x(\tau ),u(\tau )),\tau \in [t_i,t_{i+1}],x(t_i)=x_i,$$

$$\underset{X,U}{\min }$$

$${\hat{x}}_0=x_0,$$

$$F_i(x_i,u_i)=x_{i+1},$$

$$P_i{w}_i\le p_i,$$

$${\nabla }_w\mathcal{L}(w^{\ast },{\lambda }^{\ast })+P^T{\mu }^{\ast }$$

$$F(w^{\ast })$$

$$P{w}^{\ast }$$

$${\mu }^{\ast }$$

$${\mu }_j^{\ast }(P{w}^{\ast }-p{)}_j$$

$$\mathcal{L}(w,\lambda )=\sum _{i=0}^{N-1}\left(l_i(x_i,u_i)+{\lambda }_i^{\top }(F_i(w_i)-x_{i+1})\right)+l_N(x_N).$$

$$\left[\begin{array}{ccc}{H}^k& A^{k^{\top }}-E^{\top }& P_A^{\top }\\ A^k-E& & \\ P_A& & \end{array}\right]\left[\begin{array}{c}\Delta w^k\\ \Delta {\lambda }^k\\ \Delta {\mu }_A^k\end{array}\right]=-\left[\begin{array}{c}g(w^k,{\lambda }^k)\\ F(w^k)-X^k\\ P_A{w}^k-p_A\end{array}\right],$$

$$\underset{\Delta W}{\min }$$

$$\Delta x_0={\hat{x}}_0-x_0^k,$$

$$a_i^k+A_i^k\Delta w_i=\Delta x_{i+1},$$

$$P_i\Delta w_i\le p_i^k,$$

$$h_i^k:={\nabla }_{w_i}l(w_i^k)+{\left(\frac{\partial F_i}{\partial w_i}(w_i^k)-A_i^k\right)}^{\top }{\lambda }_i^k,$$

$$w^{k+1}=w^k+\Delta w^k\text{ and }{\lambda }^{k+1}={\lambda }^k+\Delta {\lambda }^k={\lambda }_{QP}^{k+1},$$

$$\text{Adjoint Condition (AC):}{\sigma }_i^{k^{\top }}{A}_i^{k+1}$$

$$\text{Forward Condition (FC):}{A}_i^{k+1}{s}_i^k$$

$$A_i^{k+1}=A_i^k+{\alpha }_i^k\left(y_i^k-A_i^k{s}_i^k\right)\left({\gamma }_i^{k^{\top }}-{\sigma }_i^{k^{\top }}{A}_i^k\right),$$

$${\alpha }_{\mathrm{A},i}^k=\frac{1}{{\sigma }_i^{k^{\top }}(y_i^k-A_i^k{s}_i^k)},$$

$${\alpha }_{\mathrm{F},i}^k=\frac{1}{({\gamma }_i^{k^{\top }}-{\sigma }_i^{k^{\top }}{A}_i^k)s_i^k},$$

$$\left|({\gamma }_i^{k^{\top }}-{\sigma }_i^{k^{\top }}{A}_i^k)s_i^k\right|\ge c_1\Vert {\sigma }_i^k\Vert \Vert y_i^k-A_i^k{s}_i^k\Vert ,$$

$$\left|{\sigma }_i^{k^{\top }}(y_i^k-A_i^k{s}_i^k)\right|\ge c_1\Vert s_i^k\Vert \Vert {\gamma }_i^k-A_i^{k^{\top }}{\sigma }_i^k\Vert ,$$

$$J^k:=\left[\begin{array}{cc}{N}^{\top }{H}^k& N^{\top }{A}^{k^{\top }}\\ A^k& \\ P_A& \end{array}\right]\approx \frac{\partial \mathcal{F}}{\partial y}(y^k),\text{where}\mathcal{F}(y):=\left[\begin{array}{c}{N}^{\top }{\nabla }_w\mathcal{L}(w,\lambda )\\ F(w)-X\\ P_{A}w-p_A\end{array}\right],$$

$$w^{k+1}=w^k+\Delta w^k\text{ and }{\lambda }^{k+1}={\lambda }^k+\Delta {\lambda }^k={\lambda }_{QP}^{k+1},$$

$$\Vert J^{{k+1}^{-1}}\left(J^k-\frac{\partial \mathcal{F}}{\partial y}(y^k+t\Delta y^k)\right)\Delta y^k\Vert \le \kappa \Vert \Delta y^k\Vert ,\forall t\in [0,1].$$

$$\underset{k\to \infty }{\lim }\left[\begin{array}{cc}{N}^{\top }{H}^{k}N& N^{\top }{A}^{k^{\top }}\\ A^{k}N& 0\end{array}\right]=\left[\begin{array}{cc}{N}^{\top }{\nabla }_{ww}^2\mathcal{L}(w^{\star },{\lambda }^{\star })N& N^{\top }\frac{\partial F}{\partial w}(w^{\star }{)}^{\top }\\ \frac{\partial F}{\partial w}(w^{\star })N& 0\end{array}\right],$$

$$\underset{k\to \infty }{\lim }\Vert \left(A_i^k-\frac{\partial F_i}{\partial w}(w_i^{\ast })\right)N_i\Vert =0,\forall i=0,\dots ,N-1,$$

$$S_{N_i}^{k_i}=\left[\begin{array}{ccc}\frac{s_{N,i}^{k_i^1}}{\Vert s_{N,i}^{k_i^1}\Vert }& \dots & \frac{s_{N,i}^{k_i^{q_i}}}{\Vert s_{N,i}^{k_i^{q_i}}\Vert }\end{array}\right]$$

$$\Vert y_i^k-A_i^l{s}_i^k\Vert$$

$$\Vert {\gamma }_i^k-A_i^{l^{\top }}{\sigma }_i^k\Vert$$

$$\Vert y_i^k-A_i^{l+1}{s}_i^k\Vert$$

Full text

{amscode}

49K15; 49M37; 90C53; 65K05

1 Introduction

Optimization based control and estimation techniques have attracted an increasing attention over the past decades. They allow a model-based design framework, in which the system dynamics, performance metrics and constraints can directly be taken into account. Receding horizon techniques such as model predictive control (MPC) and moving horizon estimation (MHE) have been studied extensively because of their desirable properties [30] and these optimization-based techniques have already been applied in a wide range of applications [21]. One of the main practical challenges in implementing such an optimization-based predictive control or estimation scheme, lies in the ability to solve the corresponding nonlinear and generally non-convex optimal control problem (OCP) under strict timing constraints and typically on embedded hardware with limited computational capabilities and available memory.

Let us consider the following continuous-time formulation of the optimal control problem that needs to be solved at each sampling instant

[TABLE]

where $T$ denotes the control horizon length, $x(t)\in\mathbb{R}^{n_{\mathrm{x}}}$ denotes the differential states and $u(t)\in\mathbb{R}^{n_{\mathrm{u}}}$ are the control inputs. The function $\ell(\cdot)$ defines the stage cost and the nonlinear system dynamics are formulated as an implicit system of ordinary differential equations (ODE) in (1c), which could additionally be extended with implicit algebraic equations. A common assumption is that the resulting system of differential-algebraic equations (DAE) is of index $1$ [7]. The optimization problem is parametric, since it depends on the state estimate $\hat{x}_{0}$ at the current sampling instant, through the initial value condition in (1b). The path constraints are defined by the function $p(\cdot)$ in Eq. (1d) and, for simplicity of notation, they are further assumed to be affine. Note that a similar problem as in (1) needs to be solved for optimization-based state and parameter estimation, without the given initial state value.

In direct optimal control methods, one forms a discrete-time approximation of the continuous-time OCP in (1) based on an appropriate parameterization of the state and control trajectories over the time horizon $t\in[0,T]$, resulting in a tractable nonlinear program (NLP) that needs to be solved. Popular examples of this approach include the direct multiple shooting method [8] and direct collocation [3, 5]. Note that these techniques often need to rely on implicit integration methods in order to deal with stiff and/or implicit systems of differential or differential-algebraic equations [36]. The resulting constrained optimization problem can be handled by standard Newton-type algorithms such as interior point methods [38] and sequential quadratic programming (SQP) [10] techniques for nonlinear optimization [31].

Quasi-Newton optimization methods are generally popular for solving such a constrained NLP. They result in computationally efficient Newton-type methods that solve the first order necessary conditions of optimality, i.e., the Karush-Kuhn-Tucker (KKT) conditions, without evaluating the complete Hessian of the Lagrangian and/or even without evaluating the Jacobian of the constraints [31]. Instead, quasi-Newton methods are based on low-rank update formulas for the Hessian and Jacobian matrix approximations [14]. Popular examples of this approach include the Broyden-Fletcher-Goldfarb-Shanno (BFGS) [11] and the symmetric rank-one (SR1) update formula [13] for approximating the Hessian of the Lagrangian. Similarly, quasi-Newton methods can be used for approximating Jacobian matrices, e.g., of the constraint functions, such as the good and bad Broyden methods [12] as well as the more recently proposed two-sided rank-one (TR1) update formula [25].

For the purpose of real-time predictive control and estimation, continuation-based online algorithms have been proposed that aim at further reducing the computational effort by exploiting the fact that a sequence of closely related parametric optimization problems is solved [7, 17]. One popular technique consists of the real-time iteration (RTI) algorithm that performs a single SQP iteration per time step, in combination with a sufficiently high sampling rate and a prediction-based warm starting in order to allow for closed-loop stability of the system [16]. The RTI algorithm can be implemented efficiently based on (fixed-step) integration schemes with tailored sensitivity propagation for discretization and linearization of the system dynamics [36] in combination with structure-exploiting quadratic programming solvers [21]. In addition, a lifted algorithm implementation has been proposed in [37] to directly embed the iterative procedure of implicit integration schemes, e.g., collocation methods, within a Newton-type optimization framework for optimal control.

Unlike standard inequality constrained optimization, nonlinear optimal control problems typically result in a particular sparsity structure in the Hessian of the Lagrangian and in the Jacobian matrix for the equality constraints. In direct optimal control methods, the objective function is typically separable resulting in a block-diagonal Hessian matrix. This property has been exploited in partitioned quasi-Newton methods that approximate and update each of the Hessian block matrices separately, as proposed and studied in [24, 23, 29]. On the other hand, the Jacobian matrix corresponding to the discretized system dynamics has a block bidiagonal sparsity structure, because of the stage wise coupling of the optimization variables at subsequent time steps of the control horizon. For this purpose, the present article analyzes a novel tailored quasi-Newton method for optimal control using a partitioned or block-structured TR1-based Jacobian update formula. This adjoint-based SQP method for nonlinear optimal control, based on a Gauss-Newton Hessian approximation in combination with inexact Jacobian matrices, was proposed recently in [28].

1.1 Contributions and outline

This paper provides a complete presentation of the block-TR1 based SQP method for nonlinear optimal control, including a detailed discussion of the lifted collocation type implementation, extending earlier work of the same authors in [28]. Unlike the latter publication, a convergence analysis of this novel quasi-Newton type optimization algorithm is provided. More specifically, we prove convergence of the block-structured quasi-Newton Jacobian approximations to the exact Jacobian matrix within the null space of the active inequality constraints. Based on this result, under mild conditions, convergence of the overall inexact SQP method can be guaranteed. Locally linear or superlinear convergence rates can be shown, respectively, when using a Gauss-Newton or quasi-Newton based Hessian approximation scheme. In addition, it is shown how this convergence analysis extends to our lifted collocation implementation that avoids any matrix factorization or matrix-matrix operations. These convergence analysis results as well as the computational performance of the optimization algorithms are illustrated numerically for two simulation case studies of nonlinear MPC.

The paper is organized as follows. Section 2 briefly introduces the direct multiple shooting based OCP problem formulation as well as the proposed adjoint-based inexact SQP method that is based on block-wise TR1 Jacobian updates. Section 3 presents the detailed convergence analysis for the optimization method and contains the main theoretical results of the present paper. A numerically efficient implementation of the block-TR1 update formula in combination with a lifted Newton-type method for direct optimal control with implicit integration schemes such as, e.g., collocation methods, is then proposed and analyzed in Section 4. Finally, Section 5 presents numerical results of the NMPC case studies and Section 6 concludes the paper.

2 Block-wise TR1 based Sequential Quadratic Programming

A popular approach for direct optimal control is based on direct multiple shooting [8] that performs a time discretization, based on a numerical integration scheme [26] to solve the following initial value problem

[TABLE]

on each of $N$ shooting intervals that are defined by a grid of consecutive time points $t_{i}$ for $i=0,\ldots,N$. For the sake of simplicity, we consider here an equidistant grid over the control horizon, i.e., $t_{i+1}-t_{i}=\frac{T}{N}$, and a piecewise constant control parametrization $u(\tau)=u_{i}$ for $\tau\in[t_{i},t_{i+1})$ in (2). An explicit fixed-step integration scheme defines the discrete-time system dynamics $x_{i+1}=F_{i}(x_{i},u_{i})$ for the shooting interval $[t_{i},t_{i+1}]$. For example, this can correspond to the popular Runge-Kutta method of order $4$ (RK4) as defined in [26]. Based on the explicit discretization scheme, the resulting block-structured optimal control problem reads as

[TABLE]

where the affine path constraints (3d) have been imposed on each of the shooting nodes and the compact notation $w_{i}:=(x_{i},u_{i})$ for $i=0,\ldots,N-1$ and $w_{N}:=x_{N}$ is defined. Note that the optimization variables for the problem in (3) are directly the state $X=[x_{0}^{\top},\ldots,x_{N}^{\top}]^{\top}$ and control trajectory $U=[u_{0}^{\top},\ldots,u_{N-1}^{\top}]^{\top}$.

2.1 SQP algorithm with inexact Jacobians

For a local minimum $w^{*}$ of the NLP in (3), for which the linear independence constraint qualification (LICQ) holds, there must exist a unique set of multiplier values $\lambda^{*}$, $\mu^{*}$ such that the following Karush-Kuhn-Tucker (KKT) conditions are satisfied

[TABLE]

where $F(\cdot)$ and $P$ are appropriate block-wise concatenations of the equality and inequality constraints, respectively, in (3c) and (3d) and $n_{\mathrm{p}}$ denotes the total number of inequality constraints. Here, we also lumped the initial condition constraint as part of the matrix $P$ since we can represent a linear equality constraint as two linear inequality constraints. Lastly, $\mathcal{L}(w,\lambda)$ denotes the ‘truncated Lagrangian’, omitting inequality constraints, and is therefore given by

[TABLE]

Given the set of indices $\mathcal{A}$ for the inequality constraints that are active at the local minimum, the KKT system reduces to a nonlinear system of equations that can be solved directly by a Newton-type method. In particular, we are interested in a quasi-Newton algorithm where we will approximate $\nabla_{ww}^{2}\mathcal{L}(w^{k},\lambda^{k})$ by a matrix $H^{k}$ and $\frac{\partial F}{\partial w}(w^{k})$ by a matrix $A^{k}$. Namely, we solve the following linearized system

[TABLE]

where $g(w^{k},\lambda^{k})=\nabla_{w}\mathcal{L}(w^{k},\lambda^{k})+\mu^{k^{\top}}_{A}(P_{A}\,w^{k}-p_{A})$ at each Newton-type iteration $k$. Note that the matrix $P_{A}$ is defined as the part of $P$ that corresponds to the inequality constraints (4c) in the active set $\mathcal{A}$, and $E$ denotes the constant matrix corresponding to the right-hand side of the equality constraints in (4b).

In order to efficiently solve the inequality constrained OCP in (3), let us consider the adjoint-based SQP algorithm with Gauss-Newton type Hessian approximation and inexact Jacobian information as introduced originally in [7, 39] for fast nonlinear MPC. Each SQP iteration solves a convex QP subproblem

[TABLE]

where notation $\Delta W=[\Delta w_{0}^{\top},\ldots,\Delta w_{N}^{\top}]^{\top}$ is used to denote the deviation variables $\Delta w_{i}:=w_{i}-w_{i}^{k}$, given the current solution guess $X^{k}$, $U^{k}$ for the state and control trajectories at iteration $k$ of the adjoint-based SQP method. The function $p(\cdot)$ that defines the path constraint (1d) was assumed to be affine and $p_{i}^{k}:=p_{i}-P_{i}\,w_{i}^{k}$. Note that tracking formulations for nonlinear MPC typically include a stage cost that is defined by a (nonlinear) least squares term $l_{i}(x_{i},u_{i})\;=\;\frac{1}{2}\|R(x_{i},u_{i})\|_{2}^{2}$ for $i=0,\ldots,N$. The generalized Gauss-Newton (GGN) method from [6] uses the block-structured Hessian approximation $H_{i}^{k}:=\nabla R(w_{i}^{k})\nabla R(w_{i}^{k})^{\top}\approx\nabla_{w_{i}w_{i}}^{2}{\mathcal{L}}(\cdot)$.

The matrix $A_{i}^{k}\approx\frac{\partial F_{i}}{\partial w_{i}}(w_{i}^{k})$ denotes the Jacobian approximation and $a_{i}^{k}:=F_{i}(w_{i}^{k})-x_{i+1}^{k}$ for the discrete-time system dynamics in Eq. (7c). For real-time NMPC, such a Jacobian approximation can be obtained by reusing information from a previous NLP solution [7, 39]. The gradient term in the objective (7a) reads as

[TABLE]

for $i=0,\ldots,N-1$, in which $\lambda_{i}^{k}$ denotes the current value of the Lagrange multipliers for the nonlinear continuity constraints in (3c). Note that the linearized KKT conditions in (6) correspond to the KKT optimality conditions for the QP in (7), for a fixed active set $\mathcal{A}$. In addition, each QP subproblem is convex because $H^{k}\succeq 0$, e.g., for the Gauss-Newton Hessian approximation. A full-step inexact SQP method will sequentially solve each QP subproblem (7) and perform the following updates:

[TABLE]

where $\lambda^{k+1}_{QP}$ denote the Lagrange multiplier values for Eq. (7c) at the QP solution. We do not need to perform explicit updates for the Lagrange multipliers associated with the inequality constraints, because they are assumed to be affine, hence not impacting any computation on the QP formulation in (7).

2.2 Dynamic block-wise TR1 Jacobian updates

At each SQP iteration, we perform the block-wise two-sided rank-one (TR1) Jacobian update, as proposed recently in [28]. Following the work in [25], given current Jacobian approximations $A^{k}_{i}$ for $i=0,\ldots,N-1$, we would like that each updated approximation matrix $A^{{k+1}}_{i}$ satisfies the following two secant conditions

[TABLE]

where we define the adjoint vector $\gamma_{i}^{k}=\frac{\partial F_{i}}{\partial w_{i}}(w^{{k+1}}_{i})^{\top}\sigma_{i}^{k}$, given $\sigma_{i}^{k^{\top}}=(\lambda^{{k+1}}_{i}-\lambda^{k}_{i})^{\top}$, and the difference in function evaluations $y_{i}^{k}=F(w^{{k+1}}_{i})-F(w^{k}_{i})$. Note that $\lambda^{{k+1}}_{i}$ and $\lambda^{k}_{i}$, respectively, denote the new and old Lagrange multipliers for the linearized equality constraints in Eq. (7c). Similarly, $w_{i}^{k}:=(x_{i}^{k},u_{i}^{k})$ and $w_{i}^{{k+1}}:=w_{i}^{k}+\Delta w_{i}^{k}$ denote, respectively, the old and new primal variables, such that $s_{i}^{k}:=w_{i}^{{k+1}}-w_{i}^{k}$. Note that the gradient $\gamma_{i}^{k}=\frac{\partial F_{i}}{\partial w_{i}}(w^{{k+1}}_{i})^{\top}\sigma_{i}^{k}$ can be computed efficiently using the backward or adjoint mode of algorithmic differentiation (AD), e.g., see [22].

The proposed block-wise TR1 update formula then reads as follows

[TABLE]

for $i=0,\ldots,N-1$ and where $\alpha_{i}^{k}$ is a scalar that will be defined further. Aside from the case where the function $F(\cdot)$ is affine, the two conditions in Eq. (10) are not consistent with each other and they can therefore generally not both be satisfied by the updated matrix $A^{{k+1}}_{i}$ at each iteration. Thus, similar to the standard TR1 update in [25], the block-wise update will only be able to satisfy one or the other. In the adjoint variant of the update, the scaling value is defined as

[TABLE]

such that the adjoint condition in (10) is satisfied exactly and the forward condition holds up to some accuracy. Similarly, this value reads as follows for the forward variant

[TABLE]

where the forward condition is satisfied exactly. It is interesting to note that, since we apply the block-wise TR1 update from (11) for each shooting interval $i=0,\ldots,N-1$, the resulting update for the complete constraint Jacobian matrix of the QP in (7) corresponds to a rank-$N$ update.

As in [25], we impose a skipping condition in order to avoid a potential blow-up of the block-wise TR1 update when the denominator of the scaling factor becomes small or even zero. For our purposes, the skipping condition itself depends on the type of formula that is used. We update the block matrix $A_{i}^{k}$ only if the following holds

[TABLE]

with $c_{1}\in(0,1)$ if $\alpha_{i}^{k}=\alpha_{\mathrm{F},i}^{k}$ in the forward TR1 update, and

[TABLE]

with $c_{1}\in(0,1)$ if $\alpha_{i}^{k}=\alpha_{\mathrm{A},i}^{k}$ in the adjoint TR1 update. In addition, in order to consistently choose either the forward or adjoint Jacobian update formula, we propose a more dynamic variant of the algorithm that picks either $\alpha_{\mathrm{F},i}^{k}$ or $\alpha_{\mathrm{A},i}^{k}$ for each block matrix at any given iteration. It may not be clear what is the best approach to select which type of update is to be executed for each block matrix at a given iteration. However, in the next section, we prove the local convergence properties of the algorithm under any arbitrary sequence of updates that satisfy the skipping conditions in (14) and (15) for each block $i$ at every iteration $k$.

The complete adjoint-based SQP method that uses parallelizable block-wise TR1 Jacobian updates is summarized in Algorithm 1. Note that, for simplicity, the SQP algorithm is presented as a full-step method without any globalization or step-length selection strategies to ensure convergence to a local minimum [31]. This is also further motivated by the use of online algorithms for real-time nonlinear MPC as discussed in [17].

3 Convergence Results for Block-wise TR1-based SQP Method

For the convergence analysis of sequential quadratic programming, it is standard to rely on a result that the active set, i.e., the set of active inequality constraints in the QP subproblems is stable in a neighbourhood around a local minimizer of the nonlinear program [31]. This allows us to study the local convergence properties of the block-TR1 based SQP method under the assumption that the active set has already been fixed, resulting, locally, in an equality constrained problem.

3.1 Stability of the active set and local convergence

Let us start by briefly repeating the result from [20] on the stability of the active set in the QP subproblems near the NLP solution and the corresponding conditions on local convergence properties for an adjoint-based SQP method with inexact Jacobians.

Theorem 3.1**.**

(Stability of active set and local convergence) Let the NLP solution vectors $w^{*}$, $\lambda^{*}$ be given and assume that:

1. (i)

at $w^{*}$ LICQ holds, and there exist Lagrange multiplier values $\mu^{*}$ such that $(w^{*},\lambda^{*},\mu^{*})$ satisfies the KKT conditions in (4). 2. (ii)

at $w^{*}$ strict complementarity holds, i.e., the multipliers $\mu_{A}^{*}$ of the active inequalities $P_{A}w^{*}=p_{A}$ satisfy $\mu_{A}^{*}>0$, where $P_{A}$ is a matrix consisting of all rows of $P$ that correspond to the active inequalities at the NLP solution. 3. (iii)

there are two sequences of uniformly bounded matrices $(A^{k},H^{k})$, each $H^{k}$ positive semidefinite on the null space of $A^{k}$, such that the sequence of matrices

[TABLE]

*is uniformly bounded and invertible with a uniformly bounded inverse. Here, $N$ is a null space matrix with appropriate dimensions with orthonormal column vectors such that $N^{\top}N=\mathbb{1}$ and $P_{A}\,N=0$. * 4. (iv)

there is a sequence of iterates $y^{k}:=(w^{k},\lambda^{k})$ generated according to

[TABLE]

where $\Delta w^{k}$ is the primal solution of the QP subproblem in (7) and $\lambda^{k+1}_{QP}$ denote the Lagrange multipliers corresponding to the equality constraints (7c). Each iteration can be written in compact form as $y^{k+1}=y^{k}-J^{k^{-1}}\mathcal{F}(y^{k})$. 5. (v)

there exists $\kappa<1$ such that, for all $k\in\mathbb{N}$, it can be guaranteed that

[TABLE]

Then, there exists a neighbourhood $\bar{\mathcal{N}}$ of $(w^{*},\lambda^{*})$ such that for all initial guesses $(w^{0},\lambda^{0})\in\bar{\mathcal{N}}$ the sequence $(w^{k},\lambda^{k})$ converges q-linearly towards $(w^{*},\lambda^{*})$ with rate $\kappa$, and the solution of each QP (7) has the same active set as $w^{*}$.

In addition to the latter result that guarantees a q-linear local convergence rate in a neighbourhood of the NLP solution, the following theorem states a condition under which q-superlinear local convergence can be obtained instead.

Theorem 3.2**.**

(Superlinear convergence) If the equality

[TABLE]

holds in addition to the assumptions of Theorem 3.1, then the local convergence rate is q-superlinear instead.

The proofs for both Theorem 3.1 and 3.2 can be found in [20] for an adjoint-based SQP method with inexact Jacobians that matches our problem formulation.

3.2 Convergence of the block-wise TR1 Jacobian updates

Theorem 3.1 holds for a general class of constraint Jacobian and Hessian approximation matrices $(A^{k},H^{k})$. Therefore, we have to show that our block-wise TR1 updates produce a sequence of block-structured matrices that converge to the exact Jacobian, which is itself block-structured, projected onto the null space of the active inequality constraint matrix $P_{A}$. Namely, defining a null space matrix $N$ as in Theorem 3.1, we need to prove that the following holds

[TABLE]

where $N_{i}$ is the projection of the null space matrix $N$ in the variable space corresponding to block $i$. The only non-zero entries that are inexact in the Jacobian approximation matrix $A^{k}$ are those corresponding to the block-TR1 matrices $A^{k}_{i}$, $i=0,\ldots,N-1$.

Assumption 3.3**.**

Let us make the following assumptions:

• (AS1)

The Lagrangian function is twice continuously differentiable.

• (AS2)

The function $\nabla_{w}F(w)$ is Lipschitz continuous, i.e., there exists a constant $c_{3}$ such that $\|\nabla_{w}F(w_{1})-\nabla_{w}F(w_{2})\|\leq c_{3}\|w_{1}-w_{2}\|$, for any $w_{1}$, $w_{2}$.

• (AS3)

*Let $\{(w^{k},\lambda^{k})\}$ be a sequence of iterates generated by our block-TR1 based SQP method in Algorithm 1, with a corresponding sequence of update parameters $\{\alpha^{k}_{i}\}$, while satisfying the skipping criterion in eqs. (14)-(15). *

• (AS4)

The SQP iterates $\{(w^{k},\lambda^{k})\}$ converge to a limit point $(w^{*},\lambda^{*})$.

• (AS5)

There is $k_{0}$ such that the active set is stable for all iterates $k\geq k_{0}$.

• (AS6)

For each block $i$, the sequence of projections of $\{s^{k}\}$ on the subspace associated with block $i$, namely $\{s^{k}_{i}\}$ is uniformly linearly independent in the projected null space $N_{i}$, i.e., there exist $c_{4}>0$ and $l\geq q_{i}$ such that, for each $k_{i}\geq k_{0}$, there exist $q_{i}$ distinct indices $k^{j}_{i}$ with $k_{i}\leq k^{1}_{i}<\ldots<k^{q_{i}}_{i}\leq k_{i}+l$, $s^{k^{j}_{i}}_{N,i}\in\mathbb{R}^{q_{i}}$, $s^{k^{j}_{i}}_{i}=N_{i}s^{k^{j}_{i}}_{N,i}$, $j=1,\ldots,q_{i}$ and the minimum singular value $\sigma_{min}(S^{k_{i}}_{N_{i}})$ of the matrix

[TABLE]

is bounded below by $c_{4}$, i.e., $\sigma_{min}(S^{k_{i}}_{N_{i}})\geq c_{4}$.

Note that the assumptions $(AS1)$-$(AS5)$ are relatively mild and quite standard in Newton-type convergence analysis of SQP methods [20]. Especially, condition $(AS5)$ holds due to the local stability result in Theorem 3.1 for the active set near the NLP solution. Even though $(AS6)$ seems relatively strong, a very similar assumption is made in existing convergence results for quasi-Newton type matrix update schemes [13, 20]. Here, we only require uniform linear independence inside each block $i$. We proceed now to prove the convergence of the quasi-Newton block-structured constraint Jacobian approximation matrices, using ideas from [13] and [20]. We start by first showing an intermediate result in the following lemma.

Lemma 3.4**.**

Given $(AS1)$-$(AS3)$ in Assumption 3.3, then the following holds for each Jacobian block matrix approximation

[TABLE]

where $i=0,\ldots,N-1$ and $\eta^{l,k}_{i}=\max\{\|w^{p}_{i}-w^{s}_{i}\|\text{ }|\text{ }k\leq s\leq p\leq l\}$ is defined.

Proof.

Our proof follows closely the proof of Lemma 4.1 in [20] but extended to our block-structured method and generalized to include both the forward and adjoint TR1 Jacobian update formulas.

Step 1: Eq. (20a) based on forward Jacobian update

We start by showing the result of Eq. (20a) when $\alpha_{i}=\alpha_{\mathrm{F},i}$. The proof is by induction for each block $i=0,\ldots,N-1$. For $l=k+1$, we know that $y^{k}_{i}-A^{l}_{i}s^{k}_{i}=0$ based on the forward update. Assume that the result in (20a) holds for all $j=k+1,\ldots,l$. Then, we have the following

[TABLE]

where $\tau^{l}_{i}=\gamma^{l}_{i}-A^{l^{\top}}_{i}\sigma^{l}_{i}$ and $\rho^{l}_{i}=y^{l}_{i}-A^{l}_{i}s^{l}_{i}$ such that $\alpha_{\mathrm{F},i}^{l}=\frac{1}{(\tau^{l}_{i},s^{l}_{i})}$. Then, using the result in Eq. (20a), we can write

[TABLE]

Using the mean-value theorem, it follows that

[TABLE]

based on the Lipschitz continuity in $(AS2)$. From the skipping condition in (14), we know that $\left|(\tau^{l}_{i},s^{l}_{i})\right|\geq c_{1}\left\|\sigma_{i}^{l}\right\|\,\left\|\rho^{l}_{i}\right\|$. In addition, given that $\eta^{l,k}_{i}\leq\eta^{l+1,k}_{i}$, we obtain

[TABLE]

Step 2: Eq. (20a) based on adjoint Jacobian update

Let us continue this proof by induction for Eq. (20a), based on the adjoint Jacobian update formula. First, we derive the following error bound for the adjoint Jacobian update in case $l=k+1$ in Eq. (20a)

[TABLE]

From the skipping conditions in (14)-(15), we know that $\left|(\sigma^{k}_{i},\rho^{k}_{i})\right|\geq c_{1}\left\|s^{k}_{i}\right\|\,\left\|\tau^{k}_{i}\right\|$ and $\left\|\rho^{k}_{i}\right\|\leq\frac{\left|(\tau^{k}_{i},s^{k}_{i})\right|}{c_{1}\left\|\sigma_{i}^{k}\right\|}\leq\frac{\left\|\tau^{k}_{i}\right\|\left\|s^{k}_{i}\right\|}{c_{1}\left\|\sigma_{i}^{k}\right\|}$ holds. We can use these lower and upper bounds to rewrite the latter expression as

[TABLE]

where we additionally used the result

[TABLE]

Note that $\eta^{k+1,k}_{i}=\left\|s^{k}_{i}\right\|$ such that Eq. (20a) holds in case $l=k+1$. Assume that the result in (20a) holds for all $j=k+1,\ldots,l$. Then, we have the following

[TABLE]

for the adjoint Jacobian update formula in which $\alpha_{\mathrm{A},i}=\frac{1}{(\sigma_{i}^{l},\rho^{l}_{i})}$. From the skipping conditions in (14)-(15), we know that $\left|(\sigma^{l}_{i},\rho^{l}_{i})\right|\geq c_{1}\left\|s^{l}_{i}\right\|\,\left\|\tau^{l}_{i}\right\|$ and $\left\|\rho^{l}_{i}\right\|\leq\frac{\left|(\tau^{l}_{i},s^{l}_{i})\right|}{c_{1}\left\|\sigma_{i}^{l}\right\|}\leq\frac{\left\|\tau^{l}_{i}\right\|\left\|s^{l}_{i}\right\|}{c_{1}\left\|\sigma_{i}^{l}\right\|}$ holds such that $\frac{\left\|\sigma^{l}_{i}\right\|\left\|\rho^{l}_{i}\right\|}{\left|(\sigma^{l}_{i},\rho^{l}_{i})\right|}\leq\frac{1}{c_{1}}\frac{\left\|\sigma^{l}_{i}\right\|\left\|\rho^{l}_{i}\right\|}{\left\|s^{l}_{i}\right\|\,\left\|\tau^{l}_{i}\right\|}\leq\frac{1}{c_{1}^{2}}$. In addition, given eqs. (22) and (23) and given that $\eta^{l,k}_{i}\leq\eta^{l+1,k}_{i}$, we obtain

[TABLE]

Note that the induction proof of step 1 and 2 implies that Eq. (20a) additionally holds when switching between the forward and adjoint Jacobian update formulas. A similar induction-based proof can be used to show the result of Eq. (20b) for the dynamic block-TR1 Jacobian updates. ∎

Now, we present the resulting theorem on the convergence of the Jacobian approximation for the block-wise TR1 scheme under any sequence of decision rules that select the adjoint or forward updates at every iteration $k$ and for each block $i=0,\ldots,N-1$.

Theorem 3.5**.**

Given $(AS1)$-$(AS6)$ in Assumption 3.3, then the following holds for each Jacobian block matrix $i=0,\ldots,N-1$

[TABLE]

such that the following holds for the complete Jacobian approximation

[TABLE]

Proof.

Based on the inequality $\|w^{p}_{i}-w^{s}_{i}\|\leq\|w^{p}_{i}-w^{*}_{i}\|+\|w^{s}_{i}-w^{*}_{i}\|$ and using the definition $\eta^{l,k}_{i}=\max\{\|w^{p}_{i}-w^{s}_{i}\|\text{ }|\text{ }k\leq s\leq p\leq l\}$, one obtains

[TABLE]

for $l\geq q_{i}$ and $q_{i}$ is defined as in Assumption 3.3. In addition, the following holds

[TABLE]

at an iteration $j$, where $k\leq j\leq k+l$, regardless of whether the forward or adjoint Jacobian update formula has been used. Moreover, from Lemma 3.4, we have that

[TABLE]

We use the triangle inequality to obtain

[TABLE]

which holds for a sequence of indices $j=k_{i}^{1},\ldots,k_{i}^{q_{i}}$. Then, we use the linear independence condition $(AS6)$ in Assumption 3.3 that guarantees both existence of the inverse $(S^{k_{i}}_{N_{i}})^{-1}$ and the upper bound $\|(S^{k_{i}}_{N_{i}})^{-1}\|\leq 1/c_{4}$, such that

[TABLE]

where c_{5}=\frac{c_{3}}{c_{4}}\Bigg{(}\frac{2}{c_{1}}\left(\frac{2}{c_{1}^{2}}+1\right)^{l+1}+1\Bigg{)}\sqrt{q_{i}} has been defined. Lastly, the result in Eq. (30) follows from the fact that assumption $(AS4)$ implies that $\nu^{k}_{i}$ tends to zero. Note that this asymptotic result holds regardless of which Jacobian update (adjoint or forward TR1 formula) is performed for each block $i=0,\ldots,N-1$. The same convergence result then holds for the complete Jacobian matrix in (31), based on separability of the active inequality constraints and of the nonlinear constraint functions.

∎

3.3 Local rate of linear convergence for Gauss-Newton based SQP

One iteration of the adjoint-based Gauss-Newton SQP method solves the linear system in Eq. (6), which can be written in the following compact form

[TABLE]

where $\mathcal{F}(\cdot)$ denotes the KKT optimality conditions in the right-hand side of Eq. (6). Let us define regularity for a local minimizer $z^{\star}:=(w^{\star},\lambda^{\star},\mu^{\star})$ of the NLP, given a particular set of active inequality constraints. For this purpose, we rely on the linear independence constraint qualification (LICQ) and the second order sufficient conditions (SOSC) for optimality, of which the latter requires that the Hessian of the Lagrangian is strictly positive definite in the directions of the critical cone [31].

Definition 3.6**.**

A minimizer of an equality constrained NLP is called a regular KKT point, if both LICQ and SOSC are satisfied at this KKT point.

The convergence of this Newton-type optimization method then follows the classical and well-known local contraction theorem from [9, 15, 20, 32]. We use a particular version of this theorem from [19, 34], providing sufficient and necessary conditions for the existence of a neighbourhood of the solution where the Newton-type iteration converges locally. Let $\rho(P)$ denote the spectral radius, i.e., the maximum absolute value of the eigenvalues for the square matrix $P$.

Theorem 3.7** (Local Newton-type contraction [19]).**

We consider the twice continuously differentiable function $\mathcal{F}(z)$ from Eq. (6) and the regular KKT point $\mathcal{F}(z^{\star})=0$ from Definition 3.6. We then apply the Newton-type iteration in Eq. (37), where $\tilde{J}_{\mathrm{IN}}(z)\approx J(z)$ is additionally assumed to be continuously differentiable and invertible in a neighbourhood of the solution. If all eigenvalues of the iteration matrix have a modulus smaller than one, i.e., if the spectral radius satisfies

[TABLE]

then this fixed point $z^{\star}$ is asymptotically stable. Additionally, the iterates $z^{k}$ converge linearly to the KKT point $z^{\star}$ with the asymptotic contraction rate $\kappa^{\star}$ when initialized sufficiently close. On the other hand, the fixed point $z^{\star}$ is unstable if $\kappa^{\star}>1$.

A proof for Theorem 3.7 can be found in [19, 36], based on nonlinear systems theory. Using this result, let us define the linear contraction rate for a Gauss-Newton method with exact Jacobian information

[TABLE]

at the local solution point $z^{\star}:=(w^{\star},\lambda^{\star},\mu^{\star})$ of the KKT conditions. In what follows, we show that the local contraction rate for the block-TR1 Gauss-Newton SQP method

[TABLE]

coincides with the exact Jacobian based linear convergence rate in (39). The following result states that the eigenspectrum of the iteration matrix $\tilde{J}_{\mathrm{IN}}(z^{\star})^{-1}J(z^{\star})-\mathbb{1}$ at the solution point $z^{\star}:=(w^{\star},\lambda^{\star},\mu^{\star})$ coincides with the eigenspectrum of the iteration matrix $\tilde{J}_{\mathrm{GN}}(z^{\star})^{-1}J(z^{\star})-\mathbb{1}$, using the notation $\sigma(P)$ to denote the spectrum, i.e., the set of eigenvalues for a matrix $P$.

Lemma 3.8**.**

For a regular KKT point $z^{\star}:=(w^{\star},\lambda^{\star},\mu^{\star})$, the eigenvalues of the block-TR1 based iteration matrix $\tilde{J}_{\mathrm{IN}}(z^{\star})^{-1}J(z^{\star})-\mathbb{1}$ satisfy

[TABLE]

Proof.

Let us define the eigenvalues $s$ of the iteration matrix $\tilde{J}_{\mathrm{IN}}(z^{\star})^{-1}J(z^{\star})-\mathbb{1}$ as the zeros of

[TABLE]

which, given that the Jacobian approximation $\tilde{J}_{\mathrm{IN}}$ is invertible, this is equivalent to

[TABLE]

This block matrix then reads as

[TABLE]

The result follows from Theorem 3.5 that claims the following asymptotic result for the block-TR1 based Jacobian approximation

[TABLE]

where $N$ is a null space matrix with appropriate dimensions and orthonormal column vectors such that $N^{\top}N=\mathbb{1}$ and $P_{A}\,N=0$. We rewrite Eq. (43) as follows

[TABLE]

It can be verified that $\text{det}\left(J(z^{\star})-(s+1)\tilde{J}_{\mathrm{IN}}(z^{\star})\right)=0$ holds for $s=0$ with an algebraic multiplicity of $2\,n_{A}$ as well as for the values of $s$ that satisfy

[TABLE]

in the limit for $k\rightarrow\infty$, where the compact notation $\Delta H:=\left(\nabla_{w}^{2}{\mathcal{L}}-(s+1)H\right)$ has been used for the Gauss-Newton Hessian approximation. Therefore, the eigenvalues of the iteration matrix $\tilde{J}_{\mathrm{IN}}(z^{\star})^{-1}J(z^{\star})-\mathbb{1}$ for the proposed block-TR1 approach, evaluated at a regular KKT point, are equal to the eigenvalues of the iteration matrix $\tilde{J}_{\mathrm{GN}}(z^{\star})^{-1}J(z^{\star})-\mathbb{1}$ for the exact Jacobian based Gauss-Newton method.

∎

Corollary 3.9**.**

Based on Lemma 41, the linear contraction rate for the block-TR1 based optimization algorithm coincides with the linear contraction rate of the exact Jacobian based Gauss-Newton method $\kappa_{\mathrm{BTR1}}^{\star}=\kappa_{\mathrm{GN}}^{\star}$, when the iterates are sufficiently close to the regular KKT point $z^{\star}:=(w^{\star},\lambda^{\star},\mu^{\star})$.

3.4 Superlinear convergence for SQP with quasi-Newton Hessian updates

Even though the majority of this article is focused on the generalized Gauss-Newton method for nonlinear least squares type optimization problems that occur frequently in predictive control applications, note that superlinear convergence results can be recovered when a block-structure preserving quasi-Newton method is additionally used to approximate the Hessian of the Lagrangian. For example, let us consider the following lemma that represents a block-structured or partitioned version [24, 23] of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) [1] or the symmetric rank-one (SR1) formula [13, 29] to approximate the block-diagonal Hessian matrix.

Theorem 3.10**.**

Given $(AS1)$-$(AS6)$ in Assumption 3.3, then the following holds for each Hessian block matrix approximation

[TABLE]

$i=0,\ldots,N-1$, such that the following holds for the complete Hessian approximation

[TABLE]

Theorem 49 on the convergence of a separable quasi-Newton type Hessian approximation method in combination with our main result in Theorem 3.5 on the block-structured quasi-Newton type Jacobian update formula can be used directly to prove the following result on convergence of the reduced KKT matrix.

Theorem 3.11**.**

Given $(AS1)$-$(AS6)$ in Assumption 3.3, the following holds

[TABLE]

Based on Theorem 3.2, the above result ensures q-superlinear convergence of the SQP iterates when using a quasi-Newton method to update both the block-structured Hessian and Jacobian matrices. The proof for Theorem 50, based on the intermediate convergence results in Theorem 3.5 and 49 can be found in [20].

4 Lifted Collocation Algorithm with Block-TR1 Jacobian Updates

As mentioned earlier, implicit integration schemes are often used in direct optimal control because of their relatively high order of accuracy and their improved numerical stability properties [26]. More specifically, problem formulations based on a system of stiff and/or implicit differential or differential-algebraic equations require the use of an implicit integration scheme. Collocation methods are a popular family of implicit Runge-Kutta methods. This section presents a novel lifted collocation algorithm based on tailored block-TR1 Jacobian updates. The standard lifted collocation method with exact Jacobian information was proposed in [37] as a structure-exploiting implementation of direct collocation, even though it shows similarities to multiple shooting.

4.1 Direct collocation for nonlinear optimal control

In direct transcription methods, such as direct collocation [3, 5], the integration scheme and its intermediate variables are directly made part of the nonlinear optimization problem. In this context, where the simulation routine is defined implicitly as part of the equality constraints in the dynamic optimization problem, one typically relies on implicit integration schemes for their relatively high order of accuracy and improved numerical stability properties. The discrete-time optimal control problem can generally be written as

[TABLE]

where the additional trajectory $K=[K_{0}^{\top},\ldots,K_{N-1}^{\top}]^{\top}$ denotes the intermediate variables of the numerical integration method. These variables are defined implicitly by the equations in (51d), such that the continuity condition reads as in Eq. (51c). More specifically, the Jacobian $\frac{\partial G_{i}}{\partial K_{i}}(\cdot)$ will generally be invertible for an integration scheme applied to a well-defined set of differential equations in (1c). A popular approach of this type is better known as direct collocation [4]. It relies on a collocation method, a subclass of implicit Runge-Kutta (IRK) methods [26], to accurately discretize the continuous time dynamics. In this case, the equations in (51d) define the collocation polynomial on each control interval $i=0,\ldots,N-1$.

In a similar fashion as in Section 2, the adjoint-based SQP method can be applied directly to the direct collocation problem in (51) by solving the following convex QP subproblem at each iteration

[TABLE]

based on $c_{i}^{k}:=G_{i}(w_{i}^{k},K_{i}^{k})$ and the Jacobian approximations $D_{i}^{k}\approx\frac{\partial G_{i}}{\partial w_{i}}(w_{i}^{k},K_{i}^{k})$ and $C_{i}^{k}\approx\frac{\partial G_{i}}{\partial K_{i}}(w_{i}^{k},K_{i}^{k})$. The corresponding gradient correction reads as

[TABLE]

where $\omega_{i}^{k}$ denotes the current value of the multipliers for the nonlinear constraints in (51d) and $\lambda_{i}^{k}$ again denotes the multipliers for the continuity constraints in (51c).

4.2 Tailored structure exploitation for direct collocation

As mentioned earlier, the Jacobian matrix $\frac{\partial G_{i}}{\partial K_{i}}$ for the collocation equations needs to be invertible. Therefore, given an invertible approximation $C_{i}^{k}\approx\frac{\partial G_{i}}{\partial K_{i}}(w_{i}^{k},K_{i}^{k})$, we can rewrite the linearized expression in Eq (52f) as follows

[TABLE]

By substituting the above expression for $\Delta K_{i}$ back into the direct collocation structured QP in (52), one obtains the condensed but equivalent formulation

[TABLE]

where $d_{i}^{k}=e_{i}^{k}-B_{i}\,C^{k^{-1}}_{i}c_{i}^{k}$ is defined and the condensed gradient reads as

[TABLE]

given the original gradient correction in (53).

Note that the resulting QP formulation in Eq. (55) is of the same problem dimensions and exhibits the same sparsity as the multiple shooting structured QP subproblem in Eq. (7). Therefore, state of the art block-structured QP solvers can be used, for which an overview can be found in [21]. After solving the condensed QP in (55), the collocation variables can be obtained from the expansion step in Eq. (54). Based on the optimality conditions of the original direct collocation structured QP in (52), the corresponding Lagrange multipliers can be updated as follows

[TABLE]

where $\lambda^{{k+1}}_{i}$ denote the new values of the Lagrange multipliers for the continuity conditions in (55c) or in (52e).

4.3 Block-TR1 Jacobian update for lifted collocation

The block-TR1 update formula from Eq. (11) can be readily applied to the direct collocation equations, resulting in

[TABLE]

where the quantities $\gamma_{i}^{k^{\top}}=\sigma_{i}^{k^{\top}}\frac{\partial G_{i}}{\partial(w_{i},K_{i})}(w^{{k+1}}_{i},K^{{k+1}}_{i})$ and $\sigma_{i}^{k}=\omega^{{k+1}}_{i}-\omega^{k}_{i}$ are defined. In addition, s_{i}^{k}:=\left[\!\!\begin{array}[]{cccccccccccc}w_{i}^{{k+1}}-w_{i}^{k}\\ K_{i}^{{k+1}}-K_{i}^{k}\end{array}\!\!\right] and $y_{i}^{k}=G_{i}(w^{{k+1}}_{i},K^{{k+1}}_{i})-G_{i}(w^{k}_{i},K^{k}_{i})$ is defined. In order to use this block-TR1 update formula in combination with the lifted collocation method, one needs to be able to efficiently form the condensed QP in Eq. (55). For this purpose, we need to avoid the costly computations of the inverse matrix $C^{k^{-1}}_{i}$ as well as the matrix-matrix multiplication $C^{k^{-1}}_{i}D_{i}^{k}$. In what follows, we present a procedure to directly obtain a rank-one update formula for the inverse matrix $C^{{k+1}^{-1}}_{i}$ and for the corresponding product $E^{{k+1}}_{i}:=C^{{k+1}^{-1}}_{i}D^{{k+1}}_{i}$.

Based on the Sherman-Morrison formula, one can directly update the matrix inverse given the previous invertible approximation $C^{k^{-1}}_{i}\approx\frac{\partial G_{i}}{\partial K_{i}}^{-1}$. Let us first rewrite the block-TR1 update from Eq. (58) as follows

[TABLE]

where $\rho_{i}^{k}=y_{i}^{k}-[D^{k}_{i}\ C^{k}_{i}]s_{i}^{k}$ and $[\tau_{\mathrm{D},i}^{k^{\top}}\ \tau_{\mathrm{C},i}^{k^{\top}}]=\gamma_{i}^{k^{\top}}-\sigma_{i}^{k^{\top}}[D^{k}_{i}\ C^{k}_{i}]$. The Sherman-Morrison formula then reads as

[TABLE]

where $\beta_{i}^{k}=\frac{1}{1\,+\,\alpha_{i}^{k}\tau_{\mathrm{C},i}^{k^{\top}}C^{k^{-1}}_{i}\rho_{i}^{k}}$. Let us define $\tilde{\rho}_{i}^{k}=C^{k^{-1}}_{i}\rho_{i}^{k}$ such that we obtain the following update for the condensed Jacobian

[TABLE]

where $\tilde{\tau}_{i}^{k^{\top}}=\tau_{\mathrm{D},i}^{k^{\top}}-\beta_{i}^{k}\tau_{\mathrm{C},i}^{k^{\top}}(E^{k}_{i}+\alpha_{i}^{k}\tilde{\rho}_{i}^{k}\tau_{\mathrm{D},i}^{k^{\top}})$ has been defined. It is readily seen that the update for $E_{i}^{k}$ in Eq. (61) is a rank-one update for the condensed Jacobian matrix. As proposed in [27], corresponding low-rank update formulas for the condensed Hessian can be obtained for the special case of a pseudospectral method based on a global collocation polynomial.

4.4 Lifted collocation SQP method with block-TR1 Jacobian updates

It is important to stress that the novel block-TR1 update formula for the condensed Jacobian matrix $E^{{k+1}}_{i}=C^{{k+1}^{-1}}_{i}D^{{k+1}}_{i}$ in Eq. (61) provides an efficient manner to directly compute the rank-one update to the matrices in the condensed QP formulation of Eq. (55), without the need for a matrix factorization, inversion and without any matrix-matrix multiplications. Instead, the proposed implementation merely requires matrix-vector multiplications and outer products, resulting in a quadratic instead of cubic computational complexity with respect to the number of optimization variables within each control interval. However, this comes at the cost of a slightly increased memory footprint, since additionally the matrices $C^{-1}_{i}$ and $E_{i}$ need to be stored from one iteration to the next. The implementation of the lifted block-TR1 based SQP method for direct collocation is presented in Algorithm 2.

4.5 Convergence results for block-TR1 based lifted collocation

We observe that the TR1 Jacobian updates of the lifted collocation implementation are equivalent to the updates of the direct collocation method. More specifically, the Jacobian approximation matrices are the same at each SQP iteration, regardless of whether we perform the condensing and expansion procedure for the collocation variables in the proposed lifted implementation of Algorithm 2. Therefore, the convergence properties shown in the previous section also hold for both the standard and lifted collocation based block-TR1 SQP method.

Corollary 4.1**.**

If the assumptions of Theorem 3.1 and Assumption 3.3 hold, then the lifted collocation SQP method with block-wise TR1 Jacobian updates in Algorithm 2, with a Gauss-Newton Hessian approximation, produces iterates $\{w^{k},\lambda^{k},\mu^{k}\}$ that converge q-linearly within a neighbourhood around the KKT point $(w^{*},\lambda^{*},\mu^{*})$ of the NLP.

Proof.

It follows from the equivalence of the SQP iterations between the direct and lifted collocation formulation based on the numerical condensing and expansion of the collocation variables in Eq. (54). In particular, the direct collocation QP subproblem (52) is a special case of the QP formulation in (7), with additional intermediate variables and corresponding equations. The block-TR1 Jacobian matrix convergence results of Theorem 3.5 therefore hold for direct collocation as well as for the proposed lifted implementation in Algorithm 2.

∎

5 Numerical Case Studies of Nonlinear Model Predictive Control

In this section, we illustrate numerically how the proposed block-TR1 SQP method can be used in the context of nonlinear MPC using an algorithm implementation based on the real-time iterations (RTI), as proposed originally in [18] with exact Jacobian information. The approach is based on one block-TR1 SQP iteration per control time step, and using a continuation-based warm starting of the state and control trajectories from one time step to the next [28]. Each iteration consists of two steps:

Preparation phase: discretize and linearize the system dynamics, linearize the remaining constraint functions, and evaluate the quadratic objective approximation to build the optimal control structured QP subproblem. 2. 2.

Feedback phase: solve the QP to update the current values for all optimization variables and obtain the next control input to apply feedback to the system.

The proposed block-wise TR1 based Jacobian updates in Algorithm 1 and 2 become part of the preparation step, in order to construct the linearized continuity equations. Therefore, the feedback step remains unchanged and the Jacobian updates do not affect the computational delay between obtaining the new state estimate and applying the next control input value to the system.

We validate the closed-loop performance of these novel block-TR1 based RTI algorithms by presenting numerical simulation results for two NMPC case studies. Motivated by real embedded control applications, we present the computation times for the proposed NMPC algorithms using the ARM Cortex-A53 processor in the Raspberry Pi 3. The block-sparse QP solution in the feedback phase will be carried out by the primal active-set method, called PRESAS, that was recently presented in [35].

5.1 Nonlinear MPC for a chain of spring-connected masses

In our first case study, the control task is to return a chain of $n_{\mathrm{m}}$ masses connected with springs to its steady state, starting from a perturbed initial configuration, without hitting a wall that is placed close to the equilibrium state configuration. The mass at one end is fixed, while the control input $u(t)\in{\mathbb{R}}^{3}$ to the system is the direct force applied to the mass at the other end of the chain. The state of each free mass $x^{j}:=[p^{j^{\top}},v^{j^{\top}}]^{\top}\in{\mathbb{R}}^{6}$ consists in its position $p^{j}:=[p_{x}^{j},p_{y}^{j},p_{z}^{j}]^{\top}\in{\mathbb{R}}^{3}$ and velocity $v^{j}\in{\mathbb{R}}^{3}$ for $j=1,\ldots,n_{\mathrm{m}}-1$, such that the dynamic system can be described by the concatenated state vector $x(t)\in{\mathbb{R}}^{6(n_{\mathrm{m}}-1)}$. Similar to the work in [37], the nonlinear chain of masses can be used to validate the computational performance and scaling of an optimal control algorithm for a range of numbers of masses $n_{\mathrm{m}}$, resulting in a range of different problem dimensions. The nonlinear system dynamics and the resulting optimal control problem formulation can be found in [39].

5.1.1 Local convergence: Gauss-Newton SQP with block-TR1 Jacobian updates

We illustrate the impact of the proposed block-wise TR1 Jacobian updates on the local convergence rate of the resulting inexact adjoint-based SQP algorithm. Figure 1 shows a comparison of the convergence between different SQP variants for the solution of the nonlinear chain of masses OCP. In particular, the comparison includes the exact Jacobian-based SQP method, the standard dense TR1 update [25], and the good and bad Broyden update formulas [12]. For the proposed block-TR1 based SQP implementation, the figure illustrates both the adjoint and forward variant by using, respectively, the scaling factor in (12) and (13). The performance of the block-TR1 method is additionally illustrated for an implementation where $\alpha_{i}$ is chosen dynamically, depending on which of the two variants results in the largest denominator in order to avoid the need to skip a block-wise Jacobian update.

It is known that an exact Jacobian-based SQP method with Gauss-Newton type Hessian approximation results in locally linear convergence, for which the asymptotic contraction rate depends on the optimal residual value in the least squares type objective [31]. It can be observed in Figure 1 that all three variants of the proposed block-wise TR1 update formula result in the same asymptotic rate of convergence as for the exact Jacobian based algorithm, i.e., the rate of convergence appears to be the same close to the local solution of the NLP. Note that this confirms numerically the result of Corollary 3.9. In addition, the block-wise TR1 Jacobian updates result in a smaller total number of SQP iterations, compared to the standard dense Jacobian update formulas for the particular example in Figure 1. In the latter case, the direct application of a standard rank-one update formula destroys the block sparsity in the QP subproblems and is therefore computationally unattractive.

5.1.2 Computational timing results for block-TR1 based lifted collocation

Figure 1 illustrates the computation times of both the preparation and feedback steps of an NMPC implementation for a chain of $n_{\mathrm{m}}=2,\ldots,8$ masses, using the lifted collocation based SQP method in Algorithm 2. It can be observed that the preparation time scales quadratically with the number of states for the block-TR1 implementation, instead of the cubic computational complexity when using the exact Jacobian. More specifically, the Jacobian evaluation, the factorization and matrix-matrix multiplications are replaced by adjoint differentiation sweeps and matrix-vector operations in Algorithm 2. On the other hand, the feedback time remains essentially the same because, after the linearization and QP preparation, both approaches lead to the solution of a similarly structured QP in Eq. (7) or (55).

Table 5.1.2 provides a more detailed comparison between the exact Jacobian and the proposed block-TR1 variant of the real-time iterations for NMPC, using an ARM Cortex-A53 processor. The table shows these results for both the explicit Runge-Kutta method of order 4 (RK4) in combination with Algorithm 1 and using the implicit 4-stage Gauss-Legendre (GL4) method within Algorithm 2. The proposed block-TR1 algorithm results in a computational speedup of about factor $6-8$ for the problem linearization step. In order to obtain a relatively fair comparison, the number of integration steps for RK4 has been chosen such that the numerical accuracy is close to that of the 4-stage GL method. However, since the system dynamics for the chain of masses are non-stiff, an explicit integration scheme should instead typically perform better in terms of computational efficiency.

5.2 Nonlinear MPC for vehicle control on a snow-covered road

Our second case study considers nonlinear MPC for real-time vehicle control as motivated by automotive applications in autonomous driving. The nonlinear optimal control problem formulation is based on single-track vehicle dynamics with a Pacejka-type tire model [33]. The experimentally validated model parameters can be found in [2]. As often the case in practice, these vehicle dynamics are rather stiff such that an implicit integration scheme should preferably be used. Therefore, it forms an ideal case study for the proposed lifted collocation based RTI method of Algorithm 2. Let us perform the closed-loop NMPC simulations as presented in [33], but using the proposed block-TR1 based RTI implementation. We carried out numerical simulations for two successive double lane changes on snow-covered road conditions. The resulting closed-loop trajectories for both the exact Jacobian and the block-TR1 method are indistinguishable from each other, as illustrated in Figure 3.

The corresponding computation times on the ARM Cortex-A53 processor are illustrated in detail by Table 5.2. Because of the relatively stiff system dynamics, the proposed block-TR1 lifted collocation method from Algorithm 2 becomes attractive and additionally provides a computational speedup of about factor $3$ over the standard exact Jacobian based implementation. Note that, even though the Raspberry Pi 3 is not an embedded processor by itself, it uses an ARM core of the same type as those that are used by multiple high-end automotive microprocessors. Therefore, the proposed algorithm implementation as well as the corresponding numerical results form a motivation for real-time embedded control applications that involve a relatively large, implicit and/or stiff system of differential equations.

6 Conclusions and outlook

In this paper, we proposed a block-wise sparsity preserving two-sided rank-one (TR1) Jacobian update for an adjoint-based inexact SQP method to efficiently solve the nonlinear optimal control problems arising in NMPC. We proved local convergence for the block-structured quasi-Newton type Jacobian matrix updates. In case of a Gauss-Newton based SQP implementation, we additionally showed that the asymptotic rate of contraction remains the same. We also presented how this approach can be implemented efficiently in a tailored lifted collocation framework, in order to avoid matrix factorizations and matrix-matrix multiplications. Finally, we illustrated the local convergence properties as well as the computational complexity results numerically for two nonlinear MPC case studies. The effect of the presented contraction properties on the convergence and closed-loop stability of the block-TR1 based real-time iterations is an important topic that is part of ongoing research.