A Semismooth Predictor Corrector Method for Suboptimal Model Predictive Control

TL;DR

This paper introduces a semismooth predictor-corrector method for suboptimal model predictive control, reducing computational costs while maintaining stability and constraint satisfaction through a dynamic interconnection approach.

Contribution

It develops a novel semismooth predictor-corrector approach for suboptimal MPC and analyzes its stability and constraint enforcement using input-to-state stability and small gain conditions.

Findings

The method achieves reduced computational complexity.

Closed-loop stability is guaranteed under certain conditions.

Numerical simulations confirm the effectiveness of the approach.

Abstract

Suboptimal model predictive control is a technique that can reduce the computational cost of model predictive control (MPC) by exploiting its robustness to incomplete optimization. Instead of solving the optimal control problem exactly, this method maintains an estimate of the optimal solution and updates it at each sampling instance. The resulting controller can be viewed as a dynamic compensator which runs in parallel with the plant. This paper explores the use of the semismooth predictor-corrector method to implement suboptimal MPC. The dynamic interconnection of the combined plant-optimizer system is studied using the input-to-state stability framework and sufficient conditions for closed-loop asymptotic stability and constraint enforcement are derived using small gain arguments. Numerical simulations demonstrate the efficacy of the scheme.

Equations

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

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

$$s.t.$$

$${\xi }_0=x,c_N({\xi }_N)\le 0,$$

$$c({\xi }_i,u_i)\le 0,i=0,\dots N-1,$$

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

$$x_{k+1}=f(x_k,d_k)=f_d(x_k,k_{mpc}(x_k)+d_k),$$

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

$$z_k=T(z_{k-1},x_k),$$

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

$$||\varphi (k,x_0,\mathbf{u})||\le \text{max}\left[\beta (||x_0||,k),\gamma \left(\underset{0\le j\le k}{\max }||u_j||\right)\right],$$

$$\underset{k\to \infty }{\overline{\lim }}||\varphi (k,x_0,\mathbf{u})||\le \gamma \left(\underset{k\to \infty }{\overline{\lim }}\left[\underset{0\le j\le k}{\max }||u_j||\right]\right),$$

$${\alpha }_1(||x||)\le V(x)\le {\alpha }_2(||x||),$$

$$V(f(x,u))-V(x)\le -cV(x)+\sigma (||u||),$$

$$\gamma (s)\le {\alpha }_1^{-1}\left(\frac{\sigma (s)}{(1-c_0)c}\right),\forall s>0.$$

$$\underset{w}{min.}$$

$$s.t.$$

$$h(w,p)\le 0,$$

$${\nabla }_{w}L(w,\lambda ,v,p)=0,$$

$$g(w,p)=0,$$

$$h(w,p)\le 0,v\ge 0,v^{T}h(w,p)=0.$$

$$\bar{S}(p)=\{z=(w,\lambda ,v)|\text{eqref}eq:KKT\text{ are satisfied}\}.$$

$$\text{rank}\left[\begin{array}{c}{\nabla }_{w}g(\bar{w},\bar{p})\\ [{\nabla }_{w}h(\bar{w},\bar{p}){]}_i\end{array}\right]=m+|I_a(\bar{w},\bar{p})|,i\in I_a(\bar{w},\bar{p}),$$

$$u^T{\nabla }_w^{2}L(\bar{z},\bar{p})u>0,\forall u\in {\mathcal{K}}_{+}(\bar{w},\bar{v},\bar{p})\setminus \{0\},$$

$$S(p)\in \bar{S}(p)$$

$$\psi (a,b)=0\iff a\ge 0,b\ge 0,ab=0.$$

$$\psi (a,b)=a+b-\sqrt{a^2+b^2}.$$

$$F(z,p)=\left[\begin{array}{c}{\nabla }_{z}L(w,\lambda ,v,p)\\ g(w,p)\\ \varphi (-h(w,p),v)\end{array}\right],$$

$$V_{k-1}(p_k-p_{k-1})+B_{k-1}({\bar{z}}_k-z_{k-1})=0,$$

$$F_k({\bar{z}}_k,p_k)+{\bar{B}}_k(z_k-{\bar{z}}_k)=0,$$

$${\partial }_{z}F=\left[\begin{array}{ccc}{\nabla }_w^{2}L(z,p)& {\nabla }_{w}g(w,p{)}^T& {\nabla }_{w}h(w,p{)}^T\\ {\nabla }_{w}g(w,p)& 0& 0\\ -C{\nabla }_{w}h(w,p)& 0& D\end{array}\right],$$

