Recursively Feasible Stochastic Model Predictive Control using Indirect Feedback

TL;DR

This paper introduces a stochastic MPC method for linear systems with unbounded, correlated disturbances, ensuring recursive feasibility and providing performance bounds, demonstrated through building control examples.

Contribution

The paper proposes a recursive feasibility initialization scheme for stochastic MPC that handles correlated disturbances and guarantees chance constraint satisfaction.

Findings

Recursive feasibility is maintained in stochastic MPC with correlated disturbances.

The method achieves an average asymptotic performance bound under i.i.d. disturbance assumptions.

Illustrative examples demonstrate the approach's effectiveness in building control scenarios.

Abstract

We present a stochastic model predictive control (MPC) method for linear discrete-time systems subject to possibly unbounded and correlated additive stochastic disturbance sequences. Chance constraints are treated in analogy to robust MPC using the concept of probabilistic reachable sets for constraint tightening. We introduce an initialization of each MPC iteration which is always recursively feasibility and thereby allows that chance constraint satisfaction for the closed-loop system can readily be shown. Under an i.i.d. zero mean assumption on the additive disturbance, we furthermore provide an average asymptotic performance bound. Two examples illustrate the approach, highlighting feedback properties of the novel initialization scheme, as well as the inclusion of time-varying, correlated disturbances in a building control setting.

Tables

Table 1. TABLE I: Comparison for Double Integrator

Controller	nom	rec	df
${\overline{J}}_{cl}\left(x\left(0\right)\right)$	$10.41$	$7.76$	$7.71$
${\overline{J}}_{cl}\left(x\left(20\right)\right)$	$7.2$	$4.33$	$4.28$
$\max {\overline{n}}_v\left(k\right)$	$10.3\%$	$8.00\%$	$6.9\%$

Table 2. TABLE II: Comparison for Double Integrator with Model Mismatch

Controller	nom	rec	df	recSC
${\overline{J}}_{cl}\left(x\left(0\right)\right)$	$130.10$	$91.22$	$145.68$	$116.87$
${\overline{J}}_{cl}\left(x\left(20\right)\right)$	$135.79$	$93.24$	$156.03$	$120.42$
${\max }_k{\overline{n}}_v\left(k\right)$	$21.70\%$	$19.8\%$	$11.0\%$	$11.3\%$

Equations

$$x(k+1)$$

$$\mathrm{Pr}\left(x(k)\in \mathcal{X}|x(0)\right)\ge p_x,$$

$$\mathrm{Pr}\left(u(k)\in \mathcal{U}|x(0)\right)\ge p_u,$$

$$\{e(k){\}}_{k\in 0,\dots ,\bar{N}},$$

$$\mathrm{Pr}(e(n)\in {\mathcal{R}}_n|e(0))\ge p.$$

$$\mathrm{Pr}(e(n)\in \mathcal{R}|e(0))\ge p\forall 0\le n\le \bar{N}.$$

$$\mathcal{R}\supseteq \bigcup _{n=0}^{\bar{N}}{\mathcal{R}}_n.$$

$$e(k+1)=A_{K}e(k)+w(k).$$

$$\mathcal{R}\supseteq {\mathcal{R}}_{\bar{N}}.$$

$$E(E)=A_{0}e(0)+\bar{A}{\mu }^W,\mathrm{var}(E)=\bar{A}{\Sigma }^W{\bar{A}}^{\mathsf{T}},$$

$$A_0=\left[\begin{array}{c}{A}_K\\ A_K^2\\ ⋮\\ A_K^{\bar{N}-1}\end{array}\right],\bar{A}=\left[\begin{array}{cccc}I& & & \\ A_K& I& & \\ ⋮& & \ddots & \\ A_K^{\bar{N}-2}& A_K^{\bar{N}-3}& \dots & I\end{array}\right].$$

$${\mathcal{R}}_n^c:=\left\{e|(e-E(e(n)){)}^{\mathsf{T}}\mathrm{var}{(e(n))}^{-1}(e-E(e(n)))<\tilde{p}\right\}$$

$$\mathrm{Pr}\left(e(n)\in {\mathcal{R}}_n|e(0)\right)\ge p$$

$$\Rightarrow$$

$$x_{i+1}=A{x}_i+B{u}_i+w_i,$$

$$p(W_k)=$$

$$p\left({[w(k{)}^{\mathsf{T}},\dots ,w(k+N{)}^{\mathsf{T}}]}^{\mathsf{T}}|{[w(0{)}^{\mathsf{T}},\dots ,w(k-1{)}^{\mathsf{T}}]}^{\mathsf{T}}\right).$$

$$W$$

$$X$$

$$W_k$$

$$X_k$$

$$z_{i+1}$$

$$e_{i+1}$$

$$z_0(k)$$

$$x_0(k)$$

$$u(k)=K{e}_0(k)+v_0^{\ast }(k),$$

$$e(k+1)=(A+BK)e(k)+w(k),$$

$${\mathcal{R}}_k^x\text{ with }\mathrm{Pr}(e(k)\in {\mathcal{R}}_k^x|e(0)=0)\ge p_x,$$

$${\mathcal{R}}_k^u\text{ with }\mathrm{Pr}(Ke(k)\in {\mathcal{R}}_k^u|Ke(0)=0)\ge p_u$$

$$z_i$$

$$v_i$$

$$\mathrm{Pr}(e_i(k)\in {\mathcal{R}}_{i+k}|x(0))\ge p,$$

$$w(k)\text{ i.i.d.},E(w(k))=0,\mathrm{var}(w(k))={\Sigma }^w,$$

$$l_k(x,u)$$

$$l_f(x)$$

Full text

I Introduction

Most real world control applications are subject to uncertainty and external disturbances, which can severely deteriorate both performance and safety of the system. In model predictive control (MPC) this problem can be addressed by explicitly assessing worst-case disturbances, leading to robust MPC approaches [1]. These approaches, however, can be overly conservative, for example in cases where occasional constraint violations are permissible. Using additional information about the disturbances in form of distributions, stochastic MPC offers advantages through a less conservative treatment of constraints, as well as by improving average performance e.g. by optimizing the expected cost [2].

