A Method for Reducing the Complexity of Model Predictive Control in Robotics Applications

TL;DR

This paper introduces a trajectory parametrization method for model predictive control in robotics that reduces computational complexity, enables warm-starting, and guarantees stability in certain cases, making it suitable for resource-limited systems.

Contribution

The proposed parametrization approach simplifies MPC problems, preserves warm-starting, and provides stability guarantees without extra terminal constraints in linear systems.

Findings

Reduces decision variables in MPC problems.

Enables effective warm-starting of optimization.

Provides stability guarantees for linear systems.

Abstract

This article describes an approach for parametrizing input and state trajectories in model predictive control. The parametrization is designed to be invariant to time shifts, which enables warm-starting the successive optimization problems and reduces the computational complexity of the online optimization. It is shown that in certain cases (e.g. for linear time-invariant dynamics with input and state constraints) the parametrization leads to inherent stability and recursive feasibility guarantees without additional terminal set constraints. Due to the fact that the number of decision variables are greatly reduced through the parametrization, while the warm-starting capabilities are preserved, the approach is suitable for applications where the available computational resources (memory and CPU-power) are limited.

Tables

Table 1. Table I: Parameters of the first-principle model of the pendulum-on-a-cart system.

parameter	value	description
$m$	0.17kg	pendulum mass
$M$	0.74kg	cart mass (lumped)
$l$	0.30m	pendulum length
$k_{\text{m}}$	0.011Nm/A	torque constant
$k_{\text{n}}$	20.62V/s	speed constant
$R_{\text{m}}$	0.30$\Omega$	resistance
$r_{\text{zr}}$	0.018m	belt wheel radius

Equations

$$\inf$$

$$x(k+1)=f(k,x(k),u(k)),$$

$$g(x(k),u(k))\le 0,k\in {\mathbb{Z}}^{+},$$

$$x(0)=x_0,$$

$$\tilde{x}(k)=(I_n\otimes \tau (k){)}^{\mathsf{T}}{\eta }_x,\tilde{u}(k)=(I_m\otimes \tau (k){)}^{\mathsf{T}}{\eta }_u,$$

$$\tau (k+1{)}^{\mathsf{T}}\eta =\tau (k{)}^{\mathsf{T}}\hat{\eta }(\eta ),\forall k\in {\mathbb{Z}}^{+}.$$

$$T:=\left[\tau(k)^{\mathsf{\tiny{T}}}\right]_{k\in I}=\left(\begin{array}[]{c}\tau(k_{1})^{\mathsf{\tiny{T}}}\\ \tau(k_{2})^{\mathsf{\tiny{T}}}\\ \vdots\end{array}\right)$$

$${\left[\tau (k+1{)}^{\mathsf{T}}\right]}_{k\in I}\eta =T\hat{\eta }(\eta ),$$

$$\tau (k+1)=M\tau (k),\forall k\in {\mathbb{Z}}^{+},$$

$$M=\left(\begin{array}[]{cccc}0&0&0&\ldots\\ 1&0&0&\ldots\\ 0&1&0&\ldots\\ \vdots&\vdots&\vdots&\ddots\end{array}\right),\quad\tau(0)=(1,0,0,\dots)^{\mathsf{\tiny{T}}},$$

$$M=e^{M_{\text{c}}T_{\text{s}}},\quad M_{\text{c}}=\left(\begin{array}[]{cccc}-\nu&0&0&\ldots\\ -2\nu&-\nu&0&\ldots\\ -2\nu&-2\nu&-\nu&\ldots\\ \vdots&\vdots&\vdots&\ddots\end{array}\right),$$

$$\tau (k)\in e^{-\nu k{T}_{\text{s}}}\text{span}(1,k{T}_{\text{s}},(k{T}_{\text{s}}{)}^2,\dots ),$$

$$\tau (k)\in e^{-\nu k{T}_{\text{s}}}$$

$$\sin (2\omega k{T}_{\text{s}}),\cos (2\omega k{T}_{\text{s}}),\dots ),$$

$$\sum _{k=0}^{\infty }\delta p(k{)}^{\mathsf{T}}(x(k+1)-f(k,x(k),u(k)))=0,$$