$$({\nu }_i,{\mu }_i)\in \{\begin{array}{ll}(1+\frac{h_i}{r_i},1-\frac{v_i}{r_i}),& \text{if }(h_i,v_i)\ne 0,\\ (1-a,1-b),& \text{if }(h_i,v_i)=0,\end{array},$$

$${\partial }_{p}F=\left[\begin{array}{c}{\nabla }_{pz}L(z,p)\\ {\nabla }_{p}g(w,p)\\ -C{\nabla }_{p}h(w,p)\end{array}\right],$$

$$||{\bar{e}}_k||\le ||e_{k-1}||+c||p_k-p_{k-1}|{|}^2,$$

$$||e_k||\le \eta ||{\bar{e}}_k|{|}^2,$$

Full text

I Introduction

In Model predictive control (MPC) [1, 2] a control law is defined by the solution of a finite horizon optimal control problem (OCP). Although MPC can systematically handle nonlinearities and constraints, it can be difficult to implement in applications where computing power is insufficient for solving a constrained non-convex OCP at each sampling instance. Developments in numerical solution methods, especially for linear-quadratic MPC, have enabled the application of MPC to a wide variety of systems, see e.g., [3, Section 2.6] and references therein. However, the application of MPC to systems requiring fast sampling rates remains challenging.

Suboptimal MPC (SOMPC) is an approach for reducing the computational cost of MPC controllers. In SOMPC, instead of solving the OCP to a high precision at each sampling instance, we maintain a guess of the optimal solution and improve it each sampling instance, e.g., by shifting the control sequence or performing one of more iterations of an optimization algorithm. The difference between an ideal model predictive controller and an suboptimal model predictive controller is illustrated in Figure 1. The ideal MPC law is a static function, while the SOMPC law is a dynamic compensator which maintains an estimate of the optimal solution of the OCP as its internal state.

The paper[4] established that, in the presence of a suitable terminal set, any feasible solution of the OCP is stabilizing. The robustness properties of SOMPC were studied in [5, 6] which established sufficient conditions on the warmstart to ensure stability of the closed-loop system.

Stability of SOMPC without any terminal conditions or constraints was studied in [7]. Continuous time SOMPC schemes using gradient type optimization methods were proposed in [8, 9], which also derive sufficient conditions for the stability of the combined system based on the convergence rate of the underlying optimization method. Discrete time gradient based schemes are considered in [10, 11]. The real-time iteration scheme [12] is a well known SOMPC strategy for NMPC wherein a single quadratic program is solved per timestep. Sufficient conditions for stability of the combined plant-optimizer system were established in [13] in the absence of inequality constraints.

In many cases the generation of points which satisfy the optimality conditions of an OCP can be cast as a parameterized rootfinding problem. This is the approach taken in [14], which considers the unconstrained case and [15] which softens constraints with barriers leading to a smooth rootfinding problem. Both paper consider the robustness of the closed-loop system to disturbances caused by suboptimality; however, the treatment of the optimizer itself as a dynamic system in the loop with the plant was not pursued.

In [16] we proposed the semismooth predictor-corrector (SSPC) method which generates solutions of the OCP by tracking the roots of a parameterized nonsmooth rootfinding problem. In this work we apply SSPC to SOMPC, and derive sufficient conditions under which the combined system is asymptotically stable using input-to-state stability (ISS) [17] and small gain arguments[18]. The performance of the method is illustrated using numerical simulations.

The contributions of this paper are as follows. First, the SSPC method exhibits second order convergence properties at the cost of a single linear system solve per iteration. This compares favorably with gradient methods, which display first order convergence, and SQP type methods, which obtain second order convergence by solving quadratic programs. In addition, compared to existing work for second order SOMPC methods[13], we relax the need for a terminal equality constraint, and consider inequality constraints. Third, the proposed stability proof forgoes the rather restrictive requirement of a monotonically decreasing cost function. This is done by using small-gain arguments to determine under what conditions the interconnected plant-optimizer system is contractive. Finally, we establish sufficient conditions for constraint satisfaction.

II Problem setting and Control Strategy

Consider the following discrete-time system,

[TABLE]

where the state and control input are denoted by $x\in X\subset\mathbb{R}^{n_{x}}$, and $u\in U\subset\mathbb{R}^{n_{u}}$ respectively.

Assumption 1**.**

The system satisfies $f_{d}(0,0)=0$ and the sets $X$ and $U$ are compact and contain the origin.

We seek to control (1) using MPC, which solves an optimal control problem (OCP) of the form

[TABLE]

where $l:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}$, $V_{f}:\mathbb{R}^{n_{x}}\to\mathbb{R}$, $c:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\to\mathbb{R}^{n_{c}}$, $c_{N}:\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{cf}}$, at each sampling instance and applies the first element of the solution sequence, $u_{0}$, to the system. We can then define the ideal MPC feedback law as $k_{mpc}(x)=u_{0}^{*}$, where $(\xi^{*},u^{*})$ is a global minimizer of (2). Let