Stochastic MPC approaches can typically be divided into two classes. While randomized methods rely on the generation of suitable disturbance realizations or scenarios, analytic approximation methods reformulate the problem into a deterministic one [3]. In this paper, we consider an analytic approximation for linear discrete-time systems subject to additive disturbances which are potentially correlated in time and have unbounded support.

Introducing feedback from state measurements into the MPC in a stochastic setting leads to the question of recursive feasibility of the underlying optimization problem. Typically, this issue is addressed in one of two ways [4]. Either it is ensured using a robust constraint tightening, which is generally restricted to the case of bounded support disturbance distributions, see e.g. [5, 6, 7], or the original MPC problem is allowed to become infeasible and a suitable recovery mechanism is employed, e.g. in [8, 9, 10]. These recovery mechanisms usually come with a loss of strict guarantees on the closed-loop chance constraint satisfaction [4]. In [11], strict guarantees were recovered under a unimodality assumption on the additive noise distribution. Additionally, a recent approach guarantees recursive feasibility also with unbounded disturbances for the case of suitably discounted violation probabilities [12].

The approach presented in this paper offers strong guarantees w.r.t. closed-loop chance constraint satisfaction and guarantees recursive feasibility, even for unbounded disturbance distributions, while incorporating feedback from the measured state in the MPC problem. To the best of our knowledge, these properties have not been established in previous approaches. The results are achieved by introducing feedback of the currently measured state $x(k)$ through the cost function only, while tightened constraints are satisfied w.r.t. a nominal system state $z(k)$, ensuring feasibility. This results in the fact that the system error $e(k)=x(k)-z(k)$ evolves linearly in closed-loop, facilitating straightforward analysis of performance and chance constraint satisfaction. Related concepts and their implications for stochastic MPC have recently been discussed in [13]. We demonstrate that this form of feedback has strong influence also on the nominal state trajectory $z(k)$ and results in closed-loop performance comparable to previous stochastic MPC methods—while requiring significantly fewer assumptions on the disturbance distribution for chance constraint satisfaction. For constraint tightening, we employ techniques related to tube-based MPC [14], making use of the concept of probabilistic reachable sets (PRS) [11, 15, 16]. These PRS take a similar role as robust invariant sets in robust MPC such that they contain the error dynamics $e(k)$ with a specified probability. In this paper, we extend this concept to non-zero mean disturbance sequences correlated in time, enabling the treatment of a broad class of problems. We finally demonstrate the flexibility resulting from these properties in a building control task, for which the disturbance sequence is non-i.i.d., non-zero mean, and strongly correlated in time.

II Preliminaries

II-A Notation

We refer to quantities of the system realized in closed-loop at time $k$ using parentheses, e.g. $x(k)$ is the state measured at time step $k$, while quantities used in the MPC prediction are indexed with subscript, e.g. $x_{i}$ is the system state predicted $i$ time steps ahead. In order to specify the time at which the prediction is made, we use $x_{i}(k)$. The weighted 2-norm is $\|x\|_{P}=\sqrt{x^{\mkern-1.5mu\mathsf{T}}Px}$, and $P\succ 0$ ($P\succeq 0$) refers to a positive (semi-)definite matrix. The notation $\mathcal{A}\ominus\mathcal{B}=\left\{a\in\mathcal{A}\,\middle|\,a+b\in\mathcal{A}\ \forall b\in\mathcal{B}\right\}$ refers to the Pontryagin set difference. The distribution $\mathcal{Q}$ of a random variable $x$ is specified as $x\sim\mathcal{Q}$. The probability density of $x$ is denoted $p(x)$, conditioned on another random variable $y$ it is $p(x|y)$. Similarly, probabilities and conditional probabilities are denoted $\Pr(A)$, $\Pr(A\,|\,B)$ and the expected value and variance of $x$ w.r.t. a random variable $w$ are $\operatorname{\mathbb{E}}_{w}(x)$ and $\operatorname{\text{var}}_{w}(x)$, respectively. Two random variables $x$, $y$ that share the same distribution are equal in distribution, denoted $x\overset{d}{=}y$.

II-B Considered System

We consider a linear time-invariant (LTI) system under additive disturbances

[TABLE]

with state $x(k)\in\mathbb{R}^{n_{x}}$, inputs $u(k)\in\mathbb{R}^{n_{u}}$ and randomly distributed disturbance realizations $w(k)$ taking values in $\mathbb{R}^{n_{x}}$. The system is subject to chance constraints on states and inputs

[TABLE]

where $\mathcal{X}$ and $\mathcal{U}$ are convex sets. The probabilities are to be understood w.r.t. knowledge at time step 0, i.e. conditioned on the given initial state. Hard constraints, e.g. on the inputs, can be included in the formulation by imposing a probability of 1. In general, however, these can only be satisfied for disturbance distributions of bounded support.

We consider control problems of arbitrarily large, but finite111We choose a finite control horizon mainly to avoid technicalities. Most properties are easily carried over to an infinite-horizon control problem, e.g. by considering the limit of $\overline{N}\rightarrow\infty$ [17]., horizon $\bar{N}$ and a disturbance sequence $W={[w(0)^{\mkern-1.5mu\mathsf{T}},\ldots,w(\bar{N})^{\mkern-1.5mu\mathsf{T}}]}^{\mkern-1.5mu\mathsf{T}}$ distributed according to distribution $W\sim\mathcal{Q}^{W}$, of which at least the first two moments are known. Note that the individual disturbances $w(k)$ are therefore not necessarily independent, identically distributed (i.i.d.) or zero mean. Using a cost function $l_{k}(x(k),u(k))$ the resulting stochastic (finite-horizon) optimal control problem can be stated as {mini!}[4] {π_k} E_W ( ∑_k=0^¯N l_k(x(k), u(k)) )

