Control Synthesis of Nonlinear Sampled Switched Systems
using Euler’s Method
A. Le Coënt∗ F. De Vuyst

L. Chamoin

CMLA, CNRS & ENS Paris-Saclay

LMT, CNRS & ENS Paris-Saclay

lecoent@cmla.ens-cachan.fr

L. Fribourg
LSV, CNRS, INRIA & ENS Paris-Saclay

In this paper, we propose a symbolic control synthesis method for nonlinear sampled switched systems whose vector fields are one-sided Lipschitz. The main idea is to use an approximate model
obtained from the forward Euler method to build a guaranteed control. The benefit of this method
is that the error introduced by symbolic modeling is bounded by choosing suitable time and space
discretizations. The method is implemented in the interpreted language Octave. Several examples
of the literature are performed and the results are compared with results obtained with a previous
method based on the Runge-Kutta integration method.

1

Introduction

As said in [10], 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 ordinary differential
equations (ODEs) is fundamental (see, e.g., [2, 14]). An interesting approach appeared recently, based
on the propagation of reachable sets using guaranteed Runge-Kutta methods with adaptive step size control (see [6, 17]). In [10] such guaranteed integration methods are used in the framework of sampled
switched systems.
Given an ODE of the form ẋ(t) = f (t, x(t)), and a set of initial values X0 , a symbolic (or “set-valued”)
integration method consists in computing a sequence of approximations (tn , x̃n ) of the solution x(t; x0 ) of
the ODE with x0 ∈ X0 such that x̃n ≈ x(tn ; xn−1 ). Symbolic integration methods extend classical numerical
integration methods which correspond to the case where X0 is just a singleton {x0 }. The simplest numerical method is Euler’s method in which tn+1 = tn + h for some step-size h and x̃n+1 = x̃n + h f (tn , x̃n );
so the derivative of x at time tn , f (tn , xn ), is used as an approximation of the derivative on the whole
time interval. This method is very simple and fast, but requires small step-sizes h. 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 of the form
x̃n+1 = x̃n + hΣsi=1 bi ki where ki = f (tn + ci h, x̃n + hΣi−1
j=1 ai j k j ) for i = 2, 3, ..., s. A challenging question
is then to compute a bound on the distance between the true solution and the numerical solution, i.e.:
kx(tn ; xn−1 ) − xn k. This distance is associated to the local truncation error of the numerical method.
In [10], such a bound is computed using the Lagrange remainders of Taylor expansions. This is achieved
using affine arithmetic [27] (by application of the Banach’s fixpoint theorem and Picard-Lindelöf operator, see [23]). In the end, the Runge-Kutta based method of [10] is an elaborated method that requires
the use of affine arithmetic, Picard iteration and computation of Lagrange remainder.
∗ corresponding

author

E. Ábrahám and S. Bogomolov (Eds.): 3rd International
Workshop on Symbolic and Numerical Methods for Reachability
Analysis (SNR 2017)
EPTCS 247, 2017, pp. 18–33, doi:10.4204/EPTCS.247.2

c A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg
This work is licensed under the
Creative Commons Attribution License.

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

19

In contrast, in this paper, we use ordinary arithmetic (instead of affine arithmetic) and a basic Euler
scheme (instead of Runge-Kutta schemes). We need neither estimate Lagrange remainders nor perform
Picard iteration in combination with Taylor series. Our simple Euler-based approach is made possible
by having recourse to the notion of one-sided Lipschitz (OSL) function [12]. This allows us to bound
directly the global error, i.e. the distance between the approximate point x̃(t) computed by the Euler
scheme and the exact solution x(t) for all t ≥ 0 (see Theorem 1).
Plan. In Section 2, we give details on related work. In Section 3, we state our main result that bounds
the global error introduced by the Euler scheme in the context of systems with OSL flows. In Section 4,
we explain how to apply this result to the synthesis of symbolic control of sampled switched systems.
We give numerical experiments and results in Section 5 for five exampes of the literature, and compare
them with results obtained with the method of [10]. We give final remarks in Section 6.

2

Related work

Most of the recent work on the symbolic (or set-valued) integration of nonlinear ODEs is based on
the upper bounding of the Lagrange remainders either in the framework of Taylor series or RungeKutta schemes [2, 5, 6, 8, 9, 10, 20, 24, 26]. Sets of states are generally represented as vectors of
intervals (or “rectangles”) and are manipulated through interval arithmetic [22] or affine arithmetic [27].
Taylor expansions with Lagrange remainders are also used in the work of [2], which uses “polynomial
zonotopes” for representing sets of states in addition to interval vectors. None of these works uses the
Euler scheme nor the notion of one-sided Lipschitz constant.
In the literature on symbolic integration, the Euler scheme with OSL conditions is explored in [12,
18]. Our approach is similar but establishes an analytical result for the global error of Euler’s estimate
(see Theorem 1) rather than analyzing, in terms of complexity, the speed of convergence to zero, the
accuracy and the stability of Euler’s method.
In the control literature, OSL conditions have been recently applied to control and stabilization [1, 7],
but do not make use of Euler’s method. To our knowledge, our work applies for the first time Euler’s
scheme with OSL conditions to the symbolic control of hybrid systems.

