Control of Nonlinear Switched Systems Based on Validated
Simulation

arXiv:1611.06692v1 [cs.SY] 21 Nov 2016

Adrien Le Coënt? a , Julien Alexandre dit Sandretto b , Alexandre Chapoutot b ,
Laurent Fribourg c
a
b

CMLA, ENS Cachan, CNRS, Université Paris-Saclay, 61 av. du Président Wilson, 94235 Cachan Cedex, France

U2IS, ENSTA ParisTech, Université Paris-Saclay, 828, Boulevard des Maréchaux, 91762 Palaiseau Cedex, France
c

LSV, ENS Cachan, CNRS, Université Paris-Saclay, 61 av. du Président Wilson, 94235 Cachan Cedex, France

Abstract
We present an algorithm of control synthesis for nonlinear switched systems, based on an existing procedure of state-space
bisection and made available for nonlinear systems with the help of validated simulation. The use of validated simulation also
permits to take bounded perturbations and varying parameters into account. It is particularly interesting for safety critical
applications, such as in aeronautical, military or medical fields. The whole approach is entirely guaranteed and the induced
controllers are correct-by-design.
Key words: Nonlinear control systems, reachability, formal methods, numerical simulation, control system synthesis

1

Introduction

We focus here on switched control systems, a class of
hybrid systems recently used with success in various domains such as automotive industry and power electronics. These systems are merely described by piecewise
dynamics, periodically sampled with a given period. At
each period, the system is in one and only one mode,
decided by a control rule [14,23]. Moreover, the considered systems can switch between any two modes instantaneously. This simplification can be easily by-passed by
the addition of intermediate facticious modes.
In this paper, we consider that these modes are represented by nonlinear ODEs. In order to compute the control of a switched system, we do need the solution of differential equations. In the general case, differential equations can not be integrated formally, and a numerical in? This paper was not presented at any IFAC meeting. A
short version of this paper appeared in SNR’16 [9]. Corresponding author A. Le Coënt Tel. +33147407429.
Email addresses: adrien.le-coent@ens-cachan.fr
(Adrien Le Coënt? ), alexandre@ensta.fr (Julien
Alexandre dit Sandretto), chapoutot@ensta.fr (Alexandre
Chapoutot), fribourg@lsv.ens-cachan.fr (Laurent
Fribourg).

Preprint submitted to Automatica

tegration scheme is used to approximate the state of the
system. With the objective of computing a guaranteed
control, we base our approach on validated simulation
(also called “reachability analysis”). The guaranteed or
validated solution of ODEs using interval arithmetic is
mainly based on two kinds of methods based on: i) Taylor
series [29,30,24,11] ii) Runge-Kutta schemes [7,15,6,2].
The former is the oldest method used in interval analysis community because the expression of the bound of a
Taylor series is simple to obtain. Nevertheless, the family
of Runge-Kutta methods is very important in the field of
numerical analysis. Indeed, Runge-Kutta methods have
several interesting stability properties which make them
suitable for an important class of problems. Our tool [1]
implements Runge-Kutta based methods which prove
their efficiency at low order for short simulation (fixed
by sampling period of controller).
In the methods of symbolic analysis and control of hybrid systems, the way of representing sets of state values
and computing reachable sets for systems defined by
autonomous ordinary differential equations (ODEs), is
fundamental (see, e.g., [16,4]). Many tools using, eg. linearization or hybridization of these dynamics are now
available (e.g., Spacex [13], Flow* [8], iSAT-ODE [12]).
An interesting approach appeared recently, based on
the propagation of reachable sets using guaranteed

22 March 2018

Problem 1 (Control Synthesis Problem) Let us
consider a sampled switched system (2). Given three sets
R, S, and B, with R ∪ B ⊂ S and R ∩ B = ∅, find a rule
σ̃(x) such that, for any x(0) ∈ R

Runge-Kutta methods with adaptive step size control
(see [6,19]). An originality of the present work is to use
such guaranteed integration methods in the framework
of switched systems. This notion of guarantee of the
results is very interesting, because we are mainly interested into critical domain, such as aeronautical, military
and medical ones. Other symbolic approaches for control synthesis of switched systems include the construction of a discrete abstraction of the original system on
a grid of the state space. This can be done by computing symbolic models that are approximately bisimilar
[17] or approximately alternatingly similar [34] to the
original system. Another recent symbolic approach relies on feedback refinement relations [31]. We compare
our work with the last two approaches, which are the
closest related methods since the associated tools (respectively PESSOA [26] and SCOTS [32]) are used to
perform control synthesis on switched systems without
any stability assumptions, such as the present method.

• τ -stability 1 : x(t) returns in R infinitely often, at some
multiples of sampling time τ .
• safety: x(t) always stays in S\B.
Under the above-mentioned notation, we propose a procedure which solves this problem by constructing a law
σ̃(x), such that for all x0 ∈ R, and under the unknown
bounded perturbation d, there exists π = σ̃(x0 ) ∈ U k
for some k such that:



x(t0 + kτ ; t0 , x0 , d, π) ∈ R
∀t ∈ [t0 , t0 + kτ ], x(t; t0 , x0 , d, π) ∈ S

∀t ∈ [t0 , t0 + kτ ], x(t; t0 , x0 , d, π) ∈
/B

The paper is divided as follows. In Section 2, we introduce some preliminaries on switched systems and some
notation used in the following. In Section 3, the guaranteed integration of nonlinear ODEs is presented. In
Section 4, we present the main algorithm of state-space
bisection used for control synthesis. In Section 5, the
whole approach is tested on three examples of the literature. We give some performance tests and compare our
approach with the state-of-the-art tools in section 6. We
conclude in section 7.
2