\addConstraintx(k​+​1)= A x(k) + B u(k) + w(k) \addConstraintu(k) = π_k(x(0),w(0),…,w(k)) \addConstraintW= [w(0)^T, …, w(¯N)^T]^T∼Q^W \addConstraintPr(x(k)∈X  —  x(0))≥p_x \addConstraintPr(u(k)∈U  —  x(0)) ≥p_u   ,

for all $k=0,\ldots,\bar{N}$, in which $\{\pi_{k}\}$ is a sequence of control laws using information up to time step $k$.

This paper presents a feasible approximate solution to control problem (1) based on receding horizon or model predictive control over a shortened horizon $N\ll\bar{N}$. We will show that the receding horizon controller $\pi^{MPC}$ derived in the following sections satisfies all constraints of optimization problem (1), in particular, it satisfies closed-loop chance constraints (2).

In the following section, we recall probabilistic reachable sets (PRS) and extend the concept to non-i.i.d. disturbance sequences. These PRS form the basis for constraint tightening used in the stochastic model predictive control approach presented in Section IV.

III Probabilistic Reachable Sets

The concept of PRS for stochastic MPC was used in [11] and is related to probabilistic set invariance [18, 15, 19]. In the following, we recall definitions of PRS and consider the extension to non-i.i.d., i.e. time-varying and correlated disturbance sequences. Note that there exists rich literature on related problems, in particular for stochastic reachability of hybrid system [16] and stochastic reach-avoid problems [20, 21].

III-A Definitions

For the following definitions, consider a stochastic process

[TABLE]

where $e(k)$ takes values in $\mathbb{R}^{n_{x}}$.

Definition 1** (Probabilistic $n$-step reachable set).**

A set $\mathcal{R}_{n}$ with $0\leq n\leq\bar{N}$ is an $n$-step probabilistic reachable set ($n$-step PRS) of probability level $p$ for process (3) initialized at $e(0)$ if

[TABLE]

Definition 2** (Probabilistic reachable set).**

A set $\mathcal{R}$ is a probabilistic reachable set (PRS) of probability level $p$ for for process (3) initialized at $e(0)$ if

[TABLE]

Note that this probability bound needs to hold at all time steps individually, as opposed to holding all time steps jointly, which would define much more restrictive sets. It is implicit in the definition of PRS that a PRS is an $n$-step PRS for all $0\leq n\leq\bar{N}$. Conversely, it follows that a set $\mathcal{R}$ is a PRS of level $p$ if it satisfies

[TABLE]

III-B PRS for State

Consider now the case in which the process (3) is defined through linear dynamics under disturbance sequence $W=[w(0)^{\mkern-1.5mu\mathsf{T}},\ldots,w(\overline{N})^{\mkern-1.5mu\mathsf{T}}]\sim\mathcal{Q}$, i.e.

[TABLE]

In [11], it is shown that under the assumption of unimodal i.i.d. disturbances $w(k)$, a convex $n$-step PRS is also an $i$-step PRS for all $i\leq n$. In this case, the relation (4) simplifies to

[TABLE]

In many applications of PRS, e.g. when using PRS for constraint tightening, it is desirable to find sets that are small in a suitable sense. In the following, we present a variance-based computational method for processes of form (5), which in particular for the special case of Gaussian disturbance sequences provides tight sets, in the sense that any down-scaling of the sets would violate the probability guarantees.

III-C Variance-based PRS Computation

We present a straightforward approach for PRS computation based on mean-variance information of $E={[e(1)^{\mkern-1.5mu\mathsf{T}},\ldots,e(\bar{N})^{\mkern-1.5mu\mathsf{T}}]}^{\mkern-1.5mu\mathsf{T}}$ under a correlated disturbance sequence $W={[w(0)^{\mkern-1.5mu\mathsf{T}},\ldots,w(\bar{N})^{\mkern-1.5mu\mathsf{T}}]}^{\mkern-1.5mu\mathsf{T}}$ with $\operatorname{\mathbb{E}}(W)=\mu^{W}$, $\operatorname{\text{var}}(W)=\Sigma^{W}$. Due to the linear dynamics (5) we have for the sequence that

[TABLE]

with

[TABLE]

Using the multivariate Chebyshev inequality we find that

[TABLE]

is an $n$-step PRS of probability level $p=1-n_{x}/\tilde{p}$. The marginal expectation and variance $\operatorname{\mathbb{E}}(e(n))$, $\operatorname{\text{var}}(e(n))$ are directly available from $\operatorname{\mathbb{E}}(E)$ and $\operatorname{\text{var}}(E)$.

Since in the case of correlated disturbance sequences the i.i.d. assumption does not hold, it is not possible to construct PRS along Definition 2 using (6). Instead, one way to construct an ellipsoidal PRS $\mathcal{R}^{c}$ is to find a minimum size ellipsoid which contains all ellipses $\mathcal{R}_{n}^{c}$, thereby satisfying (4). This can be formulated as a semidefinite program, see e.g. [22].

Remark 1* (Zero mean i.i.d. disturbances).*

When $W$ is a zero mean i.i.d. disturbance sequence with $\operatorname{\text{var}}(w(k))=\Sigma^{w}$, a PRS $\mathcal{R}^{c}$ can be constructed by solving the Lyapunov equation $\Sigma_{\infty}=A_{K}\Sigma_{\infty}A_{K}^{\mkern-1.5mu\mathsf{T}}+\Sigma^{w}$ as $\mathcal{R}^{c}:=\left\{x\,\middle|\,x^{\mkern-1.5mu\mathsf{T}}\Sigma_{\infty}^{-1}x\leq\tilde{p}\right\}$ for arbitrary horizons $\bar{N}$.

Remark 2* (Gaussian distributions).*

If $W$ is normally distributed, $\tilde{p}$ can be set to $\tilde{p}=\chi^{2}_{n_{x}}(p)$, where $\chi^{2}_{n_{x}}(p)$ is the quantile function of the chi-squared distribution with $n_{x}$ degrees of freedom, resulting in significantly smaller sets.

