Global Convergence of Online Optimization for Nonlinear Model Predictive Control

Abstract:
We study a real-time iteration (RTI) scheme for solving online optimization problem appeared in nonlinear optimal control. The proposed RTI scheme modifies the existing RTI-based model predictive control (MPC) algorithm, by selecting the stepsize of each Newton step at each sampling time using a differentiable exact augmented Lagrangian. The scheme can adaptively select the penalty parameters of augmented Lagrangian on the fly, which are shown to be stabilized after certain time periods. We prove under generic assumptions that, by involving stepsize selection instead of always using a full Newton step (like what most of the existing RTIs do), the scheme converges globally: for any initial point, the KKT residuals of the subproblems converge to zero. A key step is to show that augmented Lagrangian keeps decreasing as horizon moves forward. We demonstrate the global convergence behavior of the proposed RTI scheme in a numerical experiment.

heading: 1 Introduction
text: We consider the following time-varying nonlinear optimal control problem
P(x̄0) : min x,u N−1∑ k=0 gk(xk,uk) + gN (xN ),
s.t. xk+1 = fk(xk,uk), k ∈ [N − 1], (1) x0 = x̄0,
where xk ∈ Rnx are the state variables, uk ∈ Rnu are the control variables, gk : Rnx × Rnu → R+ (gN : Rnx → R+) are the nonnegative cost functions, fk : Rnx × Rnu → Rnx are the dynamical constraint functions, x̄0 is the given initial state, and N is the horizon length. We suppose that gk, fk are twice continuously differentiable and, without loss of generality, that the origin is a steady state, i.e. fk(0,0) = 0, and is also a stationary point with zero loss, i.e. gk(0,0) = 0 and ∇gk(0,0) = 0. The well-known linear quadratic regulator (LQR) satisfies this setup and the same setup is commonly used in the literature [8]. Problem (1) is also called dynamic program in control community, which has a close relation to reinforcement learning (RL) (see [6]). The main difference is that the dynamic fk in our paper is known, which is the case in many industrial applications and model-based RL studies.
In modern applications such as energy and autonomous control, N tends to be extremely large or even infinity, which stimulates the interest of solving (1) in real time. Let fk(xk,uk) = fk(xk,uk;dk) where dk is data at stage k, then a more realistic setting is to collect dk as a data stream. That is, instead of knowing f0 all the way to fN−1 and solving (1) offline, we only know fk sequentially and have to solve (1) online. Model predictive control (MPC), also known as receding horizon
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
control, is an important feedback control technique that is used in online optimization of control problems. Given the current system state, MPC looks ahead multiple stages on the prediction horizon, solves the subproblem on the prediction horizon to optimality, and provides the first optimal control action to the system as the feedback. Then the system evolves to a new state based on the control action and the procedure is repeated. However, the classical nonlinear MPC is being regarded as a theoretical concept rather than a practical strategy. It is not as widely used as linear MPC, and is mostly investigated in slow dynamics. This is due to the intensive computations of solving the control subproblems when we have nonlinear models. Fortunately, the design of real-time iteration (RTI) and its predecessors [26, 37] have enabled RTI-based nonlinear MPC to be applied on fast dynamics and achieve good performance [48, 4, 14, 38, 15, 24, 18, 22, 17, 20, 16, 19, 39, 1].
Instead of using optimal feedback control, RTI schemes compute a cheap approximation of control action and provide back to the system as fast as possible. The existing RTI schemes perform a single, full Newton step for each subproblem, and the output iterate of the previous subproblem is used to warmly initialize the current subproblem [8, 41, 42]. Thus, RTI schemes significantly reduce the computations of each subproblem, which makes them succeed in modern industrial applications. RTI schemes exploit the fact that two successive subproblems are closely related though different, so that the iteration can achieve the convergence “on the fly” (i.e. the error goes to zero as horizon moves forward). However, the RTI-based MPC is challenging to analyze. The main difficulty lies in the fact that the prediction horizon is shifted each time. Each subproblem is appended by a new stage that has never been scanned before and may introduce perturbation into the system. Different RTI-based MPC schemes have been studied [9, 12, 23], with an improved local analysis established in [31] recently.
Our paper complements the literature on RTI-based MPC by investigating its global convergence. We emphasize that by global convergence we mean the convergence to stationary points (i.e. KKT points) for any initialization, instead of the convergence to global optimality. The latter is not achievable for nonlinear problems without strong assumptions [34]. Instead of using a full Newton step, we design an adaptive scheme to involve the stepsize selection procedure. In fact, the stepsize selection has been suggested as a way to improve the global behavior of RTI in shrinking (not moving) horizons in an old paper [25]. The authors implemented a `1 penalized merit function1 with a watchdog line search. However, the global convergence theory remains open. Different from the suggestion of [25], we use a differentiable exact augmented Lagrangian as the merit function. This function mainly has two differences to the `1 penalized merit function: (i) it depends on both primal and dual variables, so that dual variables are updated using the same stepsize as the primal variables (`1 merit function only updates the primal variables); (ii) a simple backtracking line search can effectively overcome Maratos effect [29] and accept the unit stepsize locally (`1 merit function requires a more time-consuming line search method such as watchdog or second-order correction). We conduct experiments to show the superiority of augmented Lagrangian over `1 penalized function, especially for time-varying problems. Moreover, our algorithm adaptively selects the penalty parameters of augmented Lagrangian on the fly. When the unit stepsize is accepted, the algorithm simply boils down to the standard RTI-based MPC, and the existing local theories [8, 31] apply seamlessly. We notice that the same merit function has been employed in a recent work [32] to achieve the global convergence for solving (1) in offline fashion, with a horizon decomposition strategy. In that work, the authors showed the unit stepsize is accepted even with an inexact Newton direction. Our paper is in online regime and does not investigate the local behavior.
Related work and contribution: Our work is closely related to a growing literature on online optimization via RTI. The theories of RTI have been investigated over the past decade. [10] and [8] established the local convergence of RTI in shrinking horizons and moving horizons, respectively. When we have different models, [3, 42, 44, 35, 36] all showed similar results. A different view of RTI-based MPC is via the parametric optimization, where subproblems are parameterized by a continuous time index and a single Newton step is performed at each time instant. Similar stability results in this setup are shown in [11, 45, 13]. Recently, [31] studied a lag-L RTI-based MPC (i.e. one shifts L stages each time), and proved that RTI iterates on each stage converge linearly with a rate decaying exponentially in the prediction horizon length. Their analyses rely on the sensitivity analysis of (1) in [30, 33]. An important consequence of convergence of online iteration is the dynamic regret analysis [47, 27, 7, 43]. For example, the convergence rate in [31] and the regret order in [27] are both exponential in prediction horizon length. This is due to the widely-used Lipschitz continuity assumption on the objective.
1The merit function refers to the function that is used in the line search.
However, all above theories are local results. The algorithms always take a full (Newton) step for each subproblem, without involving any stepsize selection procedure. Thus, the analyses require to have an initial point that is sufficiently close to the local solution, or require that the first Newton step is sufficiently small. Obviously, for Problem (1), when ‖x̄0‖ is small enough, the origin is a good initial point and RTI iterates stably track the origin. However, if ‖x̄0‖ is large or we may want to track a nontrivial MPC policy, the origin is not a reasonable choice for the initial point. Based on these considerations, the global convergence of RTI addressed in this paper complements the aforementioned rich local analyses. Again, by global convergence we mean that the KKT residuals of subproblems go to zero as horizon moves forward for any initialization. The known, local stability of tracking error of RTI is not addressed in this paper, but referred to [45, Theorem 4.2].
The challenge of our analysis is that: the problem varies with time instants and only a single Newton step is performed at each time instant. Such a setup contrasts to online Newton (with line search) analysis, where a single Newton step targeting the same problem is performed. A key, technical step of our analysis is to show that the exact augmented Lagrangian merit function keeps decreasing as horizons shift. To our knowledge, this is the first paper studying the global behavior of online iteration under RTI setting, via a tool of the exact augmented Lagrangian. Another contribution to nonlinear MPC is that we consider time-varying dynamics, while most of the literature considered time-invariant dynamics (e.g. [8, 21, 40]). The subproblems of time-invariant dynamics only differ in the initial state, which results in a much simpler analysis. We notice that [31] allowed time-varying dynamics, but required to use a larger lag when shifting the horizon. We allow lag one instead.
Notations: For two vectors a, b, (a; b) denotes a long column vector that stacks a and b together. For two scalars a, b, a ∨ b = max(a, b). Without specification, ‖ · ‖ denotes the `2 norm for vectors and operator norm for matrices. We also reserve the following notations. 0 denotes zero vector or matrix, and I denotes the identity matrix. Their dimensions are clear from the context. M denotes the prediction horizon length, t indexes the sampling time, and zk = (xk;uk) ∈ Rnz with nz = nx+nu denotes the state-control pair at stage k.