[TABLE]

denote the feasible set of (2). The resulting closed-loop system is

[TABLE]

where $d_{k}$ is a disturbance that represents suboptimality. We impose the following conditions on the OCP. The first ensures the existence of second derivatives used by the Newton-type optimization method described in Section IV, and the second guarantees closed-loop stability of the nominal system (see Theorem 4).

Assumption 2**.**

All functions in (2) are $\mathcal{C}^{2}$ in their arguments.

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 X\times U$. The set $X_{f}=\{x~{}|~{}c_{N}(x)\leq 0\}$ is an admissible control invariant set for (1) and $V_{f}$ is a Control Lyapunov Function for (1) such that for all $x\in X_{f}$,

[TABLE]

where $x^{+}=f_{d}(x,u,0)$ and $Z=\{(x,u)~{}|~{}c(x,u)\leq 0\}\subseteq X\times U$.

We will also make use of robust positively invariant sets when discussing constraint satisfaction.

Definition 1**.**

[19]** The set $\Omega\subseteq\mathbb{R}^{n_{x}}$ is a Robust Positively Invariant (RPI) set for system (4) with respect to $D$ if $f(x,d)\in\Omega$ for all $x\in\Omega$, and $d\in D$. If, in addition, $\Omega\subseteq\{x~{}|~{}(x,k_{mpc}(x))\in Z\}$, where $Z$ is defined in assumption 3, then $\Omega$ is called an admissible RPI set.

Now suppose that not enough computational resources are available to accurately solve (2) at each sampling instance. To reduce computational requirements we can instead approximately track solutions of (2) as the parameter $x$ in (2c) varies in time by applying a fixed number of iterations of an appropriate numerical iterative method and warmstarting each problem with the approximate solution from the previous timestep. This can lead to considerable computational savings. However, it also introduces an error $d_{k}$ between the optimal MPC control action and the one which is applied. In essence, one has to consider a new dynamical system,

[TABLE]

where $z$ represents an estimate of the solution of (2) and the function $T$ represents the iterative method which runs in parallel with the plant. This leads to an interconnected plant-optimizer system as shown in Figure 1.

In this paper we suggest a specific Newton-type method and analyze the resulting interconnected system from a systems theoretic point of view. In Section III we provide some background on the concepts used in our analysis. In Section IV we describe our proposed method in detail. In Section V we establish the ISS properties of (5) and use them to derive sufficient conditions for the stability of the interconnected system in Section VI. Finally, numerical examples are reported in Section VII.

III Background

We will make extensive use of the notion of input-to-state stability (ISS) [17]. Since we consider constrained systems it is natural to use a local notion of ISS. We will also make extensive use of the notion of an asymptotic gain111Recall 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$..

Definition 2**.**

[18]** Consider a system,

[TABLE]

and let $\phi(k,x_{0},\mathbf{u})$ be its solution at time $k$ with inputs $\mathbf{u}=\{u_{0},u_{1},...,u_{k-1}\}$ and initial condition $x_{0}$. The system is said to be locally input-to-state stable (LISS) if there exists $\varepsilon>0$, $\beta\in\mathcal{KL}$, and $\gamma\in\mathcal{K}$ such that, $\forall k\in\mathbb{Z}_{+}$,

[TABLE]

provided $||x_{0}||\leq\varepsilon$ and $\underset{0\leq j\leq k}{max}~{}||u_{j}||\leq\varepsilon$ for all $k\geq 0$.

Definition 3**.**

[17]** Consider system (6), 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}}$.

The following theorem provides tools for establishing conditions under which a system is LISS and for characterizing its asymptotic gain.