$$\sum _{k=0}^{\infty }\delta {\eta }_p^{\mathsf{T}}(I_n\otimes \tau (k))(\tilde{x}(k+1)-f(k,\tilde{x}(k),\tilde{u}(k))=0,$$

$$\sum _{k=0}^{\infty }(I_n\otimes \tau (k))(\tilde{x}(k+1)-f(k,\tilde{x}(k),\tilde{u}(k)))=0.$$

$$(I_n\otimes$$

$$((I_n\otimes M^{\mathsf{T}}-A\otimes I_s){\eta }_x-(B\otimes I_s){\eta }_u)=0,$$

$$(I_n\otimes M^{\mathsf{T}}-A\otimes I_s){\eta }_x-(B\otimes I_s){\eta }_u=0.$$

$$\tilde{x}(k+1)=A\tilde{x}(k)+B\tilde{u}(k),\forall k\in {\mathbb{Z}}^{+},$$

$$\tilde{x}(k+1)=f(k,\tilde{x}(k),\tilde{u}(k)),k\in \{0,1,\dots ,s\},$$

$$g(\tilde{x}(k),\tilde{u}(k))\le 0,\forall k\in {\mathbb{Z}}^{+}.$$

$$\tilde{z}(k):=(\tilde{x}(k),\tilde{u}(k))=(I_{n+m}\otimes \tau (k){)}^{\mathsf{T}}{\eta }_z,$$

$$\tilde{z}(k+q)=(I_{n+m}\otimes \tau (k){)}^{\mathsf{T}}(I_{n+m}\otimes M^q{)}^{\mathsf{T}}{\eta }_z.$$

$$\tilde{z}(k+s)=\sum _{j=0}^{s-1}{c}_j\tilde{z}(k+j).$$

$$\displaystyle\bar{z}(k+1)=\left(\begin{array}[]{cccc}0&1&\dots&0\\ 0&0&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ c_{0}&c_{1}&\ldots&c_{s-1}\end{array}\right)\otimes I_{n+m}\leavevmode\nobreak\ \bar{z}(k),$$

$$g(\tilde{x}(k),\tilde{u}(k))\le 0,\forall k\in \{0,1,\dots N_{\text{max}}\},$$

$$J_i:=$$

$$\text{s.t. }g(\tilde{x}(k),\tilde{u}(k))\le 0,\text{leavevmode}\forall k\in \{0,\dots ,j\}.$$

$$\underset{{\eta }_x\in {\mathbb{R}}^{ns},{\eta }_u\in {\mathbb{R}}^{ms}}{\inf }\sum _{k=0}^{\infty }l(k,\tilde{x}(k),\tilde{u}(k))\text{s.t. }$$

$$\sum _{k=0}^{\infty }(I_n\otimes \tau (k))(\tilde{x}(k+1)-f(k,\tilde{x}(k),\tilde{u}(k)))=0,$$

$$\tilde{x}(0)=x_0,g(\tilde{x}(k),\tilde{u}(k))\le 0,\forall k\in \{0,1,\dots ,N_{\text{max}}\},$$

$$\underset{{\eta }_x\in {\mathbb{R}}^{ns},{\eta }_u\in {\mathbb{R}}^{ms}}{\inf }{\eta }_x^{\mathsf{T}}(Q\otimes \bar{J}){\eta }_x+{\eta }_u^{\mathsf{T}}(R\otimes \bar{J}){\eta }_u,\text{s.t.}$$

Full text

I Introduction

Model predictive control (MPC) has become a well-known and widely used control strategy for solving challenging control problems. Unlike many other approaches, MPC addresses input and state constraints in a systematic way. It is based on repeatedly solving an optimal control problem, including the actual state as an initial condition and a prediction of the system’s evolution. This leads naturally to an implicit feedback law, providing robustness against modeling errors and disturbances, [1].

Model predictive control has been successfully applied to many robotic systems: For example, the authors of [2] present an MPC-based steering controller for an autonomous vehicle. The controller is shown to perform complex steering maneuvers in an emergency scenario that includes a double lane change on snow. A real-time path planner for an all-terrain vehicle that is relying on MPC is presented in [3], the MPC-based control architecture of an autonomous wheelchair is discussed in [4], and MPC is used to control an unmanned rotorcraft in [5].

In many applications, however, computational resources are limited, and therefore it is often desirable to simplify the resulting optimization problems that are solved at every time step. This can be done, for example, by exploiting the specific system’s structure, [6], or by a careful design of the control structure, [7]. In this work we propose a different strategy, that is based on systematically reducing the degrees of freedom by parametrizing input and state trajectories with basis functions. We will show that the proposed parametrization evolves naturally from the requirement that the resulting trajectories should be invariant to time shifts, meaning that a previously calculated input and/or state trajectory can be used to warm start the optimization at the next time step.

In addition, it will be argued that by choosing exponentially decaying basis functions the truncation of the prediction horizon can be avoided, which, combined with the time-shift property, leads to straightforward recursive feasibility and closed-loop stability guarantees.

Related works: Several strategies have been suggested to reduce the complexity of MPC. For instance, move-blocking describes the general idea of fixing the input or its derivative to be constant over several time steps, [8]. Many different flavors and variants exist, ranging from simply clamping certain inputs together to more elaborate, time-dependent strategies, such as moving window blocking, that are designed to preserve closed-loop stability and recursive feasibility guarantees, [9]. A different point of view is adopted in [10], [11], and [12], for example, where the input is parametrized with Laguerre and/or Kauz basis functions.

In [13], multiresolution analysis is used for parametrizing the input trajectory. However, the approach is mainly applicable to open-loop stable systems, where the impulse response is assumed to be negligible after a certain time horizon. For dealing with unstable systems, the proposed approach would require additional terminal constraints on the unstable modes. Similarly, the authors from [14] apply the wavelet transformation for simplifying the control laws obtained with explicit model predictive control. They show that the resulting simplified control law is everywhere feasible and quantify the suboptimality.

Our contribution is twofold: We show that requiring the parametrization to be invariant with respect to time-shifts imposes a specific structure on the basis functions used for representing the input (and/or state) trajectories. We will highlight that the proposed parametrization captures the classical input description, where no restriction on the input sequence is imposed as a special case. Likewise, our parametrization can be chosen such that Laguerre functions are obtained as a special case, as for example described in [15]. In addition, we show that choosing the basis functions to be decaying, an infinite prediction horizon can be retained, yielding inherent closed-loop stability and recursive feasibility guarantees.

The input and state parametrizations that will be introduced in the following can be viewed as approximations to the underlying constrained linear quadratic regulator problem. This point of view has been explored in the technical report [16], where various approximation results (including convergence) are discussed. Preliminary results appeared in the conference papers [17] and [18]. In [19], and [20] a similar approach has been applied to control a flying vehicle actuated by ducted fans. Unlike the earlier contributions, this article presents a general and flexible parametrization strategy for reducing the computational complexity of MPC, when controlling nonlinear systems. It adopts a discrete-time point of view, and presents previous stability and recursive feasibility results in a unified way. The discrete-time formulation enables an exact description of the sample-and-hold process, which is not possible with the earlier continuous-time results. The effectiveness of the approach is demonstrated by performing a swing up of the inverted pendulum-on-a-cart system, which requires, unlike earlier applications, to fully account for the nonlinear dynamics. In addition, a principled way of choosing the basis functions that parametrize input and state trajectories is discussed, and a comparison to a different parametrization approach presented in the literature ([10, Ch. 5]) is presented. The comparison is in terms of execution time and closed-loop performance. We also comment on the feasibility of the resulting optimization problems.

Outline: Sec. II introduces the problem formulation and the notation that is used throughout the article. The parametrization of the input and state trajectories is discussed in Sec. III. Sec. IV covers the representation of the dynamics, whereas Sec. V shows how to include input and state constraints. The resulting optimization problems are summarized in Sec. VI, and the stability results are included in Sec. VII. The approach is compared to [10, Ch. 3] in Sec. VIII and the application to a nonlinear real-world system is shown in Sec. IX. The article concludes with a summary in Sec. X.

II Problem formulation

We consider the following optimization problem as a starting point

[TABLE]

where the function $l:\mathbb{Z}^{+}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ describes the running cost, $g:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n_{\text{c}}}$ the constraints, and $f:\mathbb{Z}^{+}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ the dynamics. The set of all nonnegative integers is denoted by $\mathbb{Z}^{+}$, the real numbers by $\mathbb{R}$. The integer $n$ describes the state dimension, $m$ the input dimension, and $n_{\text{c}}$ the number of constraints. The initial condition is denoted by $x_{0}$. The above formulation encodes finite horizon problems as a special case, as the running cost can be chosen to be zero for all time indices larger than the given horizon.

III The parametrization

We will parametrize the input and state trajectory using basis functions, that is,

[TABLE]

where $\otimes$ refers to the Kronecker product, and $\eta_{x}\in\mathbb{R}^{ns}$ and $\eta_{u}\in\mathbb{R}^{ms}$ to the parameter vectors. The basis functions are captured with the vectors $\tau(k)\in\mathbb{R}^{s}$, $k\in\mathbb{Z}^{+}$, where the integer $s>0$ refers to the number of basis functions. We use the same basis functions to describe input and state trajectories, but this not necessarily needs to be the case. The choice is motivated by the fact that, as a result, linear time-invariant dynamics can be expressed by a very simple linear relationship between the parameters $\eta_{x}$ and $\eta_{u}$. Throughout the article we will denote parametrized trajectories with a tilde. We will make the following fundamental assumption:

Assumption A1) The basis functions $\tau(k)$, $k\in\mathbb{Z}^{+}$ are linearly independent.