3
3.1

Sampled switched systems with OSL conditions
Control of switched systems

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

(1)

defined for all t ≥ 0, where x(t) ∈ Rn is the state of the system, σ (·) : R+ −→ U is the switching rule.
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 constant on the time interval [(k − 1)τ, kτ) for k ≥ 1. For all j ∈ U, f j is a function from Rn to Rn .
We make the following hypothesis:
(H0)

For all j ∈ U, f j is a locally Lipschitz continuous map.

20

Control synthesis of switched systems using Euler method

As in [15], we make the assumption that the vector field f j is such that the solutions of the differential
equation (1) are defined, e.g. by assuming that the support of the vector field f j is compact. We will
denote by φσ (t; x0 ) the solution at time t of the system:
ẋ(t) = fσ (t) (x(t)),

(2)

x(0) = x0 .

Often, we will consider φσ (t; x0 ) on the interval 0 ≤ t < τ for which σ (t) is equal to a constant,
say j ∈ U. In this case, we will abbreviate φσ (t; x0 ) as φ j (t; x0 ). We will also consider φσ (t; x0 ) on
the interval 0 ≤ t < kτ where k is a positive integer, and σ (t) is equal to a constant, say jk0 , on each
interval [(k0 − 1)τ, k0 τ) with 1 ≤ k0 ≤ k; in this case, we will abbreviate φσ (t; x0 ) as φπ (t; x0 ), where π is
a sequence of k modes (or “pattern”) of the form π = j1 · j2 · · · · · jk .
We will assume that φσ is continuous at time kτ for all positive integer k. This means that there is no
“reset” at time k0 τ (1 ≤ k0 ≤ k); the value of φσ (t, x0 ) for t ∈ [(k0 − 1)τ, kτ] corresponds to the solution of
ẋ(u) = f jk0 (x(u)) for u ∈ [0, τ] with initial value φσ ((k0 − 1)τ; x0 ).
Given a “recurrence set” R ⊂ Rn and a “safety set” S ⊂ Rn which contains R (R ⊆ S), we are interested
in the synthesis of a control such that: starting from any initial point x ∈ R, the controlled trajectory
always returns to R within a bounded time while never leaving S. We suppose that sets R and S are
compact. Furthermore, we suppose that S is convex.
We denote by T a compact overapproximation of the image by φ j of S for 0 ≤ t ≤ τ and j ∈ U, i.e. T
is such that
T ⊇ {φ j (t; x0 ) | j ∈ U, 0 ≤ t ≤ τ, x0 ∈ S}.
The existence of T is guaranteed by assumption (H0). We know furthermore by (H0) that, for all j ∈ U,
there exists a constant L j > 0 such that:
k f j (y) − f j (x)k ≤ L j ky − xk

∀x, y ∈ S.

(3)

j ∈ U.

(4)

Let us define C j for all j ∈ U:
C j = sup L j k f j (x)k

for all

x∈S

We make the additional hypothesis that the mappings f j are one-sided Lipschitz (OSL) [12]. Formally:
(H1)

For all j ∈ U, there exists a constant λ j ∈ R such that
h f j (y) − f j (x), y − xi ≤ λ j ky − xk2

∀x, y ∈ T,

where h·, ·i denotes the scalar product of two vectors of Rn .
Remark 1. Constants λ j , L j and C j ( j ∈ U) can be computed using (constrained) optimization algorithms. See Section 5 for details.

3.2

Euler approximate solutions

Given an initial point x̃0 ∈ S and a mode j ∈ U, we define the following “linear approximate solution” φ̃ j (t; x̃0 ) for t on [0, τ] by:
φ̃ j (t; x̃0 ) = x̃0 + t f j (x̃0 ).
(5)

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

21

Note that formula (5) is nothing else but the explicit forward Euler scheme with “time step” t. It is thus
a consistent approximation of order 1 in t of the exact solution of (1) under the hypothesis x̃0 = x0 .
More generally, given an initial point x̃0 ∈ S and pattern π of U k , we can define a “(piecewise linear)
approximate solution” φ̃π (t; x̃0 ) of φπ at time t ∈ [0, kτ] as follows:
• φ̃π (t; x̃0 ) = t f j (x̃0 ) + x̃0 if π = j ∈ U, k = 1 and t ∈ [0, τ], and
• φ̃π (kτ + t; x̃0 ) = t f j (z̃) + z̃ with z̃ = φ̃π 0 ((k − 1)τ; x̃0 ), if k ≥ 2, t ∈ [0, τ], π = j · π 0 for some j ∈ U
and π 0 ∈ U k−1 .
We wish to synthesize a guaranteed control σ for φσ using the approximate functions φ̃π .We define
the closed ball of center x ∈ Rn and radius r > 0, denoted B(x, r), as the set {x0 ∈ Rn | kx0 − xk ≤ r}.
Given a positive real δ , we now define the expression δ j (t) which, as we will see in Theorem 1,
represents (an upper bound on) the error associated to φ̃ j (t; x̃0 ) (i.e. kφ̃ j (t; x̃0 ) − φ j (t; x0 )k).
Definition 1. Let δ be a positive constant. Let us define, for all 0 ≤ t ≤ τ, δ j (t) as follows:
• if λ j < 0:
δ j (t) =

2 λ jt

δ e

+

C2j
λ j2

!! 12

2t
2 
λ jt
t + + 2 1−e
λj λj
2

• if λ j = 0 :
1
δ j (t) = δ 2 et +C2j (−t 2 − 2t + 2(et − 1)) 2
• if λ j > 0 :
δ j (t) =

δ 2 e3λ j t +

C2j
3λ j2

!! 12


2t
2
e3λ j t − 1
−t 2 −
+
3λ j 9λ j2

Note that δ j (t) = δ for t = 0. The function δ j (·) depends implicitly on two parameters: δ ∈ R and
j ∈ U. In Section 4, we will use the notation δ j0 (·) where the parameters are denoted by δ 0 and j.
Theorem 1. Given a sampled switched system satisfying (H0-H1), consider a point x̃0 and a positive
real δ . We have, for all x0 ∈ B(x̃0 , δ ), t ∈ [0, τ] and j ∈ U:
φ j (t; x0 ) ∈ B(φ̃ j (t; x̃0 ), δ j (t)).
Proof. Consider on t ∈ [0, τ] the differential equations
dx(t)
= f j (x(t))
dt
and
d x̃(t)
= f j (x̃0 ).
dt
with initial points x0 ∈ S, x̃0 ∈ S respectively. We will abbreviate φ j (t; x0 ) (resp. φ̃ j (t; x̃0 )) as x(t)
(resp. x̃(t)). We have

d
(x(t) − x̃(t)) = f j (x(t)) − f j (x̃0 ) ,
dt

22

Control synthesis of switched systems using Euler method

then
1 d
(kx(t) − x̃(t)k2 ) =
2 dt
=

f j (x(t)) − f j (x̃0 ), x(t) − x̃(t)
f j (x(t)) − f j (x̃(t)) + f j (x̃(t)) − f j (x̃0 ), x(t) − x̃(t)

=

f j (x(t)) − f j (x̃(t)), x(t) − x̃(t) + f j (x̃(t)) − f j (x̃0 ), x(t) − x̃(t)

≤

f j (x(t)) − f j (x̃(t)), x(t) − x̃(t) + k f j (x̃(t)) − f j (x̃0 )kkx(t) − x̃(t)k.

The last expression has been obtained using the Cauchy-Schwarz inequality. Using (H1) and (3), we
have
1 d
(kx(t) − x̃(t)k2 ) ≤ λ j kx(t) − x̃(t)k2 + k f j (x̃(t)) − f j (x̃0 )k kx(t) − x̃(t)k
2 dt
≤ λ j kx(t) − x̃(t)k2 + L j kx̃(t) − x̃0 k kx(t) − x̃(t)k
≤ λ j kx(t) − x̃(t)k2 + L j t k f j (x̃0 )k kx(t) − x̃(t)k.

Using (4) and a Young inequality, we then have
1 d
(kx(t) − x̃(t)k2 ) ≤ λ j kx(t) − x̃(t)k2 +C j t kx(t) − x̃(t)k
2 dt


1
1
2
2
αkx(t) − x̃(t)k +
≤ λ j kx(t) − x̃(t)k +C j t
2
α
for all α > 0.
• In the case λ j < 0:
λ

For t > 0, we choose α > 0 such that C j tα = −λ j , i.e. α = − C jjt . It follows, for all t ∈ [0, τ]:
λj
C jt
λj
(C j t)2
1 d
(kx(t) − x̃(t)k2 ) ≤ kx(t) − x̃(t)k2 −
= kx(t) − x̃(t)k2 −
.
2 dt
2
2α
2
2λ j
We thus get:
2

0

0 2 λ jt

kx(t) − x̃(t)k ≤ kx − x̃ k e

+

C2j
λ j2

!

2t
2 
λ jt
t + + 2 1−e
.
λj λj
2

• In the case λ j > 0:
λ

For t > 0, we choose α > 0 such that C j tα = λ j , i.e. α = C jjt . It follows, for all t ∈ [0, τ]:
3λ j
C jt
3λ j
(C j t)2
1 d
(kx(t) − x̃(t)k2 ) ≤
kx(t) − x̃(t)k2 +
=
kx(t) − x̃(t)k2 +
.
2 dt
2
2α
2
2λ j
We thus get:
2

0

0 2 3λ j t

kx(t) − x̃(t)k ≤ kx − x̃ k e

+

C2j
3λ j2

!

2t
2  3λ j t
−t −
+
e −1
3λ j 9λ j2
2

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

23

• In the case λ j = 0:
For t > 0, we choose α = C1j t . It follows:
d
(kx(t) − x̃(t)k2 ) ≤ kx(t) − x̃(t)k2 +C j t 2
dt
We thus get:
kx(t) − x̃(t)k2 ≤ kx0 − x̃0 k2 et +C2j (−t 2 − 2t + 2(et − 1))
In every case, since by hypothesis x0 ∈ B(x̃0 , δ ) (i.e. kx0 − x̃0 k2 ≤ δ 2 ), we have, for all t ∈ [0, τ]:
kx(t) − x̃(t)k ≤ δ j (t).
It follows: φ j (t; x0 ) ∈ B(φ̃ j (t; x̃0 ), δ ) for t ∈ [0, τ].

Remark 2. In Theorem 1, we have supposed that the step size h used in Euler’s method was equal to the
sampling period τ of the switching system. Actually, in order to have better approximations, it is often
τ
). Such a splitting is called “sub-sampling” in
convenient to take a fraction of τ as for h (e.g., h = 10
numerical methods. See Section 5 for details.
Corollary 1. Given a sampled switched system satisfying (H0-H1), consider a point x̃0 ∈ S, a real δ > 0
and a mode j ∈ U such that:
1. B(x̃0 , δ ) ⊆ S,
2. B(φ̃ j (τ; x̃0 ), δ j (τ)) ⊆ S, and
d 2 (δ (t))

j
3.
> 0 for all t ∈ [0, τ].
dt 2
Then we have, for all x0 ∈ B(x̃0 , δ ) and t ∈ [0, τ]: φ j (t; x0 ) ∈ S.

Proof. By items 1 and 2, B(φ̃ j (t; x̃0 ), δ j (t)) for t = 0 and t = τ. Since δ j (·) is convex on [0, τ] by
item 3, and S is convex, we have B(φ̃ j (t; x̃0 ), δ j (t)) ⊆ S for all t ∈ [0, τ]. It follows from Theorem 1 that
φ j (t; x0 ) ∈ B(φ̃ j (t; x̃0 ), δ j (t)) ⊆ S for all 1 ≤ t ≤ τ.
Remark 3. Condition 3 of Corollary 1 on the convexity of δ j (·) on [0, τ] can be established again using
an optimization function (see Section 5).

4

Application to control synthesis

Consider a point x̃0 ∈ S, a positive real δ and a pattern π of length k. Let π(k0 ) denote the k0 -th element
0
(mode) of π for 1 ≤ k0 ≤ k. Let us abbreviate the k0 -th approximate point φ̃π (k0 τ; x̃0 ) as x̃πk for k0 = 1, ..., k,
0
0
and let x̃πk = x̃0 for k0 = 0. It is easy to show that x̃πk can be defined recursively for k0 = 1, ..., k, by:
0
0
0
x̃πk = x̃πk −1 + τ f j (x̃πk −1 ) with j = π(k0 ).
0
0
0
Let us now denote by δπk (an upper bound on) the error associated to x̃πk , i.e. kx̃πk − φπ (k0 τ; x0 )k.
0
Using repeatedly Theorem 1, δπk can be defined recursively as follows:
0
0
0
For k0 = 0: δπk = δ , and for 1 ≤ k0 ≤ k: δπk = δ j0 (τ) where δ 0 denotes δπk −1 , and j denotes π(k0 ).
Likewise, for 0 ≤ t ≤ kτ, let us denote by δπ (t) (an upper bound on) the global error associated to
φ̃π (t; x̃0 ) (i.e. kφ̃π (t; x̃0 ) − φπ (t; x0 )k). Using Theorem 1, δπ (t) can be defined itself as follows:

24

Control synthesis of switched systems using Euler method
• for t = 0: δπ (t) = δ ,
• for 0 < t ≤ kτ: δπ (t) = δ j0 (t 0 ) with δ 0 = δπ`−1 , j = π(`), t 0 = t − (` − 1)τ and ` = d τt e.
0

Note that, for 0 ≤ k0 ≤ k, we have: δπ (k0 τ) = δπk . We have:
Theorem 2. Given a sampled switched system satisfying (H0-H1), consider a point x̃0 ∈ S, a positive
real δ and a pattern π of length k such that, for all 1 ≤ k0 ≤ k:
0

0

1. B(x̃πk , δπk ) ⊆ S and
2.

d 2 (δ j0 (t))
dt 2

0

> 0 for all t ∈ [0, τ], with j = π(k0 ) and δ 0 = δπk −1 .

Then we have, for all x0 ∈ B(x̃0 , δ ) and t ∈ [0, kτ]: φπ (t; x0 ) ∈ S.
Proof. By induction on k using Corollary 1.
The statement of Theorem 2 is illustrated in Figure 1 for k = 2. From Theorem 2, it easily follows:

B( ~
x ,δ)
0

0
~
x

x

0

B( ~
x , δπ )
1

x

x

1

1

1
~
x

2
2
~
x

B( ~
x , δ π)
2

2

S

Figure 1: Illustration of Theorem 2.
Corollary 2. Given a switched system satisfying (H0-H1), consider a positive real δ and a finite set
of points x̃1 , . . . x̃m of S such that all the balls B(x̃i , δ ) cover R and are included into S (i.e. R ⊆
Sm
i=1 B(x̃i , δ ) ⊆ S). Suppose furthermore that, for all 1 ≤ i ≤ m, there exists a pattern πi of length ki
such that:
0

0

1. B((x̃i )kπi , δπki ) ⊆ S, for all k0 = 1, . . . , ki − 1
2. B((x̃i )kπii , δπkii ) ⊆ R.
3.

d 2 (δ j0 (t))
dt 2

0

> 0 with j = πi (k0 ) and δ 0 = δπki −1 , for all k0 ∈ {1, ..., ki } and t ∈ [0, τ].

These properties induce a control σ 1 which guarantees
1 Given an initial point x ∈ R, the induced control σ corresponds to a sequence of patterns π , π , . . . defined as follows:
i1 i2
Since x ∈ R, there exists a a point x̃i1 with 1 ≤ i1 ≤ m such that x ∈ B(x̃i1 , δ ); then using pattern πi1 , one has: φπi1 (ki1 τ; x) ∈ R.
Let x0 = φπi1 (ki1 τ; x); there exists a point x̃i2 with 1 ≤ i2 ≤ m such that x0 ∈ B(x̃i2 , δ ), etc.

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

25

• (safety): if x ∈ R, then φσ (t; x) ∈ S for all t ≥ 0, and
• (recurrence): if x ∈ R then φσ (kτ; x) ∈ R for some k ∈ {k1 , . . . , km }.
Corollary 2 gives the theoretical foundations of the following method for synthesizing σ ensuring
recurrence in R and safety in S:
• we (pre-)compute λ j , L j ,C j for all j ∈ U;
• we find m points x̃1 , . . . x̃m of S and δ > 0 such that R ⊆

Sm

i=1 B(x̃i , δ )

⊆ S;

• we find m patterns πi (i = 1, ..., m) such that conditions 1-2-3 of Corollary 2 are satisfied.
A covering of R with balls as stated in Corollary 2 is illustrated in Figure 2. The control synthesis method
based on Corollary 2 is illustrated in Figure 3 (left) together with an illustration of method of [10] (right).

~
x1

~
x2

~
x3

~
x4

~
x5

~
x6

~
x7

~
x8

~
x9

R

S

Figure 2: A set of balls covering R and contained in S.

5

Numerical experiments and results

This method has been implemented in the interpreted language Octave, and the experiments performed
on a 2.80 GHz Intel Core i7-4810MQ CPU with 8 GB of memory.
The computation of constants L j , C j , λ j ( j ∈ U) are realized with a constrained optimization algorithm. They are performed using the “sqp” function of Octave, applied on the following optimization
problems:
• Constant L j :
k f j (y) − f j (x)k
L j = max
ky − xk
x,y∈S, x6=y
• Constant C j :
C j = max L j k f j (x)k
x∈S

• Constant λ j :
λj =

h f j (y) − f j (x), y − xi
ky − xk2
x,y∈T, x6=y
max

26

Control synthesis of switched systems using Euler method

Tube i ( Z2 )
1

δ

~
x3

B( ~
x 23 , δ 2π )
3

Z1

Z2
Post i ( Z 2 ')
2

R

R
B( ~
x 13 , δ 1π )

S

3

Z3

Z4
Tube i ( Z 2 ')
2

Z2 '= Post i ( Z 2)
1

S

Figure 3: Control of ball B(x̃3 , δ ) with our method (left); control of tile Z2 with the method of [10] (right).
d 2 (δ 0 (t))

Likewise, the convexity test dtj2 > 0 can be performed similarly.
Note that in some cases, it is advantageous to use a time sub-sampling to compute the image of a
ball. Indeed, because of the exponential growth of the radius δ j (t) within time, computing a sequence of
balls can lead to smaller ball images. It is particularly advantageous when a constant λ j is negative. We
illustrate this with the example of the DC-DC converter. It has two switched modes, for which we have
λ1 = −0.014215 and λ2 = 0.142474. In the case λ j < 0, the associated formula δ j (t) has the behavior of
Figure 4 (a). In the case λ j > 0, the associated formula δ j (t) has the behavior of Figure 4 (b). In the case
λ j < 0, if the time sub-sampling is small enough, one can compute a sequence of balls with reducing
radius, which makes the synthesis easier.

(a)

(b)

Figure 4: Behavior of δ j (t) for the DC-DC converter with δ j (0) = 0.045. (a) Evolution of δ1 (t) (with
λ1 < 0); (b) Evolution of δ2 (t) (with λ2 > 0).
In the following, we give the results obtained with our Octave implementation of this Euler-based

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

27

method on 5 examples, and compare them with those given by the C++ implementation DynIBEX [25]
of the Runge-Kutta based method used in [10].

5.1

Four-room apartment

We describe a first application on a 4-room 16-switch building ventilation case study adapted from [21].
The model has been simplified in order to get constant parameters. The system is a four room apartment
subject to heat transfer between the rooms, with the external environment, the underfloor, and human
beings. The dynamics of the system is given by the following equation:


dTi
Vi −Vi*
4
4
= ∑ ai j (T j − Ti ) + δsi bi (Tsi − Ti ) + ci max 0,
(Tu − Ti ), for i = 1, ..., 4.
dt
V̄i −Vi*
j∈N * \{i}
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 indices N * = {1, 2, 3, 4, u, o, c}. We
have T0 = 30, Tc = 30, Tu = 17, δsi = 1 for i ∈ N . The (constant) parameters Tsi , Vi* , V̄i , ai j , bi , ci
are given in [21]. The control input is Vi (i ∈ N ). 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 τ = 30s. Compared simulations are given in Figure 5. On this
example, the Euler-based method works better than DynIBEX in terms of CPU time.

R
S
τ
Time subsampling
Complete control
max j=1,...,16 λ j
max j=1,...,16 C j
Number of balls/tiles
Pattern length
CPU time

Euler
DynIBEX
2
[20, 22] × [22, 24]2
[19, 23]2 × [21, 25]2
30
No
Yes
Yes
−6.30 × 10−3
4.18 × 10−6
4096
252
1
1
63 seconds
249 seconds

Table 1: Numerical results for the four-room example.

5.2

DC-DC converter

This linear example is taken from [3] and has already been treated with the state-space bisection method
in a linear framework in [13].
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:
!
 vs 
− xrll
0
xl
A1 =
B
=
1
1
0 − x1c r0 +r
0
c

28

Control synthesis of switched systems using Euler method

Figure 5: Simulation of the four-room case study with our synthesis method (left) and with the synthesis
method of [10] (right).

A2 =

.rc
0
) − x1l r0r+r
− x1l (rl + rr00+r
c
c
1 r0
1 r0
−
xc r0 +rc
xc r0 +rc

!

 vs 
B2 =

xl

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 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]. On this example, the Euler-based method fails while DynIBEX succeeds rapidly.

R
S
τ
Complete control
λ1
λ2
C1
C2
Number of balls/tiles
Pattern length
CPU time

Euler
DynIBEX
[1.55, 2.15] × [1.0, 1.4]
[1.54, 2.16] × [0.99, 1.41]
0.5
No
Yes
−0.014215
0.142474
6.7126 × 10−5
2.6229 × 10−2
x
48
x
6
x
< 1 second

Table 2: Numerical results for the DC-DC converter example.

5.3

Polynomial example

We consider the polynomial system taken from [19]:
  

ẋ1
−x2 − 1.5x1 − 0.5x13 + u1
=
.
ẋ2
x1 + u2

(6)

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

29

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 disturbances are not taken into account. The
objective is to visit infinitely often two zones R1 and R2 , without going out of a safety zone S.
R1
R2
S
τ
Time subsampling
Complete control
λ1
λ2
λ3
λ4
C1
C2
C3
C4
Number of balls/tiles
Pattern length
CPU time

Euler
DynIBEX
[−1, 0.65] × [0.75, 1.75]
[−0.5, 0.5] × [−0.75, 0.0]
[−2.0, 2.0] × [−1.5, 3.0]
0.15
τ/20
Yes
Yes
−1.5
−1.0
−1.1992 × 10−8
−5.7336 × 10−6
641.37
138.49
204.50
198.64
16 & 16
1&1
8
7
29 & 4203 seconds <0.1 & 329 seconds

Table 3: Numerical results for the polynomial example example.
For Euler and DynIBEX, the table indicates two CPU times corresponding to the reachability from
R1 to R2 and vice versa. On this example, the Euler-based method is much slower than DynIBEX.

5.4

Two-tank system

The two-tank system is a linear example taken from [16]. The system consists of two tanks and two
valves. The first valve adds to the inflow of tank 1 and the second valve is a drain valve for tank 2. There
is also a constant outflow from tank 2 caused by a pump. The system is linearized at a desired operating
point. The objective is to keep the water level in both tanks within limits using a discrete open/close
switching strategy for the valves. Let the water level of tanks 1 and 2 be given by x1 and x2 respectively.
The behavior of x1 is given by ẋ1 = −x1 − 2 when the tank 1 valve is closed, and ẋ1 = −x1 + 3 when it
is open. Likewise, x2 is driven by ẋ2 = x1 when the tank 2 valve is closed and ẋ2 = x1 − x2 − 5 when it is
open. On this example, the Euler-based method works better than DynIBEX in terms of CPU time.

5.5

Helicopter

The helicopter is a linear example taken from [11]. The problem is to control a quadrotor helicopter
toward a particular position on top of a stationary ground vehicle, while satisfying constraints on the
relative velocity. Let g be the gravitational constant, x (reps. y) the position according to x-axis (resp.
y-axis), ẋ (resp. ẏ) the velocity according to x-axis (resp. y-axis), φ the pitch command and ψ the roll
command. The possible commands for the pitch and the roll are the following: φ , ψ ∈ {−10, 0, 10}.

30

Control synthesis of switched systems using Euler method

R
S
τ
Time subsampling
Complete control
λ1
λ2
λ3
λ4
C1
C2
C3
C4
Number of balls/tiles
Pattern length
CPU time

Euler
DynIBEX
[−1.5, 2.5] × [−0.5, 1.5]
[−3, 3] × [−3, 3]
0.2
τ/10
Yes
Yes
0.20711
-0.50000
0.20711
-0.50000
11.662
28.917
13.416
32.804
64
10
6
6
58 seconds 246 seconds

Table 4: Numerical results for the two-tank example.
Since each mode corresponds to a pair (φ , ψ), there are nine switched modes. The dynamics of the
system is given by the equation:


0
0
Ẋ = 
0
0

1
0
0
0

0
0
0
0




0
0


0
 X + g sin(−φ )



1
0
0
g sin(ψ)

where X = (x ẋ y ẏ)> . Since the variables x and y are decoupled in the equations and follow the same
equations (up to the sign of the command), it suffices to study the control for x (the control for y is the
opposite). On this example again, the Euler-based method works better than DynIBEX in terms of CPU
time.

5.6

Analysis and comparison of results

Our method presents the advantage over the work of [10] that no numerical integration is required for the
control synthesis. The computations just require the evaluation of given functions f j and (global error)
functions δ j at sampling times. The synthesis is thus a priori cheap compared to the use of numerical
integration schemes (and even compared to exact integration for linear systems). However, most of the
computation time is actually taken by the search for an appropriate radius δ of the balls Bi (1 ≤ i ≤ m)
that cover R, and the search for appropriate patterns πi that make the trajectories issued from Bi return
to R.
We observe on the examples that the resulting control strategies synthesized by our method are quite
different from those obtained by the Runge-Kutta method of [10] (which uses in particular rectangular
tiles instead of balls). This may explain why the experimental results are here contrasted: Euler’s method
works better on 3 examples and worse on the 2 others. Besides the Euler method fails on one example

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

R
S
τ
Time subsampling
Complete control
λ1
λ2
λ3
C1
C2
C3
Number of balls/tiles
Pattern length
CPU time

31

Euler
DynIBEX
[−0.3, 0.3] × [−0.5, 0.5]
[−0.4, 0.4] × [−0.7, 0.7]
0.1
τ/10
Yes
Yes
0.5
0.5
0.5
1.77535
0.5
1.77535
256
35
7
7
539 seconds 1412 seconds

Table 5: Numerical results for the helicopter motion example.
(DC-DC converter) while DynIBEX succeeds on all of them. Note however that our Euler-based implementation is made of a few hundreds lines of interpreted code Octave while DynIBEX is made of around
five thousands of compiled code C++.

6

Final Remarks

We have given a new Euler-based method for controlling sampled switched systems, and compared it
with the Runge-Kutta method of [10]. The method is remarkably simple and gives already promising
results. In future work, we plan to explore the use of the backward Euler method instead of the forward
Euler method used here (cf: [4]). We plan also to give general sufficient conditions ensuring the convexity
of the error function δ j (·); this would allow us to get rid of the convexity tests that we perform so far
numerically for each pattern.
Acknowledgement. We are grateful to Antoine Girard, Jonathan Vacher, Julien Alexandre dit Sandretto and Alexandre Chapoutot for numerous helpful discussions. This work has been partially supported by Federative Institute Farman (ENS Paris-Saclay and CNRS FR3311).

References
[1] M. Abbaszadeh & H.J. Marquez (2010): Nonlinear observer design for one-sided Lipschitz systems. In: Proceedings of the American Control Conference (ACC), IEEE, pp. 799–806, doi:10.1109/ACC.2010.5530715.
[2] Matthias Althoff (2013): 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, HSCC ’13, ACM, New York, NY, USA, pp. 173–182, doi:10.1145/2461328.2461358.
[3] A Giovanni Beccuti, Georgios Papafotiou & Manfred Morari (2005): 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, IEEE, pp. 4457–4462, doi:10.1109/CDC.2005.1582864.

32

Control synthesis of switched systems using Euler method

[4] W.-J. Beyn & J. Rieger (1998): The implicit Euler scheme for one-sided Lipschitz differential inclusions.
Discr. and Cont. Dynamical Systems B(14), pp. 409–428, doi:10.3934/dcdsb.2010.14.409.
[5] Olivier Bouissou, Alexandre Chapoutot & Adel Djoudi (2013): Enclosing Temporal Evolution of Dynamical Systems Using Numerical Methods. In: NASA Formal Methods, LNCS 7871, Springer, pp. 108–123,
doi:10.1007/978-3-642-38088-4_8.
[6] Olivier Bouissou, Samuel Mimram & Alexandre Chapoutot (2012): HySon: Set-Based Simulation of Hybrid
Systems. In: Rapid System Prototyping, IEEE, doi:10.1109/RSP.2012.6380694.
[7] Xiushan Cai, Zhenyun Wang & Leipo Liu (2015): Control Design for One-side Lipschitz
Nonlinear Differential Inclusion Systems with Time-delay.
Neurocomput. 165(C), pp. 182–189,
doi:10.1016/j.neucom.2015.03.008.
[8] Xin Chen, Erika Abraham & Sriram Sankaranarayanan (2012): Taylor Model Flowpipe Construction for
Non-linear Hybrid Systems. In: IEEE 33rd Real-Time Systems Symposium, IEEE Computer Society, pp.
183–192, doi:10.1109/RTSS.2012.70.
[9] Xin Chen, Erika Ábrahám & Sriram Sankaranarayanan (2013): Flow*: An analyzer for non-linear hybrid
systems. In: Computer Aided Verification, Springer, pp. 258–263, doi:10.1007/978-3-642-39799-8_18.
[10] A. Le Coënt, J. Alexandre dit Sandretto, A. Chapoutot & L. Fribourg (2016): Control of nonlinear switched
systems based on validated simulation. In: 2016 International Workshop on Symbolic and Numerical Methods for Reachability Analysis (SNR), pp. 1–6, doi:10.1109/SNR.2016.7479377.
[11] Jerry Ding, Eugene Li, Haomiao Huang & Claire J Tomlin (2011): Reachability-based synthesis of feedback
policies for motion planning under bounded disturbances. In: Robotics and Automation (ICRA), 2011 IEEE
International Conference on, IEEE, pp. 2160–2165, doi:10.1109/ICRA.2011.5980268.
[12] Tzanko Donchev & Elza Farkhi (1998): Stability and Euler Approximation of One-sided Lipschitz Differential Inclusions. SIAM J. Control Optim. 36(2), pp. 780–796, doi:10.1137/S0363012995293694.
[13] Laurent Fribourg, Ulrich Kühne & Romain Soulat (2014): Finite controlled invariants for sampled switched
systems. Formal Methods in System Design 45(3), pp. 303–329, doi:10.1007/s10703-014-0211-2.
[14] Antoine Girard (2005): Reachability of uncertain linear systems using zonotopes. In: Hybrid Systems:
Computation and Control, Springer, pp. 291–305, doi:10.1007/978-3-540-31954-2_19.
[15] Antoine Girard, Giordano Pola & Paulo Tabuada (2010): Approximately bisimilar symbolic models for
incrementally stable switched systems. IEEE Transactions on Automatic Control 55(1), pp. 116–126,
doi:10.1007/978-3-540-78929-1_15.
[16] Ian A Hiskens (2001): Stability of limit cycles in hybrid systems. In: System Sciences, 2001. Proceedings of
the 34th Annual Hawaii International Conference on, IEEE, doi:10.1109/HICSS.2001.926280.
[17] Fabian Immler (2015): Verified reachability analysis of continuous systems. In: Tools and Algorithms for the
Construction and Analysis of Systems, Springer, pp. 37–51, doi:10.1007/978-3-662-46681-0_3.
[18] Frank Lempio (1995): Set-Valued Interpolation, Differential Inclusions, and Sensitivity in Optimization.
In: Recent Developments in Well-Posed Variational Problems, Kluwer Academic Publishers, pp. 137–169,
doi:10.1007/978-94-015-8472-2_6.
[19] Jun Liu, Necmiye Ozay, Ufuk Topcu & Richard M Murray (2013): Synthesis of reactive switching protocols from temporal logic specifications. Automatic Control, IEEE Transactions on 58(7), pp. 1771–1785,
doi:10.1109/TAC.2013.2246095.
[20] Kyoko Makino & Martin Berz (2009): Rigorous Integration of Flows and ODEs Using Taylor Models. In:
Proceedings of the 2009 Conference on Symbolic Numeric Computation, SNC ’09, ACM, New York, USA,
pp. 79–84, doi:10.1145/1577190.1577206.
[21] Pierre-Jean Meyer (2015): Invariance and symbolic control of cooperative systems for temperature regulation
in intelligent buildings. Thèse, Université Grenoble Alpes.
[22] Ramon Moore (1966): Interval Analysis. Prentice Hall.

A. Le Coënt, F. De Vuyst, L. Chamoin & L. Fribourg

33

[23] Nedialko S. Nedialkov, K. Jackson & Georges Corliss (1999): Validated solutions of initial value problems for ordinary differential equations. Appl. Math. and Comp. 105(1), pp. 21 – 68, doi:10.1016/S00963003(98)10083-8.
[24] J. Alexandre dit Sandretto & A. Chapoutot (2015): Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes. Research Report, ENSTA
ParisTech.
[25] Julien Alexandre dit Sandretto & Alexandre Chapoutot (2015): DynIbex library.
paristech.fr/ chapoutot/dynibex/.

Http://perso.ensta-

[26] Julien Alexandre dit Sandretto & Alexandre Chapoutot (2016): Validated explicit and implicit runge-kutta
methods. Reliable Computing 22, pp. 79–103.
[27] J. Stolfi & L. H. de Figueiredo (1997): Self-Validated Numerical Methods and Applications. Brazilian Mathematics Colloquium monographs, IMPA/CNPq.