III-D PRS for Input

The methods for PRS computations outlined in the previous section are suitable for bounding a state error evolving along (5), and are ultimately used to guarantee chance constraint satisfaction on the states (2a). For the satisfaction of input chance constraints (2b) we need to similarly bound input errors, which we assume to be given by a feedback policy on the state errors $e_{u}(k)=\pi(e(k))$. A possibility of computing ($n$-step) PRS for $e_{u}$ is to transform a PRS on $e(k)$, i.e.

[TABLE]

where $\pi\!\left(\mathcal{R}_{n}\right)=\left\{\pi(e)\,\middle|\,e\in\mathcal{R}_{n}\right\}$. This, however, can lead to significant conservatism, especially if $\pi$ maps to a lower dimensional space, i.e. when $n_{u}<n_{x}$, as is often the case.

A less conservative approach is therefore to first construct the distribution (or mean and variance information) of $e_{u}$ and to find PRS based on this information directly. In the case of linear feedback laws $e_{u}(k)=Ke(k)$ this is often straightforward and PRS construction can be similarly carried out along (7) as well as Remark 1 and 2 by considering $\operatorname{\mathbb{E}}(e_{u}(n))=Ke(n)$ and $\operatorname{\text{var}}(e_{u}(n))=K\operatorname{\text{var}}(e(n))K^{\mkern-1.5mu\mathsf{T}}$.

IV Stochastic MPC using Probabilistic Reachable Sets

The goal is to find an approximate solution to the optimal control problem in (1) by solving it over a shortened horizon $N$ in a receding horizon fashion, i.e. as an MPC controller. To highlight the difference between predicted quantities and the closed-loop quantities in (1) we make use of the subscript $i$, resulting in

[TABLE]

where the predictions dynamics are initialized at the currently measured state at each time step, i.e. $x_{0}(k)=x(k)$. The predicted disturbance sequence $W_{k}=[w_{0}^{\mkern-1.5mu\mathsf{T}},\ldots,w_{N}^{\mkern-1.5mu\mathsf{T}}]^{\mkern-1.5mu\mathsf{T}}$ includes the information about $W$ up to time step $k$. Assuming access to past disturbance realizations it is therefore distributed according to $\mathcal{Q}^{W_{k}}$ defined by the conditional distribution

[TABLE]

The resulting predicted state sequence $X_{k}=[x_{0},\ldots,x_{N}]$ is therefore similarly a random variable.

The differences between closed-loop distributions of

[TABLE]

and predictive distributions of

[TABLE]

are illustrated in Figure 1. The upper plot displays different samples of the disturbance sequence $W$ as well as a confidence region for the disturbance, such that realizations lie within the shaded red region with a certain probability. Additionally, the predictive distribution $W_{k}$ at time step $k$, including confidence region and samples, is shown for a particular realization of $W$. It is evident that for highly correlated sequences, the predicted distribution strongly depends on past disturbance realizations.

The lower plot illustrates the resulting state distribution $X$ from a dynamic system under disturbance $W$ and receding horizon control, including the predictive distribution $X_{k}$ for a specific realization. Note that the distribution of $X$ is typically not available, since it results from a receding horizon control law. Under a suitable restriction of predictive control policies $\pi_{i}$, however, the predictive distribution $X_{k}$ is typically available, since the predicted system is subject to simple, often linear, predictive dynamics. Of particular interest is the relation to the stated chance constraints (2), which are indicated by the dashed line. While the closed-loop dynamics satisfy the displayed chance constraints with high probability, this is not necessarily the case for the predicted distributions, as evident from the example distribution of $X_{k}$. Constraining the predictive constraint violation probabilities can therefore lead to an overly conservative controller, and often feasibility issues, especially when considering unbounded disturbances.

With these consideration we formulate a shortened horizon optimization problem with chance constraints conditioned on the initial state $x(0)$, i.e. {mini!}[4] {π_i} E_W_k ( l_f(x_N) + ∑_i=0^N-1 l_k+i(x_i,u_i) )

\addConstraintx_i+1 = A x_i + B u_i + w_i \addConstraintu_i = π_i(x_0,w_0,…,w_i) \addConstraintW_k= [w_0^T, …, w_N^T]^T∼Q^W_k \addConstraintPr(x_i∈X —  x(0) ) ≥p_x \addConstraintPr(u_i∈U —  x(0) ) ≥p_u \addConstraintx_0= x(k)   ,

for all $i=0,\ldots,N\!-\!1$, where we use a terminal cost $l_{f}(x_{N})$ to approximate the remainder of the horizon of the stochastic optimal control problem (1) and to establish asymptotic performance bounds (see Section IV-D).

Remark 3*.*

The formulation of chance constraints (IV), (IV) considers the probability conditioned on the initial state $x(0)$, i.e. it is subject to random variables $w(0),\ldots,w(k\!-\!1)$ as well as $w_{0}(k),\ldots,w_{i}(k)$ for the $i$-th step constraint. This is in contrast to predictive chance constraints conditioned on the currently measured state $x(k)$, which are often considered in stochastic MPC formulations [3, 11] and are only subject to $w_{0}(k),\ldots,w_{i}(k)$.

IV-A Nominal Dynamics and MPC Initialization

In order to reformulate the probabilistic cost and constraints, we split dynamics (8) into a nominal part $z_{i}$ and error $e_{i}=x_{i}-z_{i}$. We furthermore restrict the class of control policies in the prediction to affine feedback laws acting on the error with a predefined gain $u_{i}=Ke_{i}+v_{i}$. This results in predictive dynamics in nominal and error state given by

[TABLE]

We couple these predicted dynamics to the closed-loop system by defining

[TABLE]