Such a law permits to perform an infinite-time statedependent control. The synthesis algorithm is described
in Section 4 and involves guaranteed set based integration presented in the next section, the main underlying
tool is interval analysis [29]. To tackle this problem, we
introduce some definitions. In the following, we will often use the notation [x] ∈ IR (the set of intervals with
real bounds) where [x] = [x, x] = {x ∈ R | x 6 x 6 x}
denotes an interval. By an abuse of notation [x] will also
denote a vector of intervals, i.e., a Cartesian product of
intervals, a.k.a. a box. In the following, the sets R, S and
B are given under the form of boxes.

Switched systems

Let us consider the nonlinear switched system
ẋ(t) = fσ(t) (x(t), d(t))

(1)
Definition 1 (Initial Value Problem (IVP))
Consider an ODE with a given initial condition

n

defined for all t ≥ 0, where x(t) ∈ R is the state of
the system, σ(·) : R+ −→ U is the switching rule, and
d(t) ∈ Rm is a bounded perturbation. The finite set U =
{1, . . . , N } is the set of switching modes of the system.
We focus on sampled switched systems: given a sampling
period τ > 0, switchings will occur at times τ , 2τ , . . .
The switching rule σ(·) is thus piecewise constant, we will
consider that σ(·) is constant on the time interval [(k −
1)τ, kτ ) for k ≥ 1. We call “pattern” a finite sequence
of modes π = (i1 , i2 , . . . , ik ) ∈ U k . With such a control
input, and under a given perturbation d, we will denote
by x(t; t0 , x0 , d, π) the solution at time t of the system

x(0) ∈ X0 , d(t) ∈ [d],
(3)
with f : R+ × Rn × Rm → Rn assumed to be continuous
in t and d and globally Lipschitz in x. We assume that parameters d are bounded (used to represent a perturbation,
a modeling error, an uncertainty on measurement, . . . ).
An IVP consists in finding a function x(t) described by
the ODE (3) for all d(t) lying in [d] and for all the initial
conditions in X0 .
ẋ(t) = f (t, x(t), d(t))

with

Definition 2 Let X ⊂ Rn be a box of the state space.
Let π = (i1 , i2 , . . . , ik ) ∈ U k . The successor set of X via
π, denoted by P ostπ (X), is the (over-approximation of
the) image of X induced by application of the pattern π,

ẋ(t) = fσ(t) (x(t), d(t)),
x(t0 ) = x0 ,
∀j ∈ {1, . . . , k}, σ(t) = ij ∈ U for t ∈ [(j − 1)τ, jτ ).
(2)
We address the problem of synthesizing a statedependent switching rule σ̃(x) for (2) in order to verify
some properties. The problem is formalized as follows:

1

This definition of stability is different from the stability in
the Lyapunov sense.

2

i.e., the solution at time t = kτ of

is defined by (see [6]):

ẋ(t) = fσ(t) (x(t), d(t)),
x(0) = x0 ∈ X,
∀t ≥ 0, d(t) ∈ [d],
∀j ∈ {1, . . . , k}, σ(t) = ij ∈ U for t ∈ [(j − 1)τ, jτ ).
(4)

x(tn ; xn−1 ) − xn =


hp+1
dp+1 φ
f (p) (ξ, x(ξ; xn−1 ), d) − p+1 (η)
(p + 1)!
dt
ξ ∈]tn , tn+1 [ and η ∈]tn , tn+1 [ . (6)
The function f (n) stands for the n-th derivative of funcn
tion f w.r.t. time t that is ddtnf and h = tn+1 − tn is
n
the step-size. The
Ps function φ : R → R is defined by
φ(t) = xn + h i=1 bi ki (t) where ki (t) are defined as
Equation (5).

Definition 3 Let X ⊂ Rn be a box of the state space.
Let π = (i1 , i2 , . . . , ik ) ∈ U k . We denote by T ubeπ (X)
the union of boxes covering the trajectories of IVP (4),
which construction is detailed in Section 3.
3

The challenge to make Runge-Kutta integration schemes
safe w.r.t. the true solution of IVP is then to compute
a bound of the result of Equation (6). In other words
we have to bound the value of f (p) (ξ, x(ξ; xn−1 ), d) and
p+1
the value of ddtp+1φ (η). The latter expression is straightforward to bound because the function φ only depends
on the value of the step-size h, and so does its (p + 1)th derivative. The bound is then obtain using the affine
arithmetic [10,3].

Validated simulation

In this section, we describe our approach for validated
simulation based on Runge-Kutta methods [6,2].
A numerical integration method computes a sequence of approximations (tn , xn ) of the solution
x(t; x0 ) of the IVP defined in Equation (3) such that
xn ≈ x(tn ; xn−1 ). The simplest method is Euler’s
method in which tn+1 = tn + h for some step-size h and
xn+1 = xn + h × f (tn , xn , d); so the derivative of x at
time tn , f (tn , xn , d), is used as an approximation of the
derivative on the whole time interval to perform a linear interpolation. This method is very simple and fast,
but requires small step-sizes. More advanced methods
coming from the Runge-Kutta family use a few intermediate computations to improve the approximation of
the derivative. The general form of an explicit s-stage
Runge-Kutta formula, that is using s evaluations of f , is
xn+1 = xn + h

s
X

However, the expression f (p) (ξ, x(ξ; xn−1 ), d) is not so
easy to bound as it requires to evaluate f for a particular value of the IVP solution x(ξ; xn−1 ) at an unknown
time ξ ∈]tn , tn+1 [. The solution used is the same as the
one found in [30,7] and it requires to bound the solution of IVP on the interval [tn , tn+1 ]. This bound is usually computed using the Banach’s fixpoint theorem applied with the Picard-Lindelöf operator, see [30]. This
operator is used to compute an enclosure of the solution
[x̃] of IVP over a time interval [tn , tn+1 ], that is for all
t ∈ [tn , tn+1 ], x(t; xn−1 ) ∈ [x̃]. We can hence bound f (p)
substituting x(ξ; xn−1 ) by [x̃].

bi k i ,

i=1

For a given pattern of switched modes π = (i1 , . . . , ik ) ∈
U k of length k, we are able to compute, for j ∈ {1, .., k},
the enclosures:


k1 = f tn , xn , d ,
i−1


X
ki = f tn + ci h, xn + h
aij kj , d , i = 2, 3, . . . , s .
j=1

• [xj ] 3 x(tj );
• [x̃j ] 3 x(t), for t ∈ [(j − 1)τ, jτ ].

(5)
The coefficients ci , aij and bi fully characterize the
method. To make Runge-Kutta validated, the challenging question is how to compute a bound on the distance
between the true solution and the numerical solution,
defined by x(tn ; xn−1 ) − xn . This distance is associated
to the local truncation error (LTE) of the numerical
method.

with respect to the system of IVPs:

















To bound the LTE, we rely on order condition [18] respected by all Runge-Kutta methods. This condition
states that a method of this family is of order p iff the
p+1 first coefficients of the Taylor expansion of the solution and the Taylor expansion of the numerical methods
are equal. In consequence, LTE is proportional to the Lagrange remainders of Taylor expansions. Formally, LTE

ẋ(t) = fσ(t) (t, x(t), d(t)),
x(t0 = 0) ∈ [x0 ], d(t) ∈ [d],
σ(t) = i1 , ∀t ∈ [0, t1 ], t1 = τ
..
.






ẋ(t) = fσ(t) (t, x(t), d(t)),





x(tk−1 ) ∈ [xk−1 ], d(t) ∈ [d],




 σ(t) = i , ∀t ∈ [t
, t ], t = kτ
k

3

k−1

k

k

Thereby, the enclosure P ostπS
([x0 ]) is included in [xk ]
and T ubeπ ([x0 ]) is included in j=1,..,k [x̃j ]. This applies
for all initial states in [x0 ] and all disturbances d(t) ∈ [d].
A view of enclosures computed by the validated simulation for one solution obtained for Example 5.2 is shown
in Figure 1.

conditions, or fails when the maximal degree of decomposition D is reached. The algorithmic form of functions
Decomposition and F ind P attern is given in Figures 3
(cf. 4 and in [14,20] for the linear case).
π

π1

Post π (V 1)

X

1

Post π (R)

(a)

Tube π ( X )

(b)

V1 V2

R

R

Fig. 2. Principle of the bisection method.
Post π ( X )

Having defined the control synthesis method, we now
introduce the main result of this paper, stated as follows:

Proposition 1 The algorithm of Figure 3 with input
(R, R, S, B, D, K) outputs, when it successfully terminates, a decomposition {Vi , πi }i∈I of R which solves
Problem 1.

Fig. 1. Functions P ostπ (X) and T ubeπ (X) for the initial box
X = [−0.69, −0.64] × [1, 1.06], with a pattern π = (1, 3, 0).

4

The state-space bisection algorithm

4.1

Proof 1 Let x0 = x(t0 = 0) be an initial condition belonging to R. If the
S decomposition has terminated successfully, we have i∈I Vi = R, and x0 thus belongs to
Vi0 for some i0 ∈ I. We can thus apply the pattern πi0
associated to Vi0 . Let us denote by k0 the length of πi0 .
We have:

Principle of the algorithm

We describe here the algorithm solving the control synthesis problem (see Problem 1, Section 2). Given the input boxes R, S, B, and given two positive integers K
and D, the algorithm provides, when it succeeds, a decomposition ∆ of R of the form {Vi , πi }i∈I , with the
properties:
S

i∈I

• x(k0 τ ; 0, x0 , d, πi0 ) ∈ R
• ∀t ∈ [0, k0 τ ], x(t; 0, x0 , d, πi0 ) ∈ S
• ∀t ∈ [0, k0 τ ], x(t; 0, x0 , d, πi0 ) ∈
/B

Vi = R,

Let x1 = x(k0 τ ; 0, x0 , d, πi0 ) ∈ R be the state reached
after application of πi0 and let t1 = k0 τ . State x1 belongs
to R, it thus belongs to Vi1 for some i1 ∈ I, and we can
apply the associated pattern πi1 of length k1 , leading to:

∀i ∈ I, P ostπi (Vi ) ⊆ R,
∀i ∈ I, T ubeπi (Vi ) ⊆ S,
∀i ∈ I, T ubeπi (Vi )

T

B = ∅.

• x(t1 + k1 τ ; t1 , x1 , d, πi1 ) ∈ R
• ∀t ∈ [t1 , t1 + k1 τ ], x(t; t1 , x1 , d, πi1 ) ∈ S
• ∀t ∈ [t1 , t1 + k1 τ ], x(t; t1 , x1 , d, πi1 ) ∈
/B

The sub-boxes {Vi }i∈I are obtained by repeated bisection. At first, function Decomposition calls sub-function
F ind P attern which looks for a pattern π of length at
most K such
T that P ostπ (R) ⊆ R, T ubeπ (R) ⊆ S and
T ubeπ (R) B = ∅. If such a pattern π is found, then a
uniform control over R is found (see Figure 2(a)). Otherwise, R is divided into two sub-boxes V1 , V2 , by bisecting R w.r.t. its longest dimension. Patterns are then
searched to control these sub-boxes (see Figure 2(b)).
If for each Vi , function F ind P attern manages to get a
pattern πi of length at most K verifying
P ostπi (Vi ) ⊆ R,
T
T ubeπi (Vi ) ⊆ S and T ubeπi (Vi ) B = ∅, then it is done.
If, for some Vj , no such pattern is found, the procedure
is recursively applied to Vj . It ends with success when
every sub-box of R has a pattern verifying the latter

We can then iterate this procedure from the new state
x2 = x(t1 + k1 τ ; t1 , x1 , d, πi1 ) ∈ R. This can be repeated
infinitely, yielding a sequence of points belonging to R
x0 , x1 , x2 , . . . attained at times t0 , t1 , t2 , . . . , at which the
patterns πi0 , πi1 , πi2 , . . . are applied.
We furthermore have that all the trajectories stay
in S and never cross B: ∀t ∈ R+ , ∃k ≥ 0, t ∈
[tk , tk+1 ] and ∀t ∈ [tk , tk+1 ], x(t; tk , xk , d, πik ) ∈
S, x(t; tk , xk , d, πik ) ∈
/ B. The trajectories thus return
infinitely often in R, while always staying in S and never
crossing B.

4

Function: Decomposition(W, R, S, B, D, K)

of all the patterns of growing length (up to K).

Input: A box W , a box R, a box S, a box
B, a degree D of bisection, a length K of input
pattern
Output:h{(Vi , πi )}i , T ruei or h , F alsei

The improved function, denoted here by F ind P attern2,
exploits heuristics to prune the search tree of patterns.
The algorithmic form of F ind P attern2 is given in Figure 5. It relies on a new data structure consisting of a
list of triplets containing:

(π, b) := F ind P attern(W, R, S, B, K)
if b = T rue then
return h{(W, P at)}, T ruei
else
if D = 0 then
return h , F alsei
else
Divide equally W into (W1 , W2 )
for i = 1, 2 do

• An initial box V ⊂ Rn ,
• A current box P ostπ (V ), image of V by the pattern π,
• The associated pattern π.
For any element e of a list of this type, we denote by
e.Yinit the initial box, e.Ycurrent the current box, and
by e.Π the associated pattern. We denote by ecurrent =
takeHead(L) the element on top of a list L (this element
is removed from list L). The function putT ail(·, L) adds
an element at the end of the list L.

(∆i , bi )
:=
Decomposition(Wi , R, S, B, D − 1, K)

end for S
V
return ( i=1,2 ∆i , i=1,2 bi )
end if
end if

Let us suppose one wants to control a box X ⊆ R. The
list L of Figure 5 is used to store the intermediate computations leading to possible solutions (patterns sending X
in R while never crossing B or Rn \ S). It is initialized as
L = {(X, X, ∅)}. First, a testing of all the control modes
is performed (a set simulation starting from W during
time τ is computed for all the modes in U ). The first level
of branches is thus tested exhaustively. If a branch leads
to crossing B or Rn \ S, the branch is cut. Otherwise, either a solution is found or an intermediate state is added
to L. The next level of branches (patterns of length 2) is
then explored from branches that are not cut. And so on
iteratively. At the end, either the tree is explored up to
level K (avoiding the cut branches), or all the branches
have been cut at lower levels. List L is thus of the form
{(X, P ostπi (X), πi )}i∈IX , whereTfor each i ∈ IX we have
P ostπi (X) ⊆ S and T ubeπi (X) B = ∅. Here, IX is the
set of indexes associated to the stored intermediate solutions, |IX | is thus the number of stored intermediate
solutions for the initial box X. The number of stored intermediate solutions grows as the search tree of patterns
is explored, then decreases as solutions are validated,
branches are cut, or the maximal level K is reached.

Fig. 3. Algorithmic form of Function Decomposition.

Function: F ind P attern(W, R, S, B, K)
Input:A box W , a box R, a box S, a box B,
a length K of input pattern
Output:hπ, T ruei or h , F alsei
for i = 1 . . . K do
Π := set of input patterns of length i
while Π is non empty do
Select π in Π
Π := Π \ {π}
if P ostπ (W ) ⊆TR and T ubeπ (W ) ⊆ S
and T ubeπ (W ) B = ∅ then
return hπ, T ruei
end if
end while
end for
return h , F alsei

The storage of the intermediate solutions P ostπi (X) allows to reuse the computations already performed. Even
if the search tree of patterns is visited exhaustively, it
already allows to obtain much better computation times
than with Function F ind P attern.

Fig. 4. Algorithmic form of Function F ind P attern.

Remark 1 Note that it is possible to perform reachability from a set R1 to another set R2 by computing
Decomposition(R1 , R2 , S, B, D, K). The set R1 is thus
decomposed with the objective to send its sub-boxes into
R2 , i.e. for a sub-box V of R1 , patterns π are searched
with the objective P ostπ (V ) ⊆ R2 (see Example 5.2).
4.2

A second list, denoted by Solution in Figure 5, is used
to store the validated patterns associated to X, i.e. a
0 , where for each
list of patterns of the form {πj }j∈IX
T
0
j ∈ IX
we have P ostπj (X) ⊆ R, T ubeπj (X) B =
0
∅ and T ubeπj (X) ⊆ S. Here, IX
is the set of indexes
0
associated the the stored validated solutions, |IX
| is thus
the number of stored validated solutions for the initial
box X. The number of stored validated solutions can
only increase, and we hope that at least one solution is

The research of patterns

We propose here an improvement of the function
F ind P attern given in [9,14,20], which is a naive testing

5

found, otherwise, the initial box X is split in two subboxes.

to [1.55, 2.15] × [1.0, 1.4], while never going out of the
safety set S = [1.54, 2.16] × [0.99, 1.41].

Note that several solutions can be returned by
F ind P attern2, further optimizations could thus be
performed, such as returning the pattern minimizing a
given cost function. In practice, and in the examples
given below, we return the first validated pattern and
stop the computation as soon as it is obtained (see
commented line in Figure 5).

The decomposition was obtained in less than one second
with a maximum length of pattern set to K = 6 and
a maximum bisection depth of D = 3. A simulation is
given in Figure 6.

Compared to [14,20], this new function highly improves
the computation times, even though the complexity of
the two functions is theoretically the same, at most in
O(N K ). A comparison between functions F ind P attern
and F ind P attern2 is given in Section 6.
5

Experimentations

In this section, we apply our approach to different case
studies taken from the literature. Our solver prototype is
written in C++ and based on DynIBEX [1]. The computations times given in the following have been performed
on a 2.80 GHz Intel Core i7-4810MQ CPU with 8 GB
of memory. Note that our algorithm is mono-threaded
so all the experimentation only uses one core to perform
the computations. The results given in this section have
been obtained with Function F ind P attern2.

Fig. 6. Simulation from the initial condition (1.55, 1.4). The
box R is in plain black. The trajectory is plotted within time
for the two state variables on the left, and in the state-space
plane on the right.

5.2
5.1

A polynomial example

A linear example: boost DC-DC converter
We consider the polynomial system taken from [25]:

This linear example is taken from [5] and has already
been treated with the state-space bisection method in a
linear framework in [14].

" #
ẋ1
ẋ2

The system is a boost DC-DC converter with one switching cell. There are two switching modes depending on
the position of the switching cell. The dynamics is given
by the equation ẋ(t) = Aσ(t) x(t) + Bσ(t) with σ(t) ∈
U = {1, 2}. The two modes are given by the matrices:

A1 =

A2 =

!

− xrll

0

0

1
− x1c r0 +r
c

− x1l (rl +

r0 .rc
r0 +rc )

r0
1
xc r0 +rc

B1 =

0
− x1l r0r+r
c
0
− x1c r0r+r
c

vs
xl

0

B2 =

vs
xl

.

(7)

x1 + u2 + d2

The control inputs are given by u = (u1 , u2 ) =
Kσ(t) (x1 , x2 ), σ(t) ∈ U = {1, 2, 3, 4}, which correspond to four different state feedback controllers
K1 (x) = (0, −x22 + 2), K2 (x) = (0, −x2 ), K3 (x) =
(2, 10), K4 (x) = (−1.5, 10). We thus have four switching modes. The disturbance d = (d1 , d2 ) lies in
[−0.005, 0.005] × [−0.005, 0.005]. The objective is to
visit infinitely often two zones R1 and R2 , without going out of a safety zone S, and while never crossing a
forbidden zone B. Two decompositions are performed:

!

!

=

"
#
−x2 − 1.5x1 − 0.5x31 + u1 + d1

• a decomposition of R1 which returns {(Vi , πi )}i∈I1
with:
S
i∈I1 Vi = R1 ,
∀i ∈ I1 , P ostπi (Vi ) ⊆ R2 ,
∀i ∈ I1 , T ubeπi (Vi ) T
⊆ S,
∀i ∈ I1 , T ubeπi (Vi ) B = ∅.
• a decomposition of R2 which returns {(Vi , πi )}i∈I2
with:
S
i∈I2 Vi = R2 ,

!

0

with xc = 70, xl = 3, rc = 0.005, rl = 0.05, r0 = 1,
vs = 1. The sampling period is τ = 0.5. The parameters
are exact and there is no perturbation. We want the
state to return infinitely often to the region R, set here

6

Function: F ind P attern2(W, R, S, B, K)
Input:A box W , a box R, a box S, a box B, a length K of input pattern
Output:hπ, T ruei or h , F alsei
Solution = {∅}
L = {(W, W, ∅)}
while L =
6 ∅ do
ecurrent = takeHead(L)
for i ∈ U do
T
if P osti (ecurrent .Ycurrent )
⊆
R and T ubei (ecurrent .Ycurrent ) B
=
∅ and
T ubei (ecurrent .Ycurrent ) ⊆ S then
putTail(Solution, ecurrent .Π + i)
/*can be replaced by: “return hecurrent .Π + i, T ruei” */
else
T
if T ubei (ecurrent .Ycurrent ) B 6= ∅ or T ubei (ecurrent .Ycurrent ) * S then
discard ecurrent
end if
else
T
if T ubei (ecurrent .Ycurrent ) B = ∅ and T ubei (ecurrent .Ycurrent ) ⊆ S then
if Length(Π) + 1 < K then
putTail(L, (ecurrent .Yinit , P osti (ecurrent .Ycurrent ), ecurrent .Π + i))
end if
end if
end if
end for
end while
return h , F alsei if no solution is found, or hπ, T ruei, π being any pattern validated in Solution.
Fig. 5. Algorithmic form of Function F ind P attern2.

∀i ∈ I2 , P ostπi (Vi ) ⊆ R1 ,
∀i ∈ I2 , T ubeπi (Vi ) T
⊆ S,
∀i ∈ I2 , T ubeπi (Vi ) B = ∅.
The input boxes are the following:
R1 = [−0.5, 0.5] × [−0.75, 0.0],
R2 = [−1.0, 0.65] × [0.75, 1.75],
S = [−2.0, 2.0] × [−1.5, 3.0],
B = [0.1, 1.0] × [0.15, 0.5].
The sampling period is set to τ = 0.15. The decompositions were obtained in 2 minutes and 30 seconds with a
maximum length of pattern set to K = 12 and a maximum bisection depth of D = 5. A simulation is given in
Figure 7 in which the disturbance d is chosen randomly
in [−0.005, 0.005] × [−0.005, 0.005] at every time step.
5.3

Fig. 7. Simulation from the initial condition (0.5, −0.75). The
trajectory is plotted within time on the left, and in the state
space plane on the right. In the sate space plane, the set R1
is in plain green, R2 in plain blue, and B in plain black.

Building ventilation

We consider a building ventilation application adapted
from [27]. The system is a four room apartment subject
to heat transfer between the rooms, with the external

environment, with the underfloor, and with human beings. The dynamics of the system is given by the follow-

7

ing equation:

The decomposition was obtained in 4 minutes with a
maximum length of pattern set to K = 2 and a maximum bisection depth of D = 4. The perturbation due
to human beings has been taken into account by setting
the parameters δsi equal to the whole interval [0, 1] for
the decomposition, and the imposed perturbation for the
simulation is given Figure 8. The temperatures To and
Tc have been set to the interval [27, 30] for the decomposition, and are set to 30◦ C for the simulation. A simulation of the controller obtained with the state-space
bisection procedure is given in Figure 9, where the control objective is to stabilize the temperature in [20, 22]4
while never going out of [19, 23]4 .

X
dTi
aij (Tj − Ti ) + δsi bi (Ts4i − Ti4 )
=
dt
*
j∈N


Vi − Vi*
(Tu − Ti ). (8)
+ ci max 0,
V̄i − Vi*
The state of the system is given by the temperatures in
the rooms Ti , for i ∈ N = {1, . . . , 4}. Room i is subject
to heat exchange with different entities stated by the
indexes N * = {1, 2, 3, 4, u, o, c}.
The heat transfer between the rooms is given by the coefficients aij for i, j ∈ N 2 , and the different perturbations are the following:
• The external environment: it has an effect on room i
with the coefficient aio and the outside temperature
To , varying between 27◦ C and 30◦ C.
• The heat transfer through the ceiling: it has an effect
on room i with the coefficient aic and the ceiling temperature Tc , varying between 27◦ C and 30◦ C.
• The heat transfer with the underfloor: it is given by
the coefficient aiu and the underfloor temperature Tu ,
set to 17◦ C (Tu is constant, regulated by a PID controller).
• The perturbation induced by the presence of humans:
it is given in room i by the term δsi bi (Ts4i − Ti4 ), the
parameter δsi is equal to 1 when someone is present
in room i, 0 otherwise, and Tsi is a given identified
parameter.

Fig. 8. Perturbation (presence of humans) imposed within
time in the different rooms.

The control Vi , i ∈ N , is applied through the term
V −V *
ci max(0, V¯i −Vi* )(Tu −Ti ). A voltage Vi is applied to force
i
i
ventilation from the underfloor to room i, and the command of an underfloor fan is subject to a dry friction.
Because we work in a switched control framework, Vi can
take only discrete values, which removes the problem of
dealing with a “max” function in interval analysis. In the
experiment, V1 and V4 can take the values 0V or 3.5V,
and V2 and V3 can take the values 0V or 3V. This leads
to a system of the form (1) with σ(t) ∈ U = {1, . . . , 16},
the 16 switching modes corresponding to the different
possible combinations of voltages Vi . The sampling period is τ = 10s.

Fig. 9. Simulation from the initial condition (22, 22, 22, 22).
The objective set R is in plain black and the safety set S is
in dotted black.

The parameters Tsi , Vi* , V̄i , aij , bi , ci are given in [27]
and have been identified with a proper identification procedure detailed
in [28]. Note that here we have neglected
P
the term j∈N δdij ci,j ∗ h(Tj − Ti ) of [27], representing
the perturbation induced by the open or closed state of
the doors between the rooms. Taking a “max” function
into account with interval analysis is actually still a difficult task. However, this term could have been taken
into account with a proper regularization (smoothing).

6

Performance tests

We present a comparison of the computation times obtained with functions F ind P attern, F ind P attern2,
and with the state-of-the-art tools PESSOA [26] and
SCOTS [32].

8

Table 1
Comparison of F ind P attern and F ind P attern2.
Example

tion times). Another advantage is that we do not require
a homogeneous discretization of the state space. We can
thus tile large parts of the state-space using only few
boxes, and this often permits to consider much less symbolic states than with discretization methods, especially
in high dimensions (see [22]).

