Shrinking Horizon Model Predictive Control with Signal Temporal Logic Constraints under Stochastic Disturbances

TL;DR

This paper introduces a novel Shrinking Horizon Model Predictive Control method for discrete-time linear systems with STL constraints, effectively handling stochastic disturbances without requiring full distribution knowledge, demonstrated on HVAC systems.

Contribution

It develops a general STL-constrained control approach under stochastic disturbances that does not need precise distribution information, extending to cases with known distributions for improved performance.

Findings

Effective control synthesis for HVAC systems.

Robust STL satisfaction under stochastic disturbances.

Optimization problems with linear constraints at each step.

Abstract

We present Shrinking Horizon Model Predictive Control (SHMPC) for discrete-time linear systems with Signal Temporal Logic (STL) specification constraints under stochastic disturbances. The control objective is to maximize an optimization function under the restriction that a given STL specification is satisfied with high probability against stochastic uncertainties. We formulate a general solution, which does not require precise knowledge of the probability distributions of the (possibly dependent) stochastic disturbances; only the bounded support intervals of the density functions and moment intervals are used. For the specific case of disturbances that are independent and normally distributed, we optimize the controllers further by utilizing knowledge of the disturbance probability distributions. We show that in both cases, the control law can be obtained by solving optimization problems with linear constraints at each step. We experimentally demonstrate effectiveness of this approach by synthesizing a controller for an HVAC system.

Tables

Table 1. TABLE I: Comparison of the statistics of the fan energy consumption using different control approaches.

Computational	Fan energy	Average
Methods	consumption [kWh]	computation time [s]
Open-loop OC	$1337.016$	3.9277
RMPC	${\mu }_1=12.2216$, ${\sigma }_1=0.045{\mu }_1$	33.4891
SHMPC	${\mu }_2=2.5101$, ${\sigma }_2=0.104{\mu }_2$	19.3622

Equations

$$X(t+1)=A(t)X(t)+B(t)u(t)+W(t),X(0)=x_0,$$

$$X(\tau )=\Phi (\tau ,t)X(t)+\sum _{k=t}^{\tau -1}\Phi (\tau ,k+1)\left(B(k)u(k)+W(k)\right),$$

$$\Phi (\tau ,t)=\{\begin{array}{ll}A(\tau -1)A(\tau -2)\dots A(t)& \tau >t\ge 0\\ {\mathbb{I}}_n& \tau =t\ge 0,\end{array}$$

$$\phi ::=\top \mid \pi \mid \neg \phi \mid \phi \wedge \psi \mid \phi {\mathcal{U}}_{[a,b]}\psi$$

$$\text{If }x(t+i)=x^{\prime }(t+i)\text{ for all }i\in \{0,\dots ,n\}$$

$$\text{Then }(\xi ,t)\models \phi \text{ iff }({\xi }^{\prime },t)\models \phi .$$

$$\phi$$

$${\rho }^{\top }(\xi ,t)$$

$${\rho }^{\pi }(\xi ,t)$$

$${\rho }^{\neg \phi }(\xi ,t)$$

$${\rho }^{\phi \wedge \psi }(\xi ,t)$$

$${\rho }^{\phi {\mathcal{U}}_{[a,b]}\psi }(\xi ,t)$$

$$\underset{\tilde{u}(0:N)}{\min }\mathbb{E}\left[J(\tilde{X}(0:N+1),\tilde{u}(0:N))\right]\text{mbox}s.t.$$

$$X(t)=\Phi (t,0)x_0+\sum _{k=0}^{t-1}\Phi (t,k+1)\left(B(k)u(k)+W(k)\right),$$

$$\text{Pr}\left[{\Xi }_N(x_0,\tilde{u}(0:N),\tilde{W}(0:N))\models \phi \right]\ge 1-\delta ,$$

$$\tilde{u}(0:N)\in U^N,$$

$$J(\tilde{X}(0:N+1),\tilde{u}(0:N)):=J_{\text{robust}}(\tilde{X}(0:N+1))+J_{\text{in}}(\tilde{u}(0:N)),$$

$$J(\bar{X}(0:t:$$

$$J_{\text{robust}}(\bar{X}(0:t:N+1))+J_{\text{in}}(\bar{u}(0:t-1:N)),$$

$$\underset{\tilde{u}(t:N)}{\min }\mathbb{E}\left[J(\bar{X}(0:t:N+1),\bar{u}(0:t-1:N))\right]\text{mbox}s.t.$$

$$X(\tau )=\Phi (\tau ,t)x(t)+\sum _{k=t}^{\tau -1}\Phi (\tau ,k+1)\left(B(k)u(k)+W(k)\right),$$

$$\text{for}t\le \tau \le N$$

$$\text{Pr}\left[{\Xi }_N(x_0,\bar{u}(0:t-1:N),\bar{W}(0:t-1:N))\models \phi \right]\ge 1-{\delta }_t$$

$$\tilde{u}(t:N)\in U^{N-t},$$

$$J_{\text{robust}}(\bar{X}(0:t:N))=\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }\left\{{\eta }_{ij}+{\lambda }_{ij}\bar{W}(0:t:N)\right\},$$

$$J_{\text{robust}}(\bar{X}(0:t:N))=\underset{i\in \{1,\dots ,K\}}{\min }\underset{j\in \{1,\dots ,n_i\}}{\max }\left\{{\zeta }_{ij}+{\gamma }_{ij}\bar{W}(0:t:N)\right\},$$

$$\min$$

$$\max \left(\min (f_1,g_1),\min (f_1,g_2),\min (f_2,g_1),\min (f_2,g_2)\right).$$

$$I_{X(\tau )}=[{\bar{a}}_{\tau }+{\bar{C}}_{\tau },{\bar{b}}_{\tau }+{\bar{C}}_{\tau }],{\mathbb{M}}_{X(\tau )}=[{\bar{c}}_{\tau }+{\bar{C}}_{\tau },{\bar{d}}_{\tau }+{\bar{C}}_{\tau }]$$

$$\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }({\hat{d}}_{ij}+{\eta }_{ij}),$$

$$I_{Y_{ij}}=[{\hat{a}}_{ij}+{\eta }_{ij},{\hat{b}}_{ij}+{\eta }_{ij}],{\mathbb{M}}_{Y_{ij}}=[{\hat{c}}_{ij}+{\eta }_{ij},{\hat{d}}_{ij}+{\eta }_{ij}]$$

$$I_{J_{\text{robust}}}=$$

$$[\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }({\hat{a}}_{ij}+{\eta }_{ij}),\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }({\hat{b}}_{ij}+{\eta }_{ij}]$$

$${\mathbb{M}}_{J_{\text{robust}}}=$$

$$[\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }({\hat{c}}_{ij}+{\eta }_{ij}),\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }({\hat{d}}_{ij}+{\eta }_{ij})].$$

$${\mu }_{\tau }$$

$${\Sigma }_{\tau }$$

$$\underset{i\in \{1,\dots ,L\}}{\max }\underset{j\in \{1,\dots ,m_i\}}{\min }{Y}_{ij}\text{ or }\underset{i\in \{1,\dots ,K\}}{\min }\underset{j\in \{1,\dots ,n_i\}}{\max }{Y}_{ij}$$

