# Linear Receding Horizon Control with Probabilistic System Parameters

Raktim Bhattacharya  James Fisher  Aerospace Engineering, Texas A&M
University, USA.  Raytheon Missile Systems.

###### Abstract

In this paper we address the problem of designing receding horizon control
algorithms for linear discrete-time systems with parametric uncertainty. We do
not consider presence of stochastic forcing or process noise in the system. It
is assumed that the parametric uncertainty is probabilistic in nature with
known probability density functions. We use generalized polynomial chaos
theory to design the proposed stochastic receding horizon control algorithms.
In this framework, the stochastic problem is converted to a deterministic
problem in higher dimensional space. The performance of the proposed receding
horizon control algorithms is assessed using a linear model with two states.

##  1 Introduction

Receding horizon control (RHC), also known as model predictive control (MPC),
has been popular in the process control industry for several years Qin and
Badgwell (1996); Bemporad and Morari (1999), and recently gaining popularity
in aerospace applications, see Bhattacharya et al. (2002). It is based on the
idea of repetitive solution of an optimal control problem and updating states
with the first input of the optimal command sequence. The repetitive nature of
the algorithm results in a state dependent feedback control law. The
attractive aspect of this method is the ability to incorporate state and
control limits as constraints in the optimization formulation. When the model
is linear, the optimization problem is quadratic if the performance index is
expressed via a $\mathcal{L}_{2}$-norm, or linear if expressed via a
$\mathcal{L}_{1}/\mathcal{L}_{\infty}$-norm. Issues regarding feasibility of
online computation, stability and performance are largely understood for
linear systems and can be found in refs. Kwon (1994); Bitmead et al. (1990).
For nonlinear systems, stability of RHC methods is guaranteed by Primbs
(1999); Jadbabaie et al. (1999), by using an appropriate control Lyapunov
function . For a survey of the state-of-the-art in nonlinear receding horizon
control problems the reader is directed to Mayne et al. (2000a).

Traditional RHC laws perform best when modeling error is small. Fisher et al.
(2007) has shown that system uncertainty can lead to significant oscillatory
behavior and possibly instabilty. Furthermore, Grimm et al. (2004) showed that
in the presence of modeling uncertainty RHC strategy may not be robust with
RHC designs. Many approaches have been taken to improve robustness of RHC
strategy in the presence of unknown disturbances and bounded uncertainty, see
work of Raković et al. (2006); Lee and Yu (1997); Kouvaritakis et al. (2000);
Mayne et al. (2000b). These approaches involve the computation of a feedback
gain to ensure robustness. The difficulty with this approach is that, even for
linear systems, the problem becomes difficult to solve, as the unknown
feedback gain transforms the quadratic programming problem into a nonlinear
programming problem.

In this paper we address the problem of RHC design for linear systems with
probabilistic uncertainty in system parameters. Parametric uncertainty arises
in systems when the physics governing the system is known and the system
parameters are either not known precisely or are expected to vary in the
operational lifetime. Such uncertainty also occurs when system models are
build from experimental data using system identification techniques. As a
result of experimental measurements, the values of the parameters in the
system model have a range of variations with quantifiable likelihood of
occurrence. In either case, the range of variation of these parameters and the
likelihood of their occurrence are assumed to be known and it is desired to
design controllers that achieve specified performance for these variations.

While the area of robust RHC is not new, approaching the problem from a
stochastic standpoint is only recently receiving attention, for example van
Hessem and Bosgra (2002); Batina et al. (2002). These approaches however
suffered from either computational complexity, high degree of conservativeness
or do not address closed-loop stability. The key difficulty in stochastic RHC
is the propagation of uncertainty over the prediction horizon. More recently,
Cannon et al. (2009) avoid this difficulty by using an autonomous augmented
formulation of the prediction dynamics. Constraint satisfaction and stability
is achieved in Cannon et al. (2009) by extending ellipsoid invariance theory
to invariance with a given probability. The cost function minimized was the
expected value of a quadratic function of random state and control
trajectories. Additionally, the uncertainty in the system parameters were
assumed to have normal distribution.

This paper presents formulation of robust RHC design problems in polynomial
chaos framework, where parametric uncertainty can be governed by any
probability density function. In this approach the solution, not the dynamics,
of the random process is approximated using a series expansion. It is assumed
that the random process to be controlled has finite second moment, which is
the assumption of the polynomial chaos framework. The polynomial chaos based
approach predicts the propagation of uncertainty more accurately, is
computationally cheaper than methods based on Monte-Carlo or series
approximation of the dynamics, and is less conservative than the invariance
based methods.

The paper is organized as follows. We first present a brief introduction to
polynomial chaos and its application in transforming linear stochastic
dynamics to linear deterministic dynamics in higher dimensional state-space.
Next stability of stochastic linear dynamics in the polynomial chaos framework
is presented. This is followed by formulation of RHC design for discrete-time
stochastic linear systems. Stability of the proposed RHC algorithm is then
analyzed. The paper concludes with numerical examples that assesses the
performance of the proposed method.

##  2 Background on Polynomial Chaos

