Time Distributed Optimization for Model Predictive Control: Stability, Robustness, and Constraint Satisfaction

TL;DR

This paper introduces a systems theoretic framework for analyzing time distributed optimization in model predictive control, focusing on stability, robustness, and constraint satisfaction, and extends theoretical results for real-time iteration schemes.

Contribution

It provides a general stability analysis framework for time distributed optimization in MPC, extending existing theory for real-time iteration methods.

Findings

The framework guarantees stability and constraint satisfaction under certain conditions.

Numerical simulations validate the effectiveness of the proposed approach.

The analysis enhances understanding of robustness in time distributed optimization schemes.

Abstract

Time distributed optimization is an implementation strategy that can significantly reduce the computational burden of model predictive control by exploiting its robustness to incomplete optimization. When using this strategy, optimization iterations are distributed over time by maintaining a running solution estimate for the optimal control problem and updating it at each sampling instant. The resulting controller can be viewed as a dynamic compensator which is placed in closed-loop with the plant. This paper presents a general systems theoretic analysis framework for time distributed optimization. The coupled plant-optimizer system is analyzed using input-to-state stability concepts and sufficient conditions for stability and constraint satisfaction are derived. When applied to time distributed sequential quadratic programming, the framework significantly extends the existing theoretical analysis for the real-time iteration scheme. Numerical simulations are presented that demonstrate the effectiveness of the scheme.

Tables

Table 1. Table 1: Bicycle Model Parameters

Name	Symbol	Value
Mass	$m$	$2041kg$
Yaw Inertia	$I_{zz}$	$4964kg{m}^2$
Front, Rear CG distance	${\ell }_f,{\ell }_r$	$1.56,1.64m$
Coefficient of friction	$\mu$	$0.8$
Tire parameters	$B,C$	$12,1.285$
Lateral Area	$A$	$7.8m^2$
Air Density	$\rho$	$1.225kg/m^2$
Lateral Drag Coefficient	$C_d$	1.5
Longitudinal Velocity	$s$	$30m/s$

Equations

$$x_{k+1}=g(x_k,u_k),$$