Full text

I Introduction

We consider the control synthesis problem for stochastic discrete-time linear systems under path constraints that are expressed as temporal logic specifications and are written in signal temporal logic (STL) [23]. Our aim is to obtain a controller that robustly satisfies desired temporal properties with high probability despite stochastic disturbances, while optimizing additional control objectives. With focus on temporal properties defined on a finite path segment, we use the model predictive control (MPC) scheme [3, 8, 20, 22] with a shrinking horizon: the horizon window is fixed and not shifted at each time step of the controller synthesis problem. We start with an initial prediction horizon dependent on the temporal logic constraints, compute the optimal control sequence for the horizon, apply the first step, observe the system evolution under the stochastic disturbance, and repeat the process (decreasing the prediction horizon by 1) until the end of the original time horizon.

Our proposed setting requires solving three technical challenges in the MPC framework.

First, in addition to optimizing control and state cost, the derived controller must ensure that the system evolution satisfies chance constraints arising from the STL specifications. Previous choices of control actions can impose temporal constraints on the rest of the path. We describe an algorithm that updates the temporal constraints based on previous actions.

Second, for some temporal constraints, we may require that the system satisfies the constraints robustly: small changes to the inputs should not invalidate the temporal constraint. To ensure robust satisfaction, we use a quantitative notion of robustness for STL [10]. We augment the control objective to maximize the expected robustness of an STL specification, in addition to minimizing control and state costs, under chance constraints. Unfortunately, the resulting optimization problem is not convex.

As a third contribution, we propose a tractable approximation method for the solution of the optimization problem. We conservatively approximate the chance constraints by linear inequalities. Second, we provide a tractable procedure to compute an upper bound for the expected value of the robustness function under these linear constraints.

Recently receding horizon control with STL constraints has been studied for a variety of domains [12, 24]. In these works, the disturbance is assumed to be deterministic but from a bounded polytope, and the worst-case MPC optimization problem is solved. The control synthesis for deterministic systems with probabilistic STL specifications is studied in [25] but only a fragment of STL is considered in order to obtain a convex optimization problem. Also, the receding horizon control has been applied to deterministic systems in the presence of perception uncertainty [17]. Additionally, chance-constrained MPC has been addressed in [26] for deterministic systems, in which the underlying probability space comes only from the measurement noise. Application of chance-constrained MPC to optimal control of drinking water networks is studied in [14].

In this paper, we assume that the the disturbance signal has an arbitrary probability distribution with bounded domain and that we only know the support and the first moment interval for each component of the disturbance signal. In order to solve the optimization problem more efficiently, we transform chance constraints into their equivalent (or approximate) linear constraints. To this end, we apply the technique presented by [4], to approximate the chance constraints with an upper bound. Also, the expected value of the robustness function can be approximated by the moment intervals of the disturbance signal, and can be computed without using numerical integration.

Furthermore, as an additional case in this study, we show that if the disturbance signal is normally distributed and hence, has no bounded support, instead of truncating the distribution to obtain a finite interval of support for random variables, we can use a different approach, which is based on the quantiles of the normally distributed random variables to replace the chance constraints by linear constraints. In this case, we also show that the expected value of the robustness function can be replaced by an upper bound using the methods presented in [13].

We empirically demonstrate the effectiveness of our approach by synthesizing a controller for a Heating, Ventilation and Air Conditioning (HVAC) system. We compare our approach with open-loop optimal controller synthesis and with robust MPC [24], and show that our approach can lead to significant energy savings.

I-A Notations

We use $\mathbb{R}$ for the set of reals and $\mathbb{N}:=\{0,1,2,\ldots\}$ for the set of non-negative integers. The set $\mathbb{B}:=\{\top,\bot\}$ indicates logical true and false. For a vector $v\in\mathbb{R}^{s}$, its components are denoted by $v_{k}$, $k\in\{1,\dots,s\}$. For a sequence $\{v(t)\in\mathbb{R}^{s},\,\,t\in\mathbb{N}\}$ and $t_{1}<t_{2}$, we define $\tilde{v}(t_{1}:t_{2}):=[v(t_{1}),v(t_{1}+1),\ldots,v(t_{2}-1)]$, In this paper, all random variables are denoted by capital letters and the deterministic variables are denoted by small letters. We also use small letter $y$ to indicate observations of a random vector $Y$. For a sequence of random vectors $\{Y(t)\in\mathbb{R}^{s},\,\,t\in\mathbb{N}\}$ and $t_{1}\leq t<t_{2}$, we define $\bar{Y}(t_{1}:t:t_{2}):=[y(t_{1}),y(t_{1}+1),\ldots,y(t),Y(t+1),\ldots,Y(t_{2}-1)]$, which is a matrix containing observations of the random variable up to time $t$ augmented with its unobserved values after $t$. For a random variable $Y(t)$ denote the support interval by $I_{Y(t)}$ and the first moment111The expected value of a random variable $X$ with support $\mathcal{D}$ and the cumulative distribution function $F$ is defined as $\mathbb{E}[X]=\int_{\mathcal{D}}xdF(x)$. The expectation exists if the integral is well-defined and yields a finite value. by $\mathbb{E}[Y(t)]$.

We consider operations on intervals according to interval arithmetic: for two arbitrary intervals $[a,b]$ and $[c,d]$, and constants $\lambda,\gamma\in\mathbb{R}$, we have $[a,b]+[c,d]=[a+c,b+d]$ and $\lambda\cdot[a,b]+\gamma=[\min(\lambda a,\lambda b)+\gamma,\max(\lambda a,\lambda b)+\gamma]$.

II Discrete-Time Stochastic Linear Systems

In this paper, we consider systems in discrete-time that can be modeled by linear difference equations perturbed by stochastic disturbances. Depending on the probability distribution of the disturbance signal, we conduct our study for two cases: a) the disturbance signal has an arbitrary probability distribution with a bounded domain for which we only know the support and their first moment intervals; and b) the disturbance signal has a normal distribution. The first case can be extended to random variables with an unbounded support, such as normal or exponential random variables, by truncating their distributions. The specific form of the distribution in the second case enables us to perform a more precise analysis using properties of the normal distribution. Note that the support of a random variable $X$ with values in $\mathbb{R}^{n}$ is defined as the set $\{x\in\mathbb{R}^{N}\,|\,\text{Pr}_{X}[\mathcal{B}(x,r)]>0,\,\forall r>0\}$, where $\mathcal{B}(x,r)$ denotes the ball with center at $x$ and radius $r$; alternatively, the support can be defined as the smallest closed set $C$ such that $\text{Pr}_{X}[C]=1$.

Consider a (time-variant) discrete-time stochastic system modeled by the difference equation

[TABLE]

where $X(t)\in\mathbb{R}^{n}$ denotes the state of the system at time instant $t$, $u(t)\in\mathbb{R}^{m}$ denotes the control input at time instant $t$, and $W(t)\in\mathbb{R}^{s}$ is a vector of random variables, the components of which have either of the above mentioned probability distributions. The random vector $W(t)$ can be interpreted as the process noise or an adversarial disturbance. Matrices $A(\cdot)\in\mathbb{R}^{n\times n}$, and $B(\cdot)\in\mathbb{R}^{n\times m}$ are possibly time-dependent appropriately defined system’s matrices, and the initial state $X(0)$ is assumed to be known. We assume that $W(0),\ldots,W(t)$ are mutually independent random vectors for all time instants $t$. Note that, for any $t\in\mathbb{N}$, the state-space model (1) provides the following explicit form for $X(\tau)$, $\tau\geq t$, as a function of $X(t)$, input $u(\cdot)$, and the process noise $W(\cdot)$:

[TABLE]

where $\Phi(\cdot,\cdot)$ is the state transition matrix of model (1) defined as

[TABLE]

with $\mathbb{I}_{n}$ being the identity matrix.

For a fixed positive integer $N$, and a given time instant $t\in\mathbb{N}$, let $\tilde{u}(t:N)=[u(t),u(t+1),\ldots,u(N-1)]$ (matrix $\tilde{W}(t:N)$ is defined similarly). A finite stochastic run of system (1) for the time interval $[t:N]$ is defined as $\Xi(t:N)=X(t)X(t+1)\ldots X(N)$, which is a finite sequence of states satisfying (2). Since each state $X(\tau)$ depends on $X(t),\tilde{u}(t:N)$, and $\tilde{W}(t:N)$, we can rewrite $\Xi(t:N)$ in a more elaborative notation as $\Xi_{N}(X(t),\tilde{u}(t:N),\tilde{W}(t:N))$. Analogously, we define an infinite stochastic run $\Xi=X(t)X(t+1)X(t+2)\ldots$ as an infinite sequence of states. Stochastic runs will be used in Section III to define the system’s specifications.

III Signal Temporal Logic

An infinite run of system (1) can be considered as a signal $\xi=x(0)x(1)x(2)\dots$, which is a sequence of observed states. We consider Signal temporal logic (STL) formulas with bounded-time temporal operators defined recursively according to the grammar [23]

[TABLE]

where $\top$ is the true predicate; $\pi$ is a predicate whose truth value is determined by the sign of a function, i.e. $\pi=\{\alpha(x)\geq 0\}$ with $\alpha:\mathbb{R}^{n}\rightarrow\mathbb{R}$ being an affine function of state variables; $\psi$ is an STL formula; $\neg$ and $\land$ indicate negation and conjunction of formulas; and ${\mathcal{U}}_{[a,b]}$ is the until operator with $a,b\in\mathbb{R}_{\geq 0}$. A run $\xi$ satisfies $\varphi$ at time $t$, denoted by $(\xi,t)\models\varphi$, if the sequence $x(t)x(t+1)\ldots$ satisfies $\varphi$. Accordingly, $\xi$ satisfies $\varphi$, if $(\xi,0)\models\varphi$.

Semantics of STL formulas are defined as follows. Every run satisfies $\top$. The run $\xi$ satisfies $\neg\varphi$ if it does not satisfy $\varphi$; it satisfies $\varphi\land\psi$ if both $\varphi$ and $\psi$ hold. For a run $\xi=x(0)x(1)x(2)\ldots$ and a predicate $\pi=\{\alpha(x)\geq 0\}$, we have $(\xi,t)\models\pi$ if $\alpha(x(t))\geq 0$. Finally, $(\xi,t)\models\varphi{\mathcal{U}}_{[a,b]}\psi$ if $\varphi$ holds at every time step starting from time $t$ before $\psi$ holds, and additionally $\psi$ holds at some time instant between $a+t$ and $b+t$. Additionally, we derive the other standard operators as follows. Disjunction $\varphi\lor\psi:=\neg(\neg\varphi\land\neg\psi)$, the eventually operator as \operatorname{\rotatebox[origin={c}]{45.0}{\Box}}_{[a,b]}\varphi:=\top{\mathcal{U}}_{[a,b]}\varphi, and the always operator as \operatorname{\Box}_{[a,b]}\varphi:=\neg\operatorname{\rotatebox[origin={c}]{45.0}{\Box}}_{[a,b]}\neg\varphi.

Thus (\xi,t)\models\operatorname{\rotatebox[origin={c}]{45.0}{\Box}}_{[a,b]}\varphi if $\varphi$ holds at some time instant between $a+t$ and $b+t$ and $(\xi,t)\models\operatorname{\Box}_{[a,b]}\varphi$ if $\varphi$ holds at every time instant between $a+t$ and $b+t$.

Formula Horizon. The horizon of an STL formula $\varphi$ is the smallest $n\in\mathbb{N}$ such that the following holds for all signals $\xi=x(0)x(1)x(2)\ldots$ and $\xi^{\prime}=x^{\prime}(0)x^{\prime}(1)x^{\prime}(2)\ldots$:

[TABLE]

Thus, in order to determine whether a signal $\xi$ satisfies an STL formula $\varphi$, we can restrict our attention to the signal prefix $x(0),\ldots,x(\Delta)$ where $\Delta$ is the horizon of $\varphi$. This horizon can be upper-approximated by a bound, denoted by $\text{len}(\varphi)$, defined to be the maximum over the sums of all nested upper bounds on the temporal operators. Formally, $\text{len}(\varphi)$ is defined recursively as:

[TABLE]

where $\varphi_{1},\varphi_{2}$ and $\psi$ are STL formulas. For example, for \varphi=\square_{[0,4]}\operatorname{\rotatebox[origin={c}]{45.0}{\Box}}_{[3,6]}\pi, we have $\text{len}(\varphi)=4+6=10$. For a given STL formula $\varphi$, it is possible to verify that $\xi\models\varphi$ using only the finite run $x(0)x(1)\ldots x(N)$, where $N$ is equal to $\text{len}(\varphi)$.

STL Robustness. In contrast to the above Boolean semantics, the quantitative semantics of STL [18] assigns to each formula $\varphi$ a real-valued function $\rho^{\varphi}$ of signal $\xi$ and $t$ such that $\rho^{\varphi}(\xi,t)>0$ implies $(\xi,t)\models\varphi$. Robustness of a formula $\varphi$ with respect to a run $\xi$ at time $t$ is defined recursively as

[TABLE]

where $x(t)$ refers to signal $\xi$ at time $t$. The robustness of the derived formula \operatorname{\rotatebox[origin={c}]{45.0}{\Box}}_{[a,b]}\varphi can be worked out to be \rho^{\operatorname{\rotatebox[origin={c}]{45.0}{\Box}}_{[a,b]}\varphi}(\xi,t)=\max_{i\in[a,b]}\rho^{\varphi}(\xi,t+i); and similarly for $\operatorname{\Box}_{[a,b]}\varphi$ as $\rho^{\operatorname{\Box}_{[a,b]}\varphi}(\xi,t)=\min_{i\in[a,b]}\rho^{\varphi}(\xi,t+i)$. The robustness of an arbitrary STL formula is computed recursively on the structure of the formula according to the above definition, by propagating the values of the functions associated with each operand using min and max operators.

STL Robustness for Stochastic Runs. With focus on stochastic runs $\Xi=X(0)X(1)X(2)\dots$ and using the bound of a formula $\varphi$, the finite stochastic run $\Xi(t:t+N)=X(t)X(1)\ldots X(t+N)$ with $N=\text{len}(\varphi)$ is sufficient to study probabilistic properties of $(\Xi,t)\models\varphi$. Analogous to the definition of robustness for deterministic run, we can define stochastic robustness $\rho^{\varphi}(\Xi,t)$ of a formula $\varphi$ with respect to the run $\Xi$ for times $t$ with the stochastic robustness being dependent on $\Xi(t:t+N)$ and $\varphi$.