Recently, use of polynomial chaos to study stochastic differential equations
is gaining popularity. It is a non-sampling based method to determine
evolution of uncertainty in dynamical system, when there is probabilistic
uncertainty in the system parameters. Polynomial chaos was first introduced by
Wiener (1938) where Hermite polynomials were used to model stochastic
processes with Gaussian random variables. It can be thought of as an extension
of Volterra’s theory of nonlinear functionals Schetzen (2006) for stochastic
systems Ghanem and Spanos (1991). According to Cameron and Martin (1947) such
an expansion converges in the $\mathcal{L}_{2}$ sense for any arbitrary
stochastic process with finite second moment. This applies to most physical
systems. Xiu and Karniadakis (2002) generalized the result of Cameron-Martin
to various continuous and discrete distributions using orthogonal polynomials
from the so called Askey-scheme Askey and Wilson (1985) and demonstrated
$\mathcal{L}_{2}$ convergence in the corresponding Hilbert functional space.
This is popularly known as the generalized polynomial chaos (gPC) framework.
The gPC framework has been applied to applications including stochastic fluid
dynamics Hou et al. (2006),stochastic finite elements Ghanem and Spanos
(1991), and solid mechanics Ghanem and Red-Horse (1999). It has been shown in
Xiu and Karniadakis (2002) that gPC based methods are computationally far
superior than Monte-Carlo based methods. However, application of gPC to
control related problems has been surprisingly limited and is only recently
gaining popularity. See Prabhakar et al. (2008); Fisher and Bhattacharya
(2008a, b) for control related application of gPC theory.

###  2.1 Wiener-Askey Polynomial Chaos

Let $(\Omega,\mathcal{F},P)$ be a probability space, where $\Omega$ is the
sample space, $\mathcal{F}$ is the $\sigma$-algebra of the subsets of
$\Omega$, and $P$ is the probability measure. Let
$\Delta(\omega)=(\Delta_{1}(\omega),\cdots,\Delta_{d}(\omega)):(\Omega,\mathcal{F})\rightarrow(\mathbb{R}^{d},\mathcal{B}^{d})$
be an $\mathbb{R}^{d}$-valued continuous random variable, where
$d\in\mathbb{N}$, and $\mathcal{B}^{d}$ is the $\sigma$-algebra of Borel
subsets of $\mathbb{R}^{d}$. A general second order process
$X(\omega)\in\mathcal{L}_{2}(\Omega,\mathcal{F},P)$ can be expressed by
polynomial chaos as

$X(\omega)=\sum_{i=0}^{\infty}x_{i}\phi_{i}({\Delta}(\omega)),$ (1)

where $\omega$ is the random event and $\phi_{i}({\Delta}(\omega))$ denotes
the gPC basis of degree $p$ in terms of the random variables $\Delta(\omega)$.
The functions $\\{\phi_{i}\\}$ are a family of orthogonal basis in
$\mathcal{L}_{2}(\Omega,\mathcal{F},P)$ satisfying the relation

$\langle\phi_{i}\phi_{j}\rangle:=\int_{\mathcal{D}_{\Delta(\omega)}}{\phi_{i}\phi_{j}w(\Delta(\omega))\,d\Delta(\omega)}=h_{i}^{2}\delta_{ij}$
(2)

where $\delta_{ij}$ is the Kronecker delta, $h_{i}$ is a constant term
corresponding to
$\int_{\mathcal{D}_{\Delta}}{\phi_{i}^{2}w(\Delta)\,d\Delta}$,
$\mathcal{D}_{\Delta}$ is the domain of the random variable $\Delta(\omega)$,
and $w(\Delta)$ is a weighting function. Henceforth, we will use $\Delta$ to
represent $\Delta(\omega)$. For random variables $\Delta$ with certain
distributions, the family of orthogonal basis functions $\\{\phi_{i}\\}$ can
be chosen in such a way that its weight function has the same form as the
probability density function $f(\Delta)$. When these types of polynomials are
chosen, we have $f(\Delta)=w(\Delta)$ and

$\int_{\mathcal{D}_{\Delta}}{\phi_{i}\phi_{j}f(\Delta)\,d\Delta}=\mathbf{E}\left[\phi_{i}\phi_{j}\right]=\mathbf{E}\left[\phi_{i}^{2}\right]\delta_{ij},$
(3)

where $\mathbf{E}\left[\cdot\right]$ denotes the expectation with respect to
the probability measure $dP(\Delta(\omega))=f(\Delta(\omega))d\Delta(\omega)$
and probability density function $f(\Delta(\omega))$. The orthogonal
polynomials that are chosen are the members of the Askey-scheme of polynomials
(Askey and Wilson (1985)), which forms a complete basis in the Hilbert space
determined by their corresponding support. Table 1 summarizes the
correspondence between the choice of polynomials for a given distribution of
$\Delta$. See Xiu and Karniadakis (2002) for more details.

Random Variable $\Delta$ |  $\phi_{i}(\Delta)$ of the Wiener-Askey Scheme
---|---
Gaussian | Hermite
Uniform | Legendre
Gamma | Laguerre
Beta | Jacobi
Table 1: Correspondence between choice of polynomials and given distribution
of $\Delta(\omega)$ Xiu and Karniadakis (2002).

###  2.2 Approximation of Stochastic Linear Dynamics Using Polynomial Chaos
Expansions

Here we derive a generalized representation of the deterministic dynamics
obtained from the stochastic system by approximating the solution with
polynomial chaos expansions.

Define a linear discete-time stochastic system in the following manner

$x(k+1,\Delta)=A(\Delta)x(k,\Delta)+B(\Delta)u(k,\Delta),$ (4)