where we initialize the nominal state at time [math] to $z(0)=x(0)$. The nominal predicted state $z_{0}(k)$ is hence updated to its value predicted in the previous time step $z_{1}(k\!-\!1)$, while $x_{0}(k)$ is updated to the measured state $x(k)$, introducing feedback into the optimization problem, as further discussed in Section IV-B and demonstrated in simulation in Section V-A. Since we assume a fixed gain $K$, the optimization is then carried out over $\{v_{i}\}$ and the resulting input applied to the closed-loop system (1) is

[TABLE]

where $e_{0}(k)$ and $v^{*}_{0}(k)$ are the state error and optimal nominal control input, respectively, at time step $k$.

Defining the nominal state $z(k)=z_{0}(k)$ and the error state $e(k)=e_{0}(k)=x(k)-z(k)$ also for the closed-loop system (1) and assuming feasibility of the optimization (which will be shown in Section IV-C), we notice that due to the choice of $z_{0}(k)=z_{1}(k\!-\!1)$ and control law (11) we have

[TABLE]

meaning that the closed-loop error dynamics remain linear, even under the receding horizon controller, which constitutes a nonlinear control law.

IV-B Stochastic MPC Formulation

In order to satisfy chance constraints (IV), (IV), we propose an analytic approximation using tightened deterministic constraints on the nominal system state $z$ and input $v$. To this end, we make use of suitable PRS of the error system, which are employed similarly to an invariant error set in robust MPC. To accommodate different probability levels for input and state constraints, we make use of two reachable sets of different probability levels for each time step in the control horizon $0\leq k\leq\bar{N}$, i.e.

[TABLE]

and arrive at tightened constraints of

[TABLE]

In order to guarantee recursive feasibility of the optimization problem, we furthermore introduce a terminal set $\mathcal{Z}_{f}$ and terminal tightening PRS $\mathcal{R}_{f}$ with $\mathcal{R}_{f}\supseteq\bigcup_{k=0}^{\bar{N}}\mathcal{R}_{k}^{x}$ and $K\mathcal{R}_{f}\supseteq\bigcup_{k=0}^{\bar{N}}\mathcal{R}_{k}^{u}$, satisfying the following properties:

Assumption 1** (Terminal invariance).**

The terminal set $\mathcal{Z}_{f}\subseteq\mathcal{X}\ominus\mathcal{R}_{f}$ is positively invariant for system (9a) under the control law $v=Kz$, i.e. for all $z\in\mathcal{Z}_{f}$ we have $(A\!+\!BK)z\in\mathcal{Z}_{f}$, and $K\mathcal{Z}_{f}\subseteq\mathcal{U}\ominus K\mathcal{R}_{f}$.

Combining these ingredients, the resulting tractable stochastic MPC optimization problem is defined as follows: {mini!}[4] {v_i} E_W_k ( l_f(x_N) + ∑_i=0^N-1 l_k+i(x_i,u_i))

\addConstraintx_i+1= z_i+1 + e_i+1 \addConstraintz_i+1= Az_i + B v_i \addConstrainte_i+1= (A+BK)e_i + w_i \addConstraintW_k= [w_0^T, …, w_N^T]^T∼Q^W_k \addConstraintz_i∈X ⊖R^x_i+k \addConstraintv_i∈U ⊖R^u_i+k \addConstraintz_N∈Z_f \addConstraintx_0= x(k),   z_0 = z_1(k​-​1)  , e_0 = x_0 - z_0   ,

for all $i=0,\ldots,N\!-\!1$. Different from other robust and stochastic MPC approaches there is no direct feedback from the measured state $x(k)$ on the updated nominal state $z_{0}$. Note that feedback also on the nominal trajectory $z(k)$ is nevertheless introduced via the cost in optimization problem (1). We will show in the following that, due to the linear evolution of the closed-loop error (12), this indirect form of feedback offers benefits for the theoretical properties in terms of recursive feasibility, closed-loop chance constraint satisfaction and stability properties, as discussed in the following sections.

IV-C Recursive Feasibility and Chance Constraint Satisfaction

Since the stochastic variables in optimization problem (1) only affect the cost, recursive feasibility can be established in terms of the nominal state $z_{i}$ and input $v_{i}$. Due to the considered update law in (10) this follows standard arguments in predictive control.

Theorem 1** (Recursive feasibility).**

Consider system (1) under the control law (11) resulting from (1). If optimization problem (1) is feasible for $x(0)=z(0)$, then it is recursively feasible, i.e. it is feasible for all times $0\leq k\leq\bar{N}\!-\!N$.

Proof.

Let $V=\{v_{0}(k),\ldots,v_{N-1}(k)\}$ be the optimal solution of optimization problem (1) at time step $k$ with $Z=\{z_{0}(k),\ldots,z_{N}(k)\}$ the resulting nominal state trajectory, satisfying the terminal constraint (1) as well as constraints (1) and (1). We want to find a candidate solution $\bar{V}=\{\bar{v}_{0}(k\!+\!1),\ldots,\bar{v}_{N-1}(k\!+\!1)\}$ which similarly satisfies the terminal constraint (1) and constraints (1), (1) with $\mathcal{R}^{x}_{i+k+1}$ and $\mathcal{R}^{u}_{i+k+1}$, i.e. for the next time step $k\!+\!1$. We choose this candidate solution by shifting $V$ and applying the linear control gain $K$ in the final time step, i.e. $\bar{V}=\{v_{1}(k),\ldots,v_{N-1}(k),Kz_{N}(k)\}$. The first $N-1$ entries evidently fulfill input constraints (1) again, since $\bar{v}_{i}(k\!+\!1)=v_{i+1}(k)\in\mathcal{U}\ominus\mathcal{R}^{u}_{i+k+1}$. Due to Assumption 1 we furthermore have that the final entry $\bar{v}_{N-1}(k\!+\!1)=Kz_{N}(k)\in\mathcal{U}\ominus K\mathcal{R}_{f}$ and since $K\mathcal{R}_{f}\supseteq\mathcal{R}^{u}_{j}$ for all $j=0,\ldots,\bar{N}$ it similarly satisfies the constraint $\bar{v}_{N}(k\!+\!1)\in\mathcal{U}\ominus\mathcal{R}^{u}_{N+k+1}$ for $k\!+\!1\leq N\!-\!\bar{N}$. Since $z_{0}(k\!+\!1)=z_{1}(k)$ the resulting candidate trajectory for the nominal state is $\bar{Z}=\{z_{1}(k),\ldots,z_{N}(k),(A+BK)z_{N}(k)\}$. Due to an analogue argument this satisfies constraints (1), and $(A+BK)z_{N}\in\mathcal{Z}_{f}$ due to Assumption 1. ∎