Theorem 1**.**

[18, Lemma 2.3]** Suppose the system in Definition 2 admits a continuous local Lyapunov function $V$ such that for some $\alpha_{1},\alpha_{2}\in\mathcal{K}_{\infty}$, $c>0$ and $\sigma\in\mathcal{K}$,

[TABLE]

for all $(x,u)$ in a neighbourhood of the origin. Then the system is LISS and for any $c_{0}\in(0,1)$ its asymptotic gain $\gamma\in\mathcal{K}$ can be chosen such that

[TABLE]

Next we will impose some conditions on (2) which are needed by the SSPC method. We can compactly write (2) as a parametrized nonlinear program

[TABLE]

where $w=(\xi,u)$, $x=p$ is the parameter, $f:\mathbb{R}^{n}\times\mathbb{R}^{n_{x}}\to\mathbb{R}$, $g:\mathbb{R}^{n}\times\mathbb{R}^{n_{x}}\to\mathbb{R}^{m}$, and $h:\mathbb{R}^{n}\times\mathbb{R}^{n_{x}}\to\mathbb{R}^{q}$. The Lagrangian of (9) is defined as $L(w,\lambda,v,p)=f(w,p)+\lambda^{T}g(w,p)+v^{T}h(w,p)$ where $\lambda\in\mathbb{R}^{m}$ and $v\in\mathbb{R}^{q}$ are dual variables. Let $z=(w,\lambda,v)$, the Karush-Kuhn-Tucker (KKT) conditions for (9) are

[TABLE]

The primal-dual solution mapping of (10), which may be multivalued since (9) is not assumed convex, will be denoted by

[TABLE]

We impose some regularity conditions on the (9) to ensure that the mapping $\bar{S}(p)$ is “well behaved”. The linear independence constraint qualification (LICQ) is said to hold at a point $(\bar{z},\bar{p})$ if

[TABLE]

where $I_{a}(w,p)=\{i\in 1~{}...~{}q~{}|~{}h_{i}(w,p)=0\}$ is the index set of active constraints. Further, if a KKT point $(\bar{z},\bar{p})$ satisfying the LICQ also satisfies

[TABLE]

where $\mathcal{K}_{+}(w,v,p)=\{u\in\mathbb{R}^{n}~{}|~{}\nabla_{w}g(\bar{w},\bar{p})u=0,~{}\nabla_{w}h_{i}(\bar{w},\bar{p})u\leq 0,i\in I^{+}_{a}(\bar{w},\bar{v},\bar{p}),\nabla_{w}f(\bar{w},\bar{p})^{T}u\leq 0\}$, and $I_{a}^{+}(w,v,p)=I_{a}(w,p)\cap\{i~{}|~{}v_{i}>0\}$ then it is said to satisfy the strong second order sufficient conditions (SSOSC). Any KKT point which satisfies the SSOSC and the LICQ is a strict local minimizer of (9). Our main regularity assumption follows.

Assumption 4**.**

(Pointwise strong regularity) The LICQ and SSOSC hold at all KKT points in $\Gamma$.

Theorem 2 establishes Lipschitz continuity of primal-dual solutions of (10) and of the optimal value function. Theorem 3 shows that the solution trajectories of the OCP are isolated and can be tracked. Finally, Theorem 4 establishes the ISS properties of (2).

Theorem 2**.**

At each $(\bar{z},\bar{p})$ satisfying (10) there exists a neighbourhood $P$ of $\bar{p}$ and a constant $L_{p}(\bar{p},\bar{z})$ such that $\bar{S}(p)$ is a single valued function satisfying $||\bar{S}(p)-\bar{z}||\leq L_{p}||p-\bar{p}||,~{}\forall p\in P$.

Proof.

The LICQ and SSOSC are necessary and sufficient for strong regularity of the KKT system, see e.g., [20, Theorem 2G.8], and strong regularity implies that $S$ is locally a Lipschitz continuous function [20, Theorem 2B.1]. ∎

Theorem 3**.**

[20, Theorem 6G.1]** Suppose $p$ is prescribed as a Lipschitz continuous function of a scalar $t\geq 0$. Then the solution trajectory mapping $S(p(t))$ is comprised of isolated Lipschitz continuous trajectories.

In this paper we will concern ourselves with

[TABLE]