where $x\in\mathbb{R}^{n},u\in\mathbb{R}^{m}$. The system has probabilistic
uncertainty in the system parameters, characterized by $A(\Delta),B(\Delta)$,
which are matrix functions of random variable
$\Delta\equiv\Delta(\omega)\in\mathbb{R}^{d}$ with certain stationary
distributions. Due to the stochastic nature of $(A,B)$, the system trajectory
$x(k,\Delta)$ will also be stochastic.

By applying the Wiener-Askey gPC expansion of finite order to
$x(k,\Delta),A(\Delta)$ and $B(\Delta)$, we get the following approximations,

$\displaystyle\hat{x}(k,\Delta)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{p}x_{i}(k)\phi_{i}(\Delta),\,x_{i}(k)\in\mathbb{R}^{n}$
(5) $\displaystyle\hat{u}(k,\Delta)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{p}u_{i}(k)\phi_{i}(\Delta),\,u_{i}(k)\in\mathbb{R}^{m}$
(6) $\displaystyle\hat{A}(\Delta)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{p}A_{i}\phi_{i}(\Delta),\,A_{i}=\frac{\langle
A(\Delta),\phi_{i}(\Delta)\rangle}{\langle\phi_{i}(\Delta)^{2}\rangle}\in\mathbb{R}^{n\times
n}$ (7) $\displaystyle\hat{B}(\Delta)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{p}B_{i}\phi_{i}(\Delta),\,B_{i}=\frac{\langle
B(\Delta),\phi_{i}(\Delta)\rangle}{\langle\phi_{i}(\Delta)^{2}\rangle}\in\mathbb{R}^{n\times
m}.$ (8)

The inner product or ensemble average $\langle\cdot,\cdot\rangle$, used in the
above equations and in the rest of the paper, utilizes the weighting function
associated with the assumed probability distribution, as listed in table 1.

The number of terms $p$ is determined by the dimension $d$ of $\Delta$ and the
order $r$ of the orthogonal polynomials $\\{\phi_{k}\\}$, satisfying
$p+1=\frac{(d+r)!}{d!r!}$. The $n(p+1)$ time varying coefficients,
$\\{x_{i}(k)\\};k=0,\cdots,p$, are obtained by substituting the approximated
solution in the governing equation (eqn.(4)) and conducting Galerkin
projection on the basis functions $\\{\phi_{k}\\}_{k=0}^{p}$, to yield
$n(p+1)$ deterministic linear system of equations, which given by

$\mathbf{X}(k+1)=\mathbf{A}\mathbf{X}(k)+\mathbf{B}\mathbf{U}(k),$ (9)

where

$\displaystyle\mathbf{X}(k)$ $\displaystyle=$
$\displaystyle[x_{0}(k)^{T}\;x_{1}(k)^{T}\;\cdots x_{p}(k)^{T}]^{T},$ (10)
$\displaystyle\mathbf{U}(k)$ $\displaystyle=$
$\displaystyle[u_{0}(k)^{T}\;u_{1}(k)^{T}\;\cdots u_{p}(k)^{T}]^{T}.$ (11)

Matrices $\mathbf{A}\in\mathbb{R}^{n(p+1)\times n(p+1)}$ and
$\mathbf{B}\in\mathbb{R}^{n(p+1)\times m}$ are defined as

$\displaystyle\mathbf{A}$ $\displaystyle=$ $\displaystyle(W\otimes
I_{n})^{-1}\left[\begin{array}[]{c}H_{A}(E_{0}\otimes I_{n})\\\ \vdots\\\
H_{A}(E_{p}\otimes I_{n})\end{array}\right],$ (16) $\displaystyle\mathbf{B}$
$\displaystyle=$ $\displaystyle(W\otimes
I_{n})^{-1}\left[\begin{array}[]{c}H_{B}(E_{0}\otimes I_{m})\\\ \vdots\\\
H_{B}(E_{p}\otimes I_{m})\end{array}\right],$ (20)

where $H_{A}=\left[A_{0}\,\cdots\,A_{p}\right]$,
$H_{B}=\left[B_{0}\,\cdots\,B_{p}\right]$,
$W=diag(\langle\phi_{0}^{2}\rangle,\cdots,\langle\phi_{p}^{2}\rangle)$, and

$E_{i}=\left[\begin{array}[]{ccc}\langle\phi_{i},\phi_{0},\phi_{0}\rangle&\cdots&\langle\phi_{i},\phi_{0},\phi_{p}\rangle\\\
\vdots&&\vdots\\\
\langle\phi_{i},\phi_{p},\phi_{0}\rangle&\cdots&\langle\phi_{i},\phi_{p},\phi_{p}\rangle\end{array}\right],$

with $I_{n}$ and $I_{m}$ as the identity matrix of dimension $n\times n$ and
$m\times m$ respectively. It can be easily shown that
$\mathbf{E}\left[x(k)\right]=x_{0}(k)$, or
$\mathbf{E}\left[x(k)\right]=\left[I_{n}\;0_{n\times np}\right]\mathbf{X}(k).$

Therefore, transformation of a stochastic linear system with
$x\in\mathbb{R}^{n},u\in\mathbb{R}^{m}$, with $p^{th}$ order gPC expansion,
results in a deterministic linear system with increased dimensionality equal
to $n(p+1)$.

##  3 Stochastic Receding Horizon Control

Here we develop a RHC methodology for stochastic linear systems similar to
that developed for deterministic systems, presented by Goodwin et al. (2005).
Let $x(k,\Delta)$ be the solution of the system in eqn.(4) with control
$u(k,\Delta)$. Consider the following optimal control problem defined by,