heading: 2 Preliminaries
text: We introduce MPC subproblems formulation, the Newton step, and the augmented Lagrangian merit function in this section. At time t, the prediction horizon shifts to [t,Mt] with Mt = t+M denoting the last stage, and the corresponding subproblem Pt(x̄t) is given by
Pt(x̄t) : minxt:Mt ut:Mt−1 Mt−1∑ k=t gk(xk,uk) + gMt(xMt ,0) + µ 2 ‖xMt‖2,
s.t. xk+1 = fk(xk,uk), k ∈ [t,Mt − 1], (2) xt = x̄t.
Here, µ is a quadratic penalty parameter independent from t; x̄t is the current system state, which is given and determined by the control feedback provided to the system before t. The formulation (2) is used in [31] for local analysis, while alternative subproblem formulations are available as well [8, 23].
Given subproblem (2), RTI schemes perform a single Newton step for Pt(x̄t) and move to the next time instant. To ease notations, we let x̃t = xt:Mt and ũt = ut:Mt−1 denote the state and control variables of the t-th subproblem, respectively. We also let z̃t = zt:Mt = (zt; . . . ; zMt) denote the whole ordered state-control variables, with slightly abusing the notation zMt to let zMt = xMt . We sometimes also express z̃t = (x̃t, ũt) when we specify each component of z̃t. Then, the Lagrangian function of (2) is
Lt(z̃t, λ̃t; x̄t) = Mt−1∑ k=t gk(xk,uk) + gMt(xMt ,0) + µ 2 ‖xMt‖2
+ Mt−1∑ k=t λTk (xk+1 − fk(xk,uk)) + λ T t−1(xt − x̄t),
where λ̃t = λt−1:Mt−1 is the dual vector. Here, λk ∈ Rnx is associated to the k-th constraint of (2) for k ∈ [t,Mt − 1], and λt−1 ∈ Rnx is associated to the initial state constraint.
Given the input iterate (z̃0t , λ̃ 0 t ), the Newton direction (∆z̃t,∆λ̃t) is solved by the KKT system( Bt (Gt)T
Gt 0 )( ∆z̃t ∆λ̃t ) = − ( ∇z̃tLt ∇λ̃tL t ) , (3)
where Bt is an approximation of the Hessian Ht = ∇2z̃tL t; Gt = ∇λ̃tz̃tL t is the Jacobian matrix; and all matrices in (3) are evaluated at (z̃0t , λ̃ 0
t ). The left hand side matrix is called the KKT matrix. In what follows, Ht, Bt, Gt are always evaluated at the input iterate (z̃0t , λ̃ 0
t ) if no evaluation point is specified. Taking a closer look at (3) (cf. Appendix A), we note that the KKT matrix and ∇z̃tLt are independent from the given initial state x̄t. This is a crucial observation. In reality, the realization of the system state may lag behind the sampling time, since the system has to apply the feedback control to evolve. Thus, x̄t may not be known as soon as the input iterate (z̃0t , λ̃ 0
t ). However, a promising property of RTI is that most of computations of (3) can be prepared without knowing x̄t.
Given the Newton direction (∆z̃t,∆λ̃t) from (3), we update the iterate by involving a stepsize αt( z̃1t λ̃ 1
t
) = ( z̃0t λ̃ 0
t
) + αt ( ∆z̃t ∆λ̃t ) . (4)
In our paper, αt is selected by forcing Armijo condition hold on a differentiable exact augmented Lagrangian, defined as
Ltη(z̃t, λ̃t; x̄t) =Lt(z̃t, λ̃t; x̄t) + η1 2 ∥∥∥∇λ̃tLt(z̃t, λ̃t; x̄t)∥∥∥2 + η22 ∥∥∥∇z̃tLt(z̃t, λ̃t; x̄t)∥∥∥2 . (5) Here η = (η1, η2) are parameters of the quadratic penalties. (5) is called exact augmented Lagrangian since one can show that, for sufficiently large η1 and sufficiently small η2, a strict local minimizer of Ltη(·) is also a strict local minimizer of Problem (2) and vice versa (see [2, Proposition 4.5]). When η2 = 0, Ltη(·) reduces to standard augmented Lagrangian, while it does not enjoy exactness property anymore, and may cause issues when it is used as merit function. For simplicity, we sometimes suppress “exact” and call Ltη(·) augmented Lagrangian. Based on (5), αt is selected by making the following Armijo condition,
Ltη(z̃1t ,λ̃ 1 t ; x̄t) ≤ Ltη(z̃0t , λ̃ 0 t ; x̄t) + αtβ
( ∇z̃tLtη(z̃ 0 t , λ̃ 0 t ; x̄t)
∇λ̃tL t η(z̃ 0 t , λ̃
0 t ; x̄t)
)T ( ∆z̃t ∆λ̃t ) , (6)
hold for a prespecified parameter β ∈ (0, 1). Note that, by mean value theory, αt can always be found by backtracking line search provided the inner product in (6) is negative (cf. Theorem 4.5).