This assumption is necessary for the solutions in (1), when $x(k)$ and $u(k)$ are replaced with $\tilde{x}(k)$ and $\tilde{u}(k)$, to be unique (provided that they exist).

In MPC, (1) will be optimized (up to a given tolerance) at every time step, subject to changing initial conditions $x_{0}$. In order to facilitate computation, we would like to warmstart the optimization at the next time step with the trajectories that were computed at the previous time step. This implies that for any $\eta\in\mathbb{R}^{s}$, there exists a parameter vector $\hat{\eta}(\eta)\in\mathbb{R}^{s}$, such that

[TABLE]

Due to the fact that the basis functions are linearly independent there are $s$ time indices $I:=\{k_{1},k_{2},\dots k_{s}\}$ such that the matrix

[TABLE]

is full rank. Thus, (3) implies that

[TABLE]

which concludes that $\hat{\eta}$ is a linear function of $\eta$ (multiplication of both sides with $T^{-1}$). Combined with the fact that (3) is required to hold for all $\eta\in\mathbb{R}^{s}$, we obtain, by taking the derivative of (3) with respect to $\eta$,

[TABLE]

where the constant matrix $M\in\mathbb{R}^{s\times s}$ contains the partial derivative of $\hat{\eta}$ with respect to $\eta$. This concludes that if the basis function are required to fulfill relation (3), they will automatically satisfy (6) for some matrix $M\in\mathbb{R}^{s\times s}$.

Conversely, if the basis functions satisfy (6), for some matrix $M\in\mathbb{R}^{s\times s}$, (3) is satisfied for all $\eta\in\mathbb{R}^{s}$ with $\hat{\eta}=M^{\mathsf{\tiny{T}}}\eta$. The result is summarized by the following proposition.

Proposition III.1

The basis functions satisfy (6) for some matrix $M\in\mathbb{R}^{s\times s}$ if and only if they are invariant under time shifts in the sense of (3).

In the following we will focus on the regulation problem, and without loss of generality we assume that $x=0$, $u=0$ is an equilibrium of the dynamics. It is thus sensible to require that the basis functions should eventually decay to zero, leading to the requirement that the dynamic system (6) should be asymptotically stable, i.e. the matrix $M$ should have all eigenvalues strictly within the unit circle. Moreover, we assume that the function $g$ describing the constraints is continuous and satisfies $g(0,0)<0$. The case where the equilibrium lies right on the boundary of the constraint manifold has limited practical importance, since, due to measurement noise, the resulting MPC algorithm is likely to become infeasible, even for initial conditions that are in a neighborhood of the equilibrium. The regularity assumptions on $g$ are needed for the closed-loop stability and recursive feasibility statements in Sec. VII, where they guarantee existence of the minimizer of (1), for example.