Computation time
F ind P attern

F ind P attern2

DC-DC Converter

1609 s

<1s

Polynomial example

Time Out

150 s

Building ventilation

272 s

228 s

7

Table 2
Comparison with state-of-the-art tools.
Example

We presented a method of control synthesis for nonlinear switched systems, based on a simple state-space bisection algorithm, and on validated simulation. The approach permits to deal with stability, reachability, safety
and forbidden region constraints. Varying parameters
and perturbations can be easily taken into account with
interval analysis. The approach has been numerically
validated on several examples taken from the literature,
a linear one with constant parameters, and two nonlinear ones with varying perturbations. Our approach compares well with the state-of-the art tools SCOTS and
PESSOA.

Computation time
FP2

SCOTS

PESSOA

DC-DC Converter

<1s

43 s

760 s

Polynomial example

150 s

131 s

Unicyle [34,31]

3619 s

492 s

Conclusion

516 s

Table 1 shows a comparison of functions F ind P attern
and F ind P attern2, which shows that the new version
highly improves the computation times. We can note
that the new version is all the more efficient as the length
of the patterns increases, and as obstacles cut the research tree of patterns. This is why we observe significant improvements on the examples of the DC-DC converter and the polynomial example, and not on the building ventilation example, which only requires patterns of
length 2, and presents no obstacle.

We would like to point out that the exponential complexity of the algorithms presented here, which is inherent to guaranteed methods, is not prohibitive. Two approaches have indeed been developed to overcome this
exponential complexity. A first approach is the use of
compositionality, which permits to split the system in
two (or more) sub-systems, and to perform control synthesis on these sub-systems of lower dimensions. This
approach has been successfully applied in [22] to a system of dimension 11, and we are currently working on
applying this approach to the more general context of
contract-based design [33]. A second approach is the use
of Model Order Reduction, which allows to approximate
the full-order system (1) with a reduced-order system, of
lower dimension, on which it is possible to perform control synthesis. The bounding of the trajectory errors between the full-order and the reduced-order systems can
be taken into account, so that the induced controller is
guaranteed. This approach, described in [21], has been
successfully applied on (space-discretized) partial differential equations, leading to systems of ODEs of dimension up to 100000. The present work is a potential ground
for the application of such methods to control of nonlinear partial differential equations, with the use of proper
nonlinear model order reduction techniques.

Table 2 shows of comparison of function F ind P attern2
with state-of-the-art tools SCOTS and PESSOA. On the
example of the DC-DC converter, our algorithm manages to control the whole state-space R = [1.55, 2.15] ×
[1.0, 1.4] in less than one second, while SCOTS and PESSOA only control a part of R, and with greater computation times. Note that these computation times vary with
the number of discretization points used in both, but
even with a very fine discretization, we never managed
to control the whole box R. For the polynomial example, we manage to control the whole boxes R1 and R2 ,
such as SCOTS and in a comparable amount of time.
However, PESSOA does not support natively this kind
of nonlinear systems. We compared our method on a
last case study on which PESSOA and SCOTS perform
well (see [34,31] for details of this case study, and see
Appendix for a simulation obtained using our method).
For this case study, we have not obtained as good computations times as they have. This comes from the fact
that this example requires a high number of switched
modes, long patterns, as well as a high number of boxes
to tile the state-space. Note that for this case study we
used an automated pre-tiling of the state-space permitting to decompose the reachability problem in a sequence
of reachability problems. This is in fact the most difficult case of application of our method. This reveals that
our method is more adapted when either the number
of switched modes of the length of patterns is not high
(though it can be handled at the cost of high computa-

Acknowledgements
This work is supported by Institut Farman (project
SWITCHDESIGN), by the French National Research
Agency through the “iCODE Institute project” funded
by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02, and
by Labex DigiCosme (project ANR-11-LABEX-0045DIGICOSME).

9

References

[17] Antoine Girard, Giordano Pola, and Paulo Tabuada.
Approximately bisimilar symbolic models for incrementally
stable switched systems. IEEE Transactions on Automatic
Control, 55(1):116–126, 2010.

[1] Julien Alexandre dit Sandretto and Alexandre Chapoutot.
Dynibex
library.
http://perso.ensta-paristech.fr/ chapoutot/dynibex/, 2015.

[18] Ernst Hairer, Syvert Paul Norsett, and Grehard Wanner.
Solving Ordinary Differential Equations I: Nonstiff Problems.
Springer-Verlag, 2nd edition, 2009.

[2] Julien Alexandre dit Sandretto and Alexandre Chapoutot.
Validated Solution of Initial Value Problem for Ordinary
Differential Equations based on Explicit and Implicit RungeKutta Schemes. Research report, ENSTA ParisTech, 2015.

[19] Fabian Immler. Verified reachability analysis of continuous
systems. In Tools and Algorithms for the Construction and
Analysis of Systems, pages 37–51. Springer, 2015.

[3] Julien Alexandre dit Sandretto and Alexandre Chapoutot.
Validated explicit and implicit runge-kutta methods. Reliable
Computing, 22:79–103, 2016.

[20] Ulrich Kühne and Romain Soulat.
Minimator 1.0.
https://bitbucket.org/ukuehne/minimator/overview, 2015.

[4] Matthias Althoff. Reachability analysis of nonlinear systems
using conservative polynomialization and non-convex sets. In
Proceedings of the 16th international conference on Hybrid
systems: computation and control, pages 173–182. ACM,
2013.

[21] Adrien Le Coënt, Florian De Vuyst, Christian Rey, Ludovic
Chamoin, and Laurent Fribourg. Control of mechanical
systems using set based methods. International Journal of
Dynamics and Control, pages 1–17, 2016.
[22] Adrien Le Coënt, Laurent Fribourg, Nicolas Markey, Florian
De Vuyst, and Ludovic Chamoin. Distributed Synthesis of
State-Dependent Switching Control, pages 119–133. Springer
International Publishing, Cham, 2016.

[5] A Giovanni Beccuti, Georgios Papafotiou, and Manfred
Morari. Optimal control of the boost dc-dc converter. In
Decision and Control, 2005 and 2005 European Control
Conference. CDC-ECC’05. 44th IEEE Conference on, pages
4457–4462. IEEE, 2005.

[23] Daniel Liberzon. Switching in systems and control. Springer
Science & Business Media, 2012.

[6] Olivier Bouissou, Alexandre Chapoutot, and Adel Djoudi.
Enclosing temporal evolution of dynamical systems using
numerical methods. In NASA Formal Methods, number 7871
in LNCS, pages 108–123. Springer, 2013.

[24] Youdong Lin and Mark A. Stadtherr. Validated solutions
of initial value problems for parametric odes. Appl. Numer.
Math., 57(10):1145–1162, 2007.

[7] Olivier Bouissou and Matthieu Martel.
GRKLib: a
Guaranteed Runge Kutta Library. In Scientific Computing,
Computer Arithmetic and Validated Numerics, 2006.

[25] Jun Liu, Necmiye Ozay, Ufuk Topcu, and Richard M Murray.
Synthesis of reactive switching protocols from temporal logic
specifications. Automatic Control, IEEE Transactions on,
58(7):1771–1785, 2013.

[8] Xin Chen, Erika Ábrahám, and Sriram Sankaranarayanan.
Flow*: An analyzer for non-linear hybrid systems. In
Computer Aided Verification, pages 258–263. Springer, 2013.

[26] Manuel Mazo, Anna Davitian, and Paulo Tabuada. PESSOA:
A Tool for Embedded Controller Synthesis, pages 566–569.
Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.

[9] A. Le Coent, J. A. dit Sandretto, A. Chapoutot, and
L. Fribourg. Control of nonlinear switched systems based
on validated simulation. In 2016 International Workshop on
Symbolic and Numerical Methods for Reachability Analysis
(SNR), pages 1–6, April 2016.

[27] Pierre-Jean Meyer. Invariance and symbolic control of
cooperative systems for temperature regulation in intelligent
buildings. Theses, Université Grenoble Alpes, September
2015.
[28] Pierre-Jean Meyer, Hosein Nazarpour, Antoine Girard,
and Emmanuel Witrant. Experimental implementation of
UFAD regulation based on robust controlled invariance. In
Proceedings of the 13th European Control Conference, pages
1468–1473, 2014.

[10] L. H. de Figueiredo and J. Stolfi.
Self-Validated
Numerical Methods and Applications. Brazilian Mathematics
Colloquium monographs. IMPA/CNPq, 1997.
[11] Tomáš Dzetkulič. Rigorous integration of non-linear ordinary
differential equations in Chebyshev basis.
Numerical
Algorithms, 69(1):183–205, 2015.

[29] Ramon Moore. Interval Analysis. Prentice Hall, 1966.

[12] Andreas Eggers, Martin Fränzle, and Christian Herde. Sat
modulo ode: A direct sat approach to hybrid systems. In
Automated Technology for Verification and Analysis, pages
171–185. Springer, 2008.

[30] Nedialko S. Nedialkov, K. Jackson, and Georges Corliss.
Validated solutions of initial value problems for ordinary
differential equations. Appl. Math. and Comp., 105(1):21 –
68, 1999.

[13] Goran Frehse, Colas Le Guernic, Alexandre Donzé, Scott
Cotton, Rajarshi Ray, Olivier Lebeltel, Rodolfo Ripado,
Antoine Girard, Thao Dang, and Oded Maler. Spaceex:
Scalable verification of hybrid systems. In Computer Aided
Verification, pages 379–395. Springer, 2011.

[31] Gunther Reissig, Alexander Weber, and Matthias Rungger.
Feedback refinement relations for the synthesis of symbolic
controllers. arXiv preprint arXiv:1503.03715, 2015.
[32] Matthias Rungger and Majid Zamani. Scots: A tool for the
synthesis of symbolic controllers. In Proceedings of the 19th
International Conference on Hybrid Systems: Computation
and Control, HSCC ’16, pages 99–104, New York, NY, USA,
2016. ACM.

[14] Laurent Fribourg, Ulrich Kühne, and Romain Soulat. Finite
controlled invariants for sampled switched systems. Formal
Methods in System Design, 45(3):303–329, 2014.

[33] Alberto Sangiovanni-Vincentelli, Werner Damm, and
Roberto Passerone. Taming dr. frankenstein: Contract-based
design for cyber-physical systems. European journal of
control, 18(3):217–238, 2012.

[15] Karol Gajda, Malgorzata Jankowska, Andrzej Marciniak, and
Barbara Szyszka. A survey of interval Runge–Kutta and
multistep methods for solving the initial value problem. In
Parallel Processing and Applied Mathematics, volume 4967
of LNCS, pages 1361–1371. Springer Berlin Heidelberg, 2008.

[34] Majid Zamani, Giordano Pola, Manuel Mazo, and Paulo
Tabuada. Symbolic models for nonlinear control systems
without stability assumptions.
IEEE Transactions on
Automatic Control, 57(7):1804–1809, 2012.

[16] Antoine Girard. Reachability of uncertain linear systems
using zonotopes. In Hybrid Systems: Computation and
Control, pages 291–305. Springer, 2005.

10

8

Appendix

Fig. 10. Simulation of the unicycle example.

11