$${\mathcal{N}}_C(v)=\{\begin{array}{ll}\{y|y^T(w-v)\le 0,\forall w\in C\},& \text{if}v\in C,\\ \varnothing & \text{else},\end{array}$$

$$A\pm B=\{y|y=a\pm b,a\in A,b\in B\}.$$

$$||x(k,x_0,\mathbf{u})||\le \max \{\beta (||x_0||,k),\gamma (||\mathbf{u}||)\},$$

$$\underset{k\to \infty }{\overline{\lim }}||x(k,x_0,\mathbf{u})||\le \gamma \left(\underset{k\to \infty }{\overline{\lim }}||u_k||\right),$$

$$x_{k+1}=f_d(x_k,u_k,d_k),$$

$$\underset{\xi ,\mu }{min.}$$

$$s.t.$$

$${\xi }_0=x,{\xi }_N\in {\mathcal{X}}_f,$$

$$({\xi }_i,{\mu }_i)\in \mathcal{Z},i=0,\dots ,N-1.$$

$$\underset{u}{\min }\{V_f(x^{+})-V_f(x)+l(x,u)|(x,u)\in \mathcal{Z},x^{+}\in {\mathcal{X}}_f\}\le 0,$$

$$\Gamma =\{x\in \mathcal{X}|\text{eqrefeq:OCP is feasible}\},$$

$$\kappa (x)={\mu }_0^{\ast }(x),$$

$$x_{k+1}=f(x_k,\Delta u_k,d_k):=f_d(x_k,\kappa (x_k)+\Delta u_k,d_k),$$

$$z_k={\mathcal{T}}_{\ell }(z_{k-1},x_k),$$

$$F(z,x)+{\mathcal{N}}_K(z)\ni 0,$$

$$S(x)=\{z|F(z,x)+{\mathcal{N}}_K(z)\ni 0\},$$

$$z_{i+1}=\mathcal{T}(z_i,x,i),$$

$${\mathcal{T}}_{\ell }(z,x)=\mathcal{T}({\mathcal{T}}_{\ell -1}(z,x,\ell -1),x,\ell ),$$

$$\underset{\ell \to \infty }{\lim }{\mathcal{T}}_{\ell }(z,x)=z^{\ast }$$

$$||{\mathcal{T}}_{\ell }(z,x)-z^{\ast }||\le \eta ||{\mathcal{T}}_{\ell -1}(z,x)-z^{\ast }|{|}^q,$$

$$z_k={\mathcal{T}}_{\ell }(z_{k-1},x_k),$$

$$u_k=\Xi z_k,$$

$$\bar{s}(x)\in S(x),$$

$$e_{k+1}={\mathcal{G}}_{\ell }(e_k,x_k,\Delta x_k),$$

$$\Delta u_k=\Xi e_k,$$

$${\mathcal{G}}_{\ell }(e_k,x_k,\Delta x_k)={\mathcal{T}}_{\ell }(e_k+\bar{s}(x_k),x_k+\Delta x_k)-\bar{s}(x_k+\Delta x_k).$$

$$||e_{k+1}||\le a(\ell )||e_k||+\theta (\ell )||\Delta x_k||,$$

$$||e_k||+b||\Delta x_k||\le \epsilon ,$$

$$||z_{k+1}-\bar{s}(x)||\le {\eta }^{\alpha (\ell )}||z_k-\bar{s}(x)|{|}^{q^{\ell }},$$

$$||z_{k+1}-$$

$$||e_{k+1}||$$

$$\le {\eta }^{\alpha (\ell )}||[z_k-\bar{s}(x_k)]-[\bar{s}(x_{k+1})-\bar{s}(x_k)]|{|}^{q^{\ell }},$$

$$\le {\eta }^{\alpha (\ell )}(||e_k||+b||\Delta x_k||{)}^{q^{\ell }},$$

Full text

††thanks: This research is supported by the Toyota Research Instituite (TRI) and by the National Science Foundation through awards CMMI 1904441 and CMMI 1562209. TRI provided funds to assist the authors with their research but this article solely reflects the opinions and conclusions of its authors and not TRI or any other Toyota entity.

, ,

1 Introduction

Model Predictive Control (MPC) [22] is a widely used control technique that computes control actions by solving an Optimal Control Problem (OCP) over a finite receding horizon. MPC can systematically handle constraints and nonlinearities but is challenging to implement since it requires the solution of a constrained and potentially non-convex OCP at each sampling instant. The development of robust and efficient quadratic and convex programming solvers, see e.g., [34, 42, 11, 43], has enabled the application of linear-quadratic MPC to a wide variety of systems. However, the implementation of MPC for systems with limited onboard computing power, fast sampling rates, and/or pronounced nonlinear dynamics remains an open problem.

One approach for reducing the computational cost of MPC is time distributed optimization (TDO). TDO distributes optimizer iterations over time by exploiting the robustness of MPC to suboptimality [2, 41, 35]. Rather than accurately solving the OCP at each sampling instant, TDO maintains a guess of the optimal solution and improves it at each timestep by performing a finite number of iterations of an optimization algorithm. TDO can be interpreted as a dynamic compensator that maintains a solution estimate as an internal state, the dynamics of which are defined by the optimizer iterations. As illustrated in Figure 1, this interpretation differs from “ideal”, or “optimal” MPC which is an implicitly defined static feedback law.

There are a variety of TDO variants proposed in the literature. The stability of input constrained TDO controllers using linearly convergent optimization algorithms are studied in [18]. Unconstrained suboptimal NMPC without terminal conditions is considered in [21]. A fixed point scheme for input constrained MPC of sampled data input affine systems is proposed in [17], a gradient based dynamic programming approach is considered in [49], a proximal gradient method for linear input-constrained MPC is studied in [51], and a continuous time gradient flow based approach is described in [38]. These methods use some combination of shifting terminal control updates and first order optimization methods. In [2, 41], a generic suboptimal MPC scheme is considered and sufficient conditions on the warmstart for robust stability are derived; the optimization algorithm is not specified and its convergence is not considered. The robustness of MPC to disturbances arising from incomplete optimizations is considered in [55, 20] and conditions for complexity certification of suboptimal state constrained linear MPC are presented in [46]. However, the treatment of the optimizer itself as a dynamic system was not pursued.

An alternative to gradient based approaches are second order methods. In particular, Time Distributed Sequential Quadratic Programming (TD-SQP) methods are attractive since they can be implemented using existing Quadratic Programming (QP) solvers. The fundamental idea behind a TD-SQP based model predictive controller is to apply a finite number of SQP iterations at each sampling instant and to warmstart the iterations with the solution estimate from the previous sampling instant. A widely used variant of TD-SQP is the Real-Time Iteration (RTI) scheme [6] which uses a Gauss-Newton Hessian approximation and performs a single SQP iteration per sampling instant. The RTI scheme has been successfully applied to a variety of applications including engines [1, 56], kites[25], cranes [52], ground vehicles [15], race cars[36], distillation columns [9] and wind turbines [19]. Software for implementing the RTI scheme is provided by the ACADO toolkit [24]. Despite its widespread success, formal stability guarantees for the RTI scheme have only been provided in the absence of inequality constraints [8].

It should be noted that TDO is distinct from so-called suboptimal solution tracking, sensitivity, or running methods, e.g., [55, 54, 32, 10, 16, 27] which are tailored numerical methods for tracking the solutions of parameterized nonlinear programs/generalized equations. These methods are typically used to accelerate or replace existing nonlinear programming solvers to reduce computation times, which is different from considering the dynamic interactions between the plant and the optimizer. Some of them, e.g., [55], consider robust stability by treating suboptimality as a bounded disturbance. This differs from our approach where we treat suboptimality as the output of a dynamic system which is coupled with the closed-loop plant.

This paper begins by presenting a system theoretic framework for analyzing a broad class of TDO algorithms. Specifically, the framework applies to any MPC feedback law that is Locally Input-to-State Stable (LISS) combined with any optimization algorithm featuring a convergence rate that is at least locally $q$-linear. Any MPC formulation with proven LISS properties can be used, including nominal MPC with terminal constraints [35], robust MPC111A nominal MPC controller does not explicitly consider the presence of disturbances in the OCP formulation, unlike a robust MPC controller. [35], and MPC formulations with no terminal constraints [22, Chapter 6]. in this paper, we establish the existence of a joint region of attraction for the state and solution estimate, i.e., we show that if the initial state is sufficiently close to the origin and if the initial solution estimate is sufficiently accurate, the state will converge to the origin and the estimate will converge to the optimal solution. Moreover, we analyze the effect of performing more iterations, establish robustness properties, and show that, if the initial solution guess is within the convergence basin of the optimization method, TDO can recover the robust region of attraction of optimal MPC with a finite number of iterations.

The proposed theoretical framework is then specialized to the RTI scheme. Our analysis extends that in [8] as follows: (i) We explicitly consider inequality constraints and relax the terminal state constraint to a terminal set constraint; (ii) We explicitly consider the robustness properties of the RTI scheme by establishing LISS of the closed-loop system; (iii) We analyze the effect of the number of SQP iterations performed at each sampling instant and establish sufficient conditions for robust constraint satisfaction. We also provide a proof which extends the classical preconditioned fixed-point type analysis of Newton’s method, see [31, Section 5.4.2], to the setting of generalized equations and establish conditions under which discrete time optimal control problems with polyhedral constraints are strongly regular. The latter property is important since it is a sufficient condition for Lipschitz continuity of the optimal value function and thus for robust stability.

This paper builds upon the results in [33] which analyzes the stability of MPC implemented using a suboptimal semismooth predictor-corrector (SSPC) method. Specifically, we generalize the previous results for suboptimal SSPC to a wide class of optimizers which are at least q-linearly convergent. Moreover, we consider external disturbances in our analysis and analyze several different variants of the RTI scheme.

The layout of the paper is as follows. We review pertinent notation and concepts in Section 2 then describe the problem setting and the class of optimization algorithms we consider in Sections 3 and 4. We establish the ISS properties of the optimization algorithms and of the coupled plant-optimizer system in Sections 5 and 6. Next, we discuss the strong regularity assumption in Section 7, we discuss relevant SQP methods in Section 8, and we illustrate how they fit into our optimization framework in Section 9. Finally, we present simulation results in Section 10.

2 Preliminaries

We denote by $\mathbb{Z}_{+(+)}$ the non-negative (positive) integers and by $\mathbb{R}_{+(+)}$ the non-negative (positive) reals. For a discrete time system

[TABLE]

given an initial state $x_{0}\in\mathbb{R}^{n}$, and an input sequence $\mathbf{u}:\mathbb{Z}_{+}\to\mathbb{R}^{m}$ we denote its solution by $x(k,x_{0},\bf{u})$. For a vector, $||\cdot||$ denotes the usual Euclidean norm, for $\mathbf{u}:\mathbb{Z}_{+}\to\mathbb{R}^{m}$ we let $||\mathbf{u}||=\mathrm{sup}\{||u_{k}||:k\in\mathbb{Z}_{+}\}$. We use $\overline{\lim}$ as shorthand for $\limsup$. Recall that a function $\gamma:\mathbb{R}_{+}\to\mathbb{R}_{+}$ is said to be of class $\mathcal{K}$ if it is continuous, strictly increasing and $\gamma(0)=0$. If it is also unbounded, then $\gamma\in\mathcal{K}_{\infty}$. A function $\beta:\mathbb{R}_{+}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$ is said to be of class $\mathcal{KL}$ if $\beta(\cdot,s)\in\mathcal{K}$ for each fixed $s\geq 0$ and $\beta(r,s)\to 0$ as $s\to\infty$ for fixed $r\geq 0$. If $\gamma_{1},\gamma_{2}:\mathbb{R}\to\mathbb{R},$ we denote their composition by $\gamma_{1}\circ\gamma_{2}$. We use $I$ to denote the identity matrix, and use $\mathrm{id}:\mathbb{R}\to\mathbb{R}$ to denote the identity function. We denote the domain of a set-valued mapping $F$ by $\mathrm{dom}~{}F$. If $A$ is a matrix then $A_{i}$ is its $i$th row. If $\mathcal{I}$ is an index set, $|\mathcal{I}|$ is its cardinality and $A_{\mathcal{I}}$ denotes the row wise concatenation of $A_{i},~{}\forall i\in\mathcal{I}$. For two vectors, $(a,b)$ denotes vertical concatenation. We denote the unit ball centered at $x$ by $\mathcal{B}(x)$, it is understood that $\mathcal{B}=\mathcal{B}(0)$. If $X$ is a closed neighbourhood of the origin, we denote its radius by $\mathrm{rad~{}X}$, i.e., the largest $r>0$ such that $\{x~{}|~{}||x||\leq r\}\subseteq X$. The normal cone mapping of a closed-convex set $C$ is defined as

[TABLE]

and set addition/subtraction is defined as

[TABLE]

We make extensive use of the concept of input-to-state stability[29]. Since MPC is a constrained control technique, meaning that it is intrinsically not “global”, it is natural to consider a local variant.

Definition 1** (LISS[28]).**

A system (1) is said to be Locally Input-to-State Stable (LISS) if there exists $\epsilon>0$, $\beta\in\mathcal{KL}$, and $\gamma\in\mathcal{K}$ such that, $\forall k\in\mathbb{Z}_{+}$,

[TABLE]

provided $||x_{0}||\leq\epsilon$ and $||\mathbf{u}||\leq\epsilon$.

Definition 2** (Asymptotic gain[29]).**

Consider system (1), we say that it has an asymptotic gain if there exists some $\gamma\in\mathcal{K}$ such that

[TABLE]

for all $x_{0}\in\mathbb{R}^{n_{x}}$.

3 Problem Setting and Control Strategy

Consider the following discrete time system,

[TABLE]

where $x_{k}\in\mathcal{X}\subset\mathbb{R}^{n_{x}}$, $u_{k}\in\mathcal{U}\subset\mathbb{R}^{n_{u}}$ and $d_{k}\in\mathcal{D}\subset\mathbb{R}^{n_{d}}$ denote the state, input, and disturbance. Throughout this paper we assume full state feedback and that the following holds.

Assumption 1**.**

The function $f_{d}$ in (4) is twice continuously differentiable in its first two arguments, Lipschitz continuous in the third, and $f_{d}(0,0,0)=0$. Moreover, the sets $\mathcal{X}$, $\mathcal{U}$ and $\mathcal{D}$ are compact and contain the origin.

We wish to control (4) using MPC and thus consider an OCP of the following form,

[TABLE]

where $N\in\mathbb{Z}_{++}$ is the horizon length, $\mathcal{Z}\subseteq\mathcal{X}\times\mathcal{U}$ are the constraints, $\mathcal{X}_{f}\subseteq\mathcal{X}$ is the terminal state constraint, $\xi=(\xi_{0},\xi_{1},\dots,\xi_{N})$ is the state sequence, and $\mu=(\mu_{0},\mu_{1},\dots,\mu_{N-1})$ is the control sequence. The OCP (5) is parameterized by the measured system state $x$. We impose the following conditions on (5) to ensure that it is well posed and can be used to construct a stabilizing control law for (4).

Assumption 2**.**

All functions in (5) are twice continuously differentiable in their arguments and their second derivatives are Lipschitz continuous.

Assumption 3**.**

The stage cost satisfies $l(0,0)=0$, and there exists $\alpha_{l}\in\mathcal{K}_{\infty}$ such that $\alpha_{l}(||x||)\leq l(x,u)$ for all $(x,u)\in\mathcal{Z}$. The terminal set $\mathcal{X}_{f}$ is a subset of $\mathcal{X}$, contains the origin in its interior, and is an admissible control invariant set for (4), i.e., for all $x\in\mathcal{X}_{f}$ there exists $u$ such that $f_{d}(x,u,0)\in\mathcal{X}_{f}$ and $(x,u)\in\mathcal{Z}$. $V_{f}$ is a Control Lyapunov Function for (4) with $d=0$, such that,

[TABLE]

for all $x\in\mathcal{X}_{f}$, where $x^{+}=f_{d}(x,u,0)$.

Denote by

[TABLE]

the set of feasible parameters. Under Assumptions 1 and 2, the set $\Gamma$ is compact and, for all $x\in\Gamma$, a minimum of (5) exists [37]. The ideal/optimal MPC feedback policy is then

[TABLE]

where $\mu^{*}(x)$ is a global minimizer of (5). To address the effects of incomplete optimization, we consider the perturbed closed loop system

[TABLE]

where the control signal is corrupted by an additive disturbance that represents suboptimality, i.e., $u_{k}=\kappa(x_{k})+\Delta u_{k}$. Before stating the stability properties of (8), recall the following notion.

Definition 3** (RPI set [35]).**

Given suitable sets $\Delta\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$ and $\mathcal{D}\subseteq\mathbb{R}^{n_{d}}$, a set $\Omega\subseteq\mathbb{R}^{n_{x}}$ is a Robust Positively Invariant (RPI) set for system (8) if $f(x,\Delta u,d)\in\Omega$ for all $x\in\Omega$, $\Delta u\in\Delta\mathcal{U}$, $d\in\mathcal{D}$. In addition, if $\Omega\subseteq\{x~{}|~{}(x,\kappa(x))\in\mathcal{X}\times\mathcal{U}\}$, then $\Omega$ is called an admissible RPI set.

The following theorem summarizes the LISS properties of nominal MPC.

Theorem 1**.**

[35, Theorem 4]** Let Assumptions 1 - 3 hold and suppose that $\phi^{*}(x)$, the optimal value function for (5), is uniformly continuous. Then, the closed-loop system (8) is LISS with respect to $(\Delta u,d)$ on a non-empty RPI set $\Omega\subseteq\Gamma$. Moreover, there exist $c_{1},c_{2}>0$ such that if $\mathrm{rad~{}\Delta\mathcal{U}}\leq c_{1}$ and $\mathrm{rad~{}\mathcal{D}}\leq c_{2}$, then $\Omega$ is an admissible RPI set for (8).

In this paper we consider the situation where not enough computational resources are available to accurately solve (5) at each sampling instant. Instead we will approximately track solution trajectories of (5) as the measured state $x$ in (5c) varies over time. To track the solution trajectories, we use an appropriate iterative optimization algorithm, e.g. SQP, which is warmstarted at time instance $t_{k}$ with the approximate solution from $t_{k-1}$. In this way we construct a dynamic system,

[TABLE]

where $z$ is an estimate of the primal-dual solution of (5) and $\mathcal{T}_{\ell}$ represents a fixed number of optimizer iterations ($\mathcal{T}_{\ell}$ is formally defined in (13)). This leads to the interconnected system illustrated in Figure 1. The objective of this paper is to analyze the interconnection between the plant and the dynamic controller from a systems theoretic perspective.

Remark 1**.**

This paper focuses on a common nominal MPC formulation for the sake of clarity. However, our analysis is applicable to any MPC formulation for which it is possible to prove LISS e.g., formulations that employ exact penalty functions [5] or robust MPC formulations [35]. Moreover, note that in (7) the MPC feedback law is defined using a global optimum of (5). This requirement is not intrinsic to our analysis but rather an artifact of the specific MPC formulation. Our analysis is performed relative to a nominal “ideal” MPC feedback law. If the nominal feedback law is input-to-state stable our analysis is applicable regardless of whether the nominal feedback law is globally optimal or not. For example, one could use the dual-mode MPC formulation in [48] which does not require global optimality.

4 Optimization Strategy

In this section we describe the class of optimization algorithms considered in this paper. We start in an abstract setting to clarify which properties are essential to our analysis. Later in Sections 8 and 9 we will illustrate how SQP fits into this framework.

Suppose the first order necessary conditions for (5) can be written as a parameterized Generalized Equation (GE) of the form

[TABLE]

where $\mathcal{N}_{K}:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is the normal cone mapping222See [14] for more background on set-valued and normal cone mappings. of a closed, convex set $K\subseteq\mathbb{R}^{n}$, $F:\mathbb{R}^{n}\times\Gamma\to\mathbb{R}^{n}$ is a function, $z\in\mathbb{R}^{n}$ are the optimization variables and $x\in\Gamma$ is the parameter. Its solution mapping is

[TABLE]

which can be set valued. Because (10) are necessary conditions for (5) and a minimizer of (5) exists under Assumptions 1 and 2 [37] we have that $\mathrm{dom}~{}S=\Gamma$.

Many optimization algorithms are designed to “solve” necessary conditions. We thus associate (10) with an iterative optimization algorithm of the form

[TABLE]

where $\mathcal{T}:\mathbb{R}^{n}\times\Gamma\times\mathbb{Z}_{+}\to\mathbb{R}^{n}$. Multiple iterations of the algorithm can be represented by the action of the function $\mathcal{T}_{\ell}(z,x):\mathbb{R}^{n}\times\Gamma\to\mathbb{R}^{n}$ which is defined recursively by the sequence

[TABLE]

where $\ell\in\mathbb{Z}_{++}$ is the iteration number and $\mathcal{T}_{0}(z,x)=z$.

Remark 2**.**

The optimality conditions of (5) can be written in multiple forms depending on the choice of (10) and (12). For example, $\xi$ can either be treated as a decision variable or a function of $\mu$. As a result, the definition of $z$ is not unique and is chosen by the designer of the optimization algorithm. Specifically, the vector $z$ always includes the control sequence, but may also include the state sequence and/or the Lagrange multipliers associated with equality (dynamic) or inequality constraints.

In this paper we consider algorithms that converge at least q-linearly for a fixed parameter $x$. The following definition formalizes this notion.

Definition 4** (At least q-linear convergence).**

For any $x\in\Gamma$ and $z^{*}\in S(x)$ an optimization algorithm $\mathcal{T}$ converges to $z^{*}$ if there exists $\varepsilon>0$ such that

[TABLE]

for all $z\in\varepsilon\mathcal{B}(z^{*})$. If there exists $\eta>0$ and $q\geq 1$ such that

[TABLE]

for all $\ell>0$ and $\eta\varepsilon^{q-1}<1$ then $\mathcal{T}$ is said to converge at least q-linearly over $\Gamma$.

Remark 3**.**

The necessary conditions (10) may be satisfied at stationary points or local maxima as well as local minima. We perform a local analysis that allows us to exclude those points.

Next we consider the regularity properties of (11), we will use of the following regularity condition for generalized equations.

Definition 5** (Strong Regularity[44]).**

A set-valued mapping $\Psi:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is said to be strongly regular at $x$ for $y$ if $y\in\Psi(x)$ and there exist neighborhoods $U$ of $x$ and $V$ of $y$ such that the truncated inverse mapping $\tilde{F}^{-1}:V\mapsto F^{-1}(V)\cap U$ is single-valued, i.e., a function, and is Lipschitz continuous on $V$.

Strong regularity reduces to non-singularity of the Jacobian matrix if $\Psi$ is a continuously differentiable function. Our main regularity assumption follows. It establishes that any solution trajectories are Lipschitz continuous. In Section 7, we discuss conditions for strong regularity for specific instances of (10).

Assumption 4**.**

All points $(z,x)$ satisfying $z\in S(x)$ that correspond to minimizers are strongly regular $\forall x\in\Gamma$.

The following theorem shows that Assumption 4 ensures that the notion of tracking a local solution trajectory is well defined.

Theorem 2**.**

[13]** Let the parameter $x\in\Gamma$ be a Lipschitz continuous function of $t\geq 0$. Then each solution trajectory $z(t)\in S(x(t))$ is isolated and Lipschitz continuous.

A local optimization algorithm, such as SQP, can be used to track a specific “isolated branch” of the solution mapping. The branch is implicitly selected through the choice of initial guess supplied to the algorithm. Some MPC formulations require global optima while some do not, as discussed in Remark 1. In practice, local methods like SQP are often used regardless due to the prohibitive computational complexity of global methods.

Remark 4**.**

To summarize, an algorithm/optimality condition pair fits in our framework if:

• •

The optimality conditions can be written in the form (10) and satisfy Assumption 4.

• •

The algorithm can be written in the form (12).

• •

The algorithm is at least q-linearly convergent.

5 LISS of time-distributed optimization

Consider the application of TDO to problem (5). In a real-time setting it is only possible to perform a finite number of iterations per sampling instant which we denote by $\ell$. In this situation the optimizer can be viewed as a dynamical system of the form,

[TABLE]

where $\Xi$ is a surjective matrix that selects $\mu_{0}$ from the solution estimate, i.e., $\kappa(x)=\Xi\bar{s}(x)$ where

[TABLE]

is an isolated single valued restriction333Our analysis is performed relative to the ideal feedback law $\kappa(x)=\Xi\bar{s}(x)$. Any choice of the restriction $\bar{s}$ that renders the origin of the closed loop system (8) LISS is admissible. of $S$ (this is possible due to Theorem 2).

In this section we establish conditions under which (16) is LISS. We consider the associated error system,

[TABLE]

where $e_{k}=z_{k}-\bar{s}(x_{k})$, and $\Delta x_{k}=x_{k+1}-x_{k}$. The error system can be explicitly constructed as follows

[TABLE]

Lemma 1**.**

Consider (16) and its error system (18) and suppose that $\mathcal{T}$ is at least q-linearly convergent. Further, let Assumption 4 hold. Then, there exists $a,\theta:\mathbb{Z}_{++}\to\mathbb{R}_{++}$, such that the error satisfies

[TABLE]

subject to the restriction

[TABLE]

where $\varepsilon$ is the convergence radius in Definition 4 and $b$ is the Lipschitz constant of $\bar{s}$ over $\Gamma$. Further, $a$ and $\theta$ are monotonically decreasing with $\ell$, $\lim_{\ell\to\infty}a(\ell)=0$, $\lim_{\ell\to\infty}\theta(\ell)=0$, and $a(\ell)\in(0,1)$.

{pf*}

Proof. If $z_{k+1}=\mathcal{T}_{\ell}(z_{k},x)$ for some fixed $x$ then the error bound (15) implies that

[TABLE]

for all $z_{k}\in\varepsilon\mathcal{B}(\bar{s}(x))$, d where $\alpha(\ell)=\sum_{i=0}^{\ell-1}q^{i}$. Now consider any $x_{k+1},x_{k}\in\Gamma$ and let $z_{k+1}=\mathcal{T}_{\ell}(z_{k},x_{k+1})$. Then, applying (21) with $x=x_{k+1}$, we obtain that

[TABLE]

where we have used that $\bar{s}$ is Lipschitz on $\Gamma$ with constant $b$ by Assumption 4. Recall that $\varepsilon$ denotes the convergence radius of $\mathcal{T}$ in $\Gamma$. A sufficient condition for the restriction $z_{k}\in\varepsilon\mathcal{B}(\bar{s}(x_{k+1}))$ is then

[TABLE]

Continuing and imposing $||e_{k}||+b||\Delta x_{k}||\leq\varepsilon$,

[TABLE]

If $q=1$, then $\eta^{\alpha(\ell)}\varepsilon^{q^{\ell}-1}=\eta^{\ell}$ since $\alpha(\ell)=\sum_{i=0}^{\ell-1}q^{i}=\ell$ and $\varepsilon^{q^{\ell}-1}=\varepsilon^{0}=1$. Otherwise, note that

[TABLE]

where we used that, for $q>1$,

[TABLE]

Thus,

[TABLE]

where

[TABLE]

and $\theta(\ell)=b~{}a(\ell)$. Since $q\geq 1$ and $\eta\varepsilon^{q-1}<1$ by assumption, $a(\ell)\in(0,1)$, the functions $a$ and $\theta$ are monotonically decreasing and $a(\ell),\theta(\ell)\to 0$ as $\ell\to\infty$. ∎

The following Theorem establishes the LISS properties of the error system.

Theorem 3**.**

Consider (16) and its error system (18) and suppose that $\mathcal{T}$ is at least q-linearly convergent. Further, let Assumption 4 hold. Then, there exists $\tau:\mathbb{Z}_{++}\to\mathbb{R}_{++}$ such that the system is LISS if $||e_{0}||\leq 0.5\varepsilon$ and $||\Delta\mathbf{x}||\leq\tau(\ell)\varepsilon$, where $\tau(\ell)=0.5~{}(\sigma(\ell)+b)^{-1}$ and $\varepsilon,b$, and $\sigma$ are defined in Lemma 1. Further, the asymptotic gain of (18) is of the form

[TABLE]

where $s_{1}$ and $s_{2}$ correspond to the $\Delta x$ and $x$ inputs, respectively, and $\sigma(\ell)\to 0$ monotonically as $\ell\to\infty$.

{pf*}

Proof. Given Lemma 1, if (20) holds for all time instants leading up to $k-1$, it follows by direct computation (see e.g. [29, Example 3.4]) that

[TABLE]

where $\sigma(\ell)=\theta(\ell)/(1-a(\ell))$. To ensure that (20) holds, we first consider the case $k=0$ and note that444Recall that $a+b\leq\max(2a,2b)$ for any two scalars.

[TABLE]

Next, assuming (20) holds for $k-1$ and recalling that $a(\ell)<1$, we can apply (29) at iteration $k$ to show that

[TABLE]

thus ensuring that (20) also holds at $k$ due to the restriction $||e_{0}||\leq 0.5\varepsilon$ and $||\Delta\mathbf{x}||\leq\tau(\ell)\varepsilon$. Since (29) recursively enforces its restrictions, we obtain

[TABLE]

which directly establishes LISS with $\beta_{\ell}(s,k)=2a(\ell)^{k}s$ and $\gamma_{\ell}(||(\Delta\mathbf{x},\mathbf{x})||)=2\sigma(\ell)||\Delta\mathbf{x}||+0\cdot||\mathbf{x}||$. The remaining claims follow from the expression $\sigma(\ell)=\theta(\ell)/(1-a(\ell))$ since $a(\ell),\theta(\ell)\to 0$ monotonically as $\ell\to\infty$.∎

6 LISS properties of suboptimal MPC

Theorem 3 establishes sufficient conditions under which an at least q-linearly convergent optimizer, viewed as a dynamic system, is LISS. It also establishes that the asymptotic gain of (18) can be made arbitrarily small by increasing the number of iterations. Since the closed-loop system (8) is itself LISS, we can derive sufficient conditions under which the coupled system, as shown in Figure 2, is LISS with respect to the disturbance input $d$.

Theorem 4**.**

Consider the dynamical systems

[TABLE]

where $h(x,\Delta u,d)=f(x,\Delta u,d)-x$ and $f$ is defined in (8). Let the optimization algorithm used to construct $\mathcal{G}_{\ell}$ satisfy (15) and let Assumptions 1 - 4 hold. Then, there exists $\ell^{*}>0$ such that if $\ell\geq\ell^{*}$ the interconnected system (33) is LISS with respect to the input $d$.

{pf*}

Proof. Under Assumptions 1 - 4, Theorem 1 holds555Recall that Assumption 4 is sufficient for Lipschitz continuity of the optimal value function $\phi^{*}(x)$.. Thus system $\Sigma_{1}$ is LISS, meaning that there exist asymptotic gains $\gamma_{1},\gamma_{2}\in\mathcal{K}$ such that

[TABLE]

for suitably restricted $d_{k}\in\mathcal{D}$ and $\Delta u_{k}\in\Delta\mathcal{U}$. Let $L$ denote the Lipschitz constant of $h$, then

[TABLE]

combining this with (34) we obtain that

[TABLE]

where $\gamma_{3}=L(\gamma_{1}+\mathrm{id})$, and $\gamma_{4}=L(\gamma_{2}+\mathrm{id})$. Similarly, due to Theorem 3, there exists $\varepsilon>0$ and positive functions $\sigma$ and $\tau$ such that

[TABLE]

given $||\Delta\mathbf{x}||\leq\tau(\ell)\varepsilon$. Therefore, it follows from (33b) that

[TABLE]

Combining (38) with (36) we obtain that

[TABLE]

Thus, if the contraction property

[TABLE]

is satisfied for all $s\in[0,\mathrm{rad~{}\Delta}\mathcal{U}]$, (33) is LISS with suitable restrictions on the initial state and on the disturbance $d$, as detailed in [50, Theorem 2]. Note that, since $\Delta u=\Xi e$, we have $\mathrm{rad~{}\Delta}\mathcal{U}\leq\Xi\varepsilon$ where $\varepsilon$ is the convergence radius666If the optimizer is globally convergent then $\varepsilon$ can be chosen arbitrarily. In that scenario it may be possible to obtain a stronger result using different analysis tools. of the optimizer defined in Theorem 3. Since $\sigma(\ell)\to 0$ monotonically as $\ell\to\infty$, the existence of $\ell^{*}<\infty$ such that (40) is satisfied follows from the finiteness of $\gamma_{3}$, $\mathrm{rad~{}\Delta}\mathcal{U}$, and $||\Xi||$. ∎ Theorem 4 establishes conditions under which the interconnected plant-optimizer system is LISS. However, this result does not provide any information about the set of admissible initial conditions and does not consider constraint satisfaction. By noting that the ideal MPC feedback law admits a robust positively invariant set, we can extend our result by deriving sufficient conditions for constraint satisfaction.

Theorem 5**.**

Suppose that the assumptions of Theorem 4 hold so the interconnected system (33) is LISS. Let $\Omega$ denote the admissible RPI set in Theorem 1, let $\gamma_{\ell}(s)=2\sigma(\ell)s$ denote the asymptotic gain of (18), and let $(x_{k},e_{k})=(x(k,x_{0},\mathbf{d}),e(k,e_{0},\mathbf{d}))$ denote the closed-loop trajectory of (33) for some initial condition $(x_{0},e_{0})$ and disturbance sequence $\mathbf{d}$. Then, if the disturbances are sufficiently small, there exists $\bar{\ell}\geq\ell^{*}$ and $\delta>0$ such that, if $||e_{0}||\leq\delta$ and $x_{0}\in\Omega$, then $x_{k}\in\Omega$ for all $k\geq 0$.

{pf*}

Proof. Due to Theorem 1, given a sufficiently small disturbance set $\mathcal{D}$, there exists a neighbourhood of the origin $\Delta\mathcal{U}$ such that, if $\Delta u_{k}\in\Delta\mathcal{U},~{}\forall k\geq 0$ and $x_{0}\in\Omega$, then $x_{k}\in\Omega,~{}\forall k\geq 0$. Since $\Delta u=\Xi e$ for a surjective matrix $\Xi$, there exists $\rho>0$ such that, if $||e_{k}||\leq\rho$, then $\Delta u_{k}\in\Delta\mathcal{U}$. Given the restriction $||\Delta\mathbf{x}||\leq\tau(\ell)\varepsilon$ and $||e_{0}||\leq 0.5\varepsilon$, where $\tau$ and $\varepsilon$ are defined in Theorem 3, it follows from (32) that $||e_{k}||\leq\rho$ can be imposed by enforcing $||e_{0}||\leq 0.5\rho$ and $2\sigma(\ell)||\Delta\mathbf{x}||\leq\rho$. To enforce $2\sigma(\ell)||\Delta\mathbf{x}||\leq\rho$, we note that the set $\Omega$ is bounded [35, Theorem 4], thus implying that $\Delta x\in\Delta\Omega=\Omega-\Omega$ (see Section 2 for a definition of set subtraction) is bounded by

[TABLE]

Since $\bar{s}$ is finite and $\sigma(\ell)\to 0$ monotonically as $\ell\to\infty$, there exists $\ell_{1}$ such that $2\sigma(\ell_{1})\bar{s}\leq\rho$. Moreover, since $\bar{s}$ is finite and $\tau(\ell)\to\infty$ monotonically as $\ell\to\infty$ there exists $\ell_{2}$ such that $\bar{s}\leq\tau(\ell_{2})\varepsilon$. Thus, letting $\delta=0.5\max\{\rho,\varepsilon\}$ and $\bar{\ell}=\max(\ell^{*},\ell_{1},\ell_{2})$, it follows that the system is LISS with restrictions on the initial conditions $x_{0}\in\Omega$ and $\|e_{0}\|\leq\delta$, as well as restrictions on the external disturbance $d\in\mathcal{D}$.∎ Theorem 5 establishes that, if enough computational resources are available and the initial solution guess is sufficiently accurate, then TDO recovers the robustness properties of optimal MPC.

Remark 5**.**

The results presented in this section are quite general: as long as the MPC formulation is LISS, the solution mapping of the OCP is strongly regular, and the convergence rate of the iterative solver is at least $q$-linear, Theorems 4 and 5 prove that it is possible to achieve robust stability and constraint satisfaction by performing a limited number of solver iterations per time step. Due to the generality of the framework, however, the actual values we obtain for $\ell^{*}$ and $\bar{\ell}$ are likely to be conservative and would be ill-suited for, e.g., complexity certification as in [46], which, it should be noted, only considers the linear case. Despite this drawback, our results significantly extend the existing analysis of the RTI scheme [8] when our framework is applied to TD-SQP. Complexity certification is significantly more challenging in the nonlinear case due to the nonconvexity of the OCPs and is left to future work.

7 Conditions for Strong Regularity

The main results in this paper are all predicated upon Assumption 4, that for each parameter value $x\in\Gamma$ the solution mapping of the OCP is strongly regular. In this section we discuss some common settings and derive strong regularity conditions for each.

7.1 Closed Convex Constraint Sets

If the constraint sets $\mathcal{Z}$ and $\mathcal{X}_{f}$ in (5) are closed and convex, it is possible to write the optimality conditions without introducing dual variables for the inequality constraints. In particular, we can express (5) compactly as

[TABLE]

where $W=\mathcal{Z}\times\mathcal{Z}\ldots\times\mathcal{X}_{f}$ and $w=(\xi_{0},\mu_{0},\ldots,\xi_{N})\in\mathbb{R}^{p}$. The Lagrangian associated with (42) is

[TABLE]

where $\lambda\in\mathbb{R}^{l}$ are dual variables (sometimes “co-states”). The KKT conditions of (42) are

[TABLE]

where $z=(w,\lambda)$ and $Z=W\times\mathbb{R}^{l}$. Note that (44) can be reduced to (10) by choosing $F=\nabla_{z}\mathcal{L}$ and $K=Z$. Our framework requires that (44) be necessary for optimality. To ensure this, we impose the following constraint qualification [45, Theorem 6.14]

[TABLE]

for all $(\bar{w},\bar{\lambda})\in S(\bar{x})$. The following lemma proves that(45) holds automatically in this setting.

Lemma 2**.**

The constraint qualification (45) holds at all points $(z,x)\in\mathbb{R}^{p+l}\times\Gamma$.

{pf*}

Proof. The constraint qualification is implied by surjectivity of the matrix $\nabla_{w}g(w,x)$. Denoting $A_{i}=\nabla_{\xi}f_{d}(\xi_{i},\mu_{i},0)$, and $B_{i}=\nabla_{\mu}f_{d}(\xi_{i},\mu_{i},0)$, the surjectivity of $\nabla_{w}g(w,x)$ becomes the condition that for every $\xi=(\xi_{0},\dots,\xi_{N})$ the system

[TABLE]

has a solution. This condition clearly holds: pick an arbitrary sequence $(\nu_{0},\dots,\nu_{N-1})$ and determine $(\zeta_{0},\ldots,\zeta_{N})$ recursively.∎

Before stating the conditions for strong regularity, we recall the following second order condition (which can be monitored numerically, see e.g., [39, Section 16.2]).

Definition 6** (SOSC).**

The Second Order Sufficient Condition (SOSC) is said to hold at $\bar{z}=(\bar{w},\bar{\lambda})\in S(\bar{x})$ if

[TABLE]

7.1.1 Convex Control Constraints

If only convex control constraints are present, Theorem 6 provides sufficient conditions for strong regularity.

Theorem 6**.**

[12, Theorem 1.2]** Suppose that $\mathcal{Z}=\mathbb{R}^{n_{x}}\times\mathcal{U}$, where $\mathcal{U}$ is closed and convex, and consider any $\bar{z}\in S(\bar{x})$. If (46) holds, then $S$ is strongly regular at $(\bar{z},\bar{x})$.

As a result of Theorem 6, Assumption 4 reduces to the assumption that (46) holds at all minimizers in $\Gamma$. In this scenario, any terminal set constraints would have to be enforced through penalty functions.

7.1.2 Polyhedral State and Control Constraints

If the state and control constraints are convex polyhedra, the following theorem applies. The result was previously asserted without proof in [10, Section 3.2], we provide a proof for completeness.

Theorem 7**.**

Suppose that $W$ in (42) is polyhedral with a representation $W=\{w~{}|~{}Mw\leq h\}$. Now consider a KKT point $\bar{z}=(\bar{w},\bar{\lambda})\in S(\bar{x})$. If (46) holds, then $S$ is strongly regular at $(\bar{z},\bar{x})$.

{pf*}

Proof. Strong regularity of the nonlinear GE (44) at $(\bar{z},\bar{x})$ follows from strong regularity of its partial linearization [44]. This can be written as

[TABLE]

where $\hat{f}=\nabla_{w}\mathcal{L}(\bar{z},\bar{x})$, $R=\nabla_{w}^{2}\mathcal{L}(\bar{z},\bar{x})$, $G=\nabla_{w}g(\bar{z},\bar{x})$, $r=\hat{f}-R\bar{w}-G^{T}\bar{\lambda}$, $g=G\bar{w}$, $\Theta=\{w~{}|~{}Gw=g,~{}Mw\leq h\}$ and $C=\Theta\times\mathbb{R}^{l}$. Equation (47) is an affine GE of the form,

[TABLE]

to which we apply [14, Theorem 2E.6] to establish strong regularity of the mapping $A+\mathcal{N}_{C}$. This requires

[TABLE]

where $\mathcal{E}^{+}=\mathcal{E}-\mathcal{E}$, $\mathcal{E}^{-}=\mathcal{E}\cap-\mathcal{E}$, and

[TABLE]

is the critical cone777See [14, Section 2E] for more details on critical cones. We’ve simplified the expression for $\mathcal{E}$ using [14, Theorem 2E.3] and (44). of $C$ at $\bar{z}$. Next, note that $\mathcal{E}\subseteq\mathcal{E}^{+}\subset\ker G\times\mathbb{R}^{l}$ thus, by the second order condition, $y^{T}Ry=y^{T}\nabla_{w}^{2}\mathcal{L}(\bar{z},x)y>0$ for all $y\in\ker G$. Thus

[TABLE]

which implies that $z^{T}Az\leq 0\implies z=0$ for all $z\in\mathcal{E}^{+}$. As a result, (49) is satisfied and (44) is strongly regular. ∎ Thus, as in the case of convex input constraints, Assumption 4 reduces to the condition that (46) holds at all minimizers in $\Gamma$.

7.2 Nonlinear Inequality Constraints

If the constraint sets in (5) can be expressed in the form $\mathcal{Z}=\{(\xi,\mu)~{}|~{}c(\xi,\mu)\leq 0\}$ and $\mathcal{X}_{f}=\{\xi~{}|~{}c_{f}(\xi)\leq 0\}$ for suitable twice continuously differentiable functions $c:\mathbb{R}^{n_{x}+n_{u}}\to\mathbb{R}^{n_{c}}$ and $c_{f}:\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{cf}}$, then (5) can be written compactly as the following Nonlinear Program (NLP),