Summarizing, we thus impose the following requirements on the basis functions

• A1)

The basis functions $\tau(k),k\in\mathbb{Z}^{+}$ are linearly independent.

• A2)

The basis functions $\tau(k),k\in\mathbb{Z}^{+}$ satisfy the relation (6) for some matrix $M\in\mathbb{R}^{s\times s}$ that has all its eigenvalues within the unit circle.

III-A Examples

There are numerous examples of basis functions that comply with Assumptions A1 and A2.

• •

The choice

[TABLE]

yields the classical approach, c.f. [1], where, in this case, the first $s$ inputs can be chosen arbitrarily. More precisely, if $\tilde{u}(k)=\tau(k)^{\mathsf{\tiny{T}}}\eta$, then $\tilde{u}(0)=\eta_{1}$, $\tilde{u}(1)=\eta_{2}$, and so on, until $\tilde{u}(s)=0$, where $\eta=(\eta_{1},\eta_{2},\dots,\eta_{s})^{\mathsf{\tiny{T}}}$.

• •

The choice

[TABLE]

$\tau(0)=\sqrt{2\nu}\leavevmode\nobreak\ (1,1,\dots)^{\mathsf{\tiny{T}}}$, yields so-called Laguerre functions, [21]. These are essentially exponentially decaying polynomials, in the following sense

[TABLE]

where $\nu$ is the decay rate ($\text{s}^{-1}$) and $T_{\text{s}}$ a sampling time.

• •

In an analogous way, the basis functions could be chosen such that

[TABLE]

where $\omega$ denotes the frequency ($\text{rad}/\text{s}$).

• •

In case the dynamics in (1) are linear time-invariant, the cost is quadratic, and constraints are absent ($g=0$), the choice $M=A+BK$, where $A$ and $B$ refer to the system dynamics, and $K$ to the infinite-horizon linear quadratic regulator gain, recovers the solutions to the infinite-horizon linear quadratic regulator problem (see [22] for the continuous-time analogue).

Any superposition of the above choices is valid as well, and is obtained by a blockdiagonal choice of the matrix $M$. In Sec. VIII we present a systematic choice for Laguerre basis functions (according to (8)) based on optimizing the closed-loop performance.

IV Dynamics

We propose to encode the dynamics using a variational formulation. We first note that the dynamics given in (1) can be restated as

[TABLE]

for all variations $\delta p(k)\in\mathbb{R}^{n}$, $k\in\mathbb{Z}^{+}$. We apply the Galerkin approach to obtain a nonlinear equality constraint that approximates the dynamics. More precisely, we require the variations to be spanned by the basis functions and insert the parametrizations for the input and state trajectories. This yields

[TABLE]

for all vectors $\delta\eta_{p}\in\mathbb{R}^{ns}$, or equivalently

[TABLE]

Note that the above equation defines a nonlinear equality constraint between the parameter vectors $\eta_{x}$ and $\eta_{u}$ that parametrize input and state trajectories. If $f$ is continuous about $x=0$, $u=0$ uniformly in $k$, then the above sum is guaranteed to exist, due to the exponential decay of the basis functions.

For gaining additional insights we will now highlight two special cases.

IV-A Linear time-invariant dynamics

In this case $f(k,x(k),u(k))=Ax(k)+Bu(k)$, where $A$ and $B$ are real matrices of appropriate size. As a result, (12) reduces to

[TABLE]

where the property (6) of the basis functions and the properties of the Kronecker product have been exploited. By the linear independence of the basis functions, the matrix $I_{n}\otimes\sum_{k=0}^{\infty}\tau(k)\tau(k)^{\mathsf{\tiny{T}}}$ has full rank,111The matrix is positive definite, as shown by a contradiction argument: If there existed a vector $\eta\neq 0$ such that $\eta^{\mathsf{\tiny{T}}}\sum_{k=0}^{\infty}\tau(k)\tau(k)^{\mathsf{\tiny{T}}}\eta=0$, this would imply that $\tau(k)^{\mathsf{\tiny{T}}}\eta=0$ $\forall k\in\mathbb{Z}^{+}$, contradicting the fact that the basis functions are linearly independent. and thus the above equation can be restated as

[TABLE]

The linear constraint (14) implies that the state and input trajectories fulfill

[TABLE]

that is, the dynamics are fulfilled exactly (as can be seen by multiplying (14) with $(I_{n}\otimes\tau(k)^{\mathsf{\tiny{T}}})$ from the left). Linear independence of the basis functions can be used to conclude the converse. We will not cover the details, since the argument parallels Sec. III, hinging on the fact that the matrix (4) is full rank. The result is summarized by the following proposition.

Proposition IV.1

Provided that the basis functions fulfill Assumptions A1 and A2, and that the dynamics are linear time-invariant, (14) and (15) are equivalent.

IV-B Basis functions chosen according to (7)

In case the basis functions are chosen according to (7), the first $s$ variations are fully decoupled. As result, (12) reduces to

[TABLE]

that is, the dynamics are fulfilled exactly for the first $s$ time steps.

V Constraints

The following section discusses the implementation of the inequality constraints,

[TABLE]

It turns out that in many cases the constraint needs only to be checked at finitely many time instances, due to the decaying nature of the basis functions. We will relate the implementation of the inequality constraints to the well-studied problem of computing the maximum output admissible set for an autonomous linear time-invariant system. To that extent we define the variable $\tilde{z}(k)$ as

[TABLE]

where $\eta_{z}:=(\eta_{x},\eta_{u})$ is obtained by stacking $\eta_{x}$ and $\eta_{u}$. For any time-shift $q\in\mathbb{Z}^{+}$, we can express $\tilde{z}(k+q)$ as