We can furthermore establish that optimization problem (1) results in a feasible solution to (IV) by noting the following.

Lemma 1**.**

Consider system (1) under the control law (11) resulting from (1). Conditioned on $x(0)$, the predicted error has the same distribution as the closed-loop error, i.e. $e_{i}(k)\overset{d}{=}e(k\!+\!i)$ for $0\leq i\leq N$.

Proof.

See Appendix. ∎

From Lemma 1 it follows directly that if $\mathcal{R}_{i+k}$ is an $i+k$ step PRS of level $p$ for the error system (12), then

[TABLE]

resulting in feasibility in (IV), i.e. $\{\pi_{i}^{*}=Ke_{i}+v_{i}^{*}\}$ obtained from (1) is a feasible solution in (IV). This directly relates to the fact that recursive feasibility due to Theorem 1 implies satisfaction of closed-loop chance constraints (2).

Theorem 2** (Chance constraint satisfaction).**

Consider system (1) under the control law (11) resulting from (1). The resulting states $x(k)$ and inputs $u(k)$ satisfy the closed-loop chance constraints (2).

Proof.

Due to the linear evolution of the closed-loop error (12) and by definition of $k$-step PRS $\mathcal{R}^{x}_{k}$ we have that $\Pr(e(k)\in\mathcal{R}^{x}_{k}\,|\,x(0))\geq p_{x}$ for all $0\leq k\leq\bar{N}$. Due to feasibility of (1) we furthermore have $z(k)=z_{0}(k)\in\mathcal{X}\ominus R^{x}_{k}$. With $x(k)=z(k)+e(k)$ it therefore holds $\Pr(x(k)\in\mathcal{X}\,|\,x(0))\geq p_{x}$. The same argument holds for the input constraints. ∎

We have therefore shown, that an MPC controller based on formulation (1) satisfies closed-loop chance constraints and is therefore a feasible solution to the stochastic optimal control problem (1).

Remark 4*.*

In (1), feedback on the nominal system is introduced through the cost only, while constraints are satisfied w.r.t. a nominal state. A straightforward extension is to consider a soft-constraint term in the cost and thereby introduce feedback from measurements also w.r.t. the constraints, which can be beneficial e.g. in cases of model mismatch while maintaining theoretical guarantees, as demonstrated in examples in Section V-A.

The presented MPC optimization problem requires the evaluation of an expected cost (IV-B). Depending on the specific form, evaluation of the expected value can be computationally expensive. For some cost functions, e.g. linear or quadratic costs, however, this can be done cheaply based on the moments of the predicted state $x_{i}$, which will be outlined in the following section, together with a resulting average asymptotic performance bound.

IV-D Quadratic Costs and Asymptotic Average Bound

Consider the special case of regularization with zero mean i.i.d. disturbances, that is

[TABLE]

for all $0\leq k\leq\bar{N}$ and a quadratic cost function

[TABLE]

with $Q\succeq 0$, $R\succ 0$. For the terminal cost we assume the following to hold.

Assumption 2** (Terminal cost weight).**

The terminal weight $P$ is chosen as the solution to the Lyapunov equation

[TABLE]

If $K$ is stabilizing, this solution always exists.

For the quadratic cost function, evaluation of the expected value can be carried out in terms of mean and variance of $x_{i}$, $u_{i}$, denoted $\mu_{i}^{x}$, $\mu^{u}_{i}$ and $\Sigma^{x}_{i}$, $\Sigma^{u}_{i}$, respectively, as

[TABLE]

Due to the linear prediction dynamics, the evolution of state and input mean and variance can be readily expressed as

[TABLE]

Note that variances $\Sigma^{x}_{i}$ and $\Sigma^{u}_{i}$ are not affected by $v$ and can therefore be neglected in an implementation of problem (1) with cost function (16). After computation of the constraint tightening, the problem is therefore of similar complexity as nominal MPC.

For cost function (16) under zero mean i.i.d. noise, we can establish an asymptotic average performance bound, based on a cost decrease in expectation.

Theorem 3** (Cost decrease).**

Consider system (1) subject to i.i.d. disturbances (15) under the control law (11) resulting from (1) with cost function (16). Let $J^{*}(x,z)$ be the optimal cost of (1), then

[TABLE]

Proof.

See Appendix. ∎

Using the cost decrease in Theorem 3 a standard argument leads to the following asymptotic cost bound.

Corollary 1** (Average asymptotic cost bound).**

Consider system (1) subject to i.i.d. disturbances (15) under the control law (11) resulting from (1) with cost functions (16). We have

[TABLE]

Proof.

See Appendix and [23, 7, 11] for similar derivations. ∎

Remark 5*.*

Note that the SMPC formulation is therefore guaranteed to provide better or equal asymptotic average cost as under linear feedback gain $K$.

V Simulation Examples

We present two simulation studies to demonstrate the proposed approach. The first illustrative example demonstrates feedback properties in a simple regulation task, while the second demonstrates the approach for time-varying correlated disturbance sequences in a building control setting.

V-A Double integrator

To demonstrate the feedback properties on the nominal system trajectory through the state update scheme (10), we consider a simple regulation problem for a double integrator under i.i.d. noise and quadratic cost with $Q=I$, $R=1$, i.e. the setup considered in Section IV-D, specifically

[TABLE]

The system is subject to constraints on the absolute velocity

[TABLE]

and we use the tube controller $K=[-0.2,-0.6]^{\mkern-1.5mu\mathsf{T}}$. For simplicity, we do not consider constraints on the inputs and make use of a constant tightening based on $\Sigma_{\infty}$ along Remark 1 and 2. Considering only the velocity state we get for all $k$

[TABLE]

We choose terminal set $\mathcal{Z}_{f}=\{[0,0]^{\mkern-1.5mu\mathsf{T}}\}$ and consider a prediction horizon of $N=30$. We compare the approach, which we call PRS-rec, to the following two variants.

PRS-nom:

Cost (IV-B) is computed for the nominal states $z_{i}$, resulting in no feedback from $x(k)$ on the nominal system trajectory.

PRS-df:

Approach with direct feedback, in which $z_{0}(k)=x(k)$ whenever feasible and $z_{0}(k)=z_{1}(k)$ otherwise, as previously presented in [11].

V-A1 Nominal Results

We simulate the system $N_{s}=1000$ times starting from initial condition $x(0)=[10,0]^{\mkern-1.5mu\mathsf{T}}$ over a horizon of $\bar{N}=100$ time steps. We compare the resulting average closed-loop cost

[TABLE]

again averaged over all trials, as well as the maximum empirical constraint violation rate in a time step $\max\bar{n}_{v}(k)$ with

[TABLE]

Additionally we report the closed-loop cost $\bar{J}_{cl}(x(20))$ starting at time step $20$ to evaluate the long term behavior without initial transient regularization. The results are given in Table I. All approaches have guaranteed closed-loop chance constraint satisfaction and it is evident that also empirically all constraints are (conservatively) fulfilled. By comparing PRS-rec and PRS-df to PRS-nom it is evident that feedback on the nominal trajectory can significantly reduce the closed-loop cost, particularly $J_{cl}(x(20))$, which Corollary 1 asymptotically guarantees to be smaller than $\operatorname{\text{tr}}(P\Sigma^{w})=7.12$, i.e. the cost under only the linear control law $K$. The example therefore exemplifies how without feedback on the nominal trajectory (PRS-nom) the controller degenerates to a linear control law over time. The considered novel form of feedback (PRS-rec) achieves similar cost to the direct feedback form (PRS-df), which, however, requires stronger assumptions for strict satisfaction of chance constraints, namely i.i.d. unimodal disturbance distributions and symmetric reachable sets [11].

V-A2 Model Mismatch and Soft Constraints

As mentioned in Remark 4, in the presented approach feedback only acts through the cost function which can have adverse effects on chance constraint satisfaction in cases of model mismatch, as is often the case in practical applications. We investigate this case by reducing the input matrix in simulation to $B_{\text{sim}}=\frac{1}{5}B$, such that the controllers are designed w.r.t. misspecified actuator gain. We repeat the regulation experiments for the aforementioned controllers, as well as a soft constrained variant

PRS-recSC:

Includes an additional cost for constraint violation in the predictive distribution. Specifically, we penalize the predicted expected value of the state $\mu^{x}_{i}$ if it lies outside the tightened constraint set $\mathcal{Z}_{i}$, i.e. such that $\mu^{x}_{i}+s_{i}\in\mathcal{Z}_{i}$. The cost on the slack variable is $c|s_{i}|$ and chosen such that the soft constraint is always satisfied, if feasible [24].

Note that the theoretical results w.r.t. constraint satisfaction for the nominal system model similarly hold for PRS-recSC.

The results of this simulation are summarized in Table II. With respect to the incurred cost, the results from this simulation are qualitatively similar to the nominal case without modeling error. PRS-nom incurs the highest cost, since no feedback acts on the nominal trajectory and the controller degenerates to the linear feedback law $Kx$. Additionally it is unable to achieve the specified level of constraint satisfaction due to the large model mismatch. The recursively feasible approach without soft constraints PRS-rec displays the lowest achieved cost, but does incur similar constraint violations as PRS-nom. For the direct feedback controller PRS-df this is not the case, since the direct feedback also acts w.r.t. to the constraints. This does, however, come at significantly higher cost. Finally, the soft constrained recursively feasible controller PRS-recSC results in similar chance constraint satisfaction, as well as good performance in terms of cost, even in this case of significant model mismatch.

In the following example we demonstrate the flexibility of the presented method by considering a building control problem, subject to non i.i.d. disturbance distributions, for which the strong guarantees in terms of closed-loop chance constraint satisfaction similarly hold.

V-B Building Control

Consider the simple building temperature control example shown in Figure 2, consisting of four rooms with individual temperatures (states) $x(k)=[T^{1}(k),T^{2}(k),T^{3}(k),T^{4}(k)]^{\top}$, combined heating/cooling units for each room (inputs) $u(k)=[u^{1}(k),u^{2}(k),u^{3}(k),u^{4}(k)]^{\top}$, and uncertain outside temperature (disturbances), represented by $w(k)$. The system dynamics are modeled using a resistance network (see e.g. [25]), in which each room is characterized by a thermal capacity, while the interactions are governed by the thermal conductance, corresponding to the resistances.

The resulting system description is given by

[TABLE]

in which $A$ captures the thermal conductances between the rooms, $B$ represents the thermal inertia with respect to heating/cooling, $B_{d}$ is the effect of the outside temperature on each room, $\bar{T}^{o}$ is given by the overall average outside temperature, which is included into the nominal system for an intuitive presentation of the results222Note that $\bar{T}^{o}$ can be similarly included in the disturbance sequence $W$, yielding equivalent results. A zero mean representation is used for increased interpretability., and $w(k)=B_{d}(T^{o}(k)-\bar{T}^{o})$ is the uncertain deviation from the mean outside temperature, see Appendix for parameters and matrices. We consider a normally distributed disturbance sequence $W=[w(0),\ldots,w(\bar{N\!}-\!1)]\sim\mathcal{(}\mu^{W},\Sigma^{W})$ which is highly correlated in time, i.e. the current outside temperature significantly influences the temperature in the following hours. For illustration, Figure 3 (top) displays samples from the conditional distributions $W_{k}$ used in prediction.

We consider regulation with respect to $T^{i}_{r}=$21.75\text{,}\mathrm{\SIUnitSymbolCelsius}$$ in each room, and a quadratic cost on the temperature deviation and a cost on the absolute value of each input, corresponding to an economic cost for heating or cooling, i.e. $l(x,u)=\|x-T_{r}\|^{2}_{Q}+\|u\|_{1}$ with $Q=550I$. We choose the tube controller feedback $K$ using an LQR design with weight matrices $Q_{LQR}=10^{5}I,~{}R_{LQR}=0.03I$ and consider for simplicity the desired set point as terminal set for each state, i.e. $\mathcal{Z}_{f}=\{21.75\!\cdot\!\mathbf{1}\}$, where $\mathbf{1}$ denotes the one vector. Comfort constraints on the room temperature and physical input constraints on the heating/cooling power are given by

[TABLE]

Simulating the system with initial condition $x(0)=22.5\!\cdot\!\mathbf{1}$ yields the closed-loop behavior for room 1 shown in Figure 3. Because of the economic cost on the heating power, the applied input $u^{1}(k)$ is close to zero for the first 10 hours until chilly outside temperatures during night cool the rooms below the desired set point. At $k=30$ hours, the reduced room temperature before nighttime requires an increased heating effort during the predicted future time steps.

Note the time-varying constraint tightening on the nominal system state $z$, which results from considering time-varying correlated disturbance sequence. This helps to significantly reduce conservatism of the MPC controller compared to time-invariant disturbance bounds. The approach is able to regulate the room temperature close to the reference point while satisfying all comfort and input constraints. The example therefore illustrates the potential applicability of the presented approach to a number of interesting engineering applications, while providing theoretical guarantees on chance constraint satisfaction.

VI Conclusion

We presented a stochastic model predictive control approach for additive correlated disturbance sequences. Using the concept of probabilistic reachable sets for constraint tightening, as well as a novel initialization scheme of the MPC, we were able to show chance constraint satisfaction for the closed-loop system. We demonstrated in simulation that this novel update scheme achieves similar or better performance to established forms of feedback in MPC, while facilitating theoretical analysis. Finally, we demonstrated the presented framework for a building control task, highlighting the applicability of the approach to a broad class of interesting engineering applications.

Acknowledgments

The authors would like to thank A. Bollinger, C. Gaehler, J. Lygeros, and R. Smith for providing the building control simulation model.

Appendix

Proof of Lemma 1.

Since $z(0)=x(0)$ we have $e(0)=0$ and due to (12) that $e(k)=\sum_{j=0}^{k-1}(A+BK)^{k-1-j}w(j)$. For the predicted errors we similarly have with $e_{0}(k)=e(k)$ that

[TABLE]

Since $[w(0)^{\mkern-1.5mu\mathsf{T}},\ldots,w(k\!+\!N)^{\mkern-1.5mu\mathsf{T}}]^{\mkern-1.5mu\mathsf{T}}$ has the same distribution as $[w(0)^{\mkern-1.5mu\mathsf{T}},\ldots,w(k-1)^{\mkern-1.5mu\mathsf{T}},w_{0}(k)^{\mkern-1.5mu\mathsf{T}},\ldots,w_{N}(k)^{\mkern-1.5mu\mathsf{T}}]^{\mkern-1.5mu\mathsf{T}}$ this is equal in distribution to

[TABLE]

∎

Proof of Theorem 3.

For notational convenience we will omit the conditioning on $x(k)$, $z(k)$. Let $J(x,z,V)$ be the cost of optimization problem (1) under input sequence $V$. We have

[TABLE]

where $\bar{V}=\{\bar{v}_{0},\ldots,\bar{v}_{N-1}\}=\{v_{1},\ldots,v_{N-1},Kz_{N}\}$ is a candidate solution with resulting nominal state trajectory $\bar{Z}=\{\bar{z}_{0},\ldots,\bar{z}_{N}\}=\{z_{1},\ldots,z_{N},(A+BK)z_{N}\}$, see also proof of Theorem 1. Note that since $w(k)$ has the same distribution as $w_{0}$, we have for the resulting initial error state

[TABLE]

such that

[TABLE]

The candidate state sequence $\bar{x}_{i}=\bar{z}_{i}+\bar{e}_{i}$ and inputs $\bar{u}_{i}=K\bar{e}_{i}+\bar{v}_{i}$ therefore result in

[TABLE]

where the last element of $\bar{U}$ follows from $\bar{u}_{N-1}=K\bar{e}_{N-1}-\bar{v}_{N-1}\overset{d}{=}Ke_{N}+Kz_{N}=Kx_{N}$.

We therefore have

[TABLE]

where the last equation follows from the fact that $x_{N}$ and $w_{N}$ are independent and $w_{N}$ is zero mean. Due to Assumption 2 this further simplifies to

[TABLE]

∎

Proof of Corollary 1.

Let $\Delta J^{*}(k)=J^{*}(x(k\!+\!1),z(k\!+\!1))-J^{*}(x(k),z(k))$. By the law of iterated expectations, we have

[TABLE]

since $\Delta J^{*}(k)$ is independent of previous $x$ and $z$ given $x(k),z(k)$. Using the cost decrease of Theorem 3 this means

[TABLE]

From this it follows that the asymptotic average is bounded from below by zero (since $J(x(0),z(0))$ is finite) and from above by

[TABLE]

and the claim follows. ∎

Building control example

The dynamics are discretized using Euler forward with step size $\Delta t=3600~{}[s]$ which yields

[TABLE]

with $C,D,H\in\mathbb{R}^{n\times n}$, where

[TABLE]

and $H_{ij}$ $[W/K]$, $h\in\mathbb{R}^{n}$ $[W/K]$ are given by

[TABLE]

and $\operatorname{\text{diag}}(C)=10^{6}\cdot[50,110,80,90]~{}[J/K]$.