[TABLE]

where $w=(\xi,\mu)\in\mathbb{R}^{p}$ are the decision variables. The Lagrangian associated with (51) is

[TABLE]

where $\lambda\in\mathbb{R}^{l}$ and $v\in\mathbb{R}^{m}$ are dual variables. Its KKT conditions [26] are

[TABLE]

where $\mathcal{N}_{+}$ is the normal cone mapping of the non-negative orthant. Comparing (53) with (10) we can identify $z=(w,\lambda,v)$, $K=\mathbb{R}^{p}\times\mathbb{R}^{l}\times\mathbb{R}_{\geq 0}^{m}$, and

[TABLE]

To ensure that (53) are necessary for optimality, as required by our framework, we need to impose a constraint qualification on (51).

Definition 7** (LICQ).**

The Linear Independence Constraint Qualification (LICQ) is said to hold at $(\bar{z},\bar{x})$ if

[TABLE]

where $\mathcal{A}(w)=\{i\in 1~{}...~{}q~{}|~{}h_{i}(w)=0\}$ is the set of active constraint indices and $l$ is the number of equality constraints.

The following Theorem summarizes necessary and sufficient conditions for strong regularity in the context of nonlinear programming.

Theorem 8**.**

[26, Prop 1.27, 1.28]** Consider a parameterized nonlinear program of the form (51) and let $S(x)$ be the solution mapping of its KKT conditions (53). A point $(\bar{z},\bar{x})$ satisfying $\bar{z}\in S(\bar{x})$ is strongly regular if it satisfies the LICQ and the strong second order sufficient condition (SSOSC)