[TABLE]

The Cayley-Hamilton theorem implies that $M^{s}$ can be expressed as $M^{s}=\sum_{j=0}^{s-1}c_{j}M^{j}$, where the constants $c_{j}\in\mathbb{R}$, $j=0,1,\dots,s-1$ are related to the characteristic polynomial of the matrix $M$. This fact can be combined with (19), enabling us to rewrite $\tilde{z}(k+s)$ as

[TABLE]

In other words, the property (6) of the basis functions implies that the trajectories $\tilde{x}(k)$, $\tilde{u}(k)$ match the trajectories of the autonomous linear time-invariant system,

[TABLE]

given the initial condition $\bar{z}(0)=(\tilde{z}(0),\dots,\tilde{z}(s-1))$, and where $\bar{z}(k):=(\tilde{z}(k),\dots,\tilde{z}(k+s-1))$. As a result, from a certain point of view, the constraint (17) describes an output admissible set of the autonomous system (25). In [23], an algorithm is derived that computes a finite dimensional representation (or a close approximation) of the output admissible set of a given autonomous linear time-invariant system. Applied to our case this yields the following result: Under favorable circumstances (see the next paragraph for a precise statement), the constraint (17) is equivalent to

[TABLE]

where $N_{\text{max}}$ can be computed by the recursive algorithm given in Alg. 1. For a general function $g$, the optimization that has to be carried out at each step of the algorithm is non-convex and there are therefore no guarantees that the constants $J_{i}$ are correctly computed. In case $g$ is affine, the optimization simplifies to a linear program that can be solved efficiently and the constants $J_{i}$ are guaranteed to be computed exactly (up to a given tolerance).

As remarked in [23] it is difficult to state general conditions on the existence of $N_{\text{max}}$ (that is, if the algorithm converges). According to [23, Thm. 4.1] an $N_{\text{max}}\in\mathbb{Z}^{+}$ is guaranteed to exist if $g$ is affine, $g(0,0)\neq 0$, and the set of all $(x,u)$ such that $g(x,u)\leq 0$ is bounded.222The observability condition, which is additionally required in [23] is met by construction.

VI Resulting Optimization

Combining the previous results, we propose to approximate (1) in the following way

[TABLE]

where $\tilde{x}(k)$ is replaced by $\tilde{x}(k)=(I_{n}\otimes\tau(k))^{\mathsf{\tiny{T}}}\eta_{x}$ and $\tilde{u}(k)$ by $\tilde{u}(k)=(I_{m}\otimes\tau(k))^{\mathsf{\tiny{T}}}\eta_{u}$ for all $k\in\mathbb{Z}^{+}$.

In the special case of linear time-invariant dynamics, a quadratic cost, and affine constraints, a quadratic program of the following form is obtained

[TABLE]

where the matrices $Q\in\mathbb{R}^{n\times n},R\in\mathbb{R}^{m\times m}$ encode the cost, $C_{\text{x}},C_{\text{u}},b$ the constraints, and $\bar{J}$ is defined as $\bar{J}:=\sum_{k=0}^{\infty}\tau(k)\tau(k)^{\mathsf{\tiny{T}}}$.

VII Closed-loop stability and recursive feasibility

Next we will discuss the question whether the resulting model predictive control algorithm, based on solving optimization problem (28) is recursively feasible and stabilizes the origin. We restrict ourselves to the special case of linear time-invariant dynamics and a time-invariant running cost. We make the assumption that the running cost $l$ is continuous, and lower and upper bounded by two quadratic functions, i.e.

[TABLE]

for all $x\in\mathbb{R}^{n}$ and $u\in\mathbb{R}^{m}$, where $\underline{\sigma}>0$, $\bar{\sigma}>0$. These assumptions ensure that the infimum in (28) is actually attained, which gives rise to (at least one) well-defined minimizer of (28), c.f. [24, p. 12, Ex.1.11]. The quadratic running cost as given in (29) satisfies these assumptions, provided that the matrix $Q$ is positive definite and $R$ is positive semi-definite. We further assume that the linear equality constraints imposing the system dynamics are regular in the sense that for any $x_{0}\in\mathbb{R}^{n}$ there exists at least one $\eta_{x}\in\mathbb{R}^{ns}$ and $\eta_{u}\in\mathbb{R}^{ms}$ that satisfy

[TABLE]

Under these assumptions we obtain the following result, which relies on the time-shift property of the basis functions. The proof can be found in App. A.

Proposition VII.1

Provided that the optimization (28) is feasible at time [math], it remains feasible for all time steps. The resulting model predictive control algorithm asymptotically stabilizes the origin in the sense of Lyapunov.

VIII Simulation Results

The next section discusses the application of the proposed MPC algorithm to a quadruple integrator system. The aim is to highlight and discuss a strategy for choosing the matrix $M$ and to compare the approach to the MPC strategy presented in [10, Ch. 3] that is likewise based on the parametrization of the input.

The dynamics of the quadruple integrator are given by $x_{\text{i}}^{(4)}(t)=u_{\text{i}}(t)$, where $x_{\text{i}}(t)$ denotes the integrator state, the superscript $(4)$ denotes the fourth derivative with respect to time, and $u_{\text{i}}(t)$ the input. The input is constrained to $u_{\text{i}}(t)\in[-0.5,0.5]$ for all $t\in[0,\infty)$. A discretization with zero-order hold and a sampling time of $20$ms yields discrete-time dynamics of the form (1) that are used as a starting point for the proposed MPC strategy. The running cost is chosen to be