which denotes the solution mapping corresponding to the global optimum 222We assume that $S$ is a function, if it is not then one could consider a restriction of $S$ thanks to Theorem 3 of (2). Local minima may cause the closed-loop system to converge to non-zero equilibrium points.

Theorem 4**.**

[19, Theorem 4]** Let Assumptions 1 - 4 hold, then the closed-loop system (4) is LISS with respect to $d$ on a robust positively invariant set $\Omega\subseteq\Gamma$ .

Remark 1**.**

Pointwise strong regularity is a common assumption in time varying optimization, e.g., [21, 22, 23, 24]. The SSOSC is generally easy to enforce through appropriate regularization [25]. If the only constraints are upper and lower control input bounds then the LICQ can be proven to hold a-priori, otherwise the problem can be reformulated, e.g., as described in [26].

IV The semismooth predictor corrector method

This section describes the SSPC method introduced in [16] which is based on mapping the KKT conditions to a parameterized rootfinding problem by replacing the complementarity conditions (10c) with nonsmooth equations. This is done using an nonlinear complementarity (NCP) function [27] $\psi:\mathbb{R}^{2}\mapsto\mathbb{R}$ which has the property that

[TABLE]

We use the Fischer-Burmeister (FB) NCP function [28]

[TABLE]

The FB function is semismooth [29], which allows us to use a nonsmooth Newton-type method based on generalized derivatives333For a function $G:\mathbb{R}^{N}\mapsto\mathbb{R}^{M}$, $\partial G(x)$ denotes Clarke’s Generalized Jacobian [30].. Following [31, 28] the SSPC method uses the FB function to map points satisfying (10) to roots of the following parameterized semismooth rootfinding problem,

[TABLE]

where $z=(w,\lambda,v)$ is the primal-dual variable and $\phi$ is the concatenation of $\psi(-h_{i}(w,p),v_{i})$ for $i=1,\ldots,q$. For any fixed $p$ the roots of $F(z,p)$ coincide with $\bar{S}(p)$ so we can obtain solution trajectories by tracking solutions of $F(x,p)=0$ as $p$ varies in time. The predictor and corrector steps are,

[TABLE]

where $V_{k-1}\in\partial_{p}F(z_{k-1},p_{k-1})$, $B_{k-1}\in\partial_{z}F(z_{k-1},p_{k-1})$, and $\bar{B}_{k}=\partial_{z}F(\bar{z}_{k},p_{k})$. The predictor solves (18a) for $\bar{z}_{k}$. The product $B_{k-1}^{-1}V_{k-1}\Delta p$ is the directional derivative of the solution mapping in the direction $\Delta p=p_{k}-p_{k-1}$. The predictor can thus be interpreted as an Euler integration step. The corrector solves (18b) for $z_{k}$ and is a single iteration of the semismooth Newton’s method. The generalized Jacobians used in (18) are

[TABLE]

where $C=\text{diag}(\nu)$ and $D=\text{diag}(\mu)$ are

[TABLE]

$h_{i}=h_{i}(w,p)$, $r_{i}=||(h_{i}(w,p),v_{i})||_{2}$, and $(a,b)$ are arbitrary scalars satisfying $||(a,b)||_{2}=1$. Any value of $(a,b)$ works, we use $a=b=2^{-1/2}$; in our experience adjusting this value does not result in any performance changes. The derivative $\partial_{p}F(z,p)$ consists of all matrices of the form

[TABLE]

where $C$ is the same matrix as in (19). All elements of $\partial_{z}F(z,\bar{p})$ are guaranteed to be invertible in a neighbourhood of any $\bar{z}\in S(\bar{p})$ [16, Proposition 1]. As detailed in [16, Theorem 2] one can establish error bounds for (18),

[TABLE]

where $e_{k}=z_{k}-S(p_{k})$, and $c,\eta>0$ are positive constants that depend on the properties of $F$. The error bound (22) is looser than the one given in [16, Theorem 2] but is algebraically cleaner. The proof is analogous to that of [16, Theorem 2].

Remark 2**.**

The error bound (22) holds provided $\Delta p_{k-1}=p_{k}-p_{k-1}\in P(p_{k-1})$ and $e_{k}\in E(p_{k})$ where $P(p),E(p)$ are set valued mappings to neighbourhoods of the origin. Since the parameter set $X$ is compact there exists sets $\mathcal{E}$ and $\mathcal{P}$ satisfying $\mathcal{E}\subseteq E(p),~{}\mathcal{P}\subseteq P(p),~{}\forall p\in\Gamma\subseteq X$.