[TABLE]

where

[TABLE]

and $\mathcal{A}^{+}(z,x)=\mathcal{A}(w)\cap\{i~{}|~{}v_{i}>0\}$. Moreover, if $\bar{w}$ is a local minimizer of (51), the LICQ and SSOSC are also necessary conditions for strong regularity.

Thus, Assumption 4 reduces to the assumption that the SSOSC and LICQ hold at all minima in $\Gamma$.

8 Sequential Quadratic Programming

Having defined under what conditions the solution mapping of the OCP is strongly regular, we investigate the convergence properties of two widely used SQP schemes to show that they can be used for TDO. To this effect, note that the OCPs (42) and (51) can both be solved using SQP. Specifically, for (51), given a solution estimate $z_{i}$, the next iterate can be computed by solving the following Quadratic Program (QP)

[TABLE]

where $B_{i}$ approximates the Hessian of the Lagrangian $\nabla_{w}^{2}L$. Specifically, if we denote the Lagrange multipliers associated with the equality and the inequality constraints by ${\pi_{i}}$ and ${\eta_{i}}$ respectively. The SQP update for (51) is $z_{i+1}=(w_{i}+\Delta w_{i},\pi_{i},\eta_{i})$. Note that (55) is fully defined by (42) or (51) except for $B_{i}$, which will depend on the specific SQP method.