Note that a general STL formula $\varphi$ consists of several Boolean and/or temporal operators. Due to the system dynamics (1), the stochastic run $\Xi(t:t+N)$ and $\rho^{\varphi}(\Xi(t:t+N),t)$ are both functions of $\tilde{W}(t:t+N)$. Therefore, $\rho^{\varphi}(\Xi(t:t+N),t)$ is a random variable since affine operators, maximization and minimization are measurable functions.

The above definition of robustness implies that, for any formula $\varphi$ and constant $\delta\in(0,1)$, a stochastic run $\Xi=X(0)X(1)X(2)\ldots$ satisfies $\varphi$ with probability greater than or equal to $1-\delta$, if the finite stochastic run $\Xi(0:N)=X(0)X(1)\ldots X(N)$ with $N\geq\text{len}(\varphi)$ satisfies $\text{Pr}\left[\rho^{\varphi}(\Xi(0:N),0)>0\right]\geq 1-\delta$.

IV Control Problem Statement

For system (1) with a given initial state $X(0)=x_{0}$, the stochastic disturbance vector $W(t)$ with a given probability distribution, STL formulas $\varphi$ and $\psi$, and some constant $N\geq\max(\text{len}(\varphi),\text{len}(\psi))$, the control problem can be defined as finding an optimal input sequence $\tilde{u}^{\ast}(0:N)=[u^{\ast}(0),\ldots,u^{\ast}(N-1)]$, that minimizes the expected value of a given objective function $J(\tilde{X}(0:N+1),\tilde{u}(0:N))$ subject to constraints on states and input variables, where $\tilde{X}(0:N+1)\!=\![X(0),X(1),\ldots,X(N)]$. This optimization problem for the time interval $0\leq t<N$ can be defined as

[TABLE]

where $\mathbb{E}[\cdot]$ denotes the expected value operator and the closed set $U^{N}\in\mathbb{R}^{mN}$ specifies the constraint set for the input variables. The chance constraints (3c) state that for a given $\delta\in(0,1)$, stochastic runs of the system should satisfy $\varphi$ with a probability greater than or equal to $1-\delta$. We consider the following objective function

[TABLE]

where the first term $J_{\text{robust}}(\tilde{X}(0:N+1)):=-\rho^{\psi}(\tilde{X}(0:N+1),0)$ represents the negative value of the robustness function on STL formula $\psi$ at time [math] that needs to be minimized; and the second term $J_{\text{in}}(\tilde{u}(0:N))$ reflects the cost on the input variables and can be defined as a linear or a quadratic function.

Note that optimization problem (3) is an open-loop optimization problem and we cannot incorporate any information related to the process noise or the states of the system.

Remark 1

The above problem formulation enables us to distinguish the following two cases: we put the robustness of a formula in the objective function if the system is required to be robust with respect to satisfying the formula; we encode the formula in the probabilistic constraint if only satisfaction of the formula is important.

IV-A Model Predictive Control

To obtain a more well-behaved control input and to include the information about the disturbances, instead of solving the optimization problem (3), we apply shrinking horizon model predictive control (SHMPC), which can be summarized as follows: at time step one, we obtain a sequence of control inputs with length $N$ (the prediction horizon) to optimize the cost function; then we only use the first component of the obtained control sequence to update the state of the system (or to observe the state in the case of having a stochastic disturbance); in the next time step, we fix the first component of the control sequence by the first component of the previously calculated optimal control sequence and hence, we only optimize for a control sequence of length $N-1$. As such, at each time step, the size of the control sequence decreases by 1. Note that in this approach, unlike the receding horizon approach, we do not shift the horizon at each time step and the end point of the prediction window is fixed. MPC allows us to incorporate the new information we obtain about the state variables and the disturbance signal, at each time step and hence, to improve the control performance comparing with the one of solving the open-loop optimization problem (3).

A natural choice for the prediction horizon $N$ in this setting with STL constraints is to set it to be greater than or equal to the bound of the formula $\varphi$, i.e., $\text{len}(\varphi)$, which was defined in the previous section. This choice provides a conservative maximum trajectory length required to make a decision about the satisfiability of the formula.

Let $\bar{X}(0:t:N+1)\!=\![x(0),\ldots,x(t),X(t+1),\ldots,X(N)]$ where $x(0),\ldots,x(t)$ are the observed states up to time $t$ and $X(\tau)$ is the random state variable at time $\tau>t$, also let $\bar{W}(0:t-1:N)\!=\![w(0),\ldots,w(t-1),W(t),W(t+1),\ldots,W(N-1)]$ such that $w(0),\ldots,w(t-1)$ are the noise realizations up to time $t-1$ and $W(\tau)$ are random vectors with given probability distributions at time $\tau\geq t$. Define $\bar{u}(0:t-1:N)=[u^{\ast}(0),\ldots,u^{\ast}(t-1),u(t),\ldots,u(N-1)]$ to be the vector of input variables such that $u^{\ast}(0),\ldots,u^{\ast}(t-1)$ are the obtained optimal control inputs up to time $t-1$ and $u(t),\ldots,u(N-1)$ are the input variables that need to be determined at time $t\geq 0$.

Given STL formula $\varphi$, observations of state variables $x(0),x(1),\ldots,x(t)$, and designed control inputs $u^{\ast}(0),\ldots,u^{\ast}(t-1)$ of system (1), the stochastic SHMPC optimization problem minimizes the expected value of the cost function

[TABLE]

at each time instant $0\leq t<N$, as follows

[TABLE]

where the expected value $\mathbb{E}[\cdot]$ in (5a) is conditioned on observations $\tilde{X}(0:t+1)=[x(0),\ldots,x(t)]$ and $\delta_{t}=\delta/N$ for all $t$. Optimization variables in (5) are the control inputs $\tilde{u}(t:N)=[u(t),\ldots,u(N-1)]$. We indicate the argument of minimum by $\tilde{u}_{opt}(t:N)=[u_{opt}(t),\ldots,u_{opt}(N-1)]$.

The complete procedure of obtaining an optimal control sequence using SHMPC is presented in Algorithm 1. Lines 3 to 8 of this algorithm specify the inputs and the parameters used in the algorithm and line 20 specifies the output. In line 10, the SHMPC optimization procedure starts for each time step $t\in[0,N-1]$. In line 11, we solve the optimization problem (5) to obtain an optimal control sequence for time instance $t$. In lines 12 to 16, we check whether the obtained solution satisfies the STL specifications or not; if yes, assign the first component of the obtained input sequence to $u^{\ast}(t)$, and if not, the optimization procedure will be terminated. Finally, in line 17, we apply $u^{\ast}(t)$ to the system (1) and observe the states at time instant $t$.

We show in the following theorem that in Algorithm 1, by using the shrinking horizon technique, the specific choice of $\delta_{t}$, and keeping track of the control inputs and observed states, the closed-loop system satisfies the STL specification $\varphi$ with probability greater than or equal to $1-\delta$.

Theorem 2

Given a constant $\delta\in(0,1)$ and an STL formula $\varphi$, if the optimization problems in Algorithm 1 are all feasible, the computed optimal control sequence $\tilde{u}^{\ast}(0:N)=[u^{\ast}(0),\ldots,u^{\ast}(N-1)]$ ensures that the closed-loop satisfy $\varphi$ with probability greater than or equal to $1-\delta$.

*Proof: * Considering the chance constraint (5c) the probability that a trajectory of the system violates $\varphi$ at time step $t$ is at most $\delta_{t}$. This implies that the probability of violating $\varphi$ in the time interval $t=0,\ldots,N-1$ is at most $\sum_{t=0}^{N-1}\delta_{t}=\sum_{t=0}^{N-1}\delta/N=\delta$, which proves that the optimal control sequence $\tilde{u}^{\ast}=[u^{\ast}(0),\ldots,u^{\ast}(N-1)]$ obtained using Algorithm 1 results in trajectories that satisfy $\varphi$ with probability greater than or equal to $1-\delta.$

Note that in practice, if at each time step a feasible solution is not found, by using the previous control value, i.e., by setting $u^{\ast}(t)=u^{\ast}(t-1)$, we can give the controller a chance to retry in the next time step after observing the next state.

Remark 3

The choice of $\delta_{t}=\delta/N$ is completely arbitrary. In general, the positive constants $\delta_{t}$ can be picked freely with the condition that $\sum_{t=0}^{N-1}\delta_{t}=\delta$.

Computation of the solution of the optimization problem (5) requires addressing two main challenges: a) the objective function (5a) depends on the optimization variables $\tilde{u}(t:N)$ and on random variables $\tilde{W}(t:N)$, thus we have to compute the expected value as a function of these variables; and b) the feasible set of the optimization restricted by the chance constraint (5c) is in general difficult to characterize. We propose approximation methods in Sections V and VI to respectively address these two challenges.

V Approximating the objective function

To solve the optimization problem (5), one needs to calculate the expected value of the objective function. One way to do this is via numerical integration methods [7]. However, numerical integration is in general both cumbersome and time-consuming. For example, the method of approximating the density function of the disturbance with piecewise polynomial functions defined on polyhedral sets [5, 19] suffers from scalability issues on top of the induced approximation error. Therefore, in this section, we discuss an efficient method that computes an upper bound for the expected value of the objective function and then, minimize this upper bound instead.

We discuss computation of such upper bounds for both cases of process noise with arbitrary probability distribution and with normal distribution in Sections V-A and V-B, respectively. For this purpose, we first provide a canonical form for the robustness function of a STL formula $\psi$, which is the mix-max or max-min of random variables. This result is inspired by [9], in which the authors provide such canonical forms for max-min-plus-scaling functions.

Theorem 4

For a given STL formula $\psi$, the robustness function $\rho^{\psi}(\Xi(0:N),0)$, and hence the function $J_{\text{robust}}(\bar{X}(0:t:N))$, can be written into a max-min canonical form

[TABLE]

and into a min-max canonical form

[TABLE]

for some integers $K,L,n_{1},\ldots,n_{K},m_{1},\ldots,m_{L}$, where $\lambda_{ij}$ and $\gamma_{ij}$ are weighting vectors and $\eta_{ij}$ and $\zeta_{ij}$ are affine functions of $\bar{u}(0:t:N)$ and $x_{0}$.

*Proof: * The proof is inductive on the structure of $\psi$ and uses the explicit form of the states in (2) utilizing the identities $-\max(f_{1},f_{2})=\min(-f_{1},-f_{2})$ and

[TABLE]

for functions $f_{1},f_{2},g_{1}$, and $g_{2}$.

V-A Arbitrary probability distributions with bounded support

Suppose the elements of the stochastic vector $W(t)$, i.e., $W_{k}(t),\;k\in\{1,\dots,n\}$ have arbitrary probability distribution with known bounded support $I_{W_{k}(t)}=[a_{k},b_{k}]$ and its first moment $\mathbb{E}[W_{k}(t)]$ belongs to the interval $\mathbb{M}_{{W_{k}(t)}}=[c_{k},d_{k}]$, with known quantities $a_{k},b_{k},c_{k},d_{k}\in\mathbb{R}$. Under this assumption, the explicit form of $X(\cdot)$ in (2) implies that, for the observed value of $X(t)$ as $x(t)$, $X(\tau)$ is a random variable with the following interval of support and the first moment interval

[TABLE]

where $\bar{C}_{\tau}=\Phi(\tau,t)x(t)+\sum_{k=t}^{\tau-1}\Phi(\tau,k+1)B(k)u(k)$, and $\bar{a}_{\tau},\bar{b}_{\tau},\bar{c}_{\tau}$ and $\bar{d}_{\tau}$ are weighted sum of $a_{k},b_{k},c_{k},d_{k},\;k\in\mathbb{N}$, obtained by using interval arithmetics mentioned in Section I-A.

The objective function in (5) can be written as $\mathbb{E}\left[J_{\text{robust}}(\bar{X}(0:t:N+1))\right]+J_{\text{in}}(\bar{u}(0:t-1:N)))$ and that $J_{\text{robust}}(\bar{X}(0:t:N+1))=-\rho^{\psi}(\bar{X}(0:t:N+1),0)$. Recall that $\bar{X}(0:t:N+1)=[x(0),\ldots,x(t),X(t+1),\ldots,X(N)]$ with observed states $x(0),\ldots,x(t)$ of system (1) and random states $X(\tau),\,\,\tau>t$. The following theorem shows how we can compute an upper bound for $\mathbb{E}[J_{\text{robust}}(\bar{X}(0:t:N+1))]$ based on the canonical form of $J_{\text{robust}}$.

Theorem 5

For a given STL formula $\psi$, $\mathbb{E}\left[J_{\text{robust}}(\bar{X}(0:t:N+1))\right]$ can be upper bounded by

[TABLE]

where the constants $\eta_{ij}$, $i\in\{1,\dots,L\},j\in\{1,\dots,m_{i}\}$, are affine functions of $\bar{u}(0:t-1:N)$ and $x(0)$, and $\hat{d}_{ij}$ are weighted sum of $w(0),\ldots,w(t-1)$ and $c_{k},d_{k}$ for $k=t,\ldots,N-1$.

*Proof: * With focus on the canonical form (6), let $Y_{ij}=\eta_{ij}+\lambda_{ij}\bar{W}(0:t:N)$. Considering the support and moment interval of the components of $W(\tau),\tau=t,\ldots,N-1$, each random variable $Y_{ij}$ has the following support and moment interval (similar to (8))

[TABLE]

where the constants $\hat{a}_{ij},\hat{b}_{ij},\hat{c}_{ij},\hat{d}_{ij}$, $i\in\{1,\dots,L\},j\in\{1,\dots,m_{i}\}$, are weighted sum of $w(0),\ldots,w(t-1)$ and $a_{k},b_{k},c_{k},d_{k}$ for $k=t,\ldots,N-1$, using interval arithmetic (cf. Section I-A). Accordingly, $J_{\text{robust}}$ is a random variable with the following support and moment intervals,

[TABLE]

Hence, as we are minimizing the cost function in (5), we can utilize the upper bound $\max_{i\in\{1,\dots,L\}}\min_{j\in\{1,\dots,m_{i}\}}(\hat{d}_{ij}+\eta_{ij})$ for $\mathbb{E}\left[J_{\text{robust}}(\bar{X}(0:t:N+1))\right]$.

Note that the approximation methodology of Theorem 5 is applicable also to the min-max canonical form (7).

By replacing the expected objective function by its upper bound given in Theorem 5, and by replacing the probabilistic constraints by their equivalent linear approximation (as is discussed in Section VI), the optimization problem (5) can be then recast as a mixed integer linear programming (MILP) problem, which can be solved using the available MILP solvers [2, 21].

V-B Normal distribution

The upper bound on the objective function provided in the previous section does not apply to process noises with unbounded support, but knowing the distribution of the process noise provides more information about the statistics of the runs of the system. In this section we address process noises with normal distribution separately due the their wide use in engineering applications.

Suppose that for any $t\in\mathbb{N}$, $W(t)$ is normally distributed with mean $\mathbb{E}[W(t)]=0$ and covariance matrix $\Sigma_{W(t)}$, i.e., $W(t)\sim\mathcal{N}(0,\Sigma_{W(t)})$. The explicit form of $X(\tau)$ in (2) and the fact that normal distribution is closed under affine transformations result in normal distribution for $X(\tau)$, $\tau\in\mathbb{N}$. Its expected value and covariance matrix with an observed value $x(t)$ of $X(t)$ are

[TABLE]

respectively, for $\tau\geq t\geq 0$.

In this section we use the canonical representation of $J_{\text{robust}}(\bar{X}(0:t:N+1))$ in Theorem 4, which states that $J_{\text{robust}}$ (for fixed $\bar{u}(0:t:N)$ and $x_{0}$) can be written in either of the forms

[TABLE]

with $Y_{ij}=\eta_{ij}+\lambda_{ij}\bar{W}(0:t-1:N)$ being affine functions of the process noise, thus normally distributed random variables (similar to $X(\tau)$ explained above). With focus on these canonical representations for $J_{\text{robust}}$ we employ Proposition 6 to show how to approximate $\mathbb{E}\left[J_{\text{robust}}\right]$ using higher order moments of $W(t)\sim\mathcal{N}(0,\Sigma)$. This proposition, also used in [13], follows due to the relation between the infinity norm and the $p$-norm of a vector and Jensen’s inequality.

Proposition 6

Consider random variables $Z_{i}$ for $i\in\{1,\dots,s\}$ and let $p$ be an even integer. Then

[TABLE]

Founded on Proposition 6, next theorem shows how we can upper bound $\mathbb{E}\left[J_{\text{robust}}\right]$ using the higher order moments of $Y_{ij}$.

Theorem 7

Considering the canonical forms in (11) for $J_{\text{robust}}$ as a function of random variables $Y_{ij}$, $\mathbb{E}\left[J_{\text{robust}}\right]$ can be upper bounded by

[TABLE]

*Proof: * For random variables $Y_{ij},\;i\in{\{1,\ldots,L\}},\;j\in{\{1,\ldots,{m_{i}}\}}$, and for a positive even integer $p$, the following inequality holds,

[TABLE]

where in $(i)$ we used the upper bound obtained in Proposition 6; in $(ii)$ we used the fact that $\min_{k\in{\{1,\ldots,r\}}}(\alpha_{k})=-\max_{k\in{\{1,\ldots,r\}}}(-\alpha_{k})$; In $(iii)$ we use again the inequality in Proposition 6. Moreover, for $i\in{\{1,\ldots,K\}},\;j\in{\{1,\ldots,n_{i}\}}$, the following inequality holds,

[TABLE]

where we apply Jensen’s inequality to the concave function $\min(\cdot)$ to get $(i)$. The inequality of Proposition 6 gives $(ii)$.

Note that random variables $Y_{ij}$ are normally distributed in both (12) and (13). Higher order moments of normally distributed random variables can be computed analytically in a closed form as a function of the first two moments, i.e., using its mean and variance. More specifically, for a normally distributed random variable $Z$ with mean $\mu$ and variance $\sigma^{2}$, the $p$-th moment has a closed form as

[TABLE]

where $i$ is the imaginary unit and

[TABLE]

is the $p$-th Hermite polynomial [1, Chapter 22 and 26]. We use (14) to compute higher order moments of normal random variables with $p$ being even integers. Note that the right-hand side of (14) is in fact real because $H_{p}(z)$ contains only even powers of $z$ when $p$ is even.

In the next section we discuss how to cope with the second challenge of characterizing the feasible set of the optimization restricted by the chance constraint (5c).

VI Under Approximation of Chance Constraints

In this section, we discuss methods for computing conservative lower approximations of the chance constraints in (5c) as linear constraints. For the sake of compact notation, we indicate the stochastic run $\Xi(0:N)=X(0)X(1)\ldots X(N)$ only by $\Xi_{N}$ without declaring its dependency on the state, input, and disturbance variables. Recall the chance constraint (5c) as $\text{Pr}\left[\!(\Xi_{N},t)\!\models\!\varphi\!\right]\geq 1-\delta_{t}.$ In order to transform this constraint to linear inequalities, we first show in the following theorem, that this constraint can be transformed into similar inequalities but $\varphi$ being an atomic predicate. Then in Sections VI-A and VI-B, we discuss how to transform the resulting constraints with atomic predicates into linear inequalities for the cases of arbitrary random variables with known bounded support and moment interval and of normally distributed random variables.

Theorem 8

for any formula $\varphi$ and a constant $\vartheta\in(0,1)$, constraints of the forms

[TABLE]

can be transformed into similar constraints with $\varphi$ being an atomic predicate using the structure of $\varphi$.

*Proof: * The proof is inductive on the structure of the formula $\varphi$ as discussed in the following three cases.

Case I: $\varphi=\neg\varphi_{1}$ we have the following equivalences

[TABLE]

Case II: $\varphi=\varphi_{1}\wedge\varphi_{2}$ we obtain the following inequalities by using the fact that for possibly joint events $\mathcal{A}$ and $\mathcal{B}$, it holds that $\text{Pr}[\mathcal{A}\wedge\mathcal{B}]\geq\vartheta\Leftrightarrow\text{Pr}(\neg\mathcal{A}\vee\neg\mathcal{B})\leq 1-\vartheta$ and $\text{Pr}(\mathcal{A}\vee\mathcal{B})\leq\text{Pr}[\mathcal{A}]+\text{Pr}[\mathcal{B}]$.

[TABLE]

Note that in the last line of (17), we assume that the probability of the two events are upper bounded by the same value, i.e., $(1-\vartheta)/2$. However, this can be replaced by any two other probabilities $\delta_{1}$ and $\delta_{2}$ such that $\delta_{1}+\delta_{2}=1-\vartheta$. Now consider the second possibility:

[TABLE]

where the last line of (18) is due to the fact that the events are disjoint. Assuming that the probabilities of these two events are lower bounded by the same values, i.e., $(1-\vartheta)/2$, we have the inequalities

[TABLE]