[TABLE]

The penality on the input is deliberately chosen to be small in order generate control signals that are likely to hit the constraints. The dynamics and the running cost are motivated by the fact that the resulting optimization problem (28) reduces to a convex quadratic program. As a result, the optimization algorithm that solves (28) (respectively (29)) is guaranteed to return a solution close to the global minimum, irrespective of the initial guess. This enables a rigorous study with randomized initial conditions. Moreover, the approach presented in [10] is based on linear time-invariant dynamics and a quadratic running cost, which is thus required for a comparison.

The following strategy is used for choosing the matrix $M$: We start with a parametrization given by (8) and change the value of $\nu$ from 0.5s*-1* to 2s*-1* in steps of 0.1s*-1*. For each value of $\nu$, we compute the averaged closed-loop cost that is achieved when running the MPC algorithm in closed-loop for 2000 steps (which amounts to 40s), starting at 100 randomized initial conditions. The initial conditions are sampled from a uniform distribution over $[-0.5,0.5]^{4}$ and the value $N_{\text{max}}$ is computed for each different choice of $M$ according to Alg. 1. For the choice $s=8$, the resulting averaged closed-loop cost as a function of $\nu$ is shown in Fig. 1. If $\nu$ is increased above 1.8s*-1* not all initial conditions are feasible. The graph reveals that the averaged closed-loop cost is a non-convex function that has multiple local minima. Due to the lack of robustness as $\nu$ increases above 1.8s*-1*, the choice $\nu=0.8$s*-1* is more favorable than $\nu=1.8$s*-1*, although the averaged closed-loop cost is slightly lower for $\nu=1.8$s*-1*. The resulting closed-loop trajectory starting from $x_{\text{i}}(0)=\dot{x}_{\text{i}}(0)=\ddot{x}_{\text{i}}(0)=\dddot{x}_{\text{i}}(0)=0.5$ is depicted in Fig. 2, where the bang-bang behavior of the input is clearly visible.

Next, the proposed approach is compared to the approach presented in [10, Ch. 3], where the MATLAB scripts given in [10, p. 109, p. 110, p. 123] are used for the implementation. The approach from [10, Ch. 3] is based on augmenting the dynamics, to parametrize the input increments (instead of the input) with Laguerre basis functions, and to eliminate the state trajectory. The value $N_{\text{max}}$, computed according to Alg. 1, is used as the length of the prediction horizon that is required for the approach from [10, Ch. 3]. In both cases, the resulting quadratic program is solved with the optimization routine qpOASES-v.3.2.1, [25] (with the default options). Fig. 3 compares the closed-loop performance, averaged over the 100 randomized initial conditions with the average (over the initial conditions) of the maximum (over a single simulation) execution time, when increasing the number of basis functions $s$ from 8 to 12 (for the approach presented herein) and from 5 to 12 (for the approach presented in [10, Ch. 3]). For values below $s=8$ not all trajectories are stabilized by the approach presented herein and for values below $s=5$ the performance of the approach presented in [10, Ch. 3] exceeds 5. The approach presented in [10, Ch. 3] is based on Laguerre basis functions, and therefore the same basis functions are used for the approach presented herein. The decay rate $\nu$ is set to $1\text{s}^{-1}$. It can be concluded that approach presented herein achieves a smaller cost at similar execution times. Note that in the approach from [10, Ch. 3] no terminal state-constraint and terminal cost is added, and hence closed-loop stability is not guaranteed, which contrasts the approach presented herein (see Prop. VII.1). It is observed that the approach from [10, Ch. 3] is still able to stabilize all initial conditions even for very low values of $s$ (we conducted experiments reducing $s$ up to 2). In contrast, the optimization problem (29) is infeasible for some initial conditions for such low values of $s$ (values $s$ below 8). This needs, however, to be put in perspective with the additional recursive feasibility and closed-loop guarantees that are inherent to (29), which restrict the region where (29) is feasible.

In addition, results that are obtained with a tailored optimization routine exploiting the structure of (29) are also included in Fig. 3. The source code is available in [26].

IX Experimental results

The next section discusses the application of the proposed MPC algorithm to a pendulum-on-a-cart system. The aim is to swing up the pendulum from its hanging equilibrium to its upright position, while accounting for the input and state constraints. A detailed description of the experimental setup can be found in App. B.

The cost function is chosen to be time-invariant and quadratic,

[TABLE]

where relatively large penalties are given to the cart’s position $x_{\text{c}}$, the pendulum angle $\varphi$, and the input $u$. Note that the physical units are omitted for simplicity in (31). The cart is attached to a rail of 0.9m length, and the input (voltage applied to the motor driving the cart) is limited to $\pm 24$V, hence $x_{\text{c}}\in[-0.45\text{m},0.45\text{m}]$ and $u\in[-24\text{V},24\text{V}]$.

The basis functions are chosen to be a combination of the standard parametrization according to (7) and Laguerre functions according to (8), that is,

[TABLE]