SQP applied to (42) is similar. The SQP update becomes $z_{i+1}=(w_{i}+\Delta w_{i},\pi_{i})$, i.e., the inequality duals are removed from the iteration, $B_{i}$ must approximate $\nabla_{z}^{2}\mathcal{L}$ instead of $\nabla_{w}^{2}L$ and (55c) becomes $M\Delta w_{i}+Mw_{i}-h\leq 0$, where $M$ and $h$ satisfy $W=\{w~{}|~{}Mw\leq h\}$.

Remark 6**.**

The convex control constraints setting (Section 7.1.1) technically allows for non-polyhedral convex constraints. In this case (55c) would need to be replaced with a convex constraint of the form $w_{i}+\Delta w_{i}\in W$ but otherwise no changes are necessary.

To provide a unified formulation, we exploit that SQP can be seen as a Newton-type process for solving GEs of the form

[TABLE]

where $z\in\mathbb{R}^{n}$, $x\in\Gamma$, $F:\mathbb{R}^{n}\times\Gamma\to\mathbb{R}^{n}$ is continuously differentiable and $\mathcal{N}_{K}:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is the normal cone mapping for a closed, convex set $K\subseteq\mathbb{R}^{n}$. Newton’s method applied to (56) is

[TABLE]

where the sequence $\{H_{i}\}$ approximates $\nabla_{z}F(z_{i},x)$. Referring to the QP subproblem (55), we note that