V ISS properties of the SSPC method

In this paper we will consider a variant of the SSPC method where $\ell\in\mathbb{Z}_{>0}$ corrector steps are taken. We can view this process as the following dynamic system,

[TABLE]

where $H$ is the matrix which selects the control input from the primal-dual solution, i.e., $u_{k}=HS(x_{k})=k_{mpc}(x_{k}),$ so that the output $u_{k}$ approximates $k_{mpc}(x_{k})$. For our ISS analysis we will work with the associated error system444The dynamic equation of the error system has been shifted forward by one time instance to bring it into the form typically used in ISS analyses.,

[TABLE]

where $\Delta x_{k}=x_{k+1}-x_{k}$. The error system obeys the following property

[TABLE]

which is obtained from (22) with the corrector applied $\ell$ times. The main result of this section is that (24) is LISS with $\Delta x$ as an input and that, under certain conditions, its ISS gain approaches 0 as $\ell\to\infty$. We begin with a technical lemma which is proven in the appendix.

Lemma 1**.**

For any scalars $a,b\geq 0,k\in\mathbb{N}$ the following holds: $(a+b)^{2^{k}}\leq 2^{2^{k-1}}(a^{2^{k}}+b^{2^{k}})$.

Theorem 5**.**

Consider the SSPC error system (24) and let Assumptions 2 and 4 hold. Then there exists a monotonically increasing function $\varepsilon_{1}(\ell)>0$, $\varepsilon_{2}>0$, and $\rho\in\mathcal{K}_{\infty}$ such that (24) is LISS if $||e_{0}||<\text{min}(\varepsilon_{1}(\ell),\varepsilon_{2})$. $2\eta\cdot\text{min}(\varepsilon_{1}(\ell),\varepsilon_{2})<1$, and $||\Delta\mathbf{x}_{k}||=\underset{0\leq j\leq k}{max}~{}||\Delta x_{j}||\leq\rho^{-1}(\text{min}(\varepsilon_{1}(\ell),\varepsilon_{2}))$ $\forall k\in\mathbb{Z}_{+}$. In addition, there exists $\varepsilon_{3}$ such that if $||\Delta\mathbf{x}_{k}||<\varepsilon_{3},~{}\forall k\in\mathbb{Z}_{+}$ then there exists $\gamma_{2}\in\mathcal{KL}$ such that the ISS gain of (24) satisfies $\gamma(s)\leq\gamma_{2}(s,\ell)$.

Proof.

We begin with (26); using Lemma 1 and performing some algebraic manipulations we obtain,

[TABLE]

where $e^{+},e$, and $\Delta x$ are shorthand for $e_{k+1},e_{k}$, and $\Delta x_{k}$. Consider the candidate ISS Lyapunov function $||\cdot||\in\mathcal{K}_{\infty}$. After further algebraic manipulations, we have that

[TABLE]

[TABLE]

Consider the function $\alpha(r,\ell)$, by analyzing the equation

[TABLE]

we see that, for a fixed $\ell$, $\alpha$ is increasing and thus of class $\mathcal{K}$ on the domain $\mathcal{D}=[0,\bar{\varepsilon}(\ell))$ where

[TABLE]

For any fixed $\ell>0$ the term

[TABLE]

is positive on the interval $[0,(2\eta)^{\frac{1}{2^{\ell}}-1})\subseteq\mathcal{D}$ and thus satisfies the inequality

[TABLE]

Pick an arbitrary scalar $\tau\in(0,\bar{\varepsilon})$, then

[TABLE]