heading: 3 An Adaptive RTI-based MPC Scheme
text: We now set the stage to propose our algorithm. The algorithm inherits the flavor of RTI schemes: it performs a single Newton step for each subproblem. Differently, it involves a stepsize selection procedure based on Ltη(·) instead of using a full step (i.e. αt = 1). The algorithm adaptively selects the parameter η on the fly (i.e. it automatically selects suitable η for each subproblem). However, since µ is required in the problem formulation, the adaptivity on µ may require to solve the KKT system multiple times for each subproblem. Thus, the gain from the adaptivity on µ for our procedure is limited. We hence prefer to tune µ manually in practice. Our experiments show that our algorithm is robust in µ (see Figure 1).
Let us now illustrate how we transit from the triple (z̃0t , λ̃ 0 t ; x̄t) to (z̃ 0 t+1, λ̃ 0
t+1; x̄t+1). To further ease notations, we let Lt,1η , Lt,0η , Lt,1, Lt,0 denote the function values of Ltη(·) and Lt(·) at the output iterate (z̃1t , λ̃ 1 t ) and the input iterate (z̃ 0 t , λ̃ 0 t ), respectively. In addition, when we specify the iterate on each stage, we let x̃Idt = x Id t:Mt,t for Id = 0 or 1 (similar for z̃Idt , λ̃ Id t , ∆z̃t, and ∆λ̃t). For example, x1t+1,t denotes the output iterate of the t-th subproblem at stage t+ 1.
Given the input iterate (z̃0t , λ̃ 0
t ), we perform four steps: (1) We computeBt,Gt,∇z̃tLt,0, and∇λ̃tL t,0 except the first row. After we observe x̄t, we compute
Algorithm 1 An Adaptive RTI-based MPC Scheme
1: Input: initial iterate (z̃00, λ̃ 0 0), initial state x̄0, and parameters µ > 0, ρ > 1, η 0 = (η01 , η 0 2) =
(µ2, 1), and β ∈ (0, 1); 2: for t = 0, 1, 2, . . . do 3: Compute Bt, Gt,∇z̃tLt,0,∇λ̃tL t,0, and solve the KKT system (3) to obtain (∆z̃t,∆λ̃t);
4: while ( ∇z̃tL t,0 ηt
∇λ̃tL t,0 ηt
)T ( ∆z̃t ∆λ̃t ) > −η t 2 4 ∥∥∥∥(∇z̃tLt,0∇λ̃tLt,0 )∥∥∥∥2 do
5: Let ηt1 = ρ 2ηt1, η t 2 = η t 2/ρ; 6: end while 7: Select αt by backtracking line search to make (6) hold, and update iterates by (4); 8: Define (z̃0t+1, λ̃ 0 t+1) and x̄t+1 according to (8), and let η t+1 = ηt; 9: end for
the first row of∇λ̃tL t,0. Then, we compute the Newton direction (∆z̃t,∆λ̃t) by solving (3). (2) We select ηt = (ηt1, ηt2) such that( ∇z̃tL t,0 ηt
∇λ̃tL t,0 ηt
)T ( ∆z̃t ∆λ̃t ) ≤ −η t 2 4 ∥∥∥∥(∇z̃tLt,0∇λ̃tLt,0 )∥∥∥∥2 . (7)
As shown in Algorithm 1, we can use a While loop to increase ηt1 and decrease η t 2 until (7) is satisfied. (3) Using ηt from the last step, we select αt by backtracking line search to make Armijo condition (6) hold. We then obtain the output iterate (z̃1t , λ̃ 1
t ) by (4). (4) We transit to the next horizon by letting
x̃0t+1 =x 0 t+1:Mt+1,t+1 := (x1t+1,t;x 1 t+2,t; . . . ;x 1 Mt,t;0),
ũ0t+1 =u 0 t+1:Mt+1−1,t+1 := (u 1 t+1,t;u 1 t+2,t; . . . ;u 1 Mt−1,t;0),
λ̃ 0
t+1 =λ 0 t:Mt+1−1,t+1 := (λ 1 t,t;λ 1 t+1,t; . . . ;λ 1 Mt−1,t;0),
x̄t+1 :=ft(x 1 t,t,u 1 t,t).
(8)
The above four steps are displayed in Algorithm 1. If we remove Step (2) and use αt = 1 in Step (3), it reduces to standard RTI schemes [8, 41, 42, 31]. However, as shown later, the difference that we select αt based on augmented Lagrangian Ltηt(·) enables Algorithm 1 to converge globally. Remark 3.1. We discuss the overhead of the algorithm. For each subproblem, the computational complexity is dominated by solving the KKT system. The KKT matrix is a square matrix with dimensions 2(M + 1)nx + Mnu. Directly inverting the matrix results in O((M(nx + nu))3) flops. Fortunately, we can exploit the sparse structure of the KKT matrix and apply sparse LDLT factorization as did in Section C of [41], which results in O(M(nx + nu)3) flops instead. So the computations are linear in terms of the prediction horizon length M .