[TABLE]

for (51)” for (42) simply discard the third row and column. Thus, the the sequence $\{H_{i}\}$ is fully determined by the Hessian approximation sequence $\{B_{i}\}$.

Remark 7**.**

In this paper we only consider “undamped” Newton methods, which are intrinsically local methods. More sophisticated implementations may include various type of regularization and/or globalization techniques such as trust regions or linesearches to enlarge the methods region of attraction. Nevertheless, undamped Newton methods are commonly used in practice, especially in the context of the RTI scheme, and the tools we develop in this paper are applicable to locally convergent algorithms. We leave the application of our tools to globalized SQP methods to future work and refer readers to e.g., [39] or [26], for more detailed treatments of SQP methods.

8.1 The Josephy-Newton (JN) method

Using the exact Hessian of the Lagrangian results in the Josephy-Newton method. The following theorem summarizes the convergence properties of the JN method applied to (56).

Theorem 9**.**

[26, Theorem 3.2]** Let $z^{*}\in\bar{S}(x)$ for some fixed $x$ and suppose that Assumption 2 holds and $(z^{*},x)$ is strongly regular. Let the sequence $\{z_{i}\}$ be generated by repeatedly solving

[TABLE]

Then, there exists $\bar{\eta}=\bar{\eta}(x)>0$ and $\bar{\epsilon}=\bar{\epsilon}(x)>0$ satisfying $\bar{\eta}\bar{\epsilon}<1$, such that, if $z_{0}\in\bar{\epsilon}\mathcal{B}(z^{*})$, then $\{z_{i}\}$ is unique and converges to $z^{*}$ q-quadratically, i.e.,