where $\bar{a}=a(\bar{\varepsilon}-\tau,\ell)$. Its clear that $\sigma\in\mathcal{K}_{\infty}$ for any fixed $\ell$ thus the dissipation inequlity in Theorem 1 holds provided $||e_{k}||<\bar{\varepsilon}(\ell)$. Applying Theorem 1 we obtain that, for any $c_{0}\in(0,1)$, the ISS gain can be chosen to satisfy $\gamma(s,\ell)\leq c_{1}\sigma(s,\ell)$, where $c_{1}^{-1}=(1-c_{0})\bar{a}$. Since (24) is LISS for $||e_{k}||$ sufficiently small, we can recursively enforce that $||e_{k}||<\bar{\varepsilon}(\ell)~{}\forall k\in\mathbb{Z}_{+}$ by restricting $||e_{0}||\leq\bar{\varepsilon}(\ell)$ and $\bar{a}^{-1}\sigma(||\Delta\mathbf{x}_{k}||,\ell)\leq\bar{\varepsilon}(\ell)$ [17, Remark 3.7]. Similarly, to enforce that $(e_{k},\Delta x_{k})\in\mathcal{E}\times\mathcal{P},~{}\forall k\in\mathbb{Z}_{+}$, (see Remark 2) we also restrict the initial condition and input to satisfy $||e_{0}||\leq\varepsilon_{2}$ and $\bar{a}^{-1}\sigma(||\Delta\mathbf{x}_{k}||,\ell)\leq\varepsilon_{2}$ where $\varepsilon_{2}$ is chosen small enough so that a ball of radius $\varepsilon_{2}$ is contained within $\mathcal{E}$ and $\mathcal{P}$. This is possible since they are neighbourhoods of the origin. In addition, if $2c\eta r<1$ then $\sigma(r,\ell)\in\mathcal{KL}$. Letting $\varepsilon_{1}(\ell)=\bar{\varepsilon}(\ell)$, $\varepsilon_{3}=\frac{1}{2c\eta}$, and $\rho=\bar{a}^{-1}\sigma$ completes the proof. ∎

VI ISS properties of suboptimal MPC

Theorem 5 illustrates that, under some conditions, the SSPC method, viewed as a dynamic system driven by parameter changes, is LISS and has an asymptotic gain which can be made arbitrarily small by performing more iterations. Since the ideal closed loop system (4) is LISS, we can treat the sub-optimality error as a disturbance and derive sufficient conditions for the stability of the interconnection between the SSPC and the plant using small gain arguments.

Theorem 6**.**

Consider the interconnected dynamic systems

[TABLE]

where $f(x,d)=f_{d}(x,k_{mpc}(x)+d)$, $f_{d}$ is defined in (1), $h(x,d)=f(x,d)-x$, and $G_{\ell}$ is defined in (24). Let Assumptions 1 - 4 and the assumptions of Theorem 5 hold. Then there exists $\ell^{*}>0$ such that, if $\ell\geq\ell^{*}$, the origin is an asymptotically stable equilibrium point for (38) whose region of attraction satisfies $\mathcal{R}\subset\Gamma\times\mathcal{E}$.

Proof.

Under Assumptions 1-4 $\Sigma_{1}$ is LISS by Theorem 4 and thus admits an asymptotic gain from $d$ to $x$ [17, Lemma 3.8]. Since the output equation is Lipschitz continuous this implies that there exists $\gamma_{1}\in\mathcal{K}$ such that

[TABLE]

Similarly by Theorem 5 we have that $\Sigma_{2}$ is LISS and that there exists $\gamma_{2}\in\mathcal{KL}$ such that

[TABLE]

Thus using small gain arguments, see e.g., [18], we have that (38) is contractive in a neighbourhood of the origin provided $\gamma_{1}\circ\gamma_{2}(s,\ell)<s,~{}\forall s\geq 0$. Since $\gamma_{2}\in\mathcal{KL}(s,\ell)$ it can be made arbitrarily small by letting $\ell\to\infty$. It follows from the finiteness of $\gamma_{1}$ that there exists $\ell^{*}$ such that the small gain condition is satisfied. To conclude the proof we define $\mathcal{R}$ as the set of initial conditions under which $(x_{k},e_{k})\in\Gamma\times\mathcal{E}$ for all $k\in\mathbb{Z}_{+}$. ∎

Theorem 6 establishes asymptotic stability of the interconnected plant-optimizer system but doesn’t consider constraint satisfaction. Since the closed-loop system under ideal feedback law is LISS on an admissible RPI set we can derive sufficient conditions for constraint satisfaction that are summarized in the following theorem.

Theorem 7**.**

Suppose that the assumptions of Theorem 6 hold so the interconnected system (38) is LISS. Let $\Omega$ denote the admissible RPI set in Theorem 4, let $\gamma_{2}(s,\ell)\in\mathcal{KL}$ upper bound $\gamma(s)$, the asymptotic gain of (38b), and let $(x_{k},e_{k})$ denote the closed-loop trajectory of (38) for some initial condition $(x_{0},e_{0})$. Then 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$.

Proof.