which are in the form of inequalities discussed previously. Note that Equations (17) to (19) discuss the case of having conjunction of two STL formulas. The results can be easily extended to conjunction of $n$ STL formulas by replacing $(1-\vartheta)/2$ with $(1-\vartheta)/n$.

Case III: $\varphi=\varphi_{1}\mathcal{U}_{[a,b]}\varphi_{2}$ The satisfaction $(\Xi_{N},t)\models\varphi_{1}\mathcal{U}_{[a,b]}\varphi_{2}$ is equivalent to $\bigvee_{j=t+a}^{t+b}\psi_{j}$ with disjoint events

[TABLE]

Thus $\text{Pr}\left[(\Xi_{N},t)\models\varphi_{1}\mathcal{U}_{[a,b]}\varphi_{2}\right]\geq\vartheta$ is equivalent to $\sum_{j=t+a}^{t+b}\text{Pr}[\psi_{j}]\geq\vartheta$. Assuming the probabilities of events are lower bounded by the same values, we have $\text{Pr}[\psi_{j}]\geq\vartheta/(b-a+1)$ for $j=a+t,\dots,b+t$, which again can be reduced as in Case II.

The second possible probabilistic constraint in Case III can be obtained as

[TABLE]

which can be again reduced as in Case II. Here also, we used the fact that $\psi_{j}$ consists of disjoint events and we assume that he probabilities of events are lower bounded by the same value, i.e., by $\vartheta/(b-a+1)$, for $j=a+t,\dots,b+t$.

So far we have shown how to inductively reduce the chance constraint (5c) to inequalities of the form (16) with atomic predicates. In the rest of this section we discuss their corresponding linear inequalities for the two types of probability distributions considered in this paper.

VI-A Arbitrary probability distributions with bounded support

To transform the chance constraints into linear constraints in the case of having random variables with arbitrary probability distributions, we apply an approximation method based on the upper bound proposed by [4]. Let $Z_{1},\ldots,Z_{n}$ be random variables with interval of bounded support $[a_{i},b_{i}]$ and let $\mathbb{E}[Z_{1}],\ldots,\mathbb{E}[Z_{n}]$ denote their expected values belonging to the moment intervals $\mathbb{M}_{i}$ for $i=1,\ldots,n$. Define $Z=\sum_{i=1}^{n}Z_{i}$ and $\mathbb{E}(Z)=\sum_{i=1}^{n}\mathbb{E}[Z_{i}]$. Using Chernoff-Hoeffding inequality, the following upper bound exists for any $\varsigma\geq 0$ [16]

[TABLE]

where $\nu>0$ is a constant. If $Z_{1},\ldots,Z_{n}$ are dependent, then the inequality applies with a constant $\nu=\chi(\hat{G})/2$, where $\hat{G}$ denotes the indirected dependency graph of $Z_{1},\ldots,Z_{n}$ and $\chi(\hat{G})$ is the chromatic number of the graph $\hat{G}$ defined as the minimum number of colors required to color $\hat{G}$. For the independent case, $\chi(\hat{G})=1$. The expression for the right tail probability is derived identically. For more details, the reader is referred to [4].

Consider the chance constraints (16) with $\varphi=\{\alpha\geq 0\}$. Since $\alpha$ is an affine function of random state variables, it is a random variable itself with the following interval of support and moment interval

[TABLE]

where for $t=0,\ldots,N$, we have $\tilde{a}_{t},\tilde{b}_{t},\tilde{c}_{t}$ and $\tilde{d}_{t}$ are weighted sum of $\bar{a}_{t},\bar{b}_{t},\bar{c}_{t},\bar{d}_{t}$ related to the interval of support and moment interval of random variables $X(t)$ (cf. (8)), and $\tilde{C}_{t}$ is a linear expression of input variables.

Let $\varsigma=\mathbb{E}\left[\alpha(X(t))\right]$; we can directly use (22) as

[TABLE]

Note that since $\delta_{t}\in(0,1)$, we have $\log(\delta_{t})<0$; hence, by multiplying both sides of the inequality by -1 in line 5 of (24), the expression $-\log(\delta_{t})\cdot\sum_{t=1}^{N}(\tilde{b}_{t}-\tilde{a}_{t})^{2}$ becomes a positive number, and hence, its square root is a real number. Note also that the last inequality is due to the fact that $\varsigma\geq 0$. Hence, we can replace $\varsigma$ in the last inequality of (24) by the lower bound of its moment interval in (23), i.e., with $\tilde{c}_{t}+\tilde{C}_{t}$, which is a linear expression in the input variables.

Consequently, in this case, the chance constraint in (5) can be replaced by

[TABLE]

For the second type of probabilistic inequality (cf. (16)), we can again use (22) for the right tail probability; hence we have

[TABLE]

and then following the same steps as in (24), we obtain the same linear expression for the chance constant as in (25) by only replacing $\delta_{t}$ by $1-\delta_{t}$ in the related expressions.

VI-B Normal distribution

To transform the chance constraints into linear constraints in the case of having normally distributed random variables, we use the quantile of the normal distribution. By definition, for a normally distributed random variable $x$ with mean $\mu$ and standard deviation $\sigma$,

[TABLE]

where $F^{-1}$ denotes the inverse of the cumulative distribution function or the quantile function and $\phi^{-1}$ is the inverse of the error function of a normally distributed random variable.

Recall the chance constraints (16) with $\varphi=\{\alpha\geq 0\}$. Since $\alpha$ is an affine function of normally distributed state variables, it is also normally distributed with appropriately defined mean $\mu_{t}$ and variance $\sigma_{t}^{2}$. Hence, we can directly use (27) and (28) as

[TABLE]

Therefore, the chance constraint can be replaced by the equivalent linear constraint (29) or (30), depending on the type of constraint we have.

VII Experimental Results

We now apply our controller synthesis approach to the room temperature control in a building. The details of the thermal model can be found in [15, 24], and is briefly explained here for clarity. The temperature of room $r$ is denoted by $T_{r}$ and the wall and the temperature of the wall between the room and its surrounding $j$ (e.g. other rooms) are denoted by $w_{j}$ and $T_{w_{j}}$, respectively. Dynamics of the temperature of wall $w_{j}$ and room $r$ can be written as [15]

[TABLE]

where $C_{j}^{w},\alpha_{j}$ and $A_{w_{j}}$ are heat capacity, a radiative heat absorption coefficient, and the area of $w_{j}$, respectively. The total thermal resistance between the centerline of wall $j$ and the side of the wall on which node $k$ is located is denoted by $R_{j_{k}}$. The radiative heat flux density on $w_{j}$ is denoted by $Q_{\text{rad}j}$, the set of all neighboring nodes to $w_{j}$ is denoted by $\mathcal{N}_{w_{j}}$, and $r_{j}$ is a wall identifier, which equals 0 for internal walls and 1 for peripheral walls, where $j$ is the outside node. The temperature, heat capacity and air mass flow into room $r$ are denoted by $T_{r},C_{j}^{r}$ and $\dot{m}_{r}$, respectively; $c_{a}$ is the specific heat capacity of air, and $T_{s}$ is the temperature of the supply air to room $r$, $w$ is a window identifier, which equals 0 if none of the walls surrounding room $r$ have windows, and 1 if at least one of them does, $\tau_{w}$ is the transmissivity of the glass of the window in room $r$, $A_{\text{win}}$ is the total area of the windows on walls surrounding room $r$, $Q_{\text{rad}}$ is the radiative heat flux density per unit area radiated to room $r$, and $\dot{Q}_{\text{int}r}$ is the internal heat generation in room $r$. Finally, $\mathcal{N}_{r}$ is the set of neighboring room nodes for room $r$. Further details on this thermal model can be found in [15].