[TABLE]

In general, we cannot expect $\nabla_{w}^{2}L$ to be positive semidefinite even in the vicinity of a solution. This may make solving the QP subproblems difficult and is a well known issue in the SQP literature. A detailed discussion is outside the scope of this paper, we refer interested readers to e.g., [26, 39, 4]. We will however briefly discuss two possible solutions. The first is to use an Augmented Lagrangian Hessian, i.e., to use $B=\nabla_{w}^{2}L+\rho(\nabla_{w}g)(\nabla_{w}g)^{T}$ for some $\rho>0$. If the second order sufficient conditions hold, then $B$ will be convex if $\rho$ is sufficiently large [26, Section 4.2]. This will shift the multipliers associated with $g$, see [26, Section 4.2] for details on how to recover the original multipliers. This method may not always be numerically efficient because it can negatively impact the sparsity of $B$. The second is to use a reduced Hessian approach, see e.g., [47, 39], which maintains a basis for the null space of the active constraints and solves the QPs on this reduced space. The Hessian projected onto the reduced space (the “reduced Hessian”) is guaranteed to be positive definite in the vicinity of a solution if a second order condition holds.

8.2 The Gauss-Newton (GN) method

The Gauss-Newton method is applicable when the objective function has the form $\phi(w)=||r(w)||^{2}_{2}$ for some residual function $r$. The Hessian of the Lagrangian is then approximated by

[TABLE]

For example if $\phi(w)=x^{T}Qx+u^{T}Ru$ the GN Hessian approximation is $B=\texttt{blkdiag}(Q,R)$. The GN method has the advantage that the Hessian approximation is guaranteed to be positive semidefinite, so the QP subproblems can be solved reliably. Because of this, the GN method is widely used in practice, see e.g., [24, 52, 19, 20, 7, 1]. The GN approximation error satisfies

[TABLE]

so the approximation error depends on the size of the residuals and on the second derivative $g$ which is related to the nonlinearity of the dynamics. We show in Theorem 10 that it is important to approximate $\nabla_{w}^{2}L(z^{*},x)$ where $z^{*}\in\bar{S}(x)$. The following theorem establishes sufficient conditions for q-linear convergence of the GN method by extending the classical fixed-point type analysis of Newton’s method, see [31, Section 5.4.2]. The nearest analysis we found in the literature is [10, Theorem 3.5] which considers a path tracking problem rather than a fixed one.

Theorem 10**.**

Fix some parameter $x\in\Gamma$, let $z^{*}\in S(x)$ and suppose that Assumptions 2 and 4 hold. Consider a sequence $\{z_{i}\}$ generated by repeatedly solving (57). Further, define $e_{i}=z_{i}-z^{*}$ and suppose that there exist $\bar{\delta}=\bar{\delta}(x)>0$ such that $||H_{i}-\nabla F(z^{*},x)||\leq\bar{\delta}$ for all $i\geq 0$. If the mapping

[TABLE]

is strongly regular for all $i\geq 0$, i.e., $J_{i}^{-1}$ is a Lipschitz continuous function with Lipschitz constant $M>0$, and $\bar{\delta}M<1$, then there exists $\bar{\epsilon}=\bar{\epsilon}(x)>0,$ and $L>0$ such that if $z_{0}\in\bar{\epsilon}\mathcal{B}(z^{*})$, then $\{z_{i}\}$ is unique, converges to $z^{*}$ q-linearly, and

[TABLE]

where $\bar{\eta}=\bar{\eta}(x)=M(\bar{\delta}+L\bar{\epsilon})$.

{pf*}

Proof. A solution, $z^{*}\in\bar{S}(x)$, exists for every $x\in\Gamma$ thanks to Assumption 4; from this point forward we will suppress the dependencies on $x$ in the subsequent expressions. The GN method can be written as

[TABLE]

where $G_{i}(z)=H_{i}z-F(z)$; note that $z^{*}=T_{i}(z^{*})$ for any choice of $\{H_{i}\}$. First consider

[TABLE]

[TABLE]

Since $\nabla F$ is Lipschitz (Assumption 2) the fundamental theorem of calculus implies that there exist $L,\epsilon_{1}>0$ such that

[TABLE]

for all $z_{i}\in\epsilon_{1}\mathcal{B}(z^{*})$, so, taking norms, we obtain that

[TABLE]

By assumption the mapping $J_{i}^{-1}$ is Lipschitz continuous so $\Delta T_{i}^{*}=||T_{i}(z)-T_{i}(z^{*})||$ satisfies

[TABLE]

for all $z_{i}\in\epsilon_{1}\mathcal{B}(z^{*})$. Now consider the update equation

[TABLE]

where we have used that $z^{*}=T_{i}(z^{*})$. Since $J_{i}$ is strongly regular, $T_{i}$ is a function and $\{z_{i}\}$ is unique. Using (66) we have

[TABLE]

Since $M\bar{\delta}<1$ by assumption, it is possible to pick $\bar{\epsilon}\in(0,\epsilon_{1})$ such that $\bar{\eta}=M(\bar{\delta}+L\bar{\epsilon})<1$. Then $\{z_{i}\}$ converges q-linearly to $z^{*}$ if $z_{0}\in\bar{\epsilon}\mathcal{B}(z^{*})$, i.e.,

[TABLE]

Theorem 10 requires that $H_{i}$ be a sufficiently good approximation of $\nabla_{z}F(z^{*})$ and that the GN subproblems be strongly regular. A sufficient condition for strong regularity is that the QP (55) satisfies the LICQ and SSOSC (Theorem 8). In practice, strong regularity can be achieved by a judicious choice of $H_{i}$. For example, if the Hessian approximation is convex then it is possible to guarantee strong regularity of the subproblems by adding a regularization term, i.e., $H\leftarrow H+\delta I$ for some small $\delta>0$. Then the mapping $H+\delta I+\mathcal{N}_{K}$ is strongly monotone, which implies strong regularity [14, Theorem 2F.6].

Remark 8**.**

Theorem 10 just requires the Hessian approximation be sufficiently good. One can conceive of useful approximation schemes other than the GN approximation, e.g., $B_{i}=\nabla_{w}^{2}\phi(w_{i})$ when $\phi$ is convex or $B_{i}=\nabla_{w}^{2}L(\bar{z},x)$ for some fixed $\bar{z}$ sufficiently close to $z^{*}$.

9 Time-distributed SQP

In this section we demonstrate that the methods described in Sections 7 and 8 satisfy the condition of Remark 4 and can therefore be used within the framework presented in Section 4.

Remark 9**.**

Note that TD-SQP using the GN Hessian approximation and with $\ell=1$ corresponds to the RTI scheme [6]. As such, when specialized to the RTI scheme, Theorems 4 and 5 are a significant extension of the existing analysis [8] which does not consider inequality constraints.

Strong Regularity Assumption: As detailed in Section 7.1, in the presence of convex constraint sets Assumption 4 can be reduced to the following:

Assumption 5**.**

The second order sufficient condition (46) holds at all minimizers in $\Gamma$.

As detailed in Section 7.2, in the nonlinear inequalities setting, Assumption 4 can instead be ensured under the following:

Assumption 6**.**

The linear independence constraint qualification (see Definition 7) and strong second order sufficient condition (see Theorem 8) hold at all minimizers in $\Gamma$.

Algorithm Definition: Both SQP methods described in Section 8 are instances of the following iterative process

[TABLE]

for specific choices of $z$, $F$, and $K$. Thus, in both cases the optimization mapping (12) can be written as

[TABLE]

Convergence Rate: If the exact Hessian is used then Theorem 9 applies and the method is at least q-linearly convergent with $q=2$. If the GN Hessian approximation is used then Theorem 10 applies under some additional assumptions regarding the accuracy of the Hessian approximation, and the method is at least q-linearly convergent with $q=1$. In both cases the definition of q-linear convergence requires that there be a uniform convergence constant $\eta$ and convergence radius $\varepsilon$ over $\Gamma$. Under the assumption that the functions $\bar{\eta}(x)$ and $\bar{\epsilon}(x)$ in Theorems 9 and 10 are upper and lower semicontinuous, respectively, these can be defined as $\varepsilon=\inf_{x\in\Gamma}\bar{\epsilon}(x)$ and $\eta=\sup_{x\in\Gamma}\bar{\eta}(x)$. Thus, SQP fits into the framework in Section 4 and can be used for time distributed optimization.

10 A Numerical Example

Figure 3 illustrates a bicycle model of a sedan. We only consider the lateral portion of the dynamics; the longitudinal velocity $s$ is assumed constant. The states and control inputs are,

[TABLE]

where $y$ is the lateral position, $\nu$ is the lateral component of velocity, $\psi$ is the yaw angle, $\omega$ is the yaw rate, $\delta_{f}$ is the front steering angle, and $\delta_{r}$ is the rear steering angle.

The equations of motion are

[TABLE]

where

[TABLE]

The tire forces are described by a Pacejka model.This model is a modified version of the one presented in [53] and roughly represents a 2017 BMW 740i. The vehicle is disturbed by normally distributed wind gusts $d$ with a mean velocity of $15~{}m/s$ and standard deviation of $5~{}m/s$. We obtain a discrete time model using a forward Euler integration scheme with a sampling period of $t_{s}=0.04s$ leading to a discrete time model of the form $x_{k+1}=f_{d}(x_{k},u_{k},d_{k})$. The model parameters are summarized in Table 1888SI units are used and all angles are in radians unless otherwise noted..

The control objective is to perform a lane change maneuver. This can be achieved by stabilizing the origin which is chosen to coincide with the center of the target lane. The vehicle begins in the neighboring lane at the initial condition $x_{0}=\left[-3.7~{}~{}0~{}~{}0~{}~{}0~{}~{}0~{}~{}0\right]^{T}$. The OCP is

[TABLE]

where $f_{d}$ is the discrete time model of the sedan. The vehicle is subject to state constraints which keep the vehicle on the road and restrict its yaw and steering angles. The state constraints on $y,\psi,v$ and $\omega$ are softened using $L_{1}$ exact penalty functions which are implemented using slack variables in order to satisfy our smoothness assumptions. The upper and lower bounds are

[TABLE]

and the weighting matrices are $Q=I_{6\times 6}$, and $R=I_{2\times 2}$. The terminal weight is obtained by solving the discrete time algebraic Riccati equation using the linearization about the origin. The matrices encoding the terminal set, $A_{f}$ and $b_{f}$, are computed using the MPT3 toolbox [23]. The natural residual

[TABLE]

is an error bound [40], i.e., it upper and lower bounds $||z-z^{*}(x)||$, where $z^{*}(x)\in S(x)$, and is commonly used as an easily computable surrogate for the error. We use it throughout this section to measure $||z-z^{*}(x)||$.

Figure 4 compares the RTI scheme [8], i.e., a TD-SQP scheme using the GN Hessian approximation with $\ell=1$, with an LQR controller and the optimal MPC feedback law999All simulations were carried out in MATLAB 2017b on a 2015 Macbook Pro with 16GB of RAM and a 2.8GHz i7 processor. We solved quadratic programs using quadprog. The optimal MPC feedback law was computed using fmincon with default settings. CASADI [3] was used to compute analytic derivatives which were supplied to the optimization routines.. The RTI feedback law successfully stabilizes the origin of the plant-optimizer system and outperforms the LQR controller. The state error and the optimization residual both converge to a ball about the origin, demonstrating the expected robustness due to the LISS properties of the combined system (Theorem 4). The closed-loop trajectories generated by the RTI controller are nearly indistinguishable from those from the optimal feedback law but are an order of magnitude cheaper to compute. The RTI scheme took $0.067s$ on average and $0.75s$ in the worst case vs. $0.65s$ and $3.2s$ for the optimal feedback law. Closed-loop responses using the RTI controller for 15 different initial position and yaw angle combinations, with all other states are initialized to zero, are shown in Figure 5.

Figure 6 compares the GN and JN methods with $\ell=1$ and $\ell=2$. In the bottom plot of Figure 6 note that if $\ell=2$ iterations are performed, the yaw angle constraint is satisfied exactly, even in the presence of disturbances, as predicted by Theorem 5. Also, note that the residuals of the computational subsystem converge faster, for a given number of iterations, if the JN method is used instead of the GN method. This is as expected since the convergence rate of the SQP algorithm is faster when the exact Hessian is used.

11 Conclusions

In this paper we presented a general framework for the stability analysis of model predictive controllers implemented using time-distributed optimization. When specialized to Sequential Quadratic Programming, our result extends the existing stability analysis of the RTI scheme by explicitly considering inequality constraints, analyzing the effect of performing additional SQP iterations, considering a wider class of Hessian approximations, and proving local input-to-state stability of the closed-loop system. Future work includes analyzing the effect of the sampling rate, applying our framework to globalized SQP methods, and developing numerical methods for estimating the the asymptotic gain functions used in the analysis.