By Theorem 4, there exists a neighbourhood $D$ of the origin such that, if $d_{k}\in D,~{}\forall k\geq 0,$ and $x_{0}\in\Omega$, then $x_{k}\in\Omega,~{}\forall k\geq 0$. Since $d=He$ for a fixed matrix $H$, this implies the existence of $\rho>0$ such that if $||e_{k}||\leq\rho$ then $d_{k}\in D$. Further, (38) is LISS, so we can enforce $||e_{k}||\leq\rho$ by noting that,

[TABLE]

by the definition of LISS. Thus we must impose that

[TABLE]

where $\beta_{0}(\cdot)=\beta(\cdot,0)\in\mathcal{K}$, and

[TABLE]

for all $\Delta x\in\Delta\Omega=\Omega-\Omega$. As proven in [19, Theorem 4], the set $\Omega$ is bounded, implying that there exists

[TABLE]

Since $\bar{s}$ is finite and $\gamma_{2}\in\mathcal{KL}$ there must exist some $\ell_{1}$ such that $\gamma_{2}(\bar{s},\ell_{1})\leq\rho$. Letting $\delta=\beta_{0}^{-1}(\rho)$ and $\bar{\ell}=\max(\ell^{*},\ell_{1})$ completes the proof. ∎

Theorem 7 establishes that, if enough computational resources are available and the initial solution guess is sufficiently accurate, then constraints are guaranteed to be satisfied.

VII Numerical examples

In this section we illustrate our theoretical results with a numerical example. The attitude dynamics of a rigid spacecraft are given by the Euler equations,

[TABLE]

where $\omega\in\mathbb{R}^{3}$ is the vector of angular velocities expressed in a body fixed frame, $\theta$ is the vector or 3-2-1 Euler angles, $x=[\omega^{T}~{}\theta^{T}]^{T}$ is the state vector, $J=diag(918,920,1365)$, is the inertia matrix, $u\in\mathbb{R}^{3}$ are external control moments and

[TABLE]

The control objective is to drive the system from $x(0)=[0~{}~{}0~{}~{}0~{}~{}15^{\circ}~{}~{}30^{\circ}~{}~{}-20^{\circ}]$ to the origin. We discretize the dynamics using explicit Euler integration with a timestep of $\tau=3~{}s$. The system is placed in closed-loop with a suboptimal MPC controller implemented using SSPC. The optimal control problem (2) to be solved at time $t_{k}$ is

[TABLE]

where $f_{p}(x,u)=x+\tau f_{c}(x,u)$, and $N=30$ is the prediction horizon. We obtain $P$ as the solution of the discrete time algebraic Riccati equation with the dynamics linearized about the origin555$Q=50~{}diag([10~{}10~{}10~{}1~{}1~{}1])$, $R=0.1I_{3\times 3}$. The terminal control invariant set $X_{f}=\{x~{}|~{}A_{f}x\leq b_{f}\}$ is computed using the MPT3 toolbox [32], derivatives are computed using CASADI[33] and the solution estimate is initialized at the origin 666The initial guess has been purposefully chosen to be poor, the initial residual is on the order of 1000..

Closed-loop simulation results are shown in Figure 2. Only one corrector iteration per timestep is needed for stability, however constraints are not satisfied. When $\ell$ is increased to 2 both the state and control constraints are satisfied as predicted by Theorem 7. Figure 3 shows the KKT residual $||F(z,x)||$, which upper and lower bounds the error [31], and the closed-loop cost function. The maximum wall clock execution time777Performed using default settings on a 2015 Macbook pro with a 2.7 GHz i7 CPU and 16 GB of RAM running MATLAB 2017b. for SSPC, implemented in native MATLAB code, for the $\ell=2$ case was $0.0371s$, compared to $1.41s$ for fmincon and $0.9701s$ for ipopt.

VIII Conclusions

In this paper we introduced a suboptimal MPC method based on the SSPC algorithm. We established conditions under which the SSPC algorithm, viewed as a dynamic system, is LISS. We also establish sufficient conditions for stability of the combined SSPC-plant system using small gain theorem based arguments and sufficient conditions for constraint satisfaction. Numerical simulations show that stability can be achieved with only one corrector iteration per time step even in the presence of a large initial estimate error. Constraint enforcement requires either more iterations or a better initial estimate. Future work includes an investigation of the effect of the sampling period, and the design of robust MPC controllers tailored for implementation using SSPC.