$\displaystyle V_{N}^{*}=\min\,V_{N}(\\{x(k+1,\Delta)\\},\\{u(k,\Delta)\\})$
(21) $\displaystyle{\rm subject\,to:}$ $\displaystyle
x(k+1,\Delta)=A(\Delta)x(k,\Delta)+B(\Delta)u(k,\Delta),$ (22)
$\displaystyle\textrm{Initial Condition: }x(0,\Delta);$ (23)
$\displaystyle\mu(u(k,\Delta))\in\mathbb{U}\subset\mathbb{R}^{m},$ (24)
$\displaystyle\mu(x(k,\Delta))\in\mathbb{X}\subset\mathbb{R}^{n},$ (25)
$\displaystyle\mu(x(N,\Delta))\in\mathbb{X}_{f}\subset\mathbb{X},$ (26)

for $k=0,\cdots,N-1$; where $N$ is the horizon length, $\mathbb{U}$ and
$\mathbb{X}$ are feasible sets for $u(k,\Delta)$ and $x(k,\Delta)$ with
respect to control and state constraints. $\mu(\cdot)$ represents moments
based constraints on state and control.The set $\mathbb{X}_{f}$ is a terminal
constraint set. The cost function $V_{N}$ is given by

$\begin{split}&V_{N}=\sum_{k=1}^{N}\mathbf{E}\left[x^{T}(k,\Delta)Qx(k,\Delta)+\right.\\\
&\left.u^{T}(k-1,\Delta)Ru(k-1,\Delta)\right]+C_{f}(x(N),\Delta),\end{split}$
(27)

where $C_{f}(x(N),\Delta)$ is a terminal cost function, and $Q=Q^{T}>0$,
$R=R^{T}>0$ are matrices with appropriate dimensions.

###  3.1 Control Structure

Here we consider the control structure,

$u(k,\Delta)=\bar{u}(k)+K(k)\left(x(k,\Delta)-\mathbf{E}\left[x(k,\Delta)\right]\right),$
(28)

where $\bar{u}(k)$, and $K(k)$ are unknown deterministic quantities. This is
similar to that proposed by Primbs et al.Primbs and Sung (2009) and enables us
to regulate the mean trajectory using open loop control and deviations about
the mean using a state-feedback control.

In terms of gPC coefficients, the system dynamics in eqn.(4) with the first
control structure is given by eqn.(9). The system dynamics in term of the gPC
expansions, with the second control structure is given by

$\mathbf{X}(k+1)=(\mathbf{A}+\mathbf{B}(\mathbf{M}\otimes
K(k)))\mathbf{X}(k)+\mathbf{B}\bar{U}(k),$ (29)

where $\bar{U}(k)=[1\,\,0_{1\times p}]^{T}\otimes\bar{u}(k)$ and
$\mathbf{M}=\left[\begin{array}[]{cc}0&0_{1\times p}\\\ 0_{p\times
1}&I_{p\times p}\end{array}\right]$.

###  3.2 Cost Functions

Here we derive the cost function in eqn.(27) is derived in terms of the gPC
coefficients $\mathbf{X}$ and $\mathbf{U}$. For scalar $x$, the quantity
$\mathbf{E}\left[x^{2}\right]$ in terms of its gPC expansions is given by

$\mathbf{E}\left[x^{2}\right]=\sum_{i=0}^{p}\sum_{j=0}^{p}x_{i}x_{j}\int_{\mathcal{D}_{\Delta}}\phi_{i}\phi_{j}fd\Delta=\mathbf{x}^{T}W\mathbf{x},$
(30)

where $\mathcal{D}_{\Delta}$ is the domain of $\Delta$ , $x_{i}$ are the gPC
expansions of $x$, $f\equiv f(\Delta)$ is the probability distribution of
$\Delta$. Here we use the notation $\mathbf{x}$ to represent the gPC state
vector for scalar $x$. The expression $\mathbf{E}\left[x^{2}\right]$ can be
generalized for $x\in\mathbb{R}^{n}$ where $\mathbf{E}\left[x^{T}Qx\right]$ is
given by

$\mathbf{E}\left[x^{T}Qx\right]=\mathbf{X}^{T}(W\otimes Q)\mathbf{X}.$ (31)

The expression for the cost function in eqn.(27), in terms of gPC states and
control is

$\begin{split}&V_{N}=\sum_{k=0}^{N-1}[\mathbf{X}^{T}(k)\bar{Q}\mathbf{X}(k)+\\\
&(\bar{U}^{T}(k)+\mathbf{X}^{T}(k)(\mathbf{M}\otimes
K^{T}(k)))\bar{R}(\bar{U}(k)+\\\ &(\mathbf{M}\otimes
K(k)))\mathbf{X}(k))]+C_{f}(x(N),\Delta),\end{split}$ (32)

where $\bar{Q}=W\otimes Q$ and $\bar{R}=W\otimes R$.

In deterministic RHC, the terminal cost is the cost-to-go from the terminal
state to the origin by the local controller Goodwin et al. (2005). In the
stochastic setting, a local controller can be synthesized using methods
presented in our previous work Fisher and Bhattacharya (2008a). The cost-to-go
from a given stochastic state variable $x(N,\Delta)$ can then be written as

$C_{f}(x(N),\Delta)=\mathbf{X}^{T}(N)P\mathbf{X}(N),$ (33)

