1

Model Predictive Control for Signal Temporal
Logic Specifications

arXiv:1703.09563v1 [] 28 Mar 2017

Vasumathi Raman1 , Alexandre Donzé2 , Mehdi Maasoumy3 ,
Richard M. Murray1 , Alberto Sangiovanni-Vincentelli2 and Sanjit A. Seshia2

Abstract—We present a mathematical programming-based
method for model predictive control of cyber-physical systems
subject to signal temporal logic (STL) specifications. We describe
the use of STL to specify a wide range of properties of these
systems, including safety, response and bounded liveness. For
synthesis, we encode STL specifications as mixed integer-linear
constraints on the system variables in the optimization problem
at each step of a receding horizon control framework. We prove
correctness of our algorithms, and present experimental results
for controller synthesis for building energy and climate control.
Index Terms—formal synthesis, timed logics, model predictive
control, cyberphysical systems

I. I NTRODUCTION
Controlling a cyber-physical system (CPS) involves handling complex interactions between computing components
and their physical environment, and often necessitates hierarchies of controllers. Typically at the highest level, a supervisory controller is responsible for making high-level decisions,
while at the lowest level traditional control laws such as PID
control are used. In general, the design of these different
controllers is done mostly in isolation at each level, and
their combination is implemented ad hoc. As the complexity
of these systems grows, reasoning about the correctness of
interactions between the various layers of control becomes
increasingly challenging, begging automation.
Formal methods is the subfield of computer science concerned with verification and synthesis, i.e., automatic and
rigorous design of digital systems. It provides mathematical
formalisms for specifying behaviors and algorithms for verification and synthesis of a system against properties specified
within these formalisms. Methods for synthesis of correctby-construction discrete supervisory controllers have been
developed and successfully used for cyber-physical systems in
domains including robotics [1] and aircraft power system design [2]. However, for physical systems that require constraints
This work was supported in part by TerraSwarm, one of six centers
of STARnet, a Semiconductor Research Corporation program sponsored by
MARCO and DARPA.
1 V. Raman and R. M. Murray are with the California Institute
of Technology, Pasadena, CA 91125, USA vasu@caltech.edu,

murray@cds.caltech.edu
2 A.
Donzé, A. Sangiovanni-Vincentelli and S. A. Seshia
are
with
the
Department
of
Electrical
Engineering
and
Computer Science, UC Berkeley, Berkeley, CA 94720, USA

donze@berkeley.edu, alberto@berkeley.edu,
sseshia@eecs.berkeley.edu
3 M.

Maasoumy is with C3 Energy Inc. Redwood City, CA 94063, USA

mehdi.maasoumy@c3energy.com
Manuscript received January 14, 2016.

not just on the order of events, but on the temporal distance
between them, simulation and testing is still the method of
choice for validating properties and establishing guarantees;
the exact exhaustive verification or synthesis of such systems
is in general undecidable [3].
Model predictive control (MPC) or receding horizon control
(RHC) is based on iterative, finite horizon optimization over
a model of the plant, i.e. the system to be controlled. At
any given time t, the current plant state is observed, and
an optimal control strategy computed for some finite time
horizon in the future, [t, t + H]. An online calculation is
performed to explore trajectories originating from the current
state, and an optimal control strategy computed up to time
t + H. To provide robustness with respect to modelling errors,
only the first step of the computed optimal control strategy is
implemented. The plant state is then sampled again, and new
calculations are performed on a horizon of H starting from
the new current state. While the global optimality of such a
receding horizon approach is not ensured, it tends to do well
in practice: in addition to reducing computational complexity,
it improves the system robustness with respect to exogenous
disturbances and modeling uncertainties [4]. Another reason
MPC is particularly attractive to industry is its ability to handle
constrained dynamical systems [5].
Signal temporal logic (STL) [6] was originally developed in
order to specify and monitor the expected behavior of physical
systems, including temporal constraints between events. STL
allows the specification of properties of dense-time, realvalued signals, and the automatic generation of monitors for
testing these properties on individual simulation traces. It
has since been applied to the analysis of several types of
continuous and hybrid systems, including dynamical systems
and analog circuits, where the continuous variables represent
quantities like currents and voltages in a circuit. STL has
the advantage of naturally admitting quantitative semantics
which, in addition to the yes/no answer to the satisfaction
question, provide a real number that grades the quality of the
satisfaction or violation. Such semantics have been defined for
timed logics, including metric temporal logic (MTL) [7] and
STL [8], to assess the robustness of systems to parameter or
timing variations.
In this paper, we solve the problem of control synthesis from
STL specifications, using a receding horizon approach. We
allow the user to specify desired properties of the system using
an STL formula, and synthesize control such that the system
satisfies that specification, while using a receding horizon
approach to ensure practicality and robustness. We do so by

2

where for all k > 0, tk+1 −tk = ∆t. A run of Σd is a sequence
ξ = (x0 u0 w0 )(x1 u1 w1 )(x2 u2 w2 )...
where xk = x(tk ) ∈ X is the state of the system at
index k, and for each k ∈ N, uk = u(tk ) ∈ U , wk =
w(tk ) ∈ W and xk+1 = fd (xk , uk , wk ). We assume that
given an initial state x0 ∈ X, a control input sequence uN =
u0 u1 u2 . . . uN −1 ∈ U N , and a sequence of environment inputs
wN = w0 w1 w2 . . . wN −1 ∈ W N , the resulting horizon-N run
of a system modeled by equation (1), which we denote by
ξ(x0 , uN , wN ) = (x0 u0 w0 )(x1 u1 w1 )(x2 u2 w2 )...(xN uN wN ),
is unique. In addition, we introduce a generic cost function
J(ξ(x0 , u, w)) that maps (infinite and finite) runs to R.
B. Signal Temporal Logic

ϕ ::= π µ | ¬ψ | ϕ1 ∧ ϕ2 |

We consider a continuous-time system Σ of the form
ẋt = f (xt , ut , wt )
where xt ∈ X ⊆ (Rnc × {0, 1}nl ) are the continuous and
binary/logical states, ut ∈ U ⊆ (Rmc ×{0, 1}ml ) are the (continuous and logical) control inputs, wt ∈ W ⊆ (Rec ×{0, 1}el )
are the external environment inputs (also referred to as “disturbances”). As xt includes binary components, the above
ODE may contain discontinuities corresponding to switches in
these components. In general, such hybrid systems are more
accurately modeled using differential-algebraic equations, but
we do not dwell on this point since we approximate their
behavior with a difference equation, as follows.
Given a sampling time ∆t > 0, we assume that Σ admits a
discrete-time approximation Σd of the form
x(tk+1 ) = fd (x(tk ), u(tk ), w(tk ))

(1)

ϕ | ϕ1 U [a,b] ϕ2 ,



ξ |= ϕ
(ξ, tk ) |= π µ
(ξ, tk ) |= ¬ψ
(ξ, tk ) |= ϕ ∧ ψ
(ξ, tk ) |= [a,b] ϕ
(ξ, tk ) |= ϕ U [a,b] ψ

⇔
⇔
⇔
⇔
⇔
⇔

(ξ, t0 ) |= ϕ
µ(xk , yk , uk , wk ) > 0
¬(ξ, tk ) |= ψ)
(ξ, tk ) |= ϕ ∧ (ξ, tk ) |= ψ
∃tk0 ∈ [tk +a, tk +b], (ξ, tk0 ) |= ϕ
∃tk0 ∈ [tk +a, tk +b] s.t. (ξ, tk0 ) |= ψ
∧∀tk00 ∈ [tk , tk0 ], (ξ, tk00 ) |= ϕ.

An STL formula ϕ is bounded-time if it contains no
unbounded operators; the bound of ϕ is the maximum over the
sums of all nested upper bounds on the temporal operators, and
provides a conservative maximum trajectory length required
to decide its satisfiability. For example, for [0,10] [1,6] ϕ, a
trajectory of length N such that tN ≥ 10+6 = 16 is sufficient
to determine whether the formula is satisfiable.


A. Systems

[a,b]

where π µ is an atomic predicate X × U × W → B whose truth
value is determined by the sign of a function µ : X ×U ×W →
R and ϕ1 , ϕ2 are STL formulas. The fact that a run ξ(x0 , u, w)
satisfies an STL formula ϕ is denoted by ξ |= ϕ. Informally,
ξ |= [a,b] ϕ if ϕ holds at some time step between a and b, and
ξ |= ϕ U [a,b] ψ if ϕ holds at every time step before ψ holds,
and ψ holds at some time step between a and b. Additionally,
we define [a,b] ϕ = ¬ [a,b] (¬ϕ), so that ξ |= [a,b] ϕ if ϕ
holds at all times between a and b. Formally, the validity of a
formula ϕ with respect to the run ξ is defined inductively as
follows



II. P RELIMINARIES



We consider STL formulas defined recursively according to
the following grammar:



decomposing the STL specifications into a series of formulas
over each time horizon, such that synthesizing a controller fulfilling the formula at each horizon results in satisfaction of the
global specification. Recent work on optimal control synthesis
of aircraft load management systems [9] represented STLlike specifications as time-dependent equality and inequality
constraints, yielding a mixed integer linear program (MILP).
The MILP was then solved in an MPC framework, yielding
an optimal control policy. However, the manual transformation
of specifications into equality and inequality constraints is
cumbersome and problem-specific. As a key contribution, this
paper presents two automatically-generated MILP encodings
for such STL specifications.
Our main contribution is a pair of bounded model checkingstyle encodings [10] for STL specifications as MILP constraints on a cyber-physical system. We show how these encodings can be used to generate open-loop control signals that
satisfy finite and infinite horizon STL properties and, moreover, to generate signals that maximize quantitative (robust)
satisfaction. We provide a fragment of STL, denoted SNNSTL, such that, under reasonable assumptions on the system
dynamics, the problem of synthesizing an open-loop control
sequence such that the system satisfies a provided specification is a linear program (LP), and therefore polynomial-time
solvable. We also demonstrate how our MILP formulation of
the STL synthesis problem can be used in an MPC framework
to compute feasible and optimal controllers for cyber-physical
systems under timed specifications. We present experimental
results comparing both encodings, and two case studies: one on
a thermal model of an heating, ventilation and air-conditioning
(HVAC) system, and another in the context of regulation
services in a micro-grid. These case studies were previously
reported in [11], [12]. We show how the MPC schemes in
these examples can be framed in terms of synthesis from an
STL specification, and present simulation results to illustrate
the effectiveness of our methodology.

Remark 1. Here we have defined a semantics for STL over
discrete-time signals, which is formally equivalent to the
simpler linear temporal logic (LTL), once time and predicates
are abstracted into steps and Boolean variables, respectively.
There are several advantages of still using STL over LTL,
though. First, STL allows us to explicitly use real time in
our specifications instead of abstract integer indices, which
improves the readability relative to the original system’s behaviors. Second, although in the rest of this paper we focus
on the control of the discrete-time system Σd , our goal is to
use the resulting controller for the control of the continuous

3

system Σ. Hence the specifications should be independent
from the sampling time ∆t. Finally, note that the relationship
between the continuous-time and discrete-time semantics of
STL, depending on discretization error and sampling time, is
beyond the scope of this paper. The interested reader can refer
to [13] for further discussion on this topic.

Note that we assume that the state of the plant is fully
observable, and the environment inputs are known in advance.
Problem 2 (closed-loop). Given a horizon 0 < L < N , for
all 0 ≤ k ≤ N − L, compute u∗k = uL∗
k , the first element of
L∗ L∗
L∗
the sequence uL∗
=
u
u
.
.
.
u
k
k
k+1
k+L−1 satisfying
uL∗
k =

C. Quantitative semantics for STL

argmin J(xk , uL
k , wk , ϕ)
L
uL
K ∈U

s.t. ξ(xk , uL
k , wk ) |= ϕ

Quantitative or robust semantics for STL are defined by
providing a real-valued function ρϕ of signal ξ and time t
such that ρϕ (ξ, t) > 0 ⇒ (ξ, t) |= ϕ. We define one such
function recursively, as follows:



The closed-loop formulation corresponds to a model predictive control scheme, where the reference disturbance can
change at each iteration k.
πµ
In Sections IV and V, we present both an open-loop solution
ρ (ξ, tk )
= µ(xk , yk , uk , wk )
to Problem 1, and a solution to Problem 2 for a large class
ρ¬ψ (ξ, tk )
= −ρψ (ξ, tk )
of STL formulas. In the absence of an objective function
ρϕ1 ∧ϕ2 (ξ, tk )
= min(ρϕ1 (ξ, tk ), ρϕ2 (ξ, tk ))
ϕ1 ∨ϕ2
J on runs of the system, we maximize the robustness of
(ξ, tk )
= max(ρϕ1 (ξ, tk ), ρϕ2 (ξ, tk ))
ρ
the generated runs with respect to ϕ. A key component of
ψ
ψ
[a,b]
(ξ, tk ) = maxtk0 ∈[t+a,t+b] ρ (ξ, tk0 )
ρ
our solution is encoding the STL specifications as MILP
ϕ1 U [a,b] ϕ2
ϕ2
ρ
(ξ, tk ) = maxtk0 ∈[t+a,t+b] (min(ρ (ξ, tk0 ),
constraints, which can be combined with MILP constraints
ϕ1
mintk00 ∈[tk ,tk0 ] ρ (ξ, tk00 ))
representing the system dynamics to efficiently solve the
µ
Note that if µ(xk , yk , uk , wk ) = 0, neither ρπ (ξ, tk ) > 0 resulting state-constrained optimization problem.
µ
nor ρπ (ξ, tk ) > 0. Therefore, (ξ, t) |= ϕ 6⇒ ρϕ (ξ, t) > 0. To
µ
simplify notation, we denote ρπ by ρµ for the remainder of
IV. O PEN - LOOP C ONTROLLER S YNTHESIS
the paper. The robustness of satisfaction for an arbitrary STL
formula is computed recursively from the above semantics by
To solve Problem 1, we extend the bounded model checking
propagating the values of the functions associated with each
encoding of [14] from finite, discrete systems to dynamical
operand using min and max operators corresponding to the
systems using mixed-integer programming instead of SAT. Our
various STL operators. For example, the robust satisfaction of
presentation and notation below follow that of [15]. For openπ µ1 where µ1 (x) = x − 3 > 0 at time 0 is ρµ1 (ξ, 0) = x0 − 3.
loop controller synthesis, we will search for a trajectory of
The robust satisfaction of µ1 ∧ µ2 is the minimum of ρµ1
length N that satisfies ϕ. To admit STL formulas describing
and ρµ2 . Temporal operators are treated as conjunctions and
infinite runs, we parametrize an infinite sequence of states
disjunctions along the time axis: since we deal with discrete
using a finite sequence with a loop. Imposing this lassotime, the robustness of satisfaction of ϕ = [0,2.1] µ1 is
shaped structure renders our synthesis procedure conservative
ϕ
µ1
for infinite-state systems, in the sense that a solution may exist
ρ (x, 0) = mintk ∈[0,2.1] ρ (x, tk )
that is not found when imposing such a structure. For finite= min{x0 − 3, x1 − 3, . . . , xK − 3},
state systems, the lasso shape is without loss of generality but
where 0 ≤ t0 < t1 < . . . < tK ≤ 2.1 < tK+1 .
we must still find an appropriate trajectory length N .
The robustness score ρϕ (ξ, t) can be interpreted as how
much ξ satisfies ϕ. Its absolute value can be viewed as the
ϕ, Σd , J
signed distance of ξ from the set of trajectories satisfying or
violating ϕ, in the space of projections with respect to the
function µ that define the predicates of ϕ [7].
III. P ROBLEM S TATEMENT
We now formally state the STL control synthesis problem
and its model predictive control formulation. Given an STL
formula ϕ, a cost function of the form J(x0 , u, w, ϕ) ∈ R, an
initial state x0 ∈ X , a horizon L and a reference disturbance
signal w ∈ W N , we formulate two problems: open-loop and
closed-loop synthesis. The two scenarios are depicted as block
diagrams in Fig. 1 and Fig. 2.

Synthesis

x0 ,
w0 , . . . , wN −1

u∗ =
u∗0 u∗1 . . . u∗N −1

Plant Σ

x, y

x0 , w

Fig. 1. Open-loop problem formulation: given an STL formula φ, the discretetime plant model Σd , the cost function J, the goal is to generate a sequence of
control inputs u∗ over a horizon of N time steps. The problem is additionally
parametrized by the the initial state x0 and disturbance vector w.

Problem 1 (open-loop). Compute u∗ = u∗0 u∗1 . . . u∗N −1 where
u∗ =

argmin J(x0 , u, w, ϕ)
u∈U N

s.t. ξ(x0 , u, w) |= ϕ

The encoding of Problem 1 as an MILP consists of system
constraints, loop constraints and STL constraints, as defined
below.

4

C. Boolean encoding of STL constraints
Given a formula ϕ, we introduce a variable ztϕ , whose value
is tied to a set of mixed integer linear constraints required for
the satisfaction of ϕ at position t in the state sequence of
horizon N . In other words, ztϕ has an associated set of MILP
constraints such that ztϕ = 1 if and only if ϕ holds at position
t. We recursively generate the MILP constraints corresponding
to z0ϕ – the value of this variable determines whether a formula
ϕ holds in the initial state.
1) Predicates: The predicates are represented by constraints on system state variables. For each predicate µ ∈ P ,
we introduce binary variables ztµ ∈ {0, 1} for times t =
0, 1, ..., N . The following constraints enforce that ztµ = 1 if
and only if µ(xt ) > 0:
µ(xt )
≤
−µ(xt ) ≤

Mt ztµ − t
Mt (1 − ztµ ) − t

ϕi
Conjunction: ztψ = ∧m
i=1 zti

ztψ ≤ ztϕii , i = 1, ..., m,
Pm
ztψ ≥ 1 − m + i=1 ztϕii
ztψ ≥ ztϕii , i = 1, ..., m,
Pm
ztψ ≤ i=1 ztϕii

ϕi
Disjunction: ztψ = ∨m
i=1 zti

Given a formula ψ containing a Boolean operation, we add
new continuous variables ztψ ∈ [0, 1], and set ztψ = ¬ztµ ,
ϕi
ψ
ϕi
m
m
ztψ = ∧m
i=1 zti , and zt = ∨i=1 zti for ψ = ¬µ, ψ = ∧i=1 ϕi
m
and ψ = ∨i=1 ϕi , respectively. These constraints enforce that
ztψ = 1 if ψ holds at time t and ztψ = 0 otherwise.
3) Temporal constraints: We first present encodings for
the  and
operators. We will use these encodings to define
the encoding for the U[a,b] operator.
Always: ψ = [a,b] ϕ
N
Let aN
t = min(t + a, N ) and bt = min(t + b, N )
W
Vb̂N
bN
N
ϕ
ψ
ϕ
j
t
Define zt = ∧i=a
lj ∧ i=â
N zi ∧ (
N zi )
j=1
t
j
The logical operation ∧ on the variables ziϕ here is as defined
in Section IV-C2. Intuitively, this encoding enforces that the
formula ϕ is satisfied at every time step on the interval [a, b]
relative to time step t.

Eventually: ψ =

[a,b]

ϕ

WN
Wb̂N
bN
ϕ
ϕ
j
t
Define ztψ = ∨i=a
lj ∧ i=â
N zi ∧ (
N zi ) This encoding
j=1
t
j
enforces that the formula ϕ is satisfied at some time step on
the interval [a, b] relative to time step t.
Until: ψ = ϕ1 U[a,b] ϕ2
The bounded until operator U[a,b] can be defined in terms of
the unbounded U (inherited from LTL) as follows [18]:
ϕ1 U[a,b] ϕ2 = [0,a] ϕ1 ∧

[a,b]

ϕ2 ∧



As described above, to synthesize open-loop control for
unbounded (infinite-horizon) specifications, we parametrize
the trajectory as a lasso, i.e. constrain it to contain a loop. This
loop encoding is again inspired by the basic idea of bounded
model checking [10], which is to consider only a finite prefix
of a path when looking for a solution to an existential model
checking problem. A crucial observation is that although the
considered path prefix is finite, it can still represent an infinite
path if there is a loop back from the last state to any of the
previous states.
To enforce the existence of a loop in the finite system
trajectory, we introduce N binary variables l1 , ..., lN , which
determine where the loop forms. These are constrained such
that only one can be high at a time, and if lj = 1, then
xj−1 = xN . The following constraints enforce these requirements:
PN
•
j=1 lj = 1
• xN ≤ xj−1 + Mj (1 − lj ), j = 1, ..., N ,
• xN ≥ xj−1 + Mj (1 − lj ), j = 1, ..., N ,
where Mj are sufficiently large positive numbers, picked based
on X .



B. Loop constraints for trajectory parametrization



The first component of the set of constraints is provided by
the system model. Our approach applies to any system that
yields to a MILP formulation for model predictive control
over horizon N . The system constraints encode valid finite
(horizon-N ) trajectories for a system with form (1) – these
constraints hold if and only if the trajectory x(x0 , uN ) satisfies
(1) for t = 0, 1, ..., N . Note that this is quite general, and
accommodates any system for which the resulting constraints
and objectives form a mixed integer-linear program. An example is the smart grid regulation control system presented in
[16]. Other useful examples include mixed logical dynamical
systems such as those presented in [17]. Other cost functions
and system dynamics can also be included by using appropriate solvers.

where Mt are sufficiently large positive numbers, and t
are sufficiently small positive numbers that serve to bound
µ(xt ) away from 0. This encoding restricts the set of STL
formulas that can be encoded using our approach to those
over linear predicates, but admits arbitrary STL formulas over
such predicates.
2) Boolean operations on MILP variables: As described
in Section IV-C1, each predicate µ has an associated binary
variable ztµ which equals 1 if µ holds at time t, and 0
otherwise. In fact, by the recursive definition of our MILP
constraints on STL formulas, each operand ϕ in a Boolean
operation has a corresponding variable ztϕ which is 1 if ϕ holds
at t and 0 otherwise. Here we define Boolean operations on
these variables: these are the building blocks of our recursive
encoding. The definitions in this subsection are consistent with
those in [15].
Logical operations on variables ztψ ∈ [0, 1] are defined as
follows:
ztψ = 1 − ztϕ
Negation: ztψ = ¬ztϕ



A. Constraints on system evolution

[a,a] (ϕ1

U ϕ2 )

We will use the encoding of the unbounded U from
[10]. When encoding over infinite trajectories, this requires

5

an auxiliary encoding that prevents the pitfalls of circular
reasoning on the finite parametrization of the infinite
sequences. The interested reader is referred to [10] for
the details of the encoding. The auxiliary encoding of the
unbounded until is
hh ϕ1(U ϕ2 iit =
ztϕ2 ∨ (ztϕ1 ∧ hh ϕ1 U ϕ2 iit+1 ), t = 1, ..., N − 1
ϕ2
zN
.
With this definition in place, we define
ztϕ1

U ϕ2

ϕ1 U
= ztϕ2 ∨ (ztϕ1 ∧ zt+1

ϕ2

)

for t = 1, ..., N − 1, and
ϕ1
zN

U ϕ2

ϕ2
ϕ1
= zN
∨ (zN
∧(

N
_

(lj ∧ hh ϕ1 U ϕ2 iij ))).

j=1

[0,a] ϕ1



ϕ1 U[a,b] ϕ2

zt

∧ zt

= zt

[a,b]

ϕ2





Given this encoding of the unbounded until and the encodings
of [a,b] and [a,b] above, we can encode
∧ zt

[a,a] (ϕ1

U ϕ2 )

.

ztϕ

By induction on the structure of STL formulas ϕ,
= 1
if and only if ϕ holds on the system at time t. With this
motivation, given a specification ϕ, we add a final constraint:
z0ϕ = 1.

(2)

For a bounded horizon formula, the union of the STL
constraints, loop constraints and system constraints gives the
MILP encoding of Problem 1; this enables checking feasibility
of this set of constraints and finding a solution using an MILP
solver. Given an objective function on runs of the system, this
approach also enables finding the optimal open-loop trajectory
that satisfies the STL specification. Algorithm 1 reviews the
procedure for solving Problem 1.
Algorithm 1 Algorithm for Problem 1
1: procedure OPEN LOOP(f, x0 , w, N, ϕ, J)
2:
LOOP CONSTRAINTS ← Sec. IV-B
3:
SYSTEM CONSTRAINTS ← Sec. IV-A
4:
STL CONSTRAINTS ← Sec. IV-C2 OR Sec. IV-D
5:

u∗ ← argmin

J(x0 , u, w, ϕ)

u∈U N

s.t.

LOOP CONSTRAINTS
SYSTEM CONSTRAINTS
STL CONSTRAINTS

Return u∗
6: end procedure

D. Quantitative Encoding
The robustness of satisfaction of the STL specification, as
defined in II-C, provides a natural objective for the MILP
defined in Section IV-C, either in the absence of, or as
a complement to domain-specific objectives on runs of the
system. The robustness can be computed recursively on the
structure of the formula in conjunction with the generation of

constraints. Moreover, since max and min operations can be
expressed in an MILP formulation using additional binary variables, this does not add complexity to the encoding, although
the additional variables do make it more computationally
expensive in practice.
In this section, we sketch the MILP encoding of the
predicates and Boolean operators using the quantitative semantics; the encoding of the temporal operators builds on
these encodings, as in Section IV-C. Given a formula ϕ,
we introduce a variable rtϕ , and an associated set of MILP
constraints such that rtϕ > 0 if and only if ϕ holds at position
t. We recursively generate the MILP constraints, such that
r0ϕ determines whether a formula ϕ holds in the initial state.
Additionally, we enforce rtϕ = ρϕ (x, t).
For each predicate µ ∈ P , we now introduce variables rtµ
for time indices t = 0, 1, ..., N , and set rtµ = µ(xt ). To define
rtψ , where ψ is a Boolean formula, we inductively assume that
each operand ϕ has a corresponding variable rtϕ = ρϕ (x, t).
Then the Boolean operations are defined as:
rtψ = −rtφ
Negation: rtψ = ¬rtφ
ϕi
Conjunction: rtψ = ∧m
i=1 rti
m
X

i
pϕ
ti = 1

i=1
rtψ ≤ rtϕii , i = 1, ..., m
ψ
i
rtϕii − (1 − pϕ
ti )M ≤ rt

(3)
(4)
≤

rtϕii

+ M (1 −

i
pϕ
ti )

(5)

i
pϕ
ti

where we introduce new binary variables
for i = 1, ..., m,
and M is a sufficiently large positive number. Then equation
(3) enforces that there is one and only one j ∈ {1, ..., m}
ϕ
such that ptij = 1, equation (4) ensures that rtψ is smaller
ϕ
than all rtϕii , and equation (5) enforces that rtψ = rtjj if and
ϕj
only if ptj = 1. Together, these constraints enforce that rtψ =
mini (rtϕii ).
ϕi
Disjunction: ψ = ∨m
i=1 rti is encoded similarly to conjunc-

tion, replacing (4) with rtψ ≥ rtϕii , i = 1, ..., m. Using a similar
reasoning to that above, this enforces rtψ = maxi (rtϕii ).
The encoding for bounded temporal operators is defined as
in Section IV-C; robustness for the unbounded until is defined
using sup and inf instead of max and min, but these are
equivalent on our finite trajectory representation with discrete
time. By induction on the structure of STL formulas ϕ, this
construction yields rtϕ > 0 if and only if ϕ is satisfied at time
t. Therefore, we can replace the constraints over ztϕ in Section
IV-C by these constraints that compute the value of rtϕ , and
instead of (2), add the constraint r0ϕ > 0.
Since we consider only the discrete time semantics of
STL in this work, the Boolean encoding in Section IV-C
could be achieved by converting each formula to LTL, and
using existing encodings such as that in [15]. However, the
robustness-based encoding we presented in this section has no
natural analog for LTL. The advantage of this encoding is that
it allows us to maximize the value of r0ϕ , obtaining a trajectory
that maximizes robustness of satisfaction. Additionally, an
encoding based on robustness has the advantage of allowing
the STL constraints to be softened or hardened as necessary.

6

For example, if the original problem is infeasible, we can
allow ρϕ
0 > − for some  > 0, thereby easily modifying
the problem to allow a limited violation of the STL property.
The disadvantage is that it is more expensive to compute,
due the the additional binary variables introduced during each
Boolean operation. Additionally, including robustness as an
objective makes the cost function inherently non-convex, with
potentially many local minima, and harder to optimize. On the
other hand, the robustness constraints are more easily relaxed,
allowing us to use a simpler cost function, which can make
the problem more tractable.
E. Complexity
In general, our synthesis algorithm has the same complexity as MILPs, which are NP-hard, hence computationally
challenging when the dimensions of the problem grow. It
is nevertheless appropriate to characterize the computational
costs of our encoding and approach in terms of the number
of variables and constraints in the resulting MILP. In practice,
one measure of problem size is the number of binary variables
required to indicate the satisfaction of the predicates µ. This
depends directly on the number of predicates used in the STL
formula ϕ.
For the Boolean encoding, if P is the set of predicates
used in the formula, then O(N · |P |) binary variables are
introduced. In addition, continuous variables are introduced
during the MILP encoding of the STL formula. The number
of continuous variables used is O(N · |ϕ|), where |ϕ| is the
length (i.e. the number of operators) of the formula.
For the robustness-based encoding, O(N · |P |) continuous
variables are introduced (one per predicate per time step).
In addition, binary variables are introduced during the MILP
encoding of each operator in the STL formula. The number
of binary variables used is thus O(N · |ϕ|), where |ϕ| is the
number of operators of the formula.
Our synthesis algorithm also has polynomial runtime for the
following fragment of STL.
Definition 1 (SNN-STL). Safe Negation-Normal STL (SNNSTL) is the fragment of STL generated by the recursive
grammar
ϕ ::= π µ | ¬π µ | ϕ1 ∧ ϕ2 | [a,b] ϕ
SNN-STL has the following properties:
• All negations appear only on atomic propositions (pushed
down to the leaf nodes of the formula abstract syntax
tree).
• The only temporal operators are  (with unbounded and
bounded intervals).
• Only conjunctions are allowed, no disjunctions.
Such specifications are expressive enough to enforce, e.g.,
safety specifications in environments where the system state
is confined to a conjunction of polyhedra.
Let
OPEN LOOP NO STL(f, x0 , N, ϕ, J)
denote the procedure that is identical to Algorithm 1, except that the optimization problem in Step 5 is solved with
STL CONSTRAINTS = ∅.

ϕ, Σd , J

Controller
u∗ =
u∗0 u∗1 . . . u∗N −1

x0

Synthesis

w0 w1 . . . wN −1

u∗0

Plant Σ

x1

x0
w0
Fig. 2. Closed-loop (MPC) problem formulation. As in the open loop scenario,
a sequence of control inputs is synthesized from the specifications, dynamics
and cost function. However, at each time step, only the first computed input
is used by the plant.

Theorem 1 (Polynomial-time Synthesis for SNN-STL). Suppose that ϕ is in SNN-LTL and has linear predicates.
Then if OPEN LOOP NO STL(f, x0 , w, N, ϕ, J) is convex, so
is OPEN LOOP(f, x0 , N, ϕ, J).
Proof. The proof proceeds by induction on the structure
of the formula, showing that the constraints added by
STL CONSTRAINTS for each operator restrict the solution to
a convex set. First note that since the predicates π µ are linear,
negation of predicates preserves convexity, since ¬π µ is also
linear. Because the intersection of convex sets is convex, the
conjunction of a set of convex constraints is also convex.
Finally, since the  operator is implemented in terms of
conjunctions, the constraints imposed by  also preserve
convexity of the resulting optimization problem.
Informally, Theorem 1 states that our encoding of SNN-STL
constraints into an MILP preserves convexity in the resulting
optimization problem. The resulting optimization problem is
therefore encodable as an LP, i.e. without the use of integer
variables.
Corollary 1. Algorithm 1 is polynomial-time for SNN-STL
specifications ϕ.
V. M ODEL P REDICTIVE C ONTROL S YNTHESIS
In this section, we will describe a solution to Problem 2
by adding STL constraints to an MPC problem formulation.
At each step t of the MPC computation, we will search for
a finite trajectory of fixed horizon length H, such that the
accumulated trajectory satisfies ϕ.
A. Synthesis for bounded-time STL formulas
The length of the horizon H is chosen to be at least the
bound of formula ϕ. At time step 0, we will synthesize control
uH,0 using the open-loop formulation in Section IV, including

7

the STL constraints on the length-H trajectory, but without the
loop constraints. We will then execute only the first time step
u0H,0 . At the next step of the MPC, we will solve for uH,1 ,
while constraining the previous values of x0 , u0 in the MILP,
and the STL constraints on the trajectory up to time H. In this
manner, we will keep track of the history of states in order to
ensure that the formula is satisfied over the length-H prefix
of the trajectory, while solving for uH,t at every time step t.

Algorithm 2 MPC Algorithm for Problem 2
1: procedure MPC(f, x0 , φ =  ϕM P C , J)
2:
Let M be a large positive constant.
3:
Let H be the bound of ϕM P C .
4:
Set P0 = 0 and Pi = −M ∀0 < i ≤ H.
5:
wt ← PREDICT W(0).
6:
Compute u0 = u00 , u01 , ...., u02H−1 as:

B. Extension to unbounded formulas

7:
8:
9:
10:
11:

For certain types of unbounded formulas, we can stitch
together trajectories of length H using a receding horizon
approach, to produce an infinite computation that satisfies the
STL formula. An example of this is safety properties, i.e.
ϕ = (ϕM P C ) for bounded STL formulas ϕM P C . For such
formulas, at each step of the MPC computation, we will search
for a finite trajectory of horizon length H (determined from
ϕM P C as above) that satisfies ϕM P C .
We now describe this approach in more detail. At each step
t of the receding horizon control computation, we will employ
the open-loop approach in Section IV to find a finite trajectory
of fixed horizon length H, such that the trajectory accumulated
over time satisfies ϕ. Given a specification ϕ =  ϕM P C ,
where ϕM P C is a bounded-time formula with bound H. In
this case, we can stitch together trajectories of length H
using a receding horizon approach to produce an infinite
computation that satisfies the STL formula. At each step of the
receding horizon computation, we search for a finite trajectory
of horizon length 2H, keeping track of the past values and
robustness constraints necessary to determine satisfaction of
ϕ at every time step in the trajectory. Note that we omit the
loop constraints in this approach, because at each step we
search for a finite trajectory, rather than an infinite trajectory
with a finite parametrization.
First we define a procedure
OPEN LOOP∗ (f, x0 , w, N,  ϕM P C , J, PH , utold )
that takes additional inputs P = {P0 , P1 , ..., PH−1 } and
utold = u0 , u1 , ..., ut−1 , and is identical to Algorithm 1, except
that the optimization problem posed in Step 5 is solved without
the loop constraints, and with the added constraints:
ρϕ (f (x0 , u, w), i) > Pi ∀i ∈ [0, H − 1]
u[i...t] = utold
We then define a receding horizon control procedure as in
Algorithm 2. At each step, we are optimizing over a horizon
of 2H. We assume available a method PREDICT W(t) for
predicting the sequence of 2H environment inputs starting at
time step t.
Algorithm 2 has two phases, a transient phase (Lines 410) and a stationary phase (Lines 11-14). The transient phase
applies until an initial control sequence of length H has been
computed, and the stationary phase follows. In the transient
phase, the number of stored previous inputs (utold ) as well as
the number of time steps at which formula ϕM P C is enforced
(i.e. time steps for which Pi = 0) grows by one at each
iteration, until they both attain a maximum of H at iteration H.

u0 ← OPEN LOOP∗ (f, x0 , w0 , 2H, [0,H] ϕM P C , J, PH , ∅)
for t=1; t¡=H;t=t+1 do
Set utold = u00 , u11 , u22 , ..., ut−1
t−1 .
Set Pi = 0 for 0 ≤ i ≤ t, Pi = −M ∀t < i ≤ H.
wt ← PREDICT W(t).
Compute ut = ut0 , ut1 , ...., ut2H−1 as:
ut ← OPEN LOOP∗ (f, xt , wt , 2H, [0,H] ϕM P C , J, PH , utold )
12:
13:
14:
15:
16:

end for
while True do
t−1
t−1
t−1
.
Set utold = ut−1
1 , u2 , u3 , ..., ut
Set Pi = 0 for 0 ≤ i ≤ H.
wt ← PREDICT W(t).
ut ← OPEN LOOP∗ (f, xt , wt , 2H, [0,H] ϕM P C , J, PH , utold )

17:
18:

end while
end procedure

Every following iteration uses a window of size H for stored
previous inputs, and sets all Pi = 0. The size-H window
of previously-computed inputs advances forward one step in
time at each iteration after step H. In this manner, we keep
a record of the previously computed inputs required to ensure
satisfaction of ϕM P C up to H time steps in the past.
We now show that if Algorithm 2 does not terminate,
then the resulting infinite sequence of control inputs enforces
satisfaction of the specification φ =  ϕM P C .
Theorem 2. Let φ =  ϕM P C , and assume that u∗ is
an infinite sequence of control inputs generated by setting
u∗ [t] = ut0 , where ut = ut0 ut1 ...ut2H−1 is the control input
sequence of length 2H generated by Algorithm 2 at time t.
Then f (x0 , u∗ , w) |= ϕ.
Proof. Since H is the bound of ϕM P C , the satisfaction of
ϕM P C at time t is established by the control inputs u∗ [t :
t + H − 1]. At time t + H,
ut+H
old

=
=
=
=

ut+H
, ut+H
, ut+H
, ..., ut+H
0
1
2
t+H−1
t+H−1
t+H−1
u1
, u2
, ut+H−1
, ..., ut+H−1
3
t+H
t+1
t+2
t+H−1
utt , ut+1 , ut+2 , ..., ut+H−1
u∗ [t : t + H − 1],

and so all the inputs required to determine satisfaction of ϕ
at time t have been fixed. Moreover, if ut+H is successfully
computed, then by the correctness of Algorithm 1, ut+H
old has
H
∗
the property that f (xt , ut+H
old , w ) |= ϕM P C . Since u [t : t +
t+H
∗
H
H −1] = uold , we see that f (xt , u [t : t+h], w ) |= ϕM P C .
It follows that f (x0 , u∗ , w) |= ϕM P C .
We have therefore shown how a control input can be

8



synthesized for infinite sequences satisfying ϕ, by repeatedly
synthesizing control for sequences of length 2H. A similar
approach applies for formulas ϕM P C and ϕM P C U ψM P C ,
where ϕM P C , ψM P C are bounded-time.
Note that we assumed that PREDICT W(t) returns an exact
prediction of the disturbance signal over the next 2H time
steps. The correctness of our approach relies on this assumption. An interesting direction of future work is to relax this
requirement, demanding only an uncertain prediction of the
disturbance signal.
Remark 2. The control objective for MPC is usually to steer
the state to the origin or to an equilibrium state. Questions that
arise include those of ensuring feasibility at each time step,
closed-loop stability and near-optimal performance [19]. There
is a mature theory of stability for MPC, where the essential
ingredients are terminal costs, terminal constraint sets, and
local stabilizing controller that ensure closed-loop stability [5].
In this work, our control objective is not closed-loop stability, but satisfaction of an STL formula. We achieve this, as
detailed above, through choice of a sufficiently large prediction
horizon H. This can be compared with the manner in which
automatic satisfaction of a terminal constraint is sometimes
attained by prior choice of a sufficiently large horizon.
VI. E XPERIMENTAL C OMPARISON OF E NCODINGS
We implemented the Boolean and robust encodings using
the tools Breach [20] and YALMIP [21], and now present
results obtained with the following formulas:

ϕ1
ϕ2
ϕ3
ϕ4

#constraints
B
R
154
488
364
897
244 1282
574 1330

YALMIP
B
1.71
1.94
1.84
2.29

Time (s)
R
2.04
2.69
3.15
3.37

Solver time (s)
B
R
0.0070 0.0085
0.0115 0.0229
0.0064 0.1356
0.2167
238.6

TABLE I
B OOLEAN (B) VS ROBUST (R) ENCODINGS . YALMIP TIME REPRESENTS
THE TIME TAKEN BY THE TOOL YALMIP IN ORDER TO GENERATE THE
MILP AND S OLVER TIME IS THE TIME TAKEN BY THE SOLVER G UROBI TO
ACTUALLY SOLVE IT.

the Boolean encoding. The reason is hard to pinpoint without
a more thorough investigation, but we can already note that
solving a MILP is NP-hard, and while solvers use sophisticated
heuristics to mitigate this complexity, instances for which these
heuristics fail are bound to appear.
VII. C ASE S TUDY: B UILDING C LIMATE C ONTROL
A. Mathematical Model of a Building
Next we consider the problem of controlling building indoor
climate, using the model proposed by Maasoumy et al [22]. In
this section we present a summary of the building’s thermal
model.
As shown in Fig. 3, the building is modeled as a resistorcapacitor circuit with n nodes, m of which are rooms and
the remaining n − m are walls. We denote the temperature of
room ri by Tri . The wall and temperature of the wall between
rooms i and j are denoted by wi,j and Twi,j , respectively.
The temperature of wall wi,j and room ri are governed by the
following equations:









ϕ1
ϕ2
ϕ3
ϕ4

(1)

= [0,0.1] xt > 0.1
(1)
(2)
•
= [0,0.1] (xt > 0.1) ∧ [0,0.1] (xt < −0.5)
(1)
•
= [0,0.5] [0,0.1] (xt > 0.1)
(1)
(2)
•
= [0,0.2] (xt > 0.1 ∧ ( [0,0.1] (xt > 0.1)
(3)
∧ [0,0.1] (xt > 0.1)))
In this study, we used the trivial system x = u, where x
(1) (2) (3)
is a 3-dimensional signal (i.e. xt = xt xt xt ), so that
no constraint is generatedPfor the system dynamics, and the
N
cost function J(x, u) = k=1 kutk k1 . Note that the output
of this procedure for a formula ϕ is a signal of minimal
norm which satisfies ϕ when using the Boolean encoding and
which satisfies ϕ with a specified robustness ρϕ (x) = 0.1
for the robust encoding. For each formula we computed the
Boolean and robust encodings for an horizon N = 30 and
sampling time τ = 0.025s and report the number of constraints
generated by each encoding, the time to create the resulting
MILP with YALMIP and the time to solve it using the solver
Gurobi.1 All experiments were run on a laptop with an Intel
Core i7 2.3 GHz processor and 16 GB of memory.
A first observation is that for both encodings, most of the
time is spent creating the MILP, while solving it is done
in a fraction of a second. Also, while the robust encoding
generates 3 to 5 times more constraints, the computational
time to create and solve the corresponding MILPs is hardly
twice more. The exception is solving the MILP for ϕ4 , which
takes significantly more time for the robust encoding than for
•

Formula

1 http://www.gurobi.com/

w
Ci,j

X
dTwi,j
=
dt

k∈Nwi,j

Trk − Twi,j
+ ri,j αi,j Awi,j Qradi,j
Ri,jk
(6)

X Tk − Tr
dTri
i
Cir
=
+ ṁri ca (Tsi − Tri )+
dt
Ri,ki
k∈Nri

wi τwi Awini Qradi + Q̇inti ,

(7)

w
where Ci,j
, αi,j and Awi,j are heat capacity, a radiative heat
absorption coefficient, and the area of wi,j , respectively. Ri,jk
is the total thermal resistance between the centerline of wall
(i, j) and the side of the wall on which node k is located.
Qradi,j is the radiative heat flux density on wi,j . Nwi,j is the
set of all neighboring nodes to wi,j . ri,j is a wall identifier,
which equals 0 for internal walls and 1 for peripheral walls,
where either i or j is the outside node. Tri , Cir and ṁri are
the temperature, heat capacity and air mass flow into room i,
respectively. ca is the specific heat capacity of air, and Tsi is
the temperature of the supply air to room i. wi is a window
identifier, which equals 0 if none of the walls surrounding
room i have windows, and 1 if at least one of them does. τwi
is the transmissivity of the glass of window i, Awini is the
total area of the windows on walls surrounding room i, Qradi
is the radiative heat flux density per unit area radiated to room
i, and Q̇inti is the internal heat generation in room i. Nri is

9

Room temperature [oC]

23
Simulation
Measurement

22
21
20
19

18
12:00 am

6:00 am

12:00 pm

6:00 pm

12:00 am

6:00 am

12:00 pm
Time of day [hr]

6:00 pm

12:00 am

2

o

Disturbance [ C/hr]

3

1
0

−1
12:00 am

Fig. 4. Simulated temperature, measured temperature and unmodelled dynamics of a thermal zone in Bancroft library on UC Berkeley campus.

Fig. 3. Resistor-capacitor representation of a typical room with a window.

the set of neighboring room nodes for room i. Further details
on this thermal model can be found in [22].
The heat transfer equations for each wall and room yield
the system dynamics:
ẋt = f (xt , ut , wt ).
Here xt ∈ Rn is the state vector representing the temperature
of the nodes in the thermal network, and ut ∈ Rlm is the
input vector representing the air mass flow rate and discharge
air temperature of conditioned air into each thermal zone (with
l being the number of inputs to each thermal zone, e.g. two
for air mass flow and supply air temperature). The HVAC
system of the building considered for this study operates with
a constant supply air temperature, while air mass flow is the
time varying control input. Hence, in the following simulations
we consider supply air temperature constant and treat air mass
flow as the control signal. Vector wt stores the estimated
disturbance values, aggregating various unmodelled dynamics
such as Tout , Q̇int and Qrad , and can be estimated using
historical data [23]. yt ∈ Rm is the output vector, representing
the temperature of the thermal zones. The building model
was trained using historical data, and the result of the system
identification is shown in Fig. 4.
B. MPC for Building Climate Control
We consider the problem of controlling the above building’s
HVAC system using an MPC scheme. We adopt the MPC
formulation proposed by Maasoumy et al. [24], with the
objective of minimizing the total energy cost (in dollar value).
τ and H denote the length of each time slot and the prediction
horizon (in number of time slots) of the MPC, respectively.
Assume that the system dynamics are also discretized with

a sampling time of τ . Here we consider τ = 0.5 hr and
H = 24. At each time t, the predictive controller solves an
optimal control problem to compute ~ut =P[ut , . . . , ut+H−1 ],
H−1
and minimizes the cumulative norm of ut : k=0 kut+k k. We
assume known an occupancy function occt which is equal to
1 when the room is occupied and to 0 otherwise. The purpose
of the MPC is to maintain a comfort temperature given by
T comf whenever the room is occupied while minimizing the
cost of heating. This problem can be expressed as follows:
min
~
ut

H−1
X

kut+k k s.t.

k=0

xt+k+1 = f (xt+k , ut+k , wt+k ),
xt |= ϕ with ϕ = [0,H] ((occt > 0) ⇒ (Tt > Ttcomf )
ut+k ∈ Ut+k , k = 0, ..., H − 1
The STL formula was encoded using the robust MILP encoding and results are presented in Fig. 5. Again we observed
that creating the MILP structure was longer than solving an
instance of it (4.1s versus 0.15s). However, by using a proper
parametrization of the problem in YALMIP, the creation of the
MILP structure can be done once offline and reused online for
each step of the MPC, which makes the approach promising
and potentially applicable even for real-time applications.
VIII. C ASE S TUDY II: R EGULATION C ONTROL FOR
S MART G RID
A. Mathematical Model
The second case study we consider is the n-areas smart
grid model presented in [16] and depicted in Fig. 6. The
interconnection of power system components, including a
governor, turbine and generator in each area is shown in the
block diagram in Fig. 7. In the diagram, δPC is a control
input which acts against an increase or decrease in power
demand to regulate the system frequency ω, and δPD denotes
fluctuations in power demand, modeled as an exogenous input
(disturbance). Under steady state, we have: ω = ωo and
o
o
PM = P G = P M
, where ωo , Vto , and PM
are the nominal
values for rated frequency, terminal voltage and mechanical
power input.

10

70
Temperature
Comfort Limit

69
68
67
66
6

12

1
0.5
0

18

24

Time of day (Hours)

Occupancy
6

12

18

60

24
Air Flow

Fig. 7. Block diagram of power system and its relation to governor, turbine,
generator, and the AGC signal for each control area. More details on the
power grid model can be found in [16].

40

where νij is the transmission line stiffness coefficient, and
the state variable xi8 is the integral of xi1 .

20
0
6

12
Time of day (Hours)

18

24

Fig. 5. Room temperature control with constraints based on occupancy,
expressed in STL.

Fig. 6. Power system grid with n areas. The dynamics in each area is depicted
in Fig.7.

Next, we present the mathematical model for one area i
(note that superscripts refer to the control area, and subscripts
index states in each area).
i
i
i
i
dxi1
(−D i xi1 + δPM
− δPD
− δPtie
+ δPanc
)
=
,
dt
M
i
i
i
i
i
i
i
(x3 − x2 )
(x4 − x3 )
(xi5 − xi4 )
dx2
dx3
dx4
=
,
=
,
=
,
dt
T7i
dt
T6i
dt
T5i
i
dxi5
(PGV
− xi5 )
=
,
dt
T4i

dxi7
dt

=

(−xi7

+

i
δPC
T1i

−

dxi6
(xi7 − xi6 )
=
,
dt
T3i
xi1 /Ri )

,

dxi8
dt

(8a)
(8b)
(8c)

i

= x1

x(tk+1 ) = Ax(tk ) + B1 uanc (tk ) + B2 w(tk ).

(10)

n T
1
. . . δPanc
] are the ancillary inputs,
Where uanc = [δPanc
and the exogenous inputs (i.e. disturbances or variations in
1
n T
demands) are denoted by w = [δPD
. . . δPD
] . We propose
controller synthesis for the ancillary services, complementing
the primary control of AGC.

(8d)

i
i
i
where δPM
and PGV
are given by δPM
= K1i xi5 + K3i xi4 +
i
i i
i i
K5 x3 + K7 x2 , and PGV = (1 − T2 /T3 )xi6 + (T2 /T3 )xi7 . D
is the damping coefficient, M is the machine inertia constant,
R is the speed regulation constant, Ti ’s are time constants
for power system components, and Ki ’s are fractions of total
mechanical power outputs associated with different operating
i
parts of the turbine. δPtie
represents power transfer from area
i to other areas. In equation (8), the first state represents the
i
frequency increment, xi1 = δωi . It can be shown that Ptie
can
be obtained from
n
X
ij
δPtie
=
νij (xi8 − xj8 ),
(9)
j=1

The classical automatic generation control (AGC) implements a simple PI control to regulate the grid frequency. In a
multi-area power system, in addition to regulating frequency
within each area, the auxiliary control should maintain the
net interchange power with neighboring areas at scheduled
values [25]. This is generally accomplished by adding a tieline flow deviation to the frequency deviation in the auxiliary
feedback control loop. A suitable linear combination of the
frequency and tie-line deviations for area i, is known as
the Area Control Error (ACE): this measures the difference
between the scheduled and actual electrical generation within
a control area while taking frequency bias into account. The
i
ACE of area i is thus defined as ACE i = δPtie
+β i xi1 , and β i
is the bias coefficient of area i. The standard industry practice
is to set the bias β i at the so-called Area Frequency Response
Characteristic (AFRC), which is defined as β i = Di + 1/Ri .
The integral of ACE is used to construct the speed changer
position feedback control signal (δPCi ), i.e., δPCi = −K i xi9 ,
dxi
where K i is the feedback gain and dt9 = ACE i .
The resulting state space model can be discretized and
written in compact form as

B. MPC for Ancillary Services
We require that uanc be bounded and satisfies a maximum
ramp constraint, i.e., uanc ≤ uanc (tk ) ≤ uanc with uanc >
0 and |uanc (tk+1 ) − uanc (tk )| ≤ λ, for some λ > 0. At each
time step k, we thus solve the following problem:
min
Uanc (k)

s.t.

J(ACE, Uanc )

(11)

x(tk+1 ) = Ax(tk ) + B1 uanc (tk ) + B2 w(tk )
uanc ≤ uanc (tk+j ) ≤ uanc
|uanc (tk+j+1 ) − uanc (tk+j )| ≤ λ

where Uanc (k) = (uanc (tk ), uanc (tk+1 ), . . . , uanc (tk+H )) is
the vector of inputs from k to k + H and H is the prediction

11

0.3
1

ACE , τ=5

0.2

ACE1, τ=10

0.1
0
−0.1
−0.2
−0.3
−0.4

0

10

20

30

40

50

60

70

80

0.8

90

100

u1anc, τ=5

0.6

u1anc, τ=10

0.4
0.2
0
−0.2
−0.4
−0.6

0

10

20

30

40

50

60

70

80

90

100

Fig. 8. Comparison of ACE1 (Area Control Error of Control Area 1)
stabilization using ϕt in (12) with τ = 5 and τ = 10. The controller enforces
the stabilization delay in both cases.

horizon. All the constraints of problem (11) that depend on j
should hold for j = 0, 1, . . . , H − 1.
The cost function proposed in [16] minimizes the `2 norm
of the ACE signal in areas i = 1, . . . , n, by exploiting the
ancillary service available in each area, while taking into
account the system dynamics and constraints. We propose to
constrain the ACE signal to satisfy a specified set of STL
properties, while minimizing the ancillary service used by
each
Thus we defined J(ACE, Uanc ) = kUanc k`2 =
P2 area.
PH−1
2
i
i=1
j=0 (Uanc [k + j]) , and an STL formula ϕ which
says that whenever |ACEi | is larger than 0.01, it should
become less than 0.01 in less than τ s. More precisely we
used ϕ = (ϕt ) with
1

| < .01)
2
(|ACE
| < .01)
[0,τ ]

[0,τ ] (|ACE

 

ϕt = ¬(|ACE1 | < .01)) ⇒ (
∧ (¬(|ACE2 | < .01)) ⇒ (

(12)

We encoded this formula and added the resulting constraints
to the MPC problem as described in the previous sections, and
solved it for different values of τ . Results are shown in Fig. 8,
and demonstrate that the STL constraint is correctly enforced
in the stabilization of the ACE signal.
IX. R ELATED W ORK
Receding horizon control for temporal logic has been considered before in the context of LTL [26], where the authors
propose a reactive synthesis scheme for specifications with
GR(1) goals. The authors in [27] also propose an MPC scheme
for specifications in synthetically co-safe LTL – our approach
extends synthesis capabilities to a wider class of temporal
logic specifications. In [28], the authors consider full LTL
but use an automata-based approach, involving potentially
expensive computations of a finite state abstraction of the
system and a Buchi automaton for the specification. We
circumvent these expensive operations using a bounded model
checking (BMC) approach to synthesis. In [17], the authors

present a model predictive control scheme to stabilize mixed
logical dynamical systems on desired reference trajectories,
while fulfilling propositional logic constraints and heuristic
rules. A major contribution of this work is to extend the
constraint specification language for such systems to STL
specifications, which allow expression of complex temporal
properties including safety, liveness, and response.
Our work extends the standard BMC paradigm for finite
discrete systems [14] to STL, which accommodates continuous
systems. In BMC, discrete state sequences of a fixed length,
representing counterexamples or plans, are obtained as satisfying assignments to a Boolean satisfiability (SAT) problem.
The approach has been extended to hybrid systems, either by
computing a discrete abstraction of the system [29], [30] or by
extending SAT solvers to reason about linear inequalities [31],
[32]. Similarly, MILP encodings inspired by BMC have been
used to generate trajectories for continuous systems with LTL
specifications [15], [33], [34], and for a restricted fragment
of MTL without nested operators [35]. While we draw much
inspiration from these early efforts, ours is the first work to
consider a BMC approach to synthesis for full STL.
X. C ONCLUDING R EMARKS
The main contribution of this paper is a pair of bounded
model checking style encodings for signal temporal logic
specifications as mixed integer linear constraints. We showed
how our encodings can be used to generate control for systems
that must satisfy STL properties, and additionally to ensure
maximum robustness of satisfaction. Our formulation of the
STL synthesis problem can be used as part of existing controller synthesis frameworks to compute feasible and optimal
controllers for cyber-physical systems. We presented experimental results for controller synthesis on simplified models of
a smart micro-grid and HVAC system, and showed how the
MPC schemes in these examples can be framed in terms of
synthesis from an STL specification, with simulation results
illustrating the effectiveness of our proposed synthesis.
We have demonstrated the ability to synthesize control for
systems on both the demand and supply sides of a smart grid.
We view this as progress toward a contract-based framework
for specifying and designing components of the smart grid
and their interactions using STL specifications. Future work
includes a reactive synthesis approach to synthesizing control
inputs for systems operating in uncertain environments: we
have already demonstrated preliminary results in this direction
in [36]. We will also further explore synthesis in an MPC
framework for unbounded STL properties. As mentioned in
Section V-B, this is an easy extension of our approach for
certain types of properties. Extending this to arbitrary properties has ties to online monitoring of STL properties [37],
which is another direction of further exploration.
R EFERENCES
[1] G. E. Fainekos, A. Girard, H. Kress-Gazit, and G. J. Pappas,
“Temporal logic motion planning for dynamic robots,” Automatica,
vol. 45, no. 2, pp. 343 – 352, 2009. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S000510980800455X

12

[2] P. Nuzzo, H. Xu, N. Ozay, J. Finn, A. Sangiovanni-Vincentelli, R. Murray, A. Donze, and S. Seshia, “A contract-based methodology for aircraft
electric power system design,” Access, IEEE, vol. PP, no. 99, pp. 1–1,
2013.
[3] R. Alur, T. A. Henzinger, G. Lafferriere, and G. J. Pappas, “Discrete
abstractions of hybrid systems,” Proceedings of the IEEE, vol. 88, no. 7,
pp. 971–984, 2000.
[4] R. M. Murray, J. Hauser, A. Jadbabaie, M. B. Milam, N. Petit, W. B.
Dunbar, and R. Franz, “Online control customization via optimizationbased control,” in In Software-Enabled Control: Information Technology
for Dynamical Systems. Wiley-Interscience, 2002, pp. 149–174.
[5] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M. Scokaert,
“Constrained model predictive control: Stability and optimality,”
Automatica, vol. 36, no. 6, pp. 789–814, 2000. [Online]. Available:
http://dx.doi.org/10.1016/S0005-1098(99)00214-9
[6] O. Maler and D. Nickovic, “Monitoring temporal properties of continuous signals,” in FORMATS/FTRTFT, 2004, pp. 152–166.
[7] G. E. Fainekos and G. J. Pappas, “Robustness of temporal
logic specifications for continuous-time signals,” Theor. Comput.
Sci., vol. 410, no. 42, pp. 4262–4291, 2009. [Online]. Available:
http://dx.doi.org/10.1016/j.tcs.2009.06.021
[8] A. Donzé and O. Maler, “Robust satisfaction of temporal logic over
real-valued signals,” in FORMATS, 2010, pp. 92–106.
[9] M. Maasoumy, P. Nuzzo, F. Iandola, M. Kamgarpour, A. SangiovanniVincentelli, and C. Tomlin, “Optimal load management system for
aircraft electric power distribution,” in IEEE Conference on Decision
and Control (CDC), 2013.
[10] A. Biere, K. Heljanko, T. A. Junttila, T. Latvala, and V. Schuppan,
“Linear encodings of bounded LTL model checking,” Logical Methods
in Computer Science, vol. 2, no. 5, 2006.
[11] V. Raman, M. Maasoumy, and A. Donzé, “Model predictive control
from signal temporal logic specifications: A case study,” in Proceedings
of the 4th ACM SIGBED International Workshop on Design, Modeling,
and Evaluation of Cyber-Physical Systems, ser. CyPhy’14. New
York, NY, USA: ACM, 2014, pp. 52–55. [Online]. Available:
http://doi.acm.org/10.1145/2593458.2593472
[12] V. Raman, M. Maasoumy, A. Donzé, R. M. Murray, A. SangiovanniVincentelli, and S. A. Seshia, “Model predictive control with signal
temporal logic specifications,” in Proc. of the IEEE Conf. on Decision
and Control, 2014.
[13] G. E. Fainekos and G. J. Pappas, “Robust sampling for MITL
specifications,” in Formal Modeling and Analysis of Timed Systems, 5th
International Conference, FORMATS 2007, Salzburg, Austria, October
3-5, 2007, Proceedings, 2007, pp. 147–162. [Online]. Available:
http://dx.doi.org/10.1007/978-3-540-75454-1 12
[14] A. Biere, A. Cimatti, E. M. Clarke, and Y. Zhu, “Symbolic model
checking without BDDs,” in TACAS, 1999, pp. 193–207.
[15] E. M. Wolff, U. Topcu, and R. M. Murray, “Optimization-based
trajectory generation with linear temporal logic specifications,” in 2014
IEEE International Conference on Robotics and Automation, ICRA
2014, Hong Kong, China, May 31 - June 7, 2014, 2014, pp. 5319–5325.
[Online]. Available: http://dx.doi.org/10.1109/ICRA.2014.6907641
[16] M. Maasoumy, B. M. Sanandaji, A. Sangiovanni-Vincentelli, and
K. Poolla, “Model predictive control of regulation services from commercial buildings to the smart grid,” in IEEE American Control Conference (ACC), 2014.
[17] A. Bemporad and M. Morari, “Control of systems integrating logic,
dynamics, and constraints,” Automatica, vol. 35, no. 3, pp. 407–427,
1999. [Online]. Available: http://dx.doi.org/10.1016/S0005-1098(98)
00178-2
[18] A. Donzé, T. Ferrère, and O. Maler, “Efficient robust monitoring for
stl,” in CAV, 2013, pp. 264–279.
[19] M. Morari and J. Lee, “Model predictive control: Past, present and
future,” Computers & Chemical Engineering, vol. 23, no. 4, pp.
667–682, 1999. [Online]. Available: https://control.ee.ethz.ch/index.cgi?
page=publications;action=details;id=1641
[20] A. Donzé, “Breach, a toolbox for verification and parameter synthesis
of hybrid systems,” in CAV, 2010, pp. 167–170.
[21] J. Lfberg, “YALMIP : A toolbox for modeling and optimization in
MATLAB,” in Proceedings of the CACSD Conference, Taipei, Taiwan,
2004. [Online]. Available: http://users.isy.liu.se/johanl/yalmip
[22] M. Maasoumy Haghighi, “Controlling energy-efficient buildings in
the context of smart grid: A cyber physical system approach,”
Ph.D. dissertation, University of California, Berkeley, Dec 2013.
[Online]. Available: http://www.eecs.berkeley.edu/Pubs/TechRpts/2013/
EECS-2013-244.html

[23] M. Maasoumy and A. Sangiovanni-Vincentelli, “Total and peak energy
consumption minimization of building hvac systems using model predictive control,” Design Test of Computers, IEEE, vol. 29, no. 4, pp. 26
–35, aug. 2012.
[24] M. Maasoumy, M. Razmara, M. Shahbakhti, and A. SangiovanniVincentelli, “Selecting building predictive control based on model uncertainty,” in IEEE American Control Conference (ACC 2014), Portland,
USA, June 2014.
[25] H. Bevrani, Robust Power System Frequency Control, ser. Power Electronics and Power Systems. Springer, 2009.
[26] T. Wongpiromsarn, U. Topcu, and R. M. Murray, “Receding horizon
temporal logic planning,” IEEE Trans. Automat. Contr., vol. 57, no. 11,
pp. 2817–2830, 2012.
[27] E. A. Gol and M. Lazar, “Temporal logic model predictive control
for discrete-time systems,” in Proceedings of the 16th international
conference on Hybrid systems: computation and control, HSCC 2013,
April 8-11, 2013, Philadelphia, PA, USA, 2013, pp. 343–352. [Online].
Available: http://doi.acm.org/10.1145/2461328.2461379
[28] X. C. Ding, M. Lazar, and C. Belta, “LTL receding horizon control for
finite deterministic systems,” Automatica, vol. 50, no. 2, pp. 399–408,
2014. [Online]. Available: http://dx.doi.org/10.1016/j.automatica.2013.
11.030
[29] G. J. P. Nicoló Giorgetti and A. Bemporad, “Bounded model checking
of hybrid dynamical systems,” in Decision and Control, 2005 and 2005
European Control Conference. CDC-ECC ’05. 44th IEEE Conference
on, Dec 2005, pp. 672–677.
[30] S. Jha, B. A. Brady, and S. A. Seshia, “Symbolic reachability analysis
of lazy linear hybrid automata,” in Proc. 5th International Conference
on Formal Modeling and Analysis of Timed Systems (FORMATS), ser.
Lecture Notes in Computer Science, vol. 4763, 2007, pp. 241–256.
[31] G. Audemard, M. Bozzano, A. Cimatti, and R. Sebastiani, “Verifying
industrial hybrid systems with MathSAT,” Electronic Notes in Theoretical Computer Science, vol. 119, no. 2, pp. 17 – 32, 2005, proceedings
of the 2nd International Workshop on Bounded Model Checking
(BMC 2004) Bounded Model Checking 2004. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S1571066105000885
[32] M. Franzle and C. Herde, “Efficient proof engines for bounded
model checking of hybrid systems,” Electronic Notes in Theoretical
Computer Science, vol. 133, no. 0, pp. 119 – 137, 2005,
proceedings of the Ninth International Workshop on Formal Methods
for Industrial Critical Systems (FMICS 2004) Formal Methods
for Industrial Critical Systems 2004. [Online]. Available: http:
//www.sciencedirect.com/science/article/pii/S1571066105050279
[33] S. Karaman, R. G. Sanfelice, and E. Frazzoli, “Optimal control of mixed
logical dynamical systems with linear temporal logic specifications,” in
Decision and Control, 2008. CDC 2008. 47th IEEE Conference on, Dec
2008, pp. 2117–2122.
[34] Y. Kwon and G. Agha, “Ltlc: Linear temporal logic for control,” in
HSCC, M. Egerstedt and B. Mishra, Eds., 2008, pp. 316–329.
[35] S. Karaman and E. Frazzoli, “Vehicle routing problem with metric
temporal logic specifications,” in Proceedings of the 47th IEEE
Conference on Decision and Control, CDC 2008, December 9-11,
2008, Cancún, México, 2008, pp. 3953–3958. [Online]. Available:
http://dx.doi.org/10.1109/CDC.2008.4739366
[36] 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, 2015, pp. 239–248.
[37] J. V. Deshmukh, A. Donzé, S. Ghosh, X. Jin, G. Juniwal, and S. A.
Seshia, “Robust online monitoring of signal temporal logic,” in Runtime
Verification - 6th International Conference, RV 2015 Vienna, Austria,
September 22-25, 2015. Proceedings, 2015, pp. 55–70. [Online].
Available: http://dx.doi.org/10.1007/978-3-319-23820-3 4
Vasumathi Raman received the B.A. degree in
Computer Science and Mathematics from Wellesley
College in 2007 and the M.S. and Ph.D. degrees
in Computer Science from Cornell University in
2011 and 2013, respectively. She was a postdoctoral
scholar in the Department of Computing & Mathematical Sciences at the California Institute of Technology from 2013-2015, and is currently a Senior
Scientist at the United Technologies Research Center
in Berkeley, CA. Her research explores algorithmic
methods for designing and controlling autonomous
systems, guaranteeing correctness with respect to user-defined specifications.

13

physical systems and systems biology.

Mehdi Maasoumy received his PhD in Electrical
Engineering & Computer Sciences in 2013 from UC
Berkeley, where he was advised by Prof. Alberto
Sangiovanni-Vincentelli. He received his MSc degree in 2010 from Mechanical Engineering Department at UC Berkeley with major in Controls and
minor in Optimization with Designated Emphasis
in Energy Systems and Technology (DEEST). He
graduated form Sharif University of Technology in
2008 with a BSc degree in Mechanical Engineering.

Richard M. Murray received the B.S. degree in
Electrical Engineering from California Institute of
Technology in 1985 and the M.S. and Ph.D. degrees
in Electrical Engineering and Computer Sciences
from the University of California, Berkeley, in 1988
and 1991, respectively. He is currently the Thomas
E. and Doris Everhart Professor of Control & Dynamical Systems and Bioengineering at Caltech.
Murray’s research is in the application of feedback
and control to networked systems, with applications
in biology and autonomy. Current projects include
analysis and design biomolecular feedback circuits; specification, design and
synthesis of networked control systems; and novel architectures for control
using slow computing.

Alberto Sangiovanni Vincentelli holds the Edgar
L. and Harold H. Buttner Chair of Electrical Engineering and Computer Sciences at the University of
California at Berkeley. He has been on the Faculty of
the Department since 1976. He obtained an electrical
engineering and computer science degree (“Dottore
in Ingegneria”) summa cum laude from the Politecnico di Milano, Milano, Italy in 1971. In 19801981, he spent a year as a Visiting Scientist at the
Mathematical Sciences Department of the IBM T.J.
Watson Research Center. In 1987, he was Visiting
Professor at MIT. He has held a number of visiting professor positions at
Italian Universities, including Politecnico di Torino, Università di Roma, La
Sapienza, Università di Roma, Tor Vergata, Università di Pavia, Università
di Pisa, Scuola di Sant’Anna. He was awarded the IEEE/RSE Wolfson
James Clerk Maxwell Medal “for groundbreaking contributions that have
had an exceptional impact on the development of electronics and electrical
engineering or related fields” and the Kaufman Award of the Electronic Design
Automation Council for “pioneering contributions to EDA” . He is Honorary
Professor at Politecnico di Torino and has Honorary Doctorates from the
University of Aalborg and KTH. He helped founding Cadence and Synopsys,
the two leading companies in EDA. He is the author of over 850 papers, 18
books and 2 patents in the area of design tools and methodologies, large scale
systems, embedded systems, hybrid systems and innovation.

Alexandre Donzé is a research scientist at the
University of California, Berkeley in the department
of Electrical Engineering and Computer Science. He
received his Ph.D. degree in Mathematics and Computer Science from the University of Joseph Fourier
at Grenoble in 2007. He worked as a post-doctoral
researcher at Carnegie Mellon University in 2008,
and at Verimag in Grenoble from 2009 to 2012. His
research interests are in simulation-based design and
verification techniques using formal methods, Signal
Temporal Logic (STL) with applications to cyberSanjit A. Seshia received the B.Tech. degree in
Computer Science and Engineering from the Indian
Institute of Technology, Bombay in 1998, and the
M.S. and Ph.D. degrees in Computer Science from
Carnegie Mellon University in 2000 and 2005 respectively. He is currently an Associate Professor
in the Department of Electrical Engineering and
Computer Sciences at the University of California,
Berkeley. His research interests are in dependable
computing and computational logic, with a current
focus on applying automated formal methods to embedded and cyber-physical systems, electronic design automation, computer
security, and synthetic biology. He has served as an Associate Editor of the
IEEE Transactions on Computer-Aided Design of Integrated Circuits and
Systems. His awards and honors include a Presidential Early Career Award
for Scientists and Engineers (PECASE) from the White House, an Alfred P.
Sloan Research Fellowship, the Prof. R. Narasimhan Lecture Award, and the
School of Computer Science Distinguished Dissertation Award at Carnegie
Mellon University.