$\tau(0)=(1,0,0,\dots)^{\mathsf{\tiny{T}}}\in\mathbb{R}^{19}$, where $M_{\text{d}}\in\mathbb{R}^{12\times 12}$ is chosen according to (7), and $\hat{M}_{\text{d}}\in\mathbb{R}^{7\times 7}$ and $\hat{\tau}(0)$ are chosen according to (8) with $\nu=14$s*-1*. The choice is motivated by the fact that a typical swing-up trajectory (see Fig. 4) involves a bang-bang control input for a bit more than 0.2s (hence roughly 12 time steps) followed by smaller inputs that guide the pendulum to its upright position. These smaller inputs are well approximated by Laguerre basis functions with a relatively small time constant of $1/(14)$s. A small time constant is desirable due to the fact that the resulting value of $N_{\text{max}}$, and hence the number of inequalities in the corresponding optimization problem, is relatively small. The basis functions are chosen to be orthogonal, that is, they are normalized such that $\sum_{k=0}^{\infty}\tau(k)\tau(k)^{\mathsf{\tiny{T}}}=I_{s}$. Running algorithm Alg. 1 subject to the input (limited voltage) and state constraint (limited rail length) results in $N_{\text{max}}=30$.

A rough initial guess for the swing-up trajectories (state and input) is generated via simulation, where the input $\hat{u}(k)$ is obtained by combining step inputs with a linear quadratic regulator for catching the pendulum in upright position. The input trajectory is shown in Fig. 4 (dashed line). The trajectories are then represented using the basis functions. In case of the input trajectory, the resulting parameter vector is given by $\hat{\eta}_{u}=\sum_{k=0}^{\infty}\tau(k)\hat{u}(k)$, and the resulting trajectory $\tau(k)^{\mathsf{\tiny{T}}}\hat{\eta}_{u}$ is shown in Fig. 4 (solid line). The representation of the state trajectories are obtained analogously. The state constraint (limited rail length) is not considered for the generation of the initial guess.

The dynamics are imposed using the variational formulation given by (12), where the infinite sum is truncated after $k=150$. The truncation is motivated by the fact that the contributions for $k>42$ are negligible, as the basis function decay with a time constant of $1/\nu=0.071$s (that is, 4 samples result in a decay of roughly 67%).

The nonlinear optimization problem is solved with Ipopt, [27] (version 3.12.10) and interfaced through C++. The functions for evaluating the derivatives of the Lagrangian and the equality constraints are generated with Casadi, [28] (version 3.4.4). Ipopt is run with the default settings, a relative tolerance of 1e-6, and with a maximum number of 8 iterations (during warmstarts). The time-shift property of the basis functions is exploited for warm-starting the optimization during the swing up. The nominal open-loop swing-up trajectory (as obtained by solving the optimization at the first time instant) is shown in Fig. 5.

The resulting swing-up trajectories are shown in Fig. 6. After roughly 1s the pendulum reaches its upright position. During the whole swing-up trajectory the input and state constraints are accounted for, and the execution time of the optimization remains below 0.02s.

Note the difference between the planned open-loop trajectory and resulting closed-loop trajectory. For instance, the constraint at $x_{\text{c}}=0.45$m is reached during the open-loop trajectory, whereas in closed-loop $x_{\text{c}}=0.45$m is never reached. This discrepancy is most likely due to modeling errors. In particular, the dry friction between the cart and the rail, which is substantial, is difficult to model accurately and to compensate for. Nevertheless, the MPC controller is able to perform the swing up. The robustness of the optimization with respect to the initial guess is further discussed and quantified in App. B.

The above choice of basis functions enables a reduction of the degrees of freedom. In case a standard parametrization according to (7) is used, 250 optimization variables are needed to capture a prediction horizon of 1s, which is roughly required for performing a swing up. The parametrization allows to capture about the same time horizon with only 85 optimization variables.

X Conclusion

The article discussed a method for reducing the complexity of the optimization algorithms encountered in MPC by parametrizing input and state trajectories with basis functions. The basis functions are chosen to be invariant to time-shifts, which enables warm-starting of consecutive optimization problems. In addition, in the case of linear time-invariant dynamics and a quadratic cost, closed-loop stability and recursive feasibility results are obtained without the addition of a terminal set constraint and a terminal cost. The method is applied to the swing-up of an inverted pendulum system, where it is shown that the number of degrees of freedom can be greatly reduced with the parametrization. The resulting MPC algorithm runs with a sampling time of 20ms and is able to swing the inverted pendulum up while accounting for input and state constraints.

Appendix A Proof of Prop. VII.1

We consider the optimal cost to be a function of the initial condition $x_{0}$,

[TABLE]

for all $x_{0}\in X$, where $X$ denotes the set of all $x_{0}$ such that the above optimization is feasible. For $x_{0}\in X$, the assumptions on $l$ and $g$ ($l$ is continuous, $l(0,0)=0$, and bounded below by a quadratic function, $g$ is continuous) assert that the above minimum exists and is attained, see [24, p. 12, Ex. 1.11]. They further imply that $J$ is lower bounded,

[TABLE]

Due to the fact that the inequality constraints in (33) describe a closed set that includes $\eta_{x}=0$, $\eta_{u}=0$ in its interior ($g$ is assumed to be continuous and $g(0,0)<0$), the regularity of equality constraints and the quadratic upper bound on the running cost imply that the inequality constraints are not active for small enough $x_{0}$. According to [24, p. 16, Thm.1.17], $J(x_{0})$ is therefore continuous at zero.

We now consider the trajectory of the closed-loop system $x(k)$, $k\in\mathbb{Z}^{+}$ starting at any $x_{0}\in X$. According to Prop. IV.1, the trajectories $\tilde{x}(k)$ and $\tilde{u}(k)$ resulting from solving the optimization (33) at time [math] satisfy the dynamics and the constraints for all times. Due to the time-shift property of the basis functions, the trajectories $\tilde{x}(k+1)$ and $\tilde{u}(k+1)$ lie likewise in the span of the basis functions, and therefore, they represent feasible candidates for the optimization at the next time step, i.e. for $k=1$. Moreover, the optimal cost at the next time step ($k=1$), $J(x(1))$, is therefore bounded by