where $\mathbf{X}(N)$ are gPC states corresponding to $x(N,\Delta)$ and
$P=P^{T}>0$ is a $n(p+1)\times n(p+1)$-dimensional matrix, obtained from the
synthesis of the terminal control law Fisher and Bhattacharya (2008a). In the
current stochastic RHC literature, the terminal cost function has been defined
on the expected value of the final state Lee and Cooley (1998); de la Penad et
al. (2005); Primbs and Sung (2009); Bertsekas (2005) or using a combination of
mean and variance Darlington et al. (2000); Nagy and Braatz (2003). The
terminal cost function in eqn.(33) is more general than the terminal cost
functions used in the literature because it penalizes all the moments of the
random variable $x(N,\Delta)$, as they are functions of $\mathbf{X}(N)$. This
can be shown as follows.

To avoid tensor notation and without loss of generality, we consider
$x(k,\Delta)\in\mathbb{R}$ and let
$\mathbf{X}(k)=[x_{0}(k),x_{1}(k),\cdots,x_{p}(k)]^{T}$ be the gPC expansion
of $x(k,\Delta)$. The $p^{th}$ moment in terms of $x_{i}(k)$ are then given by

$\begin{split}&m_{p}(k)=\sum_{i_{1}=0}^{P}\cdots\sum_{{i_{p}}=0}^{P}x_{i_{1}}(k)\cdots
x_{i_{p}}(k)\int_{\mathcal{D}_{\Delta}}\phi_{i_{1}}(\Delta)\cdots\\\
&\phi_{i_{p}}(\Delta)f(\Delta)d\Delta.\end{split}$ (34)

Thus, minimizing $C_{f}(x(N),\Delta)$ in eqn.(33) minimizes all moments of
$x(N,\Delta)$. Consequently, constraining the probability density function of
$x(N,\Delta)$.

###  3.3 State and Control Constraints

In this section we present the state and control constraints for the receding
horizon policy.

####  3.3.1 Expectation Based Constraints

Here we first consider constraints of the following form,

$\displaystyle\mathbf{E}\left[x(k,\Delta)^{T}H_{x}x(k,\Delta)+G_{x}x(k,\Delta)\right]$
$\displaystyle\leq$ $\displaystyle\alpha_{i,x},$ (35)
$\displaystyle\mathbf{E}\left[u(k,\Delta)^{T}H_{u}u(k,\Delta)+G_{u}u(k,\Delta)\right]$
$\displaystyle\leq$ $\displaystyle\alpha_{i,u},$ (36)

for $k=0\ldots N$. These constraints are on the expected value of the
quadratic functions. Thus, instead of requiring that the constraint be met for
all trajectories, they instead imply that the constraints should be satisfied
on average. These constraints can be expressed in terms of the gPC states as

$\displaystyle\mathbf{X}(k)^{T}\bar{H}_{x}\mathbf{X}(k)+\bar{G}_{x}\mathbf{X}(k)$
$\displaystyle\leq$ $\displaystyle\alpha_{i,x},$ (37)
$\displaystyle\mathbf{U}(k)^{T}\bar{H}_{u}\mathbf{U}(k)+\bar{G}_{u}\mathbf{U}(k)$
$\displaystyle\leq$ $\displaystyle\alpha_{i,u},$ (38)

where $\bar{H}_{x}=W\otimes H_{x}$, $\bar{H}_{u}=W\otimes H_{u}$,
$\bar{G}_{x}=G_{x}\left[I_{n}\;0_{n\times np}\right]$, and
$\bar{G}_{u}=G_{u}\left[I_{n}\;0_{n\times np}\right]$.

####  3.3.2 Variance Based Constraints

In many practical applications, it may be desirable to constrain the second
moment of the state trajectories, either at each time step or at final time.
One means of achieving this is to use a constraint of the form

$\mathbf{Tr}\left[\mathbf{E}\left[(x(k)-\mathbf{E}\left[x(k)\right])(x(k)-\mathbf{E}\left[x(k)\right])^{T}\right]\right]\leq\alpha_{\sigma^{2}}.$
(39)

For scalar $x$, the variance $\sigma^{2}(x)$ in terms of the gPC expansions
can be shown to be

$\sigma^{2}=\mathbf{E}\left[x-\mathbf{E}\left[x\right]\right]^{2}=\mathbf{E}\left[x^{2}\right]-\mathbf{E}\left[x\right]^{2}=\mathbf{x}^{T}W\mathbf{x}-\mathbf{E}\left[x\right]^{2},$

where

$\begin{split}&\mathbf{E}\left[x\right]=\mathbf{E}\left[\sum_{i=0}^{p}x_{i}\phi_{i}\right]=\sum_{i=0}^{p}x_{i}\mathbf{E}\left[\phi_{i}\right]=\sum_{i=0}^{p}x_{i}\int_{\mathcal{D}_{\Delta}}\phi_{i}fd\Delta\\\
&=\mathbf{x}^{T}F,\end{split}$

and $F=\left[\begin{array}[]{cccc}1\;0\;\cdots\;0\end{array}\right]^{T}$.
Therefore, $\sigma^{2}$ for scalar $x$ can be written in a compact form as

$\sigma^{2}=\mathbf{x}^{T}(W-FF^{T})\mathbf{x}.$ (40)

In order to represent the covariance for $x\in\mathbb{R}^{n}$, in terms of the
gPC states, let us define $\Phi=[\phi_{0}\cdots\phi_{p+1}]^{T}$ and write
$G=\int_{\mathcal{D}_{\Delta}}\Phi\Phi^{T}fd\Delta$. Let us represent a sub-
vector of $\mathbf{X}$, defined by elements $n_{1}$ to $n_{2}$, as
$X_{n_{1}\cdots n{2}}$, where $n_{1}$ and $n_{2}$ are positive integers. Let
us next define matrix $M_{\mathbf{X}}$, with subvectors of $\mathbf{X}$, as
$M_{\mathbf{X}}=[\mathbf{X}_{1\cdots n}\;\mathbf{X}_{n+1\cdots
2n}\;\cdots\mathbf{X}_{np+1\cdots n(p+1)}]$. For $x\in\mathbb{R}^{n}$, it can
be shown that

$\mathbf{E}\left[x\right]=(F\otimes I_{n})\mathbf{X},$ (41)

and the covariance can then be shown to be

$\textbf{Cov}(x)=M_{\mathbf{X}}GM_{\mathbf{X}}^{T}-(F\otimes
I_{n})\mathbf{X}\mathbf{X}^{T}(F^{T}\otimes I_{n}).$ (42)

The trace of the covariance matrix $\textbf{Cov}(x)$ can then be written as

$\mathbf{Tr}\left[\textbf{Cov}(x)\right]=\mathbf{X}^{T}((W-FF^{T})\otimes
I_{n})\mathbf{X}.$

Therefore, a constraint of the type

$\mathbf{Tr}\left[\textbf{Cov}(x(k))\right]\leq\alpha_{\sigma^{2}}$

can be written in term of gPC states as

$\mathbf{X}^{T}Q_{\sigma^{2}}\mathbf{X}\leq\alpha_{\sigma^{2}},$ (43)

where $Q_{\sigma^{2}}=(W-FF^{T})\otimes I_{n}$.

##  4 Stability of the RHC Policy

Here we show the stability properties of the receding horizon policy when it
is applied to the system in eqn.(9). Using gPC theory we can convert the
underlying stochastic RHC formulation in $x(t,\Delta)$ and $u(t,\Delta)$ to a
deterministic RHC formulation in $\mathbf{X}(k)$ and $\mathbf{U}(k)$. The
stability of $\mathbf{X}(k)$ in an RHC setting, with suitable terminal
controller, can be proved using results by Goodwin et al. (2005), which shows
that $\lim_{k\rightarrow\infty}\mathbf{X}(k)\rightarrow 0$, when a receding
horizon policy is employed. To relate this result to the stability of
$x(k,\Delta)$, we first present the following known result in stochastic
stability in terms of the moments of $x(k,\Delta)$. For stochastic dynamical
systems in general, stability of moments is a weaker definition of stability
than the almost sure stability definition. However, the two definitions are
equivalent for linear autonomous systems (pg. 296, Khas’minskii (1969) also
pg. 349 Chen and Hsu (1995)). Here we present the definition of asymptotic
stability in the $p^{th}$ moment for discrete-time systems.

###### Definition 1

The zero equilibrium state is said to be stable in the $p^{th}$ moment if
$\forall\epsilon>0,\,\exists\delta>0$ such that

$\sup_{k\geq 0}\mathbf{E}\left[x(k,\Delta)^{p}\right]\leq\epsilon,\;\forall
x(0,\Delta):||x(0,\Delta)||\leq\delta,\forall\Delta\in\mathcal{D}_{\Delta}.$
(44)

###### Definition 2

The zero equilibrium state is said to be asymptotically stable in the $p^{th}$
moment if it is stable in $p^{th}$ moment and

$\lim_{k\rightarrow\infty}\mathbf{E}\left[x(k,\Delta)^{p}\right]=0,$ (45)

for all $x(0,\Delta)$ in the neighbourhood of the zero equilibrium.

###### Proposition 1

For the system in eqn.(4), $\lim_{k\rightarrow\infty}\mathbf{X}(k)\rightarrow
0$ is a sufficient condition for the asymptotic stability of the zero
equilibrium state, in all moments.

######  Proof 4.1.

To avoid tensor notation and without loss of generality, we consider
$x(k,\Delta)\in\mathbb{R}$ and let
$\mathbf{X}(k)=[x_{0}(k),x_{1}(k),\cdots,x_{p}(k)]^{T}$ be the gPC expansion
of $x(k,\Delta)$. The moments in terms of $x_{i}(k)$ are given by eqn.(34).
Therefore, if $\lim_{k\rightarrow\infty}\mathbf{X}(k)\rightarrow 0$
$\implies\lim_{k\rightarrow\infty}x_{i}(k)\rightarrow 0$. Consequently,
$\lim_{k\rightarrow\infty}m_{i}(k)\rightarrow 0$ for $i=1,2,\cdots$, and
eqn.(45) is satisfied. This completes the proof. $\square$

##  5 Numerical Example

Here we consider the following linear system, similar to that considered in
Primbs and Sung (2009),

$x(k+1)=(A+G(\Delta))x(k)+Bu(k)$ (46)

where

$A=\left[\begin{array}[]{cc}1.02&-0.1\\\
.1&.98\end{array}\right],\,B=\left[\begin{array}[]{c}0.1\\\
0.05\end{array}\right],\,G=\left[\begin{array}[]{cc}0.04&0\\\
0&0.04\end{array}\right]\Delta.$

