Stochastic receding horizon control with output feedback and bounded control inputs

The research on stochastic receding horizon control is broadly subdivided into two parallel lines: the first treats multiplicative noise that enters the state equations, and the second caters to additive noise. The case of multiplicative noise has been treated in some detail in [Primbs, 2007; Primbs and Sung, 2009; Cannon et al., 2009], where the authors consider also soft constraints on the state and control input. However, in this article we exclusively focus on the case of additive noise.

The approach proposed in this article stems from and generalizes the idea of affine parametrization of control policies for finite-horizon linear quadratic problems proposed in [Ben-Tal et al., 2004, 2006], utilized within the robust MPC framework in [Ben-Tal et al., 2006; Löfberg, 2003; Goulart et al., 2006] for full state feedback, and in [van Hessem and Bosgra, 2003] for output feedback with Gaussian state and measurement noise inputs. More recently, this affine approximation was utilized in [Skaf and Boyd, 2009] for both the robust deterministic and the stochastic setups in the absence of control bounds, and optimality of the affine policies in the scalar deterministic case was reported in [Bertsimas et al., 2009]. In [Bertsimas and Brown, 2007] the authors reformulate the stochastic programming problem as a deterministic one with bounded noise and solve a robust optimization problem over a finite horizon, followed by estimating the performance when the noise can take unbounded values, i.e., when the noise is unbounded, but takes high values with low probability. Similar approach was utilized in [Oldewurtel et al., 2008] as well. There are also other approaches, e.g., those employing randomized algorithms as in [Batina, 2004; Blackmore and Williams, 2007; Maciejowski et al., 2005]. Results on obtaining lower bounds on the value functions of the stochastic optimization problem have been recently reported in [Wang and Boyd, 2009].

The remainder of this article is organized as follows. In Section 2 we formulate the receding horizon control problem with all the underlying assumptions, the construction of the estimator, and the main optimization problems to be solved. In Section 3 we provide the main results pertaining to tractability of the optimization problems and mean-square boundedness of the closed-loop system. We comment on the obtained results and provide some extensions in Section 4. We provide all the needed preliminary results, derivations, and proofs in Section 5, and some numerical examples in Section 6. Finally, we conclude in Section 7.

Let $(\Omega,\mathfrak{F},\mathbb{P})$ be a general probability space, we denote the conditional expectation given the sub-$\sigma$ algebra $\mathfrak{F}^{\prime}$ of $\mathfrak{F}$ as $\mathbb{E}_{\mathfrak{F}^{\prime}}[.]$. For any random vector $s$ we let $\Sigma_{s}\coloneqq\mathbb{E}[ss^{\mathsf{T}}]$. Hereafter we let $\mathbb{N}_{+}\coloneqq\{1,2,\ldots\}$ and $\mathbb{N}\coloneqq\mathbb{N}_{+}\cup\{0\}$ be the sets of positive and nonnegative integers, respectively. We let $\mathsf{tr}\!\left(\cdot\right)$ denote the trace of a square matrix, $\left\lVert{\cdot}\right\rVert_{p}$ denote the standard $\ell_{p}$ norm, and simply $\left\lVert{.}\right\rVert$ denote the $\ell_{2}$-norm. In a Euclidean space we denote by ${\mathcal{B}}_{r}$ the closed Euclidean ball or radius $r$ centered at the origin. For any two matrices $A$ and $B$ of compatible dimensions, we denote by $\mathfrak{R}_{k}(A,B)$ the $k$-th step reachability matrix $\mathfrak{R}_{k}(A,B)\coloneqq\left[\begin{matrix}A^{k-1}B&\cdots&AB&B\end{matrix}\right]$. For any matrix $M$, we let $\sigma_{\min}$ and $\sigma_{\max}$ be its minimal and maximal singular values, respectively. We let $(M)_{i_{1}:i_{2}}$ denote the sub-matrix obtained by selecting the rows $i_{1}$ through $i_{2}$ of $M$, and let $(M)_{i}$ denote its $i$-th row.

Consider the following affine discrete-time stochastic dynamical model:

where $t\in\mathbb{N}$, $x_{t}\in\mathbb{R}^{n}$ is the state, $u_{t}\in\mathbb{R}^{m}$ is the control input, $y_{t}\in\mathbb{R}^{p}$ is the output, $w_{t}\in\mathbb{R}^{n}$ is a random process noise, $v_{t}\in\mathbb{R}^{p}$ is a random measurement noise, and $A$, $B$, and $C$ are known matrices. We posit the following main assumption:

Without any loss of generality, we also assume that $A$ is in real Jordan canonical form. Indeed, given a linear system described by system matrices $\bigl{(}\tilde{A},\tilde{B}\bigr{)}$, there exists a coordinate transformation in the state-space that brings the pair $\bigl{(}\tilde{A},\tilde{B}\bigr{)}$ to the pair $(A,B)$, where $A$ is in real Jordan form [Horn and Johnson, 1990, p. 150]. In particular, choosing a suitable ordering of the Jordan blocks, we can ensure that the pair $(A,B)$ has the form $\left(\begin{bmatrix}A_{1}&0\\ 0&A_{2}\end{bmatrix},\begin{bmatrix}B_{1}\\ B_{2}\end{bmatrix}\right)$, where $A_{1}\in\mathbb{R}^{n_{1}\times n_{1}}$ is Schur stable, and $A_{2}\in\mathbb{R}^{n_{2}\times n_{2}}$ has its eigenvalues on the unit circle. By hypothesis (i)-(i)(b) of Assumption 1, $A_{2}$ is therefore block-diagonal with elements on the diagonal being either $\pm 1$ or $2\times 2$ rotation matrices. As a consequence, $A_{2}$ is orthogonal. Moreover, since $(A,B)$ is stabilizable, the pair $(A_{2},B_{2})$ must be reachable in a number of steps $\kappa\leqslant n_{2}$ that depends on the dimension of $A_{2}$ and the structure of $(A_{2},B_{2})$, i.e., $\text{rank}(\mathfrak{R}_{\kappa}(A_{2},B_{2}))=n_{2}$. Summing up, we can start by considering that the state equation (2.1a) has the form

where $A_{1}$ is Schur stable, $A_{2}$ is orthogonal, and there exists a nonnegative integer $\kappa\leqslant n_{2}$ such that the subsystem $(A_{2},B_{2})$ is reachable in $\kappa$ steps. This integer $\kappa$ is fixed throughout the rest of the article.

Fix a prediction horizon $N\in\mathbb{N}_{+}$, and let us consider, at any time $t$, the cost $V_{t}\coloneqq V(t,\mathcal{Y}_{t},u)$ defined by:

where $\mathcal{Y}_{t}=\{y_{0},y_{1},\ldots,y_{t}\}$ is the set of outputs up to time $t$, (or more precisely the $\sigma$-algebra generated by the set of outputs,) and $Q_{k}=Q^{\mathsf{T}}_{k}>0$, $Q_{N}=Q_{N}^{\mathsf{T}}>0$, and $R_{k}=R^{\mathsf{T}}_{k}>0$ are some given symmetric matrices of appropriate dimension, $k=0,\ldots,N-1$. At each time instant $t$, we are interested in minimizing (2.3) over the class of causal output feedback strategies $\mathcal{G}$ defined as:

for some measurable functions $\mathbf{g}_{t}\coloneqq\{g_{t},$ $\cdots,g_{t+N-1}\}$ which satisfy the hard constraint (C1) on the inputs.

The cost $V_{t}$ in (2.3) is a conditional expectation given the observations up through time $t$, and as such requires the conditional density $f(x_{t}|\mathcal{Y}_{t})$ of the state given the previous and current measurements. Our choice of causal control policies in (2.4) motivates us to rewrite the cost $V_{t}$ in (2.3) as:

The propagation of $f(x_{t}|\mathcal{Y}_{t})$ in time is not a trivial task in general. In what follows we shall report an iterative method for the computation of $f(x_{t}|\mathcal{Y}_{t})$, whenever the control input is a measurable deterministic feedback from the past and present outputs. For $t,s\in\mathbb{N}_{0}$, define $\hat{x}_{t|s}=\mathbb{E}_{\mathcal{Y}_{s}}[x_{t}]$ and $P_{t|s}=\mathbb{E}_{\mathcal{Y}_{s}}[(x_{t}-\hat{x}_{t|s})(x_{t}-\hat{x}_{t|s})^{\mathsf{T}}]$.

The following result is a slight extension of the usual Kalman Filter for which a proof can be found in Appendix A. An alternative proof can also be found in [Kumar and Varaiya, 1986, p.102].

In particular, the matrix $P_{t|t}$ together with $\hat{x}_{t|t}$ characterize the conditional density $f(x_{t}|\mathcal{Y}_{t})$, which is needed in the computation of the cost (2.5). Proposition 3 states that the conditional mean and covariances of $x_{t}$ can be propagated by an iterative algorithm which resembles the Kalman filter. Since the system is generally nonlinear due to the function $g$ being nonlinear, we cannot assume that the probability distributions in the problem are Gaussian (in fact, the a priori distributions of $x_{t}$ and of $\mathcal{Y}_{t}$ are not) and the proof cannot be developed in the standard linear framework of the Kalman filter.

Hereafter, we shall denote for notational convenience by $\hat{x}_{t}$ the estimate $\hat{x}_{t|t}$, and let

which corresponds to the Jordan decomposition in (2.2).

Having discussed the iterative update of the laws of $x_{t}$ given the measurements $\mathcal{Y}_{t}$ in Section 2.1, we return to our optimization problem in (2.5). We can rewrite our optimization problem as:

where the collection of functions $\mathbf{g}_{t}$ were defined following (2.4).

The explicit solution via dynamic programming to Problem (OP1) over the class of causal feedback policies $\mathcal{G}$, is difficult if not impossible to obtain in general [Bertsekas, 2000, 2007]. One way to circumvent this difficulty is to restrict our search for feedback policies to a subclass of $\mathcal{G}$. This will result in a suboptimal solution to our problem, but may yield a tractable optimization problem. This is the track we pursue next.