[TABLE]

The right-hand side can be expressed in terms of the cost $J(x(0))$, that is,

[TABLE]

which yields, (using $\tilde{x}(0)=x(0)$ and $\tilde{u}(0)=u(0)$),

[TABLE]

It follows by induction that (33) is recursively feasible, and that, for any $k\geq 0$,

[TABLE]

Hence, the cost $J(x(k))$ forms a monotonically decreasing sequence (in $k\in\mathbb{Z}^{+}$) that is bounded below by 0 and therefore converges. In the limit as $k\rightarrow\infty$, (37) reduces to

[TABLE]

implying that $x(k)\rightarrow 0$ as $k\rightarrow\infty$. This shows that the origin is attractive in $X$.

It remains to show that the origin is stable in the sense of Lyapunov. Note that attractivity alone does not imply stability. For proving stability we pick any $\epsilon>0$. Due to the continuity of $J(x_{0})$, there exists a $\delta>0$, such that $|x_{0}|<\delta$ implies $J(x_{0})<\underline{\sigma}\epsilon^{2}$ (and $x_{0}\in X$). As argued in the previous paragraph, the cost $J(x(k))$ of any closed-loop trajectory $x(k)$ starting at a $|x_{0}|<\delta$ is monotonically decreasing. This concludes that $\underline{\sigma}|x(k)|^{2}\leq J(x(k))\leq J(x_{0})<\underline{\sigma}\epsilon^{2}$, that is, $|x(k)|<\epsilon$. ∎

Remarks

The section concludes with the two remarks that indicate further extensions and generalization of the recursive feasibility and closed-loop stability results:

The results can be extended to a running cost that includes a discount, that is,

[TABLE]

where $\rho_{x},\rho_{u}\in\mathbb{R}$ are in the interval $(0,1)$, and $Q\in\mathbb{R}^{n\times n}$ is positive definite and $R\in\mathbb{R}^{m\times m}$ is positive semi-definite.

The results can also be extended to the trajectory-tracking case. More precisely, provided that the dynamics are linear and time-invariant, and that the reference input and state trajectories are spanned by the basis functions, the arguments of Prop. VII.1 can be used to conclude recursive feasibility and asymptotic stability of the tracking error.

• 2)

If the dynamics are nonlinear, recursive feasibility and closed-loop stability are no longer inherent to the problem formulation. However, our approach could be modified to enforce the nonlinear dynamics exactly (for example by eliminating the state trajectory) and to include similar stabilizing terminal conditions as in a finite-horizon non-parametrized approach, see e.g. [29, Ch. 5]. In that case, the arguments presented in [29, Ch. 5] would apply with minor modifications. However, the computational complexity of the resulting MPC algorithm would most likely increase. In all these stability arguments, the fact that the parametrization is invariant to time-shifts is of paramount importance.

Appendix B Experimental results

B-A Hardware setup and first-principles model

Fig. 7 shows the experimental setup. The cart is actuated by an electrical motor via a transmission belt. The voltage that can be applied to the electrical motor represents the input $u$ to the system and is limited to $\pm 24$V. The cart is attached to a rail of 0.9m length. The position of the cart and the angle of pendulum are measured with encoders. An extended Kalman filter is used to estimate the state $x:=(x_{\text{c}},\dot{x}_{\text{c}},\varphi,\dot{\varphi})$. The control and estimation algorithms are run on a laptop (Intel Core i7, 2.6GHz, 8GB random access memory) that is interfacing the pendulum system through a serial communication. All the control and estimation algorithms run at 50Hz.

The following first-principles model is used

[TABLE]

where $F_{\text{u}}$ denotes the force driving the pendulum, $x_{\text{c}}\in[-0.45$m$,0.45$m$]$ the cart position, and $\varphi$ the pendulum angle ($\varphi=0$ corresponds to the upright equilibrium). The parameters are listed in Tab. I. Friction is neglected in the model used for the MPC controller, but dry friction between the cart and rail is compensated with feedforward. The parameters of the first-principles model are identified by fitting different step-responses about the lower equilibrium.

B-B Robustness with respect to the initial guess

The robustness with respect to the initial guess that starts the optimization at each time step is studied in simulations. In general, an initial guess that is close to a successful swing-up trajectory is required to start the optimization (such as the one given by Fig. 4). In order to quantify the robustness of the optimization with respect to the initial guess, we included the following artificial disturbance of the state variable, while performing (closed-loop) swing-up maneuvers in simulation,

[TABLE]

where $f(x(k),u(k))$ are the nominal dynamics and $n(k)$ is a disturbance, uniformly distributed in $[-n_{\text{max}},n_{\text{max}}]$ and independent across time. The optimization is thus warm-started with the solution obtained from the previous time step, but has no information about the perturbation $n(k)$. The resulting closed-loop cost as a function of $n_{\text{max}}$ is shown in Fig. 8. Even with a value of $n_{\text{max}}$ of $5$, which amounts to disturbances of the cart of $\pm 0.1\text{m}$ and disturbances of the pendulum angle of $\pm 5\text{deg}$ at each time step, the optimization algorithm still converges at each time step, resulting in successful (closed-loop) swing-up trajectories. The effect of the disturbances is clearly visible in the closed-loop cost, which increases by almost a factor of 2. Increasing the value of $n_{\text{max}}$ further leads to violations of the position constraint due to the limited rail length.

In our experience, it is important that the optimization is able to handle these disturbances well, as these also occur in the experimental setup due to model mismatch.