The system in consideration is open-loop unstable and the uncertainty appears
linearly in the $G$ matrix. Here, $\Delta\in[-1,1]$ and is governed by a
uniform distribution, that doesn’t change with time. Consequently, Legendre
polynomials is used for gPC approximation and polynomials up to $4^{th}$ order
are used to formulate the control. Additionally, we assume that there is no
uncertainty in the initial condition. The expectation based constraint is
imposed on $x(k,\Delta))$ as

$\mathbf{E}\left[\;[1\;\;0]x(k,\Delta)\;\right]\geq-1,$

which in terms of the gPC states, this corresponds to

$\left[\begin{array}[]{cc}1&\mathbf{0}_{1\times
2p+1}\end{array}\right]\mathbf{X}(k)\geq-1.$

The terminal controller is designed using probabilistic LQR design techniques
described by Fisher and Bhattacharya (2008a). The cost matrices used to
determine the terminal controller are

$Q=\left[\begin{array}[]{cc}2&0\\\ 0&5\end{array}\right],\,R=1.$

Figure (1) illustrates the performance of the proposed RHC policy The
resulting optimization problem is a nonlinear programming problem which has
been solved using MATLAB’s fmincon(...) function. From the figure, we see that
the constraint on the expected value of $x_{1}$ has been satisfied and the RHC
algorithm was able to stabilize the system. These plots have been obtained
using $4^{th}$ order gPC approximation of the stochastic dynamics.

##  6 Summary

In this paper we present a RHC strategy for linear discrete time systems with
probabilistic system parameters. We have used the polynomial chaos framework
to design stochastic RHC algorithms in an equivalent deterministic setting.
The controller structure has an open loop component that controls the mean
behavior of the system, and a state-feedback component that controls
deviations about the mean trajectory. This controller structure results in a
polynomial optimization problem with polynomial constraints that is solved in
the general nonlinear programming framework. Theoretical guarantees for the
stability of the proposed algorithm has also been presented. Performance of
the RHC algorithm has been assessed using a two dimensional dynamical system.

Figure 1: State trajectories with expectation constraints.