As such, the heat transfer equations for each wall and room $r$ is in the form of nonlinear differential equation. After linearization and time-discretization, the model of the system becomes in the form of dynamical equation

[TABLE]

where $X\in\mathbb{R}^{n}$ is the state vector representing the temperature of the walls and the rooms and $u\in\mathbb{R}^{m}$ is the input vector representing the air mass flow rate and discharge air temperature of conditioned air into each thermal zone. Matrices $A$ and $B$ are obtained by time discretization of dynamics of the thermal model (31)-(32) with a sampling time of $t_{s}=30$ minutes. The disturbance $W(\cdot)$ aggregates various unmodeled dynamics and the uncertainty in physical variables such as the outside temperature, internal heat generation and radiative heat flux density. The statistics of $W(\cdot)$ can be estimated using historical data [15].

In this example, we only control the temperature of one room and include the temperature of the neighboring rooms as part of the disturbance signal $W(t)$. We also assume that there is a reference for the disturbance signal, denoted by $w_{r}(t)$, and the reference is perturbed by independent and identically distributed random vectors $e(t)\sim\mathcal{N}(0,\mathbb{I}_{n})$, i.e., the disturbance is $W(t)=w_{r}(t)+e(t)$, which is normally distributed with mean $\mu_{t}=w_{r}(t)$ and identity covariance matrix $\Sigma=\mathbb{I}_{n}$.

In contrast to [24], which considers deterministic disturbances from a bounded set, we consider stochastic disturbances and we maximize the robustness of satisfaction of the STL specifications in the presence of such disturbance. Accordingly, we handle chance constraints and include the expected value of the robustness function in the objective function.

We consider a signal $\text{occ}:\mathbb{N}\rightarrow\{-1,1\}$ representing the room occupancy; $\text{occ}(t)=1$ if the room is occupied at time $t$ and $\text{occ}(t)=-1$ otherwise. This signal is assumed to be known for the entire simulation period. The MPC prediction horizon $N$ is chosen to be $24$, representing 12 hours monitoring of the room temperature. We select $\delta=0.1$ so that the obtained control input provides confidence level of $90\%$ on the satisfaction of the desired behavior. We are interested in keeping the room temperature above a reference temperature $T_{r}$ when the room is occupied; thus the specification is

[TABLE]

At each time instant $0\leq t<N$, the optimization problem (5) obtains an optimal control input $\tilde{u}_{opt}(t:N)=[u_{opt}(t),\ldots,u_{opt}(N-1)]$ that minimizes

[TABLE]

where the robustness function is defined as

[TABLE]

The chance constraint (5c) is defined with the same specification $\varphi=\psi$. We approximate $\mathbb{E}[-\rho^{\psi}(\tilde{X}(0:t:N),0)]$ using the upper bound (13) and transform the chance constraint (5c) into linear inequalities using the approach of Section VI. We also assume that inputs are bounded, i.e., for each $0\leq t<N$, we have $0\leq u(t)\leq 380$.

The simulations are done using Matlab R2014b on a 2.6 GHz Intel Core i5 processor and the optimizations are solved using fmincon solver in Matlab. we perform $n_{s}=200$ simulations in order to check the satisfiability of the STL specifications with a probability greater than or equal to $0.9$. Figure 1 shows the results of these $200$ simulations. The top plot shows the occupancy signal and the middle plot illustrates the average, minimum, and maximum of the obtained room temperatures over 12 hours as well as the room reference temperature $T_{r}$ in Fahrenheit. The controller ensures that the room temperature goes above the reference temperature $T_{r}$ once the occupancy signal is 1 and stays there as long as the room is occupied. The minimum and maximum bounds on the room temperature shows that the specifications have never been violated in these $200$ simulations. The bottom plot shows the average, minimum, and maximum of the air flow rate in $\left[\frac{\text{ft}^{3}}{\text{min}}\right]$, which indicates that the input constraint is not violated.

Note that all these $n_{s}=200$ runs result in feasible solutions, which gives a confidence bound on the feasibility of the original problem as follows. Since all the $n_{s}$ runs of the simulation are feasible, we can claim that the original problem is also feasible with probability at least $(\beta/2)^{1/n_{s}}$, $\beta\!\in\!(0,1)$, with confidence level $1\!-\!\beta$ [6]; hence, having $200$ runs being all feasible, the optimization problem (5) is also feasible with probability $0.98$ with confidence level $0.95$.

To further illustrate the performance of the proposed method, we compare our SHMPC approach with the robust MPC (RMPC) approach of [24], in which the disturbance belongs to a bounded polyhedral set. Note that RMPC approach is not directly applicable to unbounded uncertainties. Therefore, in the optimization procedure, we truncate a normally distributed disturbance in the $2\sigma$ interval such that $e(t)\in[-1,1]$. Further, we solve the RMPC optimization problem using Monte Carlo sampling.

The total fan energy consumption is proportional to the cubic of the airflow. Table I shows the total fan energy consumption and the computation time for the three approaches. For RMPC and SHMPC, we report the average and standard deviation of total energy consumption using the sum of cubes of the optimal input sequences corresponding to the $200$ simulations. Also, for these two approaches, we report the average computation time over the $200$ simulations. Comparing statistics of these two approaches is essential because of the chance constraints in SHMPC and the Monte Carlo sampling based optimization in RMPC. The energy consumption using open-loop optimal control (OC) is very high, comparing to both RMPC and SHMPC. This is due to the fact that the open-loop strategy computes the solution of optimization problem (5) only once and hence, the computation time is smaller compared to the two other methods. As a result, the input sequence has an aggressive behavior to make sure that it reacts in time to the changes happening in the system. Since RMPC is more conservative compared to SHMPC, the average energy consumption is much higher for the RMPC controller compared to the SHMPC controller: the SHMPC controller achieves a $80\%$ reduction of total energy consumption on average compared to RMPC.

VIII Conclusions

In this paper, we presented shrinking horizon model predictive control (SHMPC) for stochastic discrete-time linear systems with signal temporal logic (STL) specifications. Our aim was to obtain an optimal control sequence that guarantees the satisfaction of STL specifications with a probability greater than a certain level. By assumption, the stochastic disturbance signal had an arbitrary probability distribution with a bounded support and the only available information related to this distribution is the intervals of support and the moment intervals of each component of the disturbance signal. Using an existing approximation technique, we showed that the chance constraints could be approximated by an upper bound, which resulted in having approximate linear constraints for the chance constraints. Moreover, in the case of having the state costs and/or the robustness function related to the degree of satisfaction of the specifications by the state trajectory, their expected value can be also approximated using the moment intervals of components of the disturbance signal. As an additional case, we further considered disturbances that are normally distributed and we showed that the chance constraints in this case can be replaced by the quantile expressions which are linear in the input variables. In the end, in an example, we applied the proposed method to control a HVAC system.