heading: 4 Global Convergence Analysis
text: We conduct global convergence analysis for Algorithm 1. We show that, for any initial iterate (z̃00, λ̃ 0
0), the KKT residual of the t-th subproblem, ‖∇Lt,0‖, converges to zero as t increases. The significance of the result is that we only perform a single Newton step for each time-varying subproblem, and our result complements the local results of RTI. We should mention that, like most of deterministic methods for nonlinear problems [34], Algorithm 1 is not selective to stationary points, and may converge to saddle points. Converging to local minimum is only possible under stronger, local assumptions, which are not addressed in this paper but we refer to [31].

heading: 4.1 Assumptions
text: We require the following assumptions.
Assumption 4.1 (Compactness). We assume all iterates lie in a convex compact set that contains the origin. In particular, we assume there exist convex compact sets 0 ∈ Z ⊆ Rnx+nu and 0 ∈ Λ ⊆ Rnx , such that (z0k,t,λ 0 k,t) ∈ Z × Λ for any t ≥ 0 and any stage k of the t-th subproblem.
Assumption 4.2 (Boundedness condition). We assume there exist absolute constants Υ, δ, γRH > 0 such that (i) for any t ≥ 0, Bt satisfies zTBtz ≥ γRH‖z‖2, ∀z ∈ {z : Gtz = 0, z 6= 0}, and ‖Bt −Ht‖ ≤ δ; (ii) for any k ≥ 0 and zk ∈ Z , we have
‖∇xkfk(zk)‖ ∨ ‖∇ukfk(zk)‖ ∨ ‖∇2zkgk(zk)‖ ∨ {‖∇ 2 zk fk,l(zk)‖}nxl=1 ≤ Υ, (9) {‖∇2zk [∇zkgk(zk)]i‖, ‖∇ 2 zk [∇zkfk,l(zk)]i‖} nz,nx i=1,l=1 ≤ Υ; (10)
(iii) for any t ≥ 0, Gt(Gt)T I/Υ2. Assumption 4.3. We assume there exist absolute constants c, C > 0 such that for any t ≥ 0 (i) ‖∆λ̃Mt−1,t‖ ≤ c‖(∆z̃t; ∆λ̃t)‖ if x1Mt,t 6= 0. (ii) ‖λ 1 t−1,t‖ ≤ C‖x1t,t − x̄t‖ if x1t,t 6= x̄t.
A 4.1 and 4.2 are standard in the literature on Newton analysis [34, 2]. We note that Jacobian Gt always has full row rank (see (13) in Appendix A), so the linear independence constraint qualification (LICQ) holds. Together with LICQ, A 4.2(i) implies the KKT matrix is nonsingular [34, Lemma 16.1]. In fact, we can simply let Bt = µI , which satisfies all conditions trivially; is sufficient for our global analysis; and is indeed used in our experiments. More generally, Bt can approximate the Hessian Ht using any methods in [34]. Different from existing analyses which only require the bounded Hessian, (10) requires one more derivative—the third derivative—to be bounded. This is standard when using augmented Lagrangian Ltη(·) in the algorithm [46, 2, 32], because the Hessian of Ltη(·) requires the third derivatives of gk, fk. We mention that the third derivatives are only used in the analysis, and do not have to be computed in the implementation. (9) is directly implied by the compactness in A 4.1. A 4.2(iii) is also reasonable to hold since Gt has full row rank. For time-invariant problems where gk = g, fk = f , A 4.2(ii) and (iii) hold trivially.
A 4.3 is particularly imposed for MPC analysis; a condition with the same flavor is in [8]. At the first glance, it seems to be a restrictive assumption. Fortunately, it is guaranteed to hold when µ is sufficiently large (the threshold depends on problem parameters), and is easily checkable during the iteration. In particular, A 4.3(i) assumes that the last dual step ∆λ̃Mt−1,t is bounded by the full Newton step (∆z̃t,∆λ̃t). It holds trivially for c = 1, however, we hope c is small. One can expect this condition to be satisfied for large enough µ, since ∆λ̃Mt−1,t measures the disturbance of x1Mt,t away from x 0 Mt,t
= 0. When µ is large, we know from the quadratic penalty in (2) that x1Mt,t ≈ 0. Hence, ∆λ̃Mt−1,t becomes arbitrarily small. See more discussions in [8, Assumption 3] for local analysis of MPC. A 4.3(ii) allows any constant C as long as C ≤ µ2, so it holds as well for large enough µ. In addition, A 4.3(ii) is not required if αt = 1 is selected, which is always true when (z̃0t , λ̃ 0 t ) is close to a stationary point of the subproblem 2. This is because that, if αt = 1, then x1t,t = x̄t for any input iterate x 0 t,t.