Given a control horizon $N_{c}$ and a prediction horizon $N\;(\geqslant N_{c})$, we would like to periodically solve Problem (OP1) at times $t=0,N_{c},2N_{c},\cdots$ over the following class of affine-like control policies

where $\hat{y}_{i}=C\hat{x}_{i}$, and for any vector $z\in\mathbb{R}^{p}$, $\varphi_{i}(z)\coloneqq\bigl{[}\varphi_{i,1}(z_{1}),\ldots,$ $\varphi_{i,p}(z_{p})\bigr{]}^{\mathsf{T}}$, where $z_{j}$ is the $j$-th element of the vector $z$ and $\varphi_{i,j}:\mathbb{R}\to\mathbb{R}$ is any function with $\sup\limits_{s\in\mathbb{R}}|\varphi_{i,j}(s)|\leqslant\varphi_{\max}<\infty$. The feedback gains $\theta_{\ell,i}\in\mathbb{R}^{m\times p}$ and the affine terms $\eta_{\ell}\in\mathbb{R}^{m}$ are the decision variables. The value of $u_{\ell}$ in (POL) depends on the values of the measured outputs from the beginning of the prediction horizon at time $t$ up to time $\ell$ only, which requires a finite amount of memory. Note that we have chosen to saturate the measurements we obtain from the vectors $(y_{i}-\hat{y}_{i})$ before inserting them into the control policy. This way we do not assume that either the process noise or the measurement noise distributions are defined over a compact domain, which is a crucial assumption in robust MPC approaches [Mayne et al., 2000]. Moreover, the choice of element-wise saturation functions $\varphi_{i}(\cdot)$ is left open. As such, we can accommodate standard saturation, piecewise linear, and sigmoidal functions, to name a few.

Finally, we collect all the ingredients of the stochastic receding horizon control problem in Algorithm 1 below.

In order to state the main results in Section 3, it is convenient to utilize a more compact notation that describes the evolution of all signals over the entire prediction horizon $N$.

The evolution of the system dynamics (2.1a)-(2.1b) over a single prediction horizon $N$, starting at $t$, can be described in a lifted form as follows:

where $X_{t}=\left[\begin{matrix}\begin{smallmatrix}x_{t}\\ x_{t+1}\\ \vdots\\ x_{t+N}\end{smallmatrix}\end{matrix}\right]$, $U_{t}=\left[\begin{matrix}\begin{smallmatrix}u_{t}\\ u_{t+1}\\ \vdots\\ u_{t+N-1}\end{smallmatrix}\end{matrix}\right]$, $W_{t}=\left[\begin{matrix}\begin{smallmatrix}w_{t}\\ w_{t+1}\\ \vdots\\ w_{t+N-1}\end{smallmatrix}\end{matrix}\right]$, $Y_{t}=\left[\begin{matrix}\begin{smallmatrix}y_{t}\\ y_{t+1}\\ \vdots\\ y_{t+N}\end{smallmatrix}\end{matrix}\right]$, $V_{t}=\left[\begin{matrix}\begin{smallmatrix}v_{t}\\ v_{t+1}\\ \vdots\\ v_{t+N}\end{smallmatrix}\end{matrix}\right]$, $\mathcal{A}=\left[\begin{matrix}\begin{smallmatrix}I\\ A\\ \vdots\\ A^{N}\end{smallmatrix}\end{matrix}\right]$,

and $\mathcal{C}=\mathrm{diag}\{C,\cdots,C\}$. Let

be the usual Kalman gain and define $\Gamma_{t}\coloneqq I-K_{t}C$, and $\Phi_{t}\coloneqq\Gamma_{t}A$. Then, we can write the estimation error vector as

where $e_{t}=x_{t}-\hat{x}_{t}$, $\mathcal{F}^{e}_{t}=\tiny\left[\begin{matrix}I\\ \Phi_{t}\\ \Phi_{t+1}\cdot\Phi_{t}\\ \vdots\\ \Phi_{t+N-1}\cdot\Phi_{t+N-2}\cdots\Phi_{t}\end{matrix}\right]$,

Finally, the innovations sequence can be written as

where $\hat{Y}_{t}\coloneqq\mathcal{C}\hat{X}_{t}$. Consequently, the innovations sequence over the prediction horizon is independent of the inputs vector $U_{t}$.

The cost function (2.3) at time $t$ can be written as

where $\mathcal{Q}=\text{diag}\{Q_{0},\cdots,Q_{n}\}$ and $\mathcal{R}=\text{diag}\{R_{0},\cdots,R_{N-1}\}$. Also, the control policy (POL) at time $t$ is given by

where

and $\boldsymbol{\Theta}_{t}$ has the following lower block diagonal structure

Since the innovation vector $Y_{t}-\hat{Y}_{t}$ in (2.13) is not a function of $\boldsymbol{\eta}_{t}$ and $\boldsymbol{\Theta}_{t}$, the control inputs $U_{t}$ in (POL) remain affine in the decision variables. This fact is important to show convexity of the optimization problem, as will be seen in the next section. Finally, the constraint (C1) can be rewritten as: