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

arXiv:1705.02152v1 [] 5 May 2017

Samira S. Farahani∗ , Rupak Majumdar∗ , Vinayak Prabhu∗ , Sadegh Esmaeil Zadeh Soudjani∗

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.

I. I NTRODUCTION
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
∗ The authors are with the Max Planck Institute for Software Systems, Germany. A limited subset of the results of this paper is accepted for presentation at American Control Conference 2017 [11].

farahani,vinayak,rupak,sadegh@mpi-sws.org

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

A. Notations
We use R for the set of reals and N := {0, 1, 2, . . .} for the set
of non-negative integers. The set B := {>, ⊥} indicates logical
true and false. For a vector v ∈ Rs , its components are denoted
by vk , k ∈ {1, . . . , s}. For a sequence {v(t) ∈ Rs , t ∈ N} and
t1 < t2 , we define ṽ(t1 : t2 ) := [v(t1 ), v(t1 + 1), . . . , v(t2 − 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) ∈ Rs , t ∈ N} and t1 ≤ t < t2 , we define Ȳ (t1 : t : t2 ) :=
[y(t1 ), y(t1 + 1), . . . , y(t),Y (t + 1), . . . ,Y (t2 − 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 IY (t) and
the first moment1 by E[Y (t)].
We consider operations on intervals according to interval
arithmetic: for two arbitrary intervals [a, b] and [c, d], and
constants λ , γ ∈ R, we have [a, b] + [c, d] = [a + c, b + d] and
λ · [a, b] + γ = [min(λ a, λ b) + γ, max(λ a, λ b) + γ].
II. D ISCRETE -T IME S TOCHASTIC L INEAR S YSTEMS
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 Rn is defined
as the set {x ∈ RN | PrX [B(x, r)] > 0, ∀r > 0}, where 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 PrX [C] = 1.
Consider a (time-variant) discrete-time stochastic system
modeled by the difference equation
X(t + 1) = A(t)X(t) + B(t)u(t) +W (t), X(0) = x0 ,

(1)

where X(t) ∈ Rn denotes the state of the system at time instant
t, u(t) ∈ Rm denotes the control input at time instant t, and
W (t) ∈ Rs 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
1 The expected value of a random variable X with support D and the
R
cumulative distribution function F is defined as E[X] = D xdF(x). The
expectation exists if the integral is well-defined and yields a finite value.

as the process noise or an adversarial disturbance. Matrices
A(·) ∈ Rn×n , and B(·) ∈ Rn×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), . . . ,W (t)
are mutually independent random vectors for all time instants
t. Note that, for any t ∈ N, the state-space model (1) provides
the following explicit form for X(τ), τ ≥ t, as a function of
X(t), input u(·), and the process noise W (·):
τ−1

X(τ) = Φ(τ,t)X(t) + ∑ Φ(τ, k + 1) (B(k)u(k) +W (k)) , (2)
k=t

where Φ(·, ·) is the state transition matrix of model (1) defined
as
(
A(τ − 1)A(τ − 2) . . . A(t) τ > t ≥ 0
Φ(τ,t) =
In
τ = t ≥ 0,
with In being the identity matrix.
For a fixed positive integer N, and a given time instant t ∈ N,
let ũ(t : N) = [u(t), u(t + 1), . . . , u(N − 1)] (matrix W̃ (t : N) is
defined similarly). A finite stochastic run of system (1) for
the time interval [t : N] is defined as Ξ(t : N) = X(t)X(t +
1) . . . X(N), which is a finite sequence of states satisfying (2).
Since each state X(τ) depends on X(t), ũ(t : N), and W̃ (t : N),
we can rewrite Ξ(t : N) in a more elaborative notation as
ΞN (X(t), ũ(t : N), W̃ (t : N)). Analogously, we define an infinite
stochastic run Ξ = X(t)X(t + 1)X(t + 2) . . . as an infinite
sequence of states. Stochastic runs will be used in Section
III to define the system’s specifications.
III. S IGNAL T EMPORAL L OGIC
An infinite run of system (1) can be considered as a
signal ξ = x(0)x(1)x(2) . . . , 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]
ϕ ::= > | π | ¬ϕ | ϕ ∧ ψ | ϕU[a,b] ψ
where > is the true predicate; π is a predicate whose truth
value is determined by the sign of a function, i.e. π =
{α(x) ≥ 0} with α : Rn → R being an affine function of state
variables; ψ is an STL formula; ¬ and ∧ indicate negation and
conjunction of formulas; and U[a,b] is the until operator with
a, b ∈ R≥0 . A run ξ satisfies ϕ at time t, denoted by (ξ ,t) |= ϕ,
if the sequence x(t)x(t + 1) . . . satisfies ϕ. Accordingly, ξ
satisfies ϕ, if (ξ , 0) |= ϕ.
Semantics of STL formulas are defined as follows. Every
run satisfies >. The run ξ satisfies ¬ϕ if it does not satisfy
ϕ; it satisfies ϕ ∧ ψ if both ϕ and ψ hold. For a run
ξ = x(0)x(1)x(2) . . . and a predicate π = {α(x) ≥ 0}, we have
(ξ ,t) |= π if α(x(t)) ≥ 0. Finally, (ξ ,t) |= ϕU[a,b] ψ if ϕ holds
at every time step starting from time t before ψ holds, and
additionally ψ holds at some time instant between a + t and
b + t. Additionally, we derive the other standard operators
as follows. Disjunction ϕ ∨ ψ := ¬(¬ϕ ∧ ¬ψ), the eventually
operator as [a,b] ϕ := >U[a,b] ϕ, and the always operator as
2[a,b] ϕ := ¬ [a,b] ¬ϕ.
2 2

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.

2

Thus (ξ ,t) |= [a,b] ϕ if ϕ holds at some time instant
between a + t and b + t and (ξ ,t) |= 2[a,b] ϕ if ϕ holds at
every time instant between a + t and b + t.
Formula Horizon. The horizon of an STL formula ϕ is the
smallest n ∈ N such that the following holds for all signals
ξ = x(0)x(1)x(2) . . . and ξ 0 = x0 (0)x0 (1)x0 (2) . . .:
If x(t + i) = x0 (t + i) for all i ∈ {0, . . . , n}
Then (ξ ,t) |= ϕ iff (ξ 0 ,t) |= ϕ.
Thus, in order to determine whether a signal ξ satisfies an STL
formula ϕ, we can restrict our attention to the signal prefix
x(0), . . . , x(∆) where ∆ is the horizon of ϕ. This horizon can be
upper-approximated by a bound, denoted by len(ϕ), defined to
be the maximum over the sums of all nested upper bounds on
the temporal operators. Formally, len(ϕ) is defined recursively
as:
ϕ := > ⇒ len(ϕ) = 0,

ϕ := π ⇒ len(ϕ) = 0,

ϕ := ¬ϕ1 ⇒ len(ϕ) = len(ϕ1 ),
ϕ := ϕ1 ∧ ϕ2 ⇒ len(ϕ) = max(len(ϕ1 ), len(ϕ2 )),
ϕ := ϕ1 U[a,b] ϕ2 ⇒ len(ϕ) = b + max(len(ϕ1 ), len(ϕ2 )),

2

where ϕ1 , ϕ2 and ψ are STL formulas. For example, for ϕ =
[0,4] [3,6] π, we have len(ϕ) = 4 + 6 = 10. For a given STL
formula ϕ, it is possible to verify that ξ |= ϕ using only the
finite run x(0)x(1) . . . x(N), where N is equal to len(ϕ).
STL Robustness. In contrast to the above Boolean semantics,
the quantitative semantics of STL [18] assigns to each formula
ϕ a real-valued function ρ ϕ of signal ξ and t such that
ρ ϕ (ξ ,t) > 0 implies (ξ ,t) |= ϕ. Robustness of a formula ϕ
with respect to a run ξ at time t is defined recursively as

Note that a general STL formula ϕ consists of several
Boolean and/or temporal operators. Due to the system dynamics (1), the stochastic run Ξ(t : t +N) and ρ ϕ (Ξ(t : t +N),t) are
both functions of W̃ (t : t + N). Therefore, ρ ϕ (Ξ(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 ϕ and constant δ ∈ (0, 1), a stochastic run
Ξ = X(0)X(1)X(2) . . . satisfies ϕ with probability greater
than or equal to 1 − δ , if the finite stochastic run
Ξ(0 : N) = X(0)X(1) . . . X(N) with N ≥ len(ϕ) satisfies
Pr [ρ ϕ (Ξ(0 : N), 0) > 0] ≥ 1 − δ .
IV. C ONTROL P ROBLEM S TATEMENT
For system (1) with a given initial state X(0) = x0 , the
stochastic disturbance vector W (t) with a given probability distribution, STL formulas ϕ and ψ, and some constant N ≥ max(len(ϕ), len(ψ)), the control problem can be
defined as finding an optimal input sequence ũ∗ (0 : N) =
[u∗ (0), . . . , u∗ (N − 1)], that minimizes the expected value of
a given objective function J(X̃(0 : N + 1), ũ(0 : N)) subject to
constraints on states and input variables, where X̃(0 : N + 1)=
[X(0), X(1), . . . , X(N)]. This optimization problem for the time
interval 0 ≤ t < N can be defined as


min E J(X̃(0 : N + 1), ũ(0 : N))
s.t.
(3a)
ũ(0:N)

t−1

X(t) = Φ(t, 0)x0 + ∑ Φ(t, k + 1) (B(k)u(k) +W (k)) , (3b)
k=0


Pr ΞN (x0 , ũ(0 : N), W̃ (0 : N)) |= ϕ ≥ 1 − δ ,
(3c)
ũ(0 : N) ∈ U N ,

ρ > (ξ ,t) = +∞,

(3d)

π

ρ (ξ ,t) = α(x(t)) with π = {α(x) ≥ 0},
ρ ¬ϕ (ξ ,t) = −ρ ϕ (ξ ,t)
ρ ϕ∧ψ (ξ ,t) = min(ρ ϕ (ξ ,t), ρ ψ (ξ ,t)),

ρ ϕ U[a,b] ψ (ξ ,t) = max min(ρ ψ (ξ ,t + i), min ρ ϕ (ξ ,t + j)) ,
j∈[0,i)

i∈[a,b]

2

where x(t) refers to signal ξ at time t. The robustness
of the derived formula
[a,b] ϕ can be worked out to
2

be ρ [a,b] ϕ (ξ ,t) = maxi∈[a,b] ρ ϕ (ξ ,t + i); and similarly for
2[a,b] ϕ as ρ 2[a,b] ϕ (ξ ,t) = mini∈[a,b] ρ ϕ (ξ ,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 Ξ = X(0)X(1)X(2) . . . and using the bound
of a formula ϕ, the finite stochastic run Ξ(t : t + N) =
X(t)X(1) . . . X(t + N) with N = len(ϕ) is sufficient to study
probabilistic properties of (Ξ,t) |= ϕ. Analogous to the definition of robustness for deterministic run, we can define
stochastic robustness ρ ϕ (Ξ,t) of a formula ϕ with respect
to the run Ξ for times t with the stochastic robustness being
dependent on Ξ(t : t + N) and ϕ.

where E[·] denotes the expected value operator and the closed
set U N ∈ RmN specifies the constraint set for the input
variables. The chance constraints (3c) state that for a given
δ ∈ (0, 1), stochastic runs of the system should satisfy ϕ with
a probability greater than or equal to 1 − δ . We consider the
following objective function
J(X̃(0 : N + 1), ũ(0 : N)) := Jrobust (X̃(0 : N + 1)) + Jin (ũ(0 : N)),
(4)
where the first term Jrobust (X̃(0 : N + 1)) := −ρ ψ (X̃(0 : N +
1), 0) represents the negative value of the robustness function
on STL formula ψ at time 0 that needs to be minimized;
and the second term Jin (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.

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 ϕ, i.e., len(ϕ), 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 X̄(0 : t : N + 1) = [x(0), . . . , x(t), X(t + 1), . . . , X(N)]
where x(0), . . . , x(t) are the observed states up to time t and
X(τ) is the random state variable at time τ > t, also let W̄ (0 :
t − 1 : N) = [w(0), . . . , w(t − 1),W (t),W (t + 1), . . . ,W (N − 1)]
such that w(0), . . . , w(t − 1) are the noise realizations up to
time t − 1 and W (τ) are random vectors with given probability distributions at time τ ≥ t. Define ū(0 : t − 1 : N) =
[u∗ (0), . . . , u∗ (t − 1), u(t), . . . , u(N − 1)] to be the vector of
input variables such that u∗ (0), . . . , u∗ (t − 1) are the obtained
optimal control inputs up to time t − 1 and u(t), . . . , u(N − 1)
are the input variables that need to be determined at time t ≥ 0.
Given STL formula ϕ, observations of state
variables x(0), x(1), . . . , x(t), and designed control inputs
u∗ (0), . . . , u∗ (t − 1) of system (1), the stochastic SHMPC
optimization problem minimizes the expected value of the
cost function
J(X̄(0 : t :N + 1), ū(0 : t − 1 : N)) =
Jrobust (X̄(0 : t : N + 1)) + Jin (ū(0 : t − 1 : N)),
at each time instant 0 ≤ t < N, as follows
min E [J(X̄(0 : t : N + 1), ū(0 : t − 1 : N))] s.t.

(5a)

ũ(t:N)

τ−1

X(τ) = Φ(τ,t)x(t) + ∑ Φ(τ, k + 1) (B(k)u(k) +W (k)) ,
k=t

for t ≤ τ ≤ N

(5b)

Pr [ΞN (x0 , ū(0 : t − 1 : N), W̄ (0 : t − 1 : N)) |= ϕ] ≥ 1 − δt
(5c)
ũ(t : N) ∈ U N−t ,

(5d)

where the expected value E[·] in (5a) is conditioned on
observations X̃(0 : t + 1) = [x(0), . . . , x(t)] and δt = δ /N for
all t. Optimization variables in (5) are the control inputs
ũ(t : N) = [u(t), . . . , u(N − 1)]. We indicate the argument of
minimum by ũopt (t : N) = [uopt (t), . . . , uopt (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 ∈ [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∗ (t),
and if not, the optimization procedure will be terminated.
Finally, in line 17, we apply u∗ (t) to the system (1) and observe
the states at time instant t.
Algorithm 1
1: procedure
C HANCE -C ONSTRAINED
S TOCHASTIC
SHMPC
2: input:
3:
STL formulas ϕ and ψ and a fixed δ ∈ (0, 1)
4: parameters:
5:
N ≥ max(len(ϕ), len(ψ))
6:
probability distributions of the process noise
{W (t), t = 0, . . . , N − 1}
7:
Initial state x0
8:
δt = δ /N for t = 0, . . . , N − 1
9: closed-loop optimization problem:
10:
for t = 0; t < N; t = t + 1 do
11:
Compute ũopt (t : N) = [uopt (t), . . . , uopt (N − 1)] by
solving the optimization problem (5)
12:
if the solution of optimization problem (5) exists
then
13:
Let u∗ (t) := uopt (t);
14:
else
15:
Return Infeasible Solution and terminate the
optimization procedure;
16:
end if
17:
Apply u∗ (t) to the system and observe the value
of X(t + 1) as x(t + 1)
18:
end for
19: output:
20:
u∗ (0 : N) = [u∗ (0), . . . , u∗ (N − 1)]
21: end procedure
We show in the following theorem that in Algorithm 1, by
using the shrinking horizon technique, the specific choice of
δt , and keeping track of the control inputs and observed states,
the closed-loop system satisfies the STL specification ϕ with
probability greater than or equal to 1 − δ .
Theorem 2: Given a constant δ ∈ (0, 1) and an STL formula ϕ, if the optimization problems in Algorithm 1 are all
feasible, the computed optimal control sequence ũ∗ (0 : N) =
[u∗ (0), . . . , u∗ (N − 1)] ensures that the closed-loop satisfy ϕ
with probability greater than or equal to 1 − δ .

Proof: Considering the chance constraint (5c) the probability that a trajectory of the system violates ϕ at time step
t is at most δt . This implies that the probability of violating
N−1
ϕ in the time interval t = 0, . . . , N − 1 is at most ∑t=0
δt =
N−1
∑t=0 δ /N = δ , which proves that the optimal control sequence
ũ∗ = [u∗ (0), . . . , u∗ (N − 1)] obtained using Algorithm 1 results
in trajectories that satisfy ϕ with probability greater than or
equal to 1 − δ .
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∗ (t) = u∗ (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 δt = δ /N is completely arbitrary.
In general, the positive constants δt can be picked freely with
N−1
the condition that ∑t=0
δt = δ .
Computation of the solution of the optimization problem
(5) requires addressing two main challenges: a) the objective
function (5a) depends on the optimization variables ũ(t : N)
and on random variables 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. A PPROXIMATING 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 timeconsuming. 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 ψ, which is the mixmax 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 ψ, the robustness
function ρ ψ (Ξ(0 : N), 0), and hence the function Jrobust (X̄(0 :
t : N)), can be written into a max-min canonical form

Jrobust (X̄(0 : t : N)) = max
min
ηi j + λi jW̄ (0 : t : N) ,
i∈{1,...,L} j∈{1,...,mi }

(6)
and into a min-max canonical form
Jrobust (X̄(0 : t : N)) =

min

max

i∈{1,...,K} j∈{1,...,ni }


ζi j + γi jW̄ (0 : t : N) ,
(7)

for some integers K, L, n1 , . . . , nK , m1 , . . . , mL , where λi j and γi j
are weighting vectors and ηi j and ζi j are affine functions of
ū(0 : t : N) and x0 .
Proof: The proof is inductive on the structure of ψ and
uses the explicit form of the states in (2) utilizing the identities
− max( f1 , f2 ) = min(− f1 , − f2 ) and
min (max( f1 , f2 ), max(g1 , g2 )) =
max (min( f1 , g1 ), min( f1 , g2 ), min( f2 , g1 ), min( f2 , g2 )) .
for functions f1 , f2 , g1 , and g2 .
A. Arbitrary probability distributions with bounded support
Suppose the elements of the stochastic vector W (t), i.e.,
Wk (t), k ∈ {1, . . . , n} have arbitrary probability distribution
with known bounded support IWk (t) = [ak , bk ] and its first
moment E[Wk (t)] belongs to the interval MWk (t) = [ck , dk ], with
known quantities ak , bk , ck , dk ∈ R. Under this assumption, the
explicit form of X(·) in (2) implies that, for the observed value
of X(t) as x(t), X(τ) is a random variable with the following
interval of support and the first moment interval
IX(τ) = [āτ + C̄τ , b̄τ + C̄τ ], MX(τ) = [c̄τ + C̄τ , d¯τ + C̄τ ]

(8)

where C̄τ = Φ(τ,t)x(t) + ∑τ−1
k=t Φ(τ, k + 1)B(k)u(k), and
āτ , b̄τ , c̄τ and d¯τ are weighted sum of ak , bk , ck , dk , k ∈ N,
obtained by using interval arithmetics mentioned in Section
I-A.
The objective function in (5) can be written as
E [Jrobust (X̄(0 : t : N + 1))] + Jin (ū(0 : t − 1 : N))) and that
Jrobust (X̄(0 : t : N + 1)) = −ρ ψ (X̄(0 : t : N + 1), 0). Recall
that X̄(0 : t : N + 1) = [x(0), . . . , x(t), X(t + 1), . . . , X(N)] with
observed states x(0), . . . , x(t) of system (1) and random states
X(τ), τ > t. The following theorem shows how we can compute an upper bound for E[Jrobust (X̄(0 : t : N + 1))] based on
the canonical form of Jrobust .
Theorem 5: For a given STL formula ψ,
E [Jrobust (X̄(0 : t : N + 1))] can be upper bounded by
max

min

i∈{1,...,L} j∈{1,...,mi }

(dˆi j + ηi j ),

where the constants ηi j , i ∈ {1, . . . , L}, j ∈ {1, . . . , mi }, are
affine functions of ū(0 : t − 1 : N) and x(0), and dˆi j are
weighted sum of w(0), . . . , w(t −1) and ck , dk for k = t, . . . , N −
1.
Proof: With focus on the canonical form (6), let Yi j =
ηi j + λi jW̄ (0 : t : N). Considering the support and moment
interval of the components of W (τ), τ = t, . . . , N − 1, each
random variable Yi j has the following support and moment
interval (similar to (8))
IYi j = [âi j + ηi j , b̂i j + ηi j ], MYi j = [ĉi j + ηi j , dˆi j + ηi j ]

(9)

where the constants âi j , b̂i j , ĉi j , dˆi j , i ∈ {1, . . . , L}, j ∈
{1, . . . , mi }, are weighted sum of w(0), . . . , w(t − 1) and
ak , bk , ck , dk for k = t, . . . , N − 1, using interval arithmetic (cf.
Section I-A). Accordingly, Jrobust is a random variable with
the following support and moment intervals,
IJrobust =

[ max

min

(âi j +ηi j ), max

i∈{1,...,L} j∈{1,...,mi }

h
i
≤ E ((Z1 ) p + . . . + (Zs ) p )1/p
!1/p

(b̂i j +ηi j ]

min

i∈{1,...,L} j∈{1,...,mi }

s

MJrobust =
[ max

min

(ĉi j +ηi j ), max

i∈{1,...,L} j∈{1,...,mi }

min

(dˆi j +ηi j )].

Hence, as we are minimizing the cost function in (5), we can
utilize the upper bound maxi∈{1,...,L} min j∈{1,...,mi } (dˆi j + ηi j )
for E [Jrobust (X̄(0 : t : N + 1))].
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].
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 ∈ N, W (t) is normally distributed with
mean E[W (t)] = 0 and covariance matrix ΣW (t) , i.e., W (t) ∼
N (0, ΣW (t) ). The explicit form of X(τ) in (2) and the fact
that normal distribution is closed under affine transformations
result in normal distribution for X(τ), τ ∈ N. Its expected value
and covariance matrix with an observed value x(t) of X(t) are
τ−1

µτ = Φ(τ,t)x(t) + ∑ Φ(τ, k + 1)B(k)u(k) and
k=t

max

min

i∈{1,...,L} j∈{1,...,mi }


E

max Yi j ≤

i∈{1,...,K} j∈{1,...,ni }

min
i∈{1,...,K}

or

min

max Yi j

i∈{1,...,K} j∈{1,...,ni }

(11)

with Yi j = ηi j + λi jW̄ (0 : t − 1 : N) being affine functions of
the process noise, thus normally distributed random variables
(similar to X(τ) explained above). With focus on these canonical representations for Jrobust we employ Proposition 6 to show
how to approximate E [Jrobust ] using higher order moments of
W (t) ∼ N (0, Σ). 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 Zi for i ∈
{1, . . . , s} and let p be an even integer. Then
E [max(Z1 , . . . , Zs )] ≤ E [max(|Z1 |, . . . , |Zs |)]

j=1

Proof: For random variables Yi j , i ∈ {1, . . . , L}, j ∈
{1, . . . , mi }, and for a positive even integer p, the following
inequality holds,



 !1/p
max

min

i∈{1,...,L} j∈{1,...,mi }

p

L

(i)

E

Yi j ≤

min

∑E

j∈{1,...,mi }

i=1

Yi j
 p !1/p

L


E
∑ −

(ii)

=

i=1
L

(iii)

≤

mi

∑∑

max

j∈{1,...,mi }

−Yi j

!1/p
E[Yipj ]

,

i=1 j=1

where in (i) we used the upper bound obtained in Proposition 6; in (ii) we used the fact that mink∈{1,...,r} (αk ) =
− maxk∈{1,...,r} (−αk ); In (iii) we use again the inequality in
Proposition 6. Moreover, for i ∈ {1, . . . , K}, j ∈ {1, . . . , ni }, the
following inequality holds,



 (i)
max Yi j
E
min
max Yi j ≤ min E
j∈{1,...,ni }

i∈{1,...,K}

≤

Yi j

!
h i 1/p
p
. (13)
∑ E Yi j
ni


min

, (12)

i=1 j=1

k=t

min

∑ ∑ E[Yipj ]

Yi j ≤

(ii)

i∈{1,...,L} j∈{1,...,mi }

mi

L

E

∑ Φ(τ, k + 1)ΣW (k) Φ(τ, k + 1)T ,

respectively, for τ ≥ t ≥ 0.
In this section we use the canonical representation of
Jrobust (X̄(0 : t : N + 1)) in Theorem 4, which states that Jrobust
(for fixed ū(0 : t : N) and x0 ) can be written in either of the
forms
max

Founded on Proposition 6, next theorem shows how we can
upper bound E [Jrobust ] using the higher order moments of Yi j .
Theorem 7: Considering the canonical forms in (11) for
Jrobust as a function of random variables Yi j , E [Jrobust ] can be
upper bounded by
!1/p



i∈{1,...,K} j∈{1,...,ni }

τ−1

.

i=1

i∈{1,...,L} j∈{1,...,mi }

(10)

Στ =

∑ E [(Zi ) p ]

≤

!
h i 1/p
p
∑ E Yi j ,
ni

min
i∈{1,...,K}

j=1

where we apply Jensen’s inequality to the concave function
min(·) to get (i). The inequality of Proposition 6 gives (ii).
Note that random variables Yi j 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 µ and variance σ 2 ,
the p-th moment has a closed form as
E [Z p ] = σ p i−p H p (iµ/σ )

(14)

where i is the imaginary unit and
p/2

(−1)l z p−2l
l
l=0 2 l!(p − 2l)!

H p (z) = p! ∑

(15)

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).

⇔Pr [(ΞN ,t) |= ¬ϕ1 ∨ (ϕ1 ∧ ¬ϕ2 )] ≥ 1 − ϑ
⇔Pr [(ΞN ,t) |= ¬ϕ1 ]+Pr [(ΞN ,t) |= ϕ1 ∧ ¬ϕ2 ] ≥ 1 − ϑ , (18)
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 − ϑ )/2, we
have the inequalities

VI. U NDER A PPROXIMATION OF C HANCE C ONSTRAINTS
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 Ξ(0 : N) = X(0)X(1) . . . X(N) only
by ΞN without declaring its dependency on the state, input,
and disturbance variables. Recall the chance constraint (5c) as
Pr [(ΞN ,t) |= ϕ] ≥ 1 − δ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 ϕ 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 ϕ and a constant ϑ ∈ (0, 1),
constraints of the forms
Pr [(ΞN ,t) |= ϕ] ≥ ϑ and Pr [(ΞN ,t) |= ϕ] ≤ ϑ

(16)

can be transformed into similar constraints with ϕ being an
atomic predicate using the structure of ϕ.
Proof: The proof is inductive on the structure of the
formula ϕ as discussed in the following three cases.
Case I: ϕ = ¬ϕ1 we have the following equivalences

Pr [(ΞN ,t) |= ¬ϕ1 ] ≥ (1 − ϑ )/2,
1−ϑ
,
Pr [(ΞN ,t) |= ϕ1 ∧ ¬ϕ2 ] ≥
2

Pr [(ΞN ,t) |= ¬ϕ1 ] ≤ ϑ ⇔Pr [(ΞN ,t) 2 ϕ1 ] ≤ ϑ
⇔Pr [(ΞN ,t) |= ϕ1 ] ≥ 1 − ϑ .
Case II: ϕ = ϕ1 ∧ ϕ2 we obtain the following inequalities
by using the fact that for possibly joint events A and B,
it holds that Pr[A ∧ B] ≥ ϑ ⇔ Pr(¬A ∨ ¬B) ≤ 1 − ϑ and
Pr(A ∨ B) ≤ Pr[A ] + Pr[B].
Pr [(ΞN ,t) |= ϕ1 ∧ ϕ2 ] ≥ ϑ
⇔Pr [(ΞN ,t) |= ϕ1 ∧ (ΞN ,t) |= ϕ2 ] ≥ ϑ
⇔Pr [(ΞN ,t) 2 ϕ1 ∨ (ΞN ,t) 2 ϕ2 ] ≤ 1 − ϑ
⇐Pr [(ΞN ,t) 2 ϕ1 ] + Pr [(ΞN ,t) 2 ϕ2 ] ≤ 1 − ϑ
1−ϑ
i = 1, 2.
(17)
⇐Pr [(ΞN ,t) 2 ϕi ] ≤
2
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 − ϑ )/2. However, this can be replaced by any two other
probabilities δ1 and δ2 such that δ1 +δ2 = 1−ϑ . Now consider
the second possibility:
Pr [(ΞN ,t) |= ϕ1 ∧ ϕ2 ] ≤ ϑ
⇔Pr [(ΞN ,t) |= ¬ϕ1 ∨ ¬ϕ2 ] ≥ 1 − ϑ

(20)

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 −
ϑ )/2 with (1 − ϑ )/n.
Case III: ϕ = ϕ1 U[a,b] ϕ2 The satisfaction (ΞN ,t) |= ϕ1 U[a,b] ϕ2
W
is equivalent to t+b
j=t+a ψ j with disjoint events
t+a−1
^

j−1
^

i=t

i=a+t

(ΞN , i) |= ϕ1

ψj =

(ΞN , i) |= (ϕ1 ∧ ¬ϕ2 ) ∧ (ΞN , j) |= ϕ2 .



Thus Pr (ΞN ,t) |= ϕ1 U[a,b] ϕ2 ≥ ϑ is equivalent to
∑t+b
j=t+a Pr[ψ j ] ≥ ϑ . Assuming the probabilities of events
are lower bounded by the same values, we have
Pr[ψ j ] ≥ ϑ /(b − a + 1) for j = a + t, . . . , b + t, which
again can be reduced as in Case II.
The second possible probabilistic constraint in Case III can
be obtained as
"
#
b+t
_


Pr (ΞN ,t) |= ϕ1 U[a,b] ϕ2 ≤ ϑ ⇔ Pr
ψj ≤ ϑ
j=a+t

Pr [(ΞN ,t) |= ¬ϕ1 ] ≥ ϑ ⇔Pr [(ΞN ,t) 2 ϕ1 ] ≥ ϑ
⇔Pr [(ΞN ,t) |= ϕ1 ] ≤ 1 − ϑ ,

(19)

t+b

⇔

∑

Pr[ψ j ] ≤ ϑ

j=t+a

⇔ Pr[ψ j ] ≥ ϑ /(b − a + 1),
(21)
which can be again reduced as in Case II. Here also, we used
the fact that ψ j consists of disjoint events and we assume
that he probabilities of events are lower bounded by the same
value, i.e., by ϑ /(b − a + 1), for j = a + t, . . . , 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.
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 Z1 , . . . , Zn be random
variables with interval of bounded support [ai , bi ] and let
E[Z1 ], . . . , E[Zn ] denote their expected values belonging to the
moment intervals Mi for i = 1, . . . , n. Define Z = ∑ni=1 Zi and

E(Z) = ∑ni=1 E[Zi ]. Using Chernoff-Hoeffding inequality, the
following upper bound exists for any ς ≥ 0 [16]


−ς 2
.
(22)
Pr [Z − E[Z] ≤ −ς ] ≤ exp
ν ∑ni=1 (bi − ai )2
where ν > 0 is a constant. If Z1 , . . . , Zn are dependent, then
the inequality applies with a constant ν = χ(Ĝ)/2, where Ĝ
denotes the indirected dependency graph of Z1 , . . . , Zn and
χ(Ĝ) is the chromatic number of the graph Ĝ defined as
the minimum number of colors required to color Ĝ. For the
independent case, χ(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 ϕ = {α ≥ 0}.
Since α is an affine function of random state variables, it is a
random variable itself with the following interval of support
and moment interval
Iα(X(t)) = [ãt + C̃t , b̃t + C̃t ]
Mα(X(t)) = [c̃t + C̃t , d˜t + C̃t ]

(23)

where for t = 0, . . . , N, we have ãt , b̃t , c̃t and d˜t are weighted
sum of āt , b̄t , c̄t , d¯t related to the interval of support and
moment interval of random variables X(t) (cf. (8)), and C̃t
is a linear expression of input variables.
Let ς = E [α(X(t))]; we can directly use (22) as


⇐exp

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 µ and
standard deviation σ ,
Pr[x ≤ b] ≤ δt ⇔ F −1 (δt ) ≥ b ⇔ µ + σ φ −1 (δt ) ≥ b
Pr[x ≤ b] ≥ δt ⇔ F

−1

(δt ) ≤ b ⇔ µ + σ φ

Pr [(ΞN ,t) |= π α ] ≥ 1 − δt ⇔Pr [α(X(t)) > 0] ≥ 1 − δt

Pr [(ΞN ,t) |= π ] ≤ 1 − δt ⇔Pr [α(X(t)) > 0] ≤ 1 − δt
⇔Pr [α(X(t)) ≤ 0] ≥ δt
⇔F −1 (δt ) ≤ 0
⇔µt + σt φ −1 (δt ) ≤ 0.

(24)

Note that since δt ∈ (0, 1), we have log(δt ) < 0; hence, by
multiplying both sides of the inequality by -1 in line 5 of (24),
N
the expression − log(δt ) · ∑t=1
(b̃t − 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 ς ≥ 0. Hence,
we can replace ς in the last inequality of (24) by the lower
bound of its moment interval in (23), i.e., with c̃t + C̃t , which
is a linear expression in the input variables.
Consequently, in this case, the chance constraint in (5) can
be replaced by
s
N

(25)

t=1

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

⇐Pr [α(X(t)) ≥ 0] ≤ 1 − δt

(29)

α

t=1

Pr [(ΞN ,t) |= π ] ≤ 1 − δt

(28)

where
denotes the inverse of the cumulative distribution
function or the quantile function and φ −1 is the inverse of the
error function of a normally distributed random variable.
Recall the chance constraints (16) with ϕ = {α ≥ 0}. Since
α is an affine function of normally distributed state variables,
it is also normally distributed with appropriately defined mean
µt and variance σt2 . Hence, we can directly use (27) and (28)
as

⇔µt + σt φ −1 (δt ) ≥ 0,

N

α

(δt ) ≤ b

(27)

⇔F −1 (δt ) ≥ 0

t=1

−ν log(δt ) · ∑ (b̃t − ãt )2 .

−1

F −1

⇔ − ς 2 ≤ ν log(δt ) ∑ (b̃t − ãt )2

c̃t + C̃t ≥

(26)

B. Normal distribution

N

−ν log(δt ) ∑ (b̃t − ãt )2

≤ 1 − δt ,

⇔Pr [α(X(t)) ≤ 0] ≤ δt

⇔Pr [α(X(t)) ≤ 0] ≤ δt


−ς 2
≤ δt
⇐exp
N
ν ∑t=1
(b̃t − ãt )2
−ς 2
⇔
≤ log(δt )
N
ν ∑t=1 (b̃t − ãt )2

⇐ς ≥



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 δt by 1 − δt in the related expressions.

Pr [(ΞN ,t) |= π α ] ≥1 − δt ⇔ Pr [α(X(t)) > 0] ≥ 1 − δt

s

−ς 2
N
ν ∑t=1 (b̃t − ãt )2

(30)

Therefore, the chance constraint can be replaced by the
equivalent linear constraint (29) or (30), depending on the type
of constraint we have.
VII. E XPERIMENTAL R ESULTS
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 Tr 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 Tw j , respectively. Dynamics of the temperature of wall w j
and room r can be written as [15]
Cwj

Tr,k − Tw j
dTw j
= ∑
+ r j α j Aw j Qrad j
dt
R j,k
k∈Nw

(31)

j

dTr
Crj
dt

=

Tk − Tr
+ ṁr ca (Ts − Tr )
Rk
k∈Nr

∑

+ wi τw Awin Qrad + Q̇int

(32)

where Cwj , α j and Aw 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 jk . The radiative heat flux density on w j is denoted by
Qrad j , the set of all neighboring nodes to w j is denoted by
Nw 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 Tr ,Crj and ṁr , respectively; ca is the specific
heat capacity of air, and Ts 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, τw is the transmissivity of the
glass of the window in room r, Awin is the total area of the
windows on walls surrounding room r, Qrad is the radiative
heat flux density per unit area radiated to room r, and Q̇intr is
the internal heat generation in room r. Finally, Nr 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
X(t + 1) = AX(t) + Bu(t) +W (t),
where X ∈ Rn is the state vector representing the temperature
of the walls and the rooms and u ∈ Rm 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 ts = 30
minutes. The disturbance W (·) 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 (·) 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
wr (t), and the reference is perturbed by independent and
identically distributed random vectors e(t) ∼ N (0, In ), i.e.,
the disturbance is W (t) = wr (t) + e(t), which is normally
distributed with mean µt = wr (t) and identity covariance
matrix Σ = In .
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 occ : N → {−1, 1} representing the
room occupancy; occ(t) = 1 if the room is occupied at time
t and 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 δ = 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
Tr when the room is occupied; thus the specification is

ψ = 2[0,N] occ(t) = 1 → X(t) > Tr .
At each time instant 0 ≤ t < N, the optimization problem (5) obtains an optimal control input ũopt (t : N) =
[uopt (t), . . . , uopt (N − 1)] that minimizes
N−1

E[−ρ ψ (X̄(0 : t : N), 0)] +

∑ ||u(k)||1 ,
k=t

where the robustness function is defined as
ρ ψ (X̄(0 : t : N), 0) = min{X(τ) − Tr | τ ∈ [0, N], occ(τ) > 0}.
The chance constraint (5c) is defined with the same specification ϕ = ψ. We approximate E[−ρ ψ (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 ≤ t < N,
we have 0 ≤ u(t) ≤ 380.
The simulations are done using M ATLAB R2014b on a
2.6 GHz Intel Core i5 processor and the optimizations are
solved using fmincon solver in M ATLAB. we perform ns = 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 Tr in Fahrenheit. The controller ensures that the
room temperature goes above the reference temperature Tr
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
minimum, and maximum of the air flow
i
h 3average,
ft
, which indicates that the input constraint is not
rate in min
violated.
Note that all these ns = 200 runs result in feasible solutions,
which gives a confidence bound on the feasibility of the
original problem as follows. Since all the ns runs of the
simulation are feasible, we can claim that the original problem
is also feasible with probability at least (β /2)1/ns , β ∈ (0, 1),
with confidence level 1−β [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σ interval such
that e(t) ∈ [−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.

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.
Fig. 1. Controlling the room temperature using SHMPC in the presence of
normally distributed disturbance and STL constraints.
TABLE I
C OMPARISON OF THE STATISTICS OF THE FAN ENERGY CONSUMPTION
USING DIFFERENT CONTROL APPROACHES .
Computational
Methods
Open-loop OC
RMPC
SHMPC

Fan energy
consumption [kWh]
1337.016
µ1 = 12.2216, σ1 = 0.045µ1
µ2 = 2.5101, σ2 = 0.104µ2

Average
computation time [s]
3.9277
33.4891
19.3622

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 openloop 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. C ONCLUSIONS
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

R EFERENCES
[1] M. A. Abramowitz and I. Stegun. Handbook of Mathematical Functions. National Bureau of Standards, US Government Printing Office,
Washington DC, 1964.
[2] A. Atamtürk and M. W. P. Savelsbergh. Integer-programming software
systems. Annals of Operations Research, 140(1):67–124, November
2005.
[3] A. Bemporad, W. P. M. H. Heemels, and B. De Schutter. On hybrid
systems and closed-loop MPC systems. IEEE Transactions on Automatic
Control, 47(5):863–869, May 2002.
[4] O. Bouissou, E. Goubault, S. Putot, A. Chakarov, and S. Sankaranarayanan. Uncertainty propagation using probabilistic affine forms
and concentration of measure inequalities. In Proceedings of Tools and
Algorithms for Construction and Analysis of Systems (TACAS), volume
TBA of Lecture Notes in Computer Science, page TBA. Springer-Verlag,
2016.
[5] B. Büeler, A. Enge, and K. Fukuda. Exact volume computation for
convex polytopes: A practical study. In G. Kalai and G.M. Ziegler,
editors, Polytopes – Combinatorics and Computation, pages 131–154.
Birkäuser Verlag, Basel, Switzerland, 2000.
[6] C. Clopper and E. S. Pearson. The use of confidence or fiducial limits
illustrated in the case of the binomial. Biometrika, 26(4):404–413, 1934.
[7] P. J. Davis and P. Rabinowitz. Methods of Numerical Integration.
Academic Press, New York, 2nd edition, 1984.
[8] B. De Schutter and T. van den Boom. Model predictive control for
max-plus-linear discrete event systems. Automatica, 37(7):1049–1056,
July 2001.
[9] B. De Schutter and T. J. J. van den Boom. MPC for continuous
piecewise-affine systems. Systems & Control Letters, 52(3–4):179–192,
July 2004.
[10] G. E. Fainekos and G. J. Pappas. Robustness of temporal logic specifications for continuous-time signals. Theoretical Computer Science,
410(42):4262–4291, 2009.
[11] S. S. Farahani, R. Majumdar, V.S. Prabhu, and S. Esmaeil Zadeh
Soudjani. Shrinking horizon model predictive control with chanceconstrained signal temporal logic specifications. In To appear in the
proceedings of American Control Conference (ACC) 2017.
[12] S. S. Farahani, V. Raman, and R. M. Murray. Robust model predictive
control for signal temporal logic synthesis. In Proceedings of the IFAC
Conference on Analysis and Design of Hybrid Systems, Atlanta, Georgia,
October 2015.
[13] S. S. Farahani, T. van den Boom, H. van der Weide, and B. De Schutter.
An approximation method for computing the expected value of maxaffine expressions. European Journal of Control, 27:17–27, January
2016.
[14] J. M. Grosso, C. Ocampo-Martinez, V. Puig, and B. Joseph. Chanceconstrained model predictive control for drinking waternetworks. Journal of process control, 24:504–516, 2014.
[15] M. Maasoumy Haghighi. Controlling Energy Efficient Buildings in the
Context of Smart Grid: A Cyber Physical System Approach. PhD thesis,
University of California, Berkeley, Berkeley, California, December 2013.

[16] S. Janson. Large deviations for sums of partly dependent random
variables. Random Structures Algorithms, 24(3):234–248, 2004.
[17] S. Jha and V. Raman. Automated synthesis of safe autonomous vehicle
control under perception uncertainty. In 8th International Symposium
of NASA Formal Methods, Lecture Notes in Computer Science, pages
117–132, Minneapolis, Minnesota, June 2016.
[18] X. Jin, A. Donzé, J. V. Deshmukh, and S. A. Seshia. Mining
requirements from closed-loop control models. In Hybrid Systems:
Computation and Control, HSCC 2013, April 8-11, 2013, Philadelphia,
PA, USA, pages 43–52, 2013.
[19] J. B. Lasserre. Integration on a convex polytope. Proceedings of the
American Mathematical Society, 126(8):2433–2441, August 1998.
[20] M. Lazar, M. Heemels, S. Weiland, and A. Bemporad. Stability of hybrid
model predictive control. IEEE Transactions on Automatic Control,
15(11):1813–1818, 2006.
[21] J. T. Linderoth and T. K. Ralphs. Noncommercial software for mixedinteger linear programming. In J. Karlof, editor, Reinforcement Learning: State-Of-The-Art, CRC Press Operations Research Series, pages
253–303. 2005.
[22] J. M. Maciejowski. Predictive Control with Constraints. Prentice Hall,
Harlow, England, 2002.
[23] O. Maler and D. Nickovic. Monitoring temporal properties of continuous
signals. In In FORMATS/FTRTFT, pages 152–166, 2004.
[24] V. Raman, A. Donzé, D. Sadigh, R. M. Murray, and S. A. Seshia.
Reactive synthesis from signal temporal logic specifications. In Hybrid
Systems: Computation and Control, HSCC 2015, Seattle, WA, USA, April
14-16, 2015, pages 239–248, 2015.
[25] D. Sadigh and A. Kapoor. Safe control under uncertainty with probabilistic signal temporal logic. In Proceedings of Robotics: Science and
Systems Conference, AnnArbor, Michigan, June 2016.
[26] A. T. Schwarm and M. Nikolaou. Chance-constrained model predictive
control. AIChE Journal, 45(8):1743–1752, 1999.

Samira S. Farahani is a postdoctoral researcher at
Max Planck Institute for Software Systems. Her research interests are in control synthesis for stochastic
hybrid systems, temporal logic and formal methods,
and discrete-event systems. Prior to this position, she
was a postdoctoral researcher at CMS department,
California Institute of Technology, and at Energy
and Industry Department, Delft University of Technology. Dr. Farahani has received her PhD degree
in Systems and Control and her MSc. degree in
Applied Mathematics, both from Delft University of
Technology, the Netherlands, in 2012 and 2008, respectively. She obtained her
BSc. degree in Applied Mathematics from Sharif University of Technology,
Iran, in 2005.

Rupak Majumdar is a Scientific Director at the
Max Planck Institute for Software Systems. His research interests are in the verification and control of
reactive, real-time, hybrid, and probabilistic systems,
software verification and programming languages,
logic, and automata theory. Dr. Majumdar received
the President’s Gold Medal from IIT, Kanpur, the
Leon O. Chua award from UC Berkeley, an NSF
CAREER award, a Sloan Foundation Fellowship,
an ERC Synergy award, ”Most Influential Paper”
awards from PLDI and POPL, and several best paper
awards (including from SIGBED, EAPLS, and SIGDA). He received the
B.Tech. degree in Computer Science from the Indian Institute of Technology
at Kanpur and the Ph.D. degree in Computer Science from the University of
California at Berkeley.

Vinayak Prabhu did his PhD in the department of
Electrical Engineering and Computer Sciences at the
University of California, Berkeley (USA) in 2008,
and is currently a postdoctoral researcher at Max
Planck Institute for Software Systems. His research
interests are in modelling, analysis, verification, and
control, of reactive, real-time, hybrid, and multiagent systems.

Sadegh Esmaeil Zadeh Soudjani is a postdoctoral
researcher at the Max Planck Institute for Software
Systems, Germany. His research interests are formal synthesis, abstraction, and verification of complex dynamical systems with application in CyberPhysical Systems, particularly involving power and
energy networks, smart grids, and transportation systems. Sadegh received the B.Sc. degrees in electrical
engineering and pure mathematics, and the M.Sc.
degree in control engineering from the University of
Tehran, Tehran, Iran, in 2007 and 2009, respectively.
He received the Ph.D. degree in Systems and Control in November 2014
from the Delft Center for Systems and Control at the Delft University of
Technology, Delft, the Netherlands. His PhD thesis received the best thesis
award from the Dutch Institute of Systems and Control. Before joining
Max Planck Institute, he was a postdoctoral researcher at the department of
Computer Science, University of Oxford, United Kingdom.