heading: 4.2 Global convergence
text: We first show that Algorithm 1 is well-posed, that is, each iteration can be performed successfully. The computation of the Newton direction in Line 3 is ensured by the nonsingularity of the KKT matrix, implied by A 4.2 as discussed before. The backtracking line search in Line 7 is guaranteed by the mean value theorem [5]. It hence suffices to show that the While loop in Line 4 terminates in finite time. Theorem 4.4. Consider Algorithm 1 under Assumptions 4.1 and 4.2. Let us define thresholds τ1, τ2, τ3 > 0 by
τ1 = 25µ2Υ2
γRH , τ2 = γRH 2δ2 , τ3 = 32µ2Υ2 γ2RH .
For any t ≥ 0, if η1 ≥ τ1, η2 ≤ τ2, η1η2 ≥ τ3, then( ∇z̃tLt,0η ∇λ̃tL t,0 η )T ( ∆z̃t ∆λ̃t ) ≤ −η2 4 ∥∥∥∥(∇z̃tLt,0∇λ̃tLt,0 )∥∥∥∥2 .
2[28] showed under certain conditions that the backtracking line search would always select a unit stepsize locally when using exact augmented Lagrangian as the merit function.
Further, there exist an integer τ > 0 and positive constants η̄ = (η̄1, η̄2) satisfying
η̄1 ≤ 322ρ2µ2Υ4
γ4RH , η̄2 ≥ γRH 32ρΥ2 ,
such that ηt = η̄, ∀t ≥ τ .
Theorem 4.4 shows that the While loop indeed terminates in finite time, and more importantly, the penalty parameter ηt stabilizes after certain iterations. This result is critical for global convergence since, after τ subproblems, all the succeeding subproblems try to decrease a much similar merit function. We note that η̄1 depends on µ, so the penalty of the objective connects to the penalty of augmented Lagrangian. We mention that the constants in {τi}3i=1 and η̄1, η̄2 can be significantly improved by a careful calculation. However, we leave it for future work since it’s not main target here.
It is also worth mentioning that the stabilized parameters (η̄1, η̄2) may not ensure the equivalence between stationary points of Ltη̄(·) and KKT points of Pt(x̄t), because the stabilized values rely on particular iteration sequence. The parameters thresholds for the exactness property of Ltη ensure that any stationary points of augmented Lagrangian are equivalent to KKT points, which is a stronger guarantee. However, the stabilized parameters indeed ensure that our particular iteration converges to a KKT point, not just a stationary point of augmented Lagrangian.
We now study how the merit function decreases from one subproblem to another. Theorem 4.5. Consider Algorithm 1 under Assumptions 4.1 and 4.2. There exists a constant ᾱ > 0 such that αt ≥ ᾱ, ∀t ≥ τ where τ is from Theorem 4.4. Moreover,
Lt,1η̄ ≤ L t,0 η̄ −
η̄2ᾱβ
4 ‖∇Lt,0‖2.
Theorem 4.5 shows that Armijo condition is satisfied for sufficiently small αt, with a threshold independent of t. Moreover, after one Newton step, we can observe a sufficient decrease of the merit function. However, since we only perform one Newton step and transit to the next subproblem immediately, we cannot claim the decrease is sustained. This is indeed the main challenge of global analysis of RTI schemes. Different subproblems target different merit functions, even η is shared.
We show in the next lemma that, the transition between subproblems may increase the merit function, but not too much. Overall, the merit function Lt,0η̄ keeps decreasing. Lemma 4.6. Under Assumptions 4.1-4.3 and let κ = µγRH denote the condition number of the Hessian matrix. Suppose the constants c, C in Assumption 4.3 satisfies c . γRH/κ3 and C . µ2. Then,3
Lt+1,0η̄ ≤ L t,1 η̄ +
η̄2ᾱβ
8 ‖∇Lt,0‖2.
Combining Lemma 4.6 with Theorem 4.5, we immediately have
one-step error recursion: Lt+1,0η̄ ≤ L t,0 η̄ −
η̄2ᾱβ
8 ‖∇Lt,0‖2. (11)
Using (11), we sum over t and finally establish the global convergence in the next theorem. Theorem 4.7. Consider Algorithm 1 under Assumptions 4.1-4.3. Suppose c . γRH/κ3 and C . µ2, then for any input iterate (z̃00, λ̃ 0 0), we have ‖∇Lt,0‖ → 0 as t→∞.

heading: 5 Numerical Experiment
text: We apply Algorithm 1 on a dynamic program to demonstrate the global behavior of RTI-based MPC. We compare augmented Lagrangian merit function with `1 penalized merit function, defined as
L̄tν(z̃t; x̄t) = Mt−1∑ k=t gk(xk,uk) + gMt(xMt ,0) + µ 2 ‖xMt‖ 2
+ ν ( Mt−1∑ k=t ‖xk+1 − fk(xk,uk)‖1 + ‖xt − x̄t‖1 ) .
3“.” means the inequality holds up to a problem-independent constant.
Using L̄tν(·) as the merit function, Armijo condition is rewritten as L̄t,1ν ≤ L̄t,0ν +αtβ ·D(L̄t,0ν ,∆z̃t), where D(L̄t,0ν ,∆z̃t) is the directional derivative of L̄tν(·) at z̃0t along the direction ∆z̃t, given by
D(L̄t,0ν ,∆z̃t) = Mt−1∑ k=t ∇T gk(z0k,t)∆zk,t +∇TxMt gMt(x 0 Mt,t,0)∆xMt,t
− ν ( Mt−1∑ k=t ‖xk+1 − fk(xk,uk)‖1 + ‖xt − x̄t‖1 ) =: D1 − νD2.
To ensure D(L̄t,0ν ,∆z̃t) to be negative, we let νt = max(D1/ {(ρ− 1)D2} , νt−1) with ρ ∈ (1, 2). Thus, RTI-based MPC using `1 penalized merit function (called `1 MPC) is obtained by replacing Lines 4 and 7 of Algorithm 1 by the above corresponding steps.
Simulation setups. We consider 1D trigonometric perturbed LQR problem. In particular, we let gk(zk) = akx 2 k+bku 2 k+ck sin
2(xk) and fk = xk+uk+sin(xk). We study three cases summarized in Table 1. For each case, we perform 1000 independent runs with randomly generalized initial iterate (z̃00, λ̃ 0 0), by letting (x 0 k,0,u 0 k,0,λ 0 k,0) ∼ N (0, 25I), ∀k. We stop the iteration if either t > N −M (i.e. attains the iteration threshold) or ‖∇Lt,0‖ ≤ = 10−8 (i.e. attains the error threshold). We let Bt = µI , ∀t. For each case, we try three prediction horizon lengths, M = 5, 10, 15. We also vary µ in different cases to test the robustness of Algorithm 1 on µ. For Case 1, we let ak = bk = ck = 1, which results in a time-invariant problem. For Cases 2 and 3, we let ak = bk = ck = k and k2, which results in time-varying problems. The linear and quadratic increase of coefficients bring significant challenges for solving subproblems online. The code is implemented in Julia 1.5.4 and is publicly available (with high resolution figures) at https://github.com/senna1128/Global-RTI-MPC.
Evaluation of results. For each case, we compare the average running time between two merit functions, and compare the maximum (over 1000 independent runs) of the number of iterations. The comparisons are shown in Table 2. From Table 2, we see that two merit functions work equally well for the time-invariant problem in Case 1, while augmented Lagrangian works significantly better than `1 MPC for the time-varying problems in Cases 2 and 3. The latter fails to attain 10−8 KKT residuals within 250 iterations in both Case 2 and Case 3.
Moreover, we plot ‖∇Lt,0‖ v.s. t in log scale for all cases. To have a better visualization, among 1000 initializations, we only plot the first 5 in Figure 1, but draw an error bar plot for all initializations in Figure 2. We see from Figure 1 (a)-(f) and Figure 2 (a)-(f) that both MPC schemes achieve global convergence within 25 iterations for the time-invariant problem in Case 1. Their (linear) convergence speed is comparable as well. However, for the time-varying problems in Cases 2 and 3, Figure 1 (g)-(r) and Figure 2 (g)-(r) suggest that augmented Lagrangian MPC performs much better than `1 MPC. Although both schemes decrease KKT residuals as t increases, which validates our global convergence theory, augmented Lagrangian MPC decreases the residuals linearly in t, and below the threshold = 10−8 within a moderate number of iterations. As a comparison, `1 MPC decreases the residuals only to [10−2, 10−4] using 250 iterations. Moreover, it is clear to see that `1 MPC has a slow convergence speed locally (because the tail of ‖∇Lt,0‖ is flat), which does not appear in augmented Lagrangian MPC. Therefore, for time-varying problems, the proposed augmented Lagrangian MPC is much preferable than `1 MPC.
We also take Case 2 as an example to show our algorithm is robust to µ. We vary µ = 1, 5, 10, 100, 500, 1000 and randomly generate an initial point for each µ. We see from Figure 1 (s)-(x) that the performance of our algorithm is stable for different choices of µ.

heading: 6 Conclusion
text: We studied the global behavior of RTI-based nonlinear MPC algorithm. Different from the existing schemes, we imposed line search for each subproblem to select the stepsize, using a differentiable exact augmented Lagrangian as the merit function. We proved that the proposed algorithm enjoys global convergence: for any initial iterate, the KKT residuals of subproblems converge to zero. Our analysis complements the existing local analyses of RTI-based MPC, which require a sufficiently good initial point.
One extension is to replace the Newton step by the quasi-Newton step. This may lead to a slow local convergence, however, the global convergence should be preserved. Another important extension is to study stochastic dynamics, like [40]. We can evaluate the iteration by either regret [27] or E‖∇Lt,0‖. Then even the local stability property of RTI-based MPC is not well understood. Furthermore, studying RTI schemes for model-based RL is an attractive research direction and, to some extent, our analysis is still applicable. However, how to apply RTI on model-free RL is still not clear to our knowledge. Finally, our experiments show that RTI-based MPC converge globally even for moderate µ (e.g. µ = 1), which in some sense suggests Assumption 4.3 is not necessary. How to remove Assumption 4.3 is worth investigating.

Review summaries:
This paper gives an algorithm for nonlinear MPC based on augmented Lagrangians, and demonstrates that the iterates converge globally to a critical point of the KKT equations. The key insight is a rule for tuning the regularization parameters of the augmented Lagrangian that ensures that each algorithm step reduces the KKT residuals sufficiently to overcome any increase that would arise due to taking a step in the dynamics. As a result, after alternating between algorithm steps and dynamics steps for several iterations, all iterates produced by the algorithm will be near critical points. 

This paper studies a real-time iteration scheme for solving nonlinear MPC. The paper provides an adaptive way of Newton stepsize selection and establishes global convergence. The paper also provides abundant numerical results.

This  paper proposes an adaptive extension of a real-time iteration (RTI)  scheme used to solve online non linear optimal control problems, with  time-varying dynamics.  It introduces a parametric merit function, named exact augmented  Lagrangian, which is used to adapt the step-size of each unique Newton  steps involved to solve the succession of MPC sub-problems. Such  step-size follows an Armijo-type condition. Before each back-tracking  line search, the penalty parameters of this merit function are adjusted  according to a well-defined scheme.  Under generic assumptions, authors prove the global convergence of their  approach. That is, the KKT residual of t-th sub-problem converges to  zero as t, the time horizon, increases. Furthermore, for large t,  penalty parameters of the merit function used are no longer changed.   Finally, the paper proposes two types of numerical experiments. First,  it compares the performances (in terms of speed and precision) of their  approach with the one using a $l_1$ merit function to perform a  backtracking scheme. Secondly, it analyses the robustness of their  assumptions with respect to a quadratic penalty parameter $\mu$, which  is tuned manually in their approach.   Numerical results show indeed that this adaptive RTI has a greater  performance, especially for time-varying problems, and is robust to  $\mu$ parameter choice.