## References

  * Askey and Wilson [1985] R. Askey and J. Wilson.  Some basic hypergeometric polynomials that generalize jacobi polynomials.  _Memoirs Amer. Math. Soc._ , 319, 1985.
  * Batina et al. [2002] I. Batina, A. A. Stoorvogel, and S. Weiland.  Optimal control of linear, stochastic systems with state and input constraints.  _Proceedings of the 41st IEEE Conference on Decision and Control_ , 2:1564– 1569, 2002.
  * Bemporad and Morari [1999] A. Bemporad and M. Morari.  Robust model predictive control: A survey.  Technical report, Automatic Control Laboratory, Swiss Federal Institute of Technology (ETH), Physikstrasse 3, CH-8092 Zürich, Switzerland, www.control.ethz.ch, 1999.
  * Bertsekas [2005] D. Bertsekas.  Dynamic programming and suboptimal control: A survey from ADP to MPC.  _European Journal of Control_ , 11(4-5):310–334, 2005.
  * Bhattacharya et al. [2002] Raktim Bhattacharya, Gary J. Balas, M. Alpay Kaya, and Andy Packard.  Nonlinear receding horizon control of an f-16 aircraft.  _Journal of Guidance, Control, and Dynamics_ , 25(5):924–931, 2002.
  * Bitmead et al. [1990] R.R. Bitmead, M. Gevers, and V. Wertz.  _Adaptive Optimal Control: The Thinking Man’s GPC_.  International Series in Systems and Control Engineering. Prentice Hall, Englewood Cliffs, NJ, 1990.
  * Cameron and Martin [1947] R. H. Cameron and W. T. Martin.  The orthogonal development of non-linear functionals in series of fourier-hermite functionals.  _The Annals of Mathematics_ , 48(2):385–392, 1947\.
  * Cannon et al. [2009] M. Cannon, B. Kouvaritakis, and X. Wu.  Model predictive control for systems with stochastic multiplicative uncertainty and probabilistic constraints.  _Automatica_ , 45(1):167 – 172, 2009.
  * Chen and Hsu [1995] By G. Chen and S. H. Hsu.  _Linear Stochastic Control Systems_.  CRC Press, 1995.
  * Darlington et al. [2000] J. Darlington, CC Pantelides, B. Rustem, and BA Tanyi.  Decreasing the sensitivity of open-loop optimal solutions in decision making under uncertainty.  _European Journal of Operational Research_ , 121(2):343–362, 2000.
  * de la Penad et al. [2005] DM de la Penad, A. Bemporad, and T. Alamo.  Stochastic programming applied to model predictive control.  In _44th IEEE Conference on Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC’05_ , pages 1361–1366, 2005.
  * Fisher and Bhattacharya [2008a] J. Fisher and R. Bhattacharya.  On stochastic LQR design and polynomial chaos.  In _American Control Conference, 2008_ , pages 95–100, 2008a.
  * Fisher and Bhattacharya [2008b] J. Fisher and R. Bhattacharya.  Optimal Trajectory Generation with Probabilistic System Uncertainty Using Polynomial Chaos.  _In Press Journal of Dynamic Systems, Measurement and Control_ , 2008b.
  * Fisher et al. [2007] James Fisher, Raktim Bhattacharya, and S. R. Vadali.  Spacecraft momentum management and attitude control using a receding horizon approach.  In _Proceedings of the 2007 AIAA Guidance, Navigation, and Control Conference and Exhibit_ , Hilton Head, SC, August 2007. AIAA.
  * Ghanem and Red-Horse [1999] Roger Ghanem and John Red-Horse.  Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach.  _Phys. D_ , 133(1-4):137–144, 1999.  ISSN 0167-2789.  http://dx.doi.org/10.1016/S0167-2789(99)00102-5.
  * Ghanem and Spanos [1991] Roger G. Ghanem and Pol D. Spanos.  _Stochastic Finite Elements: A Spectral Approach_.  Springer-Verlag Inc., New York, NY, 1991.  ISBN 0-387-97456-3.
  * Goodwin et al. [2005] G.C. Goodwin, M. Seron, and J. De Dona.  _Constrained control and estimation: an optimisation approach_.  Springer, 2005.
  * Grimm et al. [2004] Gene Grimm, Michael J. Messina, Sezai E. Tuna, and Andrew R. Teel.  Examples when nonlinear model predictive control is nonrobust.  _Automatica_ , 40:1729–1738, 2004.
  * Hou et al. [2006] Thomas Y. Hou, Wuan Luo, Boris Rozovskii, and Hao-Min Zhou.  Wiener chaos expansions and numerical solutions of randomly forced equations of fluid mechanics.  _J. Comput. Phys._ , 216(2):687–706, 2006.  ISSN 0021-9991.  http://dx.doi.org/10.1016/j.jcp.2006.01.008.
  * Jadbabaie et al. [1999] A. Jadbabaie, J. Yu, and J. Hauser.  Stabilizing receding horizon control of nonlinear systems: A control lyapunov function approach.  _Proceedings of the 1999 American Control Conference_ , 3:1535–1539, 1999.
  * Khas’minskii [1969] R. Z. Khas’minskii.  _Stability of Systems of Differential Equations in the Presence of Random Disturbances (in Russian)_.  Nauka, Moscow, 1969.
  * Kouvaritakis et al. [2000] B. Kouvaritakis, J. A. Rossiter, and J. Schuurmans.  Efficient robust predictive control.  _IEEE Transactions on Automatic Control_ , 45(8):1545–1549, 2000.
  * Kwon [1994] W.H. Kwon.  Advances in predictive control: Theory and application.  _1st Asian Control Conference, Tokyo_ , 1994.
  * Lee and Yu [1997] J. H. Lee and Zhenghong Yu.  Worst-case formulations of model predictive control for systems with bounbded parameters.  _Automatica_ , 33(5):763–781, 1997.
  * Lee and Cooley [1998] JH Lee and BL Cooley.  Optimal feedback control strategies for state-space systems with stochastic parameters.  _IEEE Transactions on Automatic Control_ , 43(10):1469–1475, 1998.
  * Mayne et al. [2000a] D. Mayne, J. Rawlings, C. Rao, and P. Scokaert.  Constrained model predictive control, stability and optimality.  _Automatica_ , 36:789–814, 2000a.
  * Mayne et al. [2000b] D.Q. Mayne, J.B. Rawlings, C.V. Rao, and PO Scokaert.  Constrained model predictive control: Stability and optimality.  _AUTOMATICA-OXFORD-_ , 36:789–814, 2000b.
  * Nagy and Braatz [2003] Z.K. Nagy and R.D. Braatz.  Robust nonlinear model predictive control of batch processes.  _AIChE Journal_ , 49(7):1776–1786, 2003.
  * Prabhakar et al. [2008] A. Prabhakar, J. Fisher, and R. Bhattacharya.  Polynomial Chaos Based Analysis of Probabilistic Uncertainty in Hypersonic Flight Dynamics.  _submitted AIAA Journal of Guidance, Control, and Dynamics_ , 2008\.
  * Primbs [1999] J.A. Primbs.  _Nonlinear Optimal Control: A Receding Horizon Approach_.  PhD thesis, California Institute of Technology, Pasadena, CA, 1999.
  * Primbs and Sung [2009] J.A. Primbs and C.H. Sung.  Stochastic receding horizon control of constrained linear systems with state and control multiplicative noise.  _Automatic Control, IEEE Transactions on_ , 54(2):221–230, 2009.
  * Qin and Badgwell [1996] S.J. Qin and T. Badgwell.  An overview of industrial model predictive control technology.  _AIChE Symposium Series_ , 93:232–256, 1996.
  * Raković et al. [2006] S. V. Raković, A. R. Teel, D. Q. Mayne, and A. Astolfi.  Simple robust control invariant tubes for some classes of nonlinear discrete time systems.  In _Proceedings of the 45th IEEE Conference on Decision and Control_ , pages 6397–6402, San Diego, CA, December 2006. IEEE.
  * Schetzen [2006] M. Schetzen.  _The Volterra and Wiener Theories of Nonlinear Systems_.  Krieger Pub., 2006.
  * van Hessem and Bosgra [2002] D. H. van Hessem and O. H. Bosgra.  A conic reformulation of model predictive control including bounded and stochastic disturbances and input constraints.  In _Proceedings of the 2002 Conference on Decision and Control_ , volume 4, pages 4643–4648, Las Vegas, NV, December 2002. IEEE.
  * Wiener [1938] N. Wiener.  The homogeneous chaos.  _American Journal of Mathematics_ , 60(4):897–936, 1938.
  * Xiu and Karniadakis [2002] Dongbin Xiu and George Em Karniadakis.  The wiener–askey polynomial chaos for stochastic differential equations.  _SIAM J. Sci. Comput._ , 24(2):619–644, 2002\.  ISSN 1064-8275.