Sebti Mouelhi

Antoine Girard

Gregor Gössler

CoSyMA: A Tool for Controller Synthesis Using Multi-scale Abstractions

Keywords: Switched systems, Multi-scale abstractions, Controller synthesis, Symbolic Algorithms Systèmes commutés, Abstraction multi-échelles, Synthèse de contrôleurs, Algorithmes symboliques

In this report we introduce CoSyMA, a tool for automatic controller synthesis for incrementally stable switched systems based on multi-scale discrete abstractions. The tool accepts a description of a switched system represented by a set of differential equations and the sampling parameters used to define an approximation of the state-space on which discrete abstractions are computed. The tool generates a controller -if it exists -for the system that enforces a given safety or time-bounded reachability specification. We illustrate by examples the synthesized controllers and the significant performance gains during their computation.

Introduction

Controller synthesis for hybrid systems using discrete abstractions has become an established approach (see [START_REF] Tabuada | Verification and Control of Hybrid Systems: A Symbolic Approach[END_REF] and the references therein). In [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], it has been demonstrated that the (sampled) continuous behavior of incrementally stable [START_REF] Angeli | A lyapunov approach to incremental stability properties[END_REF][START_REF] Liberzon | Switching in systems and control[END_REF] switched systems are approximately bisimilar to some discrete abstractions. The states of these abstractions are elements of lattices that approximate the continuous state-space. Time and space sampling parameters are chosen to achieve a desired precision; the smaller the time sampling parameter, the finer the lattice used for approximating the statespace, and consequently, the larger the number of states in the abstraction. The approach detailed in [START_REF] Cámara | Synthesis of switching controllers using approximately bisimilar multiscale abstractions[END_REF][START_REF] Cámara | Safety Controller Synthesis for Switched Systems Using Multi-Scale Symbolic Models[END_REF] based on multi-scale discrete abstractions can be used to cope with this problem. These abstractions are defined over a set of embedded lattices. The finer lattices are only explored when the specification cannot be met at the coarsest level.

In this paper, we present CoSyMA (COntroller SYnthesis using Multi-scale Abstractions), a tool implementing symbolic approaches based on multi-scale abstractions to synthesize controllers for incrementally stable switched systems. CoSyMA accepts as input a switched system defined by differential equations indexed by a set of modes, time and space sampling parameters used to set an approximation of the continuous state-space, and a safety or a time-bounded reachability specification. If it exists, it computes a controller satisfying the specification. The tool is implemented using OCaml [START_REF] Leroy | The Objective Caml system release 3.12[END_REF] and it is available for download at multiscale-dcs.gforge.inria.fr.

2 Theoretical background 2.1 Incrementally stable switched systems Definition 1. A switched system is a quadruple Σ = R n , P, P, F where R n is the state-space; P = {1, . . . , m} is a finite set of modes; P is the set of piecewise constant functions from R + to P, continuous from the right and with a finite number of discontinuities on every bounded interval of R + ; F = {f 1 , . . . , f m } is a collection of smooth vector fields indexed by P. For all mode p ∈ P, f p : R n → R n is a locally Lipschitz continuous map.

A switching signal of Σ is a function p ∈ P. A piecewise C 1 function x : R + → R n is said to be a trajectory of Σ if it is continuous and there exists a switching signal p ∈ P such that, at each t ∈ R + where the function p is continuous, x is continuously differentiable and satisfies ẋ(t) = f p(t) (x(t)). We denote by x(t, x, p) the point reached at time t, starting from the state x and applying the constant switching signal p(s) = p, for all s ∈ [0, t]. A switched system is δ-GUAS (i.e. globally uniformly asymptotically incrementally stable [START_REF] Angeli | A lyapunov approach to incremental stability properties[END_REF][START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]) if all the trajectories associated with the same switching signal converge asymptotically to the same reference trajectory independently of their initial states.

Approximate bisimulation

In this section we provide a brief introduction of the notion of approximate bisimulation that relates a switched system to the specific discrete abstraction we construct. Let us consider a class of transition systems of the form T = Q, L, r, O, H, I consisting of a set of states Q; a set of labels L; a transition relation r ⊆ Q × L × Q; an output set O; an output function H : Q → O; a set of initial states I ⊆ Q. T is said to be metric if the output set O is equipped with a metric d, discrete if Q and L are finite or countable sets. For q ∈ Q and l ∈ L let succs l (q) = {q ′ ∈ Q | (q, l, q ′ ) ∈ r}. An action l ∈ L belongs to the set of enabled actions at the state q, denoted Enab(q), if succs l (q) = / 0. The transition system is said to be deterministic if for all q ∈ Q and l ∈ Enab(q), succs l (q) has only one element denoted by succ l (q). A trajectory of the transition system is a finite or infinite sequence of transitions σ = q 0 l 0 q 1 l 1 q 2 l 2 . . . , it is initialized if q 0 ∈ I. A state q ∈ Q is reachable if there exists an initialized trajectory reaching q. We denote by Φ(T) the set of all the trajectories of T.

Transition systems can described the dynamics of switched systems. Given a switched system Σ = R n , P, P, F , let T(Σ) = Q, L, r, O, H, I be the transition system where the set of states is Q = R n ; the set of labels is L = P × R + ; the transition relation is given by (x, (p, τ), x ′ ) ∈ r iff x(τ, x, p) = x ′ , i.e. the switched system Σ goes from state x to state x ′ by applying the constant mode p for a duration τ; the set of outputs is O = R n ; the observation map H is the identity map over R n ; the set of initial states is I = R n . T(Σ) is deterministic and metric when the set of outputs O = R n is equipped with the metric d(x, x ′ ) = xx ′ . The relation between the discrete abstractions and T(Σ) can be defined by an approximate bisimulation [START_REF] Girard | Approximation metrics for discrete and continuous systems[END_REF].

q 1 ∈ Q 1 q ′ 1 ∈ Q 1 T 1 : q 2 ∈ Q 2 q ′ 2 ∈ Q 2 T 2 : (q 1 , q 2 ) ∈ R (q ′ 1 , q ′ 2 ) ∈ R d(H 1 (q 1 ), H 2 (q 2 )) ≤ ε d(H 1 (q ′ 1 ), H 2 (q ′ 2 )) ≤ ε l l Figure 1: The principal of ε-approximate bisimulation relation Definition 2. Let T i = Q i , L i , r i , O i , H i , I i , for i ∈ {1,
≥ 0. A relation R ⊆ Q 1 × Q 2 is
said to be an ε-approximate bisimulation relation between T 1 and T 2 if for all (q 1 , q 2 ) ∈ R:

• d(H 1 (q 1 ), H 2 (q 2 )) ≤ ε; • ∀(q 1 , l, q ′ 1 ) ∈ r 1 , ∃(q 2 , l, q ′ 2 ) ∈ r 2 , such that (q ′ 1 , q ′ 2 ) ∈ R; • ∀(q 2 , l, q ′ 2 ) ∈ r 2 , ∃(q 1 , l, q ′ 1 ) ∈ r 1 , such that (q ′ 1 , q ′
2 ) ∈ R. T 1 and T 2 are said to be approximately bisimilar with precision ε, denoted T 1 ∼ ε T 2 , if for all q 1 ∈ I 1 , there exists q 2 ∈ I 2 , such that (q 1 , q 2 ) ∈ R, and for all q 2 ∈ I 2 , there exists q 1 ∈ I 1 , such that (q 1 , q 2 ) ∈ R.

The diagram shown in Figure 1 illustrates the principle of the ε-approximate bisimulation relation.

Multi-scale abstractions

In applications where the switching has to be fast, uniform approximately bisimilar abstractions, as defined in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF], approximate the state-space using fine lattices that results Inria in a a huge number of abstract states. In practice, fast switching is generally necessary only on a restricted part of the state space. For instance, for safety specifications, fast switching is needed only when the system gets closer to unsafe regions. In order to enable fast switching while using abstractions with a reasonable number of states, we consider discrete abstractions enabling transitions of different durations. For transitions of long duration, it is sufficient to consider abstract states on the coarse lattice. The finer ones are reached by the shorter transitions only when the specification cannot be met at the coarsest level.

Let us consider a switched system Σ whose switching is determined by a timetriggered controller with time-periods in the finite set Θ N τ = {2 -s τ | s = 0, . . . , N} that consists of dyadic fractions of a time sampling parameter τ > 0 up to some scale parameter N ∈ N. The dynamics of a switched system Σ is then described by the transition system

T N τ (Σ) = Q 1 , P × Θ N τ , r 1 , O, H 1 , I 1 where Q 1 = O = I 1 = R n , H 1 is the identity map over R n , and (x, (p, 2 -s τ), x ′ ) ∈ r 1 iff x(2 -s τ, x, p) = x ′ .
The discrete abstraction of T N τ (Σ) is defined on an approximation of

Q 1 = R n by a set of embedded lattices [R n ] s η defined by [R n ] s η = q ∈ R n q[i] = k i 2 -s+1 η √ n , k i ∈ Z, i = 1, ..., n
where s = 0, . . . , N, q[i] is the i-th coordinate of q and η > 0 is a state space discretization parameter. By simple geometrical considerations, we can check that for all x ∈ R n and s = 0, . . . , N, there exists q ∈ [R n ] s η such that xq ≤ 2 -s η. Then, we can define the abstraction of T N τ (Σ) as the transition system

T N τ,η (Σ) = Q 2 , P × Θ N τ , r 2 , O, H 2 , I 2 , where the set of states is Q 2 = [R n ] N
η ; the set of actions remains L = P × Θ N τ ; r 2 is defined such that (q, (p, 2 -s τ), q ′ ) ∈ r 2 iff q ′ = arg min m∈[R n ] s η ( x(2 -s τ, q, p)m ). The approximation principle is illustrated in Figure 2.3. The observation map H 2 is the natural inclusion map from [R n ] N η to R n ; the set of initial states is

I 2 = [R n ] 0 η . q ′ q q ′′
x(τ/2, q, p)

x(τ, q, p)

q ′ = succ 2 (q, (p, τ 2 )) = arg min m∈[R n ] 1 η ( x( τ 2 , q, p) -m ) and q ′′ = succ 2 (q, (p, τ)) = arg min m∈[R n ] 0 η ( x(τ, q, p) - m )
The resulting abstraction T N τ,η (Σ) is discrete and deterministic, its set of states and its set of actions are respectively countable and finite. For N = 0, we recover the "uniform" abstractions introduced in [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. In [START_REF] Cámara | Synthesis of switching controllers using approximately bisimilar multiscale abstractions[END_REF], it was proved that for a switched system Σ admitting a common δ-GUAS Lyapunov function, T N τ (Σ) ∼ ε T N τ,η (Σ) where the precision ε can be made arbitrarily small by reducing the state sampling parameter η.

Controller synthesis using multi-scale abstractions

Before presenting the tool details and our experimental results, we explain briefly how we use multi-scale abstractions for synthesizing safety and time-bounded reachability controllers. Let T = Q, L, r, O, H, I be a deterministic transition system, a controller for T is a map S : Q → 2 L such that for all q ∈ Q, S(q) ⊆ Enab(q). The system T controlled by S is T/S = Q, L, r S , O, H, I where the transition relation is given by (q, l, q ′ ) ∈ r S iff (l ∈ S(q) ∧ (q, l, q ′ ) ∈ r). The support of S is defined by supp(S) = {q ∈ Q | S(q) = / 0}.

Safety controller synthesis

Given a safety specification Q S ⊆ Q (obtained from a subset O S ⊆ O of safe outputs), a state q of T is controllable with respect to a safety specification Q S if q ∈ Q S and there exists an infinite trajectory σ ∈ Φ(T) starting from q and remaining in Q S . We denote the set of controllable states of T with respect to the safety specification Q S by SCont(T, Q S ). A safety controller S for T and Q S is defined such that supp(S) ⊆ SCont(T, Q S ) and for all q ∈ supp(S): (1) q ∈ Q S (safety) and (2) ∀l ∈ S(q), succ l (q) ∈ supp(S) (deadend freedom). The set SCont(T, Q S ) is computable for discrete abstractions. However, the larger the number of states, the more expensive the computation.

For that reason, we want to capitalize on multi-scale abstractions to propose an efficient algorithm for safety controller synthesis. The lazy safety synthesis problem consists in controlling a system so as to keep any trajectory starting from some initial state in I within the safe subset of states Q S , while applying at each state transitions of the longest possible duration for which safety can be guaranteed. For that purpose we define a priority relation on the set of labels L = P × Θ N τ giving priority to transitions of longer duration: for l, l ′ ∈ L with l = (p, τ),

l ′ = (p ′ , τ ′ ), l l ′ iff τ ≤ τ ′ , l ≺ l ′ iff τ < τ ′ and l ∼ = l ′ iff τ = τ ′ . Given a subset of labels L ′ ⊆ L, we define max (L ′ ) = {l ′ ∈ L ′ | ∀l ∈ L ′ , l l ′ }.
Definition 3. A maximal lazy safety (MLS) controller S : Q → 2 L for T and Q S is a safety controller such that I ∩ SCont(T, Q S ) ⊆ supp(S) and for all states q ∈ supp(S), we have: (1) if l ∈ S(q), then for any l ≺ l ′ , succ l ′ (q) / ∈ SCont(T, Q S ) (laziness), and (2) if l ∈ S(q), then for any l ∼ = l ′ , l ′ ∈ S(q) iff succ l ′ (q) ∈ SCont(T, Q S ) (maximality).

The MLS controller exists and is unique as proved in [START_REF] Cámara | Safety Controller Synthesis for Switched Systems Using Multi-Scale Symbolic Models[END_REF].

Time-bounded reachability controller synthesis

Given a transition system T = Q, L, r, O, H, I , for all transitions (q, l, q ′ ) ∈ r, let δ(l) be the time needed by T to reach q ′ from q by action l. For all finite trajectories σ = q 0 l 0 q 1 l 1 . . . l n-1 q n in Φ(T), we define its duration by ∆(σ) = δ(l 0 ) + δ(l 1 ) + • • • + δ(l n-1 ). For instance, for the transition systems T N τ (Σ) and T N τ,η (Σ), we have L = P × Θ N τ and for all l = (p, 2 -s τ) ∈ L, we have δ(l) = 2 -s τ.

To formally define time-bounded reachability controller, we define

C(T) = Q c , L, r c , O, H c , I c the transition system with clock of T where Q c = Q × R + is the set of states Q ex- tended by a clock; for all ((q, c), l, (q ′ , c ′ )) ∈ r c , (q, l, q ′ ) ∈ r and c ′ = c + δ(l); H c ((q, c)) = H(q) for all (q, c) ∈ Q c ; I c = I × {0}. The set of reachable states of C T N τ,η (Σ) is countable and defined by Q × 2 -N τN. Given a maximal time bound B ∈ R + , the state (q, c) of C(T) is controllable with respect to a time-bounded reachability specification Inria (Q S , Q T , B), where Q T ⊆ Q S if (1) q ∈ Q S and
there exists a finite trajectory σ of T starting from q, eventually reaching Q T , and remaining in Q S until reaching Q T , such that c + ∆(σ) ≤ B. The set of these states is denoted by RCont(C(T), Q S , Q T , B).

Next we define time-bounded reachability controllers using the notion of safety controllers. We start by defining the notion of stuttering (•) actions. An outgoing transition from a state q labeled by a stuttering action loops on the same state (succ • (q) = q and δ(•) = 0).

Let us now define

T Q ′ = Q, L ∪ {•}, r , O, H, I for Q ′ ⊆ Q such that (q, l, q ′ ) ∈ r iff q = q ′ if l = • and q ∈ Q ′ ; (q, l, q ′ ) ∈ r if l = • and q ∈ Q \ Q ′ .
T Q ′ is the transition system derived from T where the only actions enabled from a state in Q ′ are stuttering.

Since we are not concerned with the evolution of the system after reaching the target Q T , we will use the transition system

C(T Q T ) = Q c , L, r c , O, H c , I c rather than C(T). We easily show that RCont(C(T Q T ), Q S , Q T , B) = RCont(C(T), Q S , Q T , B) = SCont(C(T Q T ), Q S × [0, B]).
A time-bounded reachability controller for the transition system

C(T Q T ) and (Q S , Q T , B) is a safety controller for C(T Q T ) and Q S × [0, B].
: Q × R + → 2 L for C(T Q T ) and (Q S , Q T , B) is the MLS controller for C(T Q T ) and Q S × [0, B].
It is clear, based on the previous definition, that R m is unique and that Q T is reachable within the specified time bound starting from controllable initial states.

Example 2.1. We shall give a simple example to describe the principle of MLBR controller. Given a time and space sampling parameters τ and η, we consider the transition system T 1 τ,η (Σ) of a switched system Σ with a unique mode p. The set of labels is defined by L = P × Θ 1 τ where P = {p}. In Figure 3 (left), we show only a part of the transition system T 1 τ,η (Σ) where, from initial states {1, 5, 7, 9}, we take only paths that enter the set of target states Q T (the set {3, 4, 6, 8} is included in Q T ) within a time lesser or equals to B = 2τ. We suppose that all states {1, . . . , 9} are in the set of safe states Q S . The system controlled by the MLBR controller according to the specification (Q S , Q T , 2τ) is described in Figure 3 (right). The reader can easily check that R m is a time-bounded reachability controller for T 1 τ,η (Σ) and (Q S , Q T , 2τ) such that the maximality and laziness properties are satisfied.

We can use the algorithm synthesizing MLS controllers proposed in [START_REF] Cámara | Safety Controller Synthesis for Switched Systems Using Multi-Scale Symbolic Models[END_REF] to synthesize MLBR controllers. However, their computation is expensive because dealing with problems of n dimensions amounts to handle their equivalents of n + 1 dimensions by adding a clock. To avoid this constraint, we can settle for a sub-controller V of R m such that all controllable initial states of R m are also controllable by V.

Let R m be the MLBR controller for

C(T Q T ) and (Q S , Q T , B). A sub-controller V of R m is a time-bounded reachability controller for C(T Q T ) and (Q S , Q T , B) such that for all (q, c) ∈ supp(V), V((q, c)) ⊆ R m ((q, c)). The sub-controller V is complete if {i ∈ I | (i, 0) ∈ supp(R m )} = {i ∈ I | ∃0 ≤ c ≤ B | (i, c) ∈ supp(V)}.
We can now define static reachability controllers based on the previous definition. (p, τ)

(p, τ) (p, τ) (p, τ) (p, τ) (p, τ/2) (p, τ/2) (p, τ/2) (1, 0) (2, τ) (3, 2τ) (5, 0) (1, τ) (4, 1.5τ) (7, 0) (5, τ) (6, 1.5τ) (9, 0) (7, τ) (8, 1.5τ) (p, τ) (p, τ) (p, τ) (p, τ) (p, τ) (p, τ/2) (p, τ/2) (p, τ/2) • • • • Figure 3: A part T of T 1 τ,η (Σ) (above); The transition system C(T Q T )/R m where R m is the MLBR controller for T and (Q S , Q T , B) (below) Definition 5. Consider a complete sub-controller V of R m . The static reachability controller R ls V : Q → 2 L for T Q T and (Q S , Q T , B) obtained from V is the controller such that for all (q, c) ∈ supp(V), q ∈ supp(R ls V ) and for all q ∈ supp(R ls V ), R ls V (q) = V (q, c max (q)) where c max (q) = max{c ′ | (q, c ′ ) ∈ supp(V)}. (5, 0) (1, τ) (4, 1.5τ) (9, 0) (7, τ) (8, 1.5τ) (p, τ) (p, τ) (p
• (p, τ/2) (p, τ) • (p, τ/2) 9 7 8 5 1 6 4 (p, τ) • • • (p, τ/2) (p, τ/2) (p, τ/2) Figure 5: The transition system T Q T /R ls V controlled by R ls V where V is the controller of C(T Q T )/V shown in Figure4 (above); the transition system T Q T /R ls V where V = R m (below)
Example 2.3. Consider the complete sub-controller V where C(T Q T )/V is shown in Figure 4. The static reachability controller R ls V derived from V is defined by the transition system shown in Figure 5 (left). If V = R m , the correspondent static reachability controller R ls V is shown in Figure 5 (right). Theorem 1. Any trajectory σ of T Q T /R ls V starting from q ∈ supp(R ls V ) eventually reaches Q T and ∆(σ) ≤ B.

Inria

Proof. For all states q ∈ supp(R ls V ), the state (q, c), where c = max{d | (q, d) ∈ supp(V)}, is controllable with respect to (Q S , Q T , B). This means that any path σ in C(T Q T )/V starting from (q, c) eventually reach Q T such that c + ∆(σ) ≤ B. Lets take (q ′ , c ′ ) the first state reachable by σ such that c ′ < max{d | (q ′ , d) ∈ supp(V)}, if we replace the rest of σ starting from (q ′ , c ′ ) by a path

σ ′ of C(T Q T )/V starting from (q ′ , c m q ′ ) such that c m q ′ = max{d | (q ′ , d) ∈ supp(V)}, σ eventually reach Q T and c ′ + ∆(σ ′ ) < c m q ′ + ∆(σ ′ ) ≤ B.
If we apply the same operations in such way, for each state (q i , c i ) reachable by σ,

c i = max{d | (q i , d) ∈ supp(V)}, σ eventually reach Q T and c + ∆(σ) ≤ B. These paths without clocks are those of T Q T /R ls V .

Static reachability controller synthesis

In this section, we present an algorithm computing the static reachability controller for a given complete sub-controller

V for R m where R m is a MLBR controller for C(T Q T ) and (Q S , Q T , B).
Algorithm 1 guards the system at its original dimension and explore states of T Q T rather then states of C(T Q T ) to reduce the synthesis complexity caused by growing the system by the supplementary clock parameter. For each explored trajectory σ, it stores for each reachable state q by σ three clocks: (i) the current encountered clock tmp_clock(q) for q, (ii) the effective maximal encountered clock(q) such that the state (q, clock(q)) ∈ RCont(C(T), Q S , Q T , B) and reachable in C(T Q T )/V, and (iii) the minimal encountered clock cutoff (q) for q such that (q, cutoff (q)) is unreachable in C(T Q T )/V. A state q is controllable and added to the abstraction reachable states if clock(q) is set at least one time. Given a map f : X → Y, we denote by dom(f ) = X the domain of f . The algorithm uses a first-depth traversal starting from initial states q 0 ∈ I with a temporary clock tmp_clock(q 0 ) = 0. For each initialized trajectory σ = q 0 l 0 q 1 l 1 . . . l i-1 q i . . . of T Q T , we use a stack path containing tuples (-, q 0 ), (l 0 , q 1 ), . . . , (l i-1 , q i ) from the bottom to the top. Since q 0 has no predecessor in σ, the tuple corresponding to q 0 , added to path, is (-, q 0 ). The stack path is grown until reaching a state q ∈ Q T or a previously controlled state q such that tmp_clock(q) ≤ clock(q) (lines 31 to 35) if there is no tuple (l, q ′ ) such that q ′ = q was inserted before in path to avoid the treatment of cyclic paths looping on q. In other words, a state q is pushed on path and explored if clock(q) was not set before during the exploration of another initialized trajectory σ ′ or tmp_clock(q) > clock(q) because, in the latter case, we can find traces reaching the target from (q, tmp_clock(q)) in a shorter time than those starting from (q, clock(q)). For a top element (l, q), transitions of longer duration starting from q are explored first (laziness). We use for that two maps lexp and rexp associating, respectively, to each state the set of explored labels and remaining labels to explore. (-, q 0 ) := pop(initials) ; 5 if q 0 ∈ dom(clock) then //if q 0 was set controllable then it is not explored

Algorithm 1: STATIC REACHABILITY CONTROLLER SYNTHESIS Input: Time-bounded specification (Q S , Q T , B); T Q T = Q, L, r • , O, H, I where I ⊆ [R n ] 0 η ∩ Q S ; priority order . Output: An LSB reachability controller R : Q → 2 L . Data: initials: stack({-} × I); path: stack(L ∪ {-} × Q); clock,tmp_clock,cutoff : Q → N; rexp: Q → 2 L ; lexp: Q → 2 L ; → R
6 if q 0 ∈ Q T then clock(q 0 ) := 0;
7 else push(path,(-, q 0 )) ; rexp(q 0 ) := max (L); lexp(q 0 ) := / 0; tmp_clock(q 0 ) := 0;

8 while path = / 0 do 9 (l, q) := top(path) ; 10 if rexp(q) = / 0 ∧ ∃(q, l ′ , q ′ ) ∈→ R then //if q is fully explored and has successors 11 pop(path) ;

12 if q ∈ dom(clock) then clock(q) := tmp_clock(q); 13 if q = q 0 then //while q is not initial if rexp(q) = / 0 then rexp(q) := max (L \ lexp(q)); 20 choose l ′ from rexp(q); q ′ := succ l ′ (q); 21 clk := tmp_clock(q) + δ(l ′ );

22 if q ′ ∈ dom(cutoff ) ∨ cutoff (q ′ ) > clk then 23 if clk ≤ clock(q ′ ) then 24
call COPYTRANSITION((q, l ′ , q ′ ));

25 if q ′ ∈ dom(clock) ∨ clk > clock(q ′ ) then 26 if q ′ ∈ Q S ∧ clk ≤ B then //if q ′ is safe 27 if q ′ ∈ Q T then //if q ′ within the target 28 clock(q ′ ) := clk;
29 call COPYTRANSITION((q, l ′ , q ′ )); 30 else if q ′ was not pushed before in path then 31 tmp_clock(q ′ ) := clk; 32 //initialize or reinitialize labels to explore for q ′ 33 rexp(q ′ ) := max (L); lexp(q ′ ) := / 0; 34 push(path,(l ′ , q ′ ));

35

//update the set of explored labels for q 36 lexp(q) := lexp(q) ∪ {l ′ }; rexp(q) := rexp(q) \ {l ′ }; 37 else//all labels in L are explored for q without finding successors 38 cutoff (q) := tmp_clock(q); 39 pop(path);

40 return {(q, l) | ∃(q, l, q ′ ) ∈→ R } 41 end 42
A tuple (l, q) is popped from the head of path if (1) (tmp_clock(q), q) is uncontrollable (lines 38 to 40), or (2) it is controllable and labels l = (m, 2 -s τ), for all m ∈ P by incrementing s from 0 to N, are all explored until finding controllable successors for q.

In the case (2), q is set to be controllable if not the case (lines 11 to 13) and the function COPYTRANSITION with the input (p, l, q) is called. If tmp_clock(p) > clock(p), clock(p) is updated to tmp_clock(p) and all the transitions starting from p previously Inria added to → R are removed and replaced by (p, l, q). The update of tmp_clock(p) have to be made (line 4 of the function COPYTRANSITION only one time to replace transitions starting from p with the old value of clock(p) by the first computed transition starting from p with the new value of tmp_clock(p). The rest of outgoing transitions from p with the new clock are added after the update and lines 2 to 5 are not executed when COPYTRANSITION is called. The function COPYTRANSITION is also called with input (q, l ′ , q ′ ) where q ′ = succ l ′ (q) if tmp_clock(q ′ ) ≤ clock(q ′ ) (line 24) or q ′ ∈ Q T (line 28).

Procedure COPYTRANSITION Input: a transition (p, l, q). begin

1 if tmp_clock(p) > clock(p) then 2 Invariant (∃(p, l ′ , q ′ ) ∈ → R );
→ R := → R ∪{(p, l, q)} 6 end 7 ≫ 1 path → R clock 1 (3, (p, τ)) (2, (p, τ)) (1, -) 1 2 3 (p, τ) (p, τ) • 1 → 0 2 → τ 3 → 2τ 5 (4, (p, τ/2)) (1, (p, τ)) (5, -) 5 1 4
(p, τ)

• • • (p, τ/2) (p, τ/2) (p, τ/2) 1 → τ 5 → 1 7 → τ 9 → 0 4 → 1.5τ 6 → 1.5τ 8 → 1.5τ
Table 1: The LSB controller synthesis according to ≫ 1 Algorithm 1 terminates: for each reachable state q with a temporary clock c, a tuple (l, q ′ ) with q ′ = q is inserted at most once in path because if q was set controllable (clock(q) = c), it cannot be reinserted in path except in the case where c is grater than clock(q i ) (line 24). In the case where (q, c) is uncontrollable, cutoff (q) is set to c which prevent adding the state with the same temporary clock a second time in path (line 23). Also, eventually, q will be popped from path (lines 12 and 40).

From an operational view point, the complete sub-controller V ⊆ R m is selected according to the order of exploration ≫ I of initial states (q 0 , 0) ∈ I × {0}. More clearly, let's take the transition system shown in Fig. 3 (above), according to the order ≫ 1 = 1 ⊲ 5 ⊲ 7 ⊲ 9, the complete sub-controller V ≫ 1 is the MLBR controller R m shown in Figure 3 (below) (cf. Table. 1). However, if we consider the order ≫ 2 = 9 ⊲ 7 ⊲ 5 ⊲ 1, the complete sub-controller V ≫ 2 is shown in Figure 4 (cf. Table 2). Algorithm 1 iterates the stack initials containing tuples (q 0 , -) where q 0 ∈ I according to the order ≫ I from top to bottom.

≫ 2 path → R clock 9 ( 8
(p, τ)

• (p, τ/2) (p, τ) • (p, τ/2) 9 → 0 7 → τ 8 → 1.5τ 1 → τ 5 → 0 4 → 1.5τ 1
Nothing is done because tmp_clock(1) ≤ clock(1), i.e. 0 ≤ τ Table 2: The LSB controller synthesis according to ≫ 2

Tool details

In this section, we present the internal architecture of CoSyMA, its description language, and our implementation choices. CoSyMA accepts a configuration (.conf) file describing the system and the synthesis parameters. It contains in the order: the description of a switched system Σ in terms of differential equations ẋ(t) = f p(t) (x(t)); time τ and space η sampling parameters, and the scale N used to compute the finer lattice [R] N η that approximate the continuous state-space; a safety specification Q S or a timebounded reachability specification (Q S , Q T , B); the plot parameters. The grammar is detailed in Appendix A).

The tool architecture, shown in Figure 6, represents the different steps by which the tool synthesizes the controllers of the described system according to the given safety or a time-bounded reachability specification. After the parsing of the configuration file, the tool represents the vector fields f p for each switching mode p by an OCaml function of type float -> float array -> float array. It computes the successor of a state q under a label l = (p, 2 -s τ) by solving ẋ(t) = f p (x(t)) for t ∈ [0, δ(l)] and x(0) = q using the common fourth-order Runge-Kutta method. The tool synthesizes the controllers based on T N τ,η approximately bisimilar to T N τ (cf. Section 2.3). For safety specifications Q S , the tool synthesizes the MLS controller based on the algorithm presented in [START_REF] Cámara | Safety Controller Synthesis for Switched Systems Using Multi-Scale Symbolic Models[END_REF] which computes the abstraction T N τ,η on the fly.

Inria

where id is the name of the configuration. The dimension of the system is defined by the statement "Dimension" = n where n is a strictly positive natural number. The modes are defined by the statement Modes = { 'p1 , ... , 'pk };

where 'p1, ..., 'pk for k >= 1 are the different modes. The keyword "Constants" is followed by constants and their values used to define the differential equation systems.

The variable declaration are defined after the keyword "Variables". The variables may be dependent or independent of time. The variables are used to describe the functions used to represent the trajectory fields of a system x(t). After the keyword "Trajectory", we define the piecewise function x(t) representing the trajectory in terms of time-dependent variables. Many trajectories can be defined. For example, we define a trajectory x(t) by the following statement.

x = [v1,v2] ;

We mention that the special variable "@t" is used for the time parameter. Under the keyword "Coefficients", given a linear differential equation systems ẋ(t) = a(t)x(t) + b(t), we define the coefficients a and b. Two types of coefficients are allowed : timedependent matrices and time-dependent vectors.

a : Matrix (2); b : Vector (2);
The type Matrix (n) represents square matrices of dimension n and Vector (n) represents vectors of dimension n. For example, we can define the matrix a by the following statement a = [[c1 * @t, sin(2 * @t) + 1],[5, 2 * @t]; and the vector b by the following one. b = [5 * @t, c2 + #t]; where c1 and c2 are constants. The keyword "VectorFields" is used to define the differential equation system ẋ(t) = f p (x(t)) of the declared modes p. The differential equation systems may be simplified or expanded. The statement used to define a linear differential system for a given mode p is the following Inria simplified when 'p :: x_dot = a * x + b; where x is the system trajectory. The statement used to define an expanded differential equation system for a given mode q is the following nonlinear when 'q :: v1_dot = v1 + cos(v2) -@t, v2_dot = ln(v1);

where v1 and v2 are the elements of x. Under the keyword "SamplingParameters", we declare the time and space sampling parameters, and a desired finer refinement scale. time = 0.5 ; space = 1 / (40 * sqrt ( 2)); finer = 1;

CoSyMA allows also the declaration of state-spaces (see the grammar in Appendix A), and the specification of safety and time-bounded reachability problems. The command for the controller synthesis for safety problems is written as follows SafetySynthesis (x,safe,init,ds,rk4s,hash_size) ; where x is the trajectory identifier, safe is the identifier of safe state-space, init is the identifier of initial state-space, ds equals to "Enum" or "Bdd" according to the desired abstraction representation, rk4s indicates the precision (number of iterations) of the fourth-order Runge-Kutta method, and finally hash_size is the initial size of the hash table used to represent the abstraction (it is considered only when ds = "Enum"). The command time-bounded reachability problems is written as follows ReachabilitySynthesis (x,safe,init,target,b,ds,rk4s,hash_size) ; where target is the target state-space and b is a natural such that the maximal time bound B = bτ.

Finally, the tool CoSyMA also allows to generate TikZ [START_REF] Tantau | The TikZ and PGF Packages, Manual for version 2[END_REF] plots to visualize controllers and simulate them. It allows the generation of plots for systems of two, three, and four dimensions. The plot configuration for two dimensional systems is defined as follows.

Plot2D (dora,mors,modes_colors,scales_colors,x_res,y_res) ; where dora equals to "dots" or "arrows" according to the desired plot with dots (printing only states) or arrows (printing states and enabled transitions), mors equals to "modes" if we want to print states according to the modes of their outgoing transitions or "scales" if we want to print states according to the durations of their outgoing transitions. The lists modes_colors and scales_colors represents respectively the colors used for modes and durations. The reals x_res and y_res represents respectively the resolution of the x-axis and the y-axis. The plot configurations Plot2Dof3D and Plot2Dof4D are similar to Plot2D except we have to add after the parameter y_res parameter the variables of the trajectory corresponding to the x-axis and y-axis by indicating the values of the rest of variables.

Abstraction representation and manipulation

The user has the choice to represent the system abstractions either by enumerated types or boolean functions. The tool uses hash tables as an enumerated type to operate on the system abstractions. Hash tables turn out to be more efficient than search trees or other table lookup structures, especially for large numbers of entries.

Alternatively, the tool can use BDDs (Binary Decision Diagrams) [START_REF] Bryant | Graph-based algorithms for boolean function manipulation[END_REF] to represent the system abstractions. They are able to represent sets and relations compactly in memory as boolean functions. We implement BDDs using the modules Cudd.Man and Cudd.Bdd of the OCaml IDL interface MLCUDDIDL [START_REF] Jeannet | MLCUDDIDL: OCaml interface for CUDD library, version 2.2.0[END_REF] of the CUDD (CU Decision Diagram) package [START_REF] Somenzi | CUDD: CU Decision Diagram Package[END_REF]. However, using BDDs makes the controller synthesis algorithms presented above more costly than hash tables since the symbolic abstraction is constructed by enumerating the states and computing their successors using the Runge-Kutta method.

Experimental results

DC-DC Converter

As a first case study, we apply our approach to a boost DC-DC converter. It is a switched system with two modes, the two dimensional dynamics associated with both modes are affine of the form ẋ(t) = a p x(t) + b for p ∈ {1, 2} (see [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF] for numerical values). It can be shown that it is incrementally stable and thus approximately bisimilar discrete abstractions can be computed. We consider the problem of keeping the state of the system in a desired region of operation given by the safe set O S = [1.15, 1.55] × [5.45, 5.85].

We use approximately bisimilar abstractions to synthesize MLS controllers for the DC-DC converter. We compare the cost of controller synthesis for the uniform abstraction T 0 τ 1 ,η 1 for parameters τ 1 = 0.5s and η 1 = 10 -3 √ 2/4 (containing transitions of duration 0.5s) and the multi-scale abstractions T2 τ 2 ,η 2 for parameters τ 2 = 4τ 1 and η 2 = 4η 1 (containing transitions of durations in Θ 2 τ = {2s, 1s, 0.5s}). These two abstractions have the same precision. Table 3 details the experimental results obtained for the synthesis of the T 0 τ 1 ,η 1 and T 2 τ 2 ,η 2 . We can see that there is a noteworthy reduction of the time used to compute the controller using multi-scale abstractions instead of using uniform ones (up to a 86% improvement between T 0 τ 1 ,η 1 and T 2 τ 2 ,η 2 ). This is due to the fact that the size of uniform abstractions grows exponentially with higher resolutions, whereas using multi-scale abstractions are refined only when we get closer to unsafe regions (reduction of more than 91% between T 0 τ 1 ,η 1 and T 2 τ 2 ,η 2 ). Interestingly, this reduction in computation time and size does not affect the performance of the multi-scale controllers, which yield a ratio of controllable initial states1 (CR) over the safety specification comparable to that of their its uniform counterparts. It is worth emphasizing that using CoSyMA there is a remarkable reduction of the computation times compared to those reported in in [START_REF] Cámara | Safety Controller Synthesis for Switched Systems Using Multi-Scale Symbolic Models[END_REF] obtained by a prototype implementation of the algorithm. Figure 7 depicts the maximal lazy safety controller for T 2 τ 2 ,η 2 and Q S and the trace of its simulation starting from the state (1.15,5. size static reachability controllers respectively for T 0 τ 1 ,η 1 where τ 1 = 0.25 and η 1 = 10 -3 √ 2/4, and T 2 τ 2 ,η 2 where τ 2 = 4τ 1 and η 2 = 4η 1 , and the specification (Q S , Q T , B). As shown in Section 2.4, the MLBR controller for the abstraction C(T Q T ) and (Q S , Q T , B) where T equal to T 0 τ 1 ,η 1 or T 2 τ 2 ,η 2 is the maximal lazy safety controller for C(T Q T ) and the safety specification Q S × [0, B]. Its computation is costly because the problem is grown form 2 to 3 dimensions by considering the supplementary clock parameter. Synthesizing a static reachability controller rather than an MLBR controller significantly reduces complexity. In Table 4, we observe a considerable reduction of the size of the controlled abstraction of T 2 τ 2 ,η 2 of more than 91.46% compared to that of T 0 τ 1 ,η 1 with comparable controllability ratios of initial states. Also, the computation time of the controller of T 2 τ 2 ,η 2 is significantly shorter than that of T 0 τ 1 ,η 1 .

Abstractions T N τ,η N = 0, τ = 0.5s, N = 2, τ = 2s, η = 10 .9] and B = 20s. We synthesize static reachability controllers respectively for T 0 τ 1 ,η 1 where τ 1 = 0.25 and η 1 = 10 -3 √ 2/4, and T 2 τ 2 ,η 2 where τ 2 = 4τ 1 and η 2 = 4η 1 , and the specification (Q S , Q T , B). As shown in Section 2.4, the MLBR controller for the abstraction C(T Q T ) and (Q S , Q T , B) where T equal to T 0 τ 1 ,η 1 or T 2 τ 2 ,η 2 is the maximal lazy safety controller for C(T Q T ) and the safety specification Q S × [0, B]. Its computation is costly because the problem is grown form 2 to 3 dimensions by considering the supplementary clock parameter. Synthesizing static reachability controllers rather than MLBR controllers, reduces widely the synthesis complexity. In Table 4, we observe a considerable reduction of the size of the controlled abstraction of T 2 τ 2 ,η 2 of more than 91.46% compared to that of T 0 τ 1 ,η 1 with comparable controllability ratios of initial states. Also, the com-

Building Temperature Regulation

The second case study deals with temperature regulation in a circular building. Each room is equipped with a heater and at a given instant at most one heater is switched on. The temperature t i of the room i is defined by the differential equation ṫi = α(t i+1 + t i-1 -2t i )+β(t e -t i )+γ(t h -t i )u i (t) where t i-1 is the temperature of the room i-1; t i+1 the temperature of the room i + 1; t e is the temperature of the external environment of the building; t h is the temperature of the heater; α is the temperature transfer ratio between the rooms i ± 1 and the room i; β is the temperature transfer ratio between the external environment and the room i; γ is the temperature transfer ratio between the heater and the room i; u i (t) equals to 1 if the room i is heated, or 0 otherwise. Given a number n ≥ 2 of rooms, we distinguish n + 1 switching modes. For 1 ≤ i ≤ n, the mode p i represents the mode of activating the heater of room i. The mode p n+1 represents that no heater is activated. The values of α, β, γ, t e , and t h are respectively 1/20, 1/200, 1/100, 10, and 50. We will increase the system dimension to test the limits of the tool in terms of memory usage and computation time. Given the safety specification Q S = [20.0, 22.0] n for n ∈ {3, 4, 5}, we synthesize safety controllers for buildings of three, four, and five rooms. The values of τ and η are given in Table 5. By looking Inria to the results, we can see the combinatorial explosion of the size of abstractions by increasing the system dimension from 3 to 5. Also, it makes sense that the ratio of controllability of initial states decreases by increasing the number of rooms. On our machine equipped with 4GiB of RAM, synthesis fails for the 6-dimensional instance due to to memory overflow. Table 5: Comparison of experimental results for the safety synthesis for the temperature regulator system of three, four, and five dimensions Figure 9 shows the MLS controller for the transition system T 2 20,0.05 of three dimension and the safety specification Q S = [20.0, 22.0] 3 . The plots are slices of the state space in the dimensions (t 1 , t 2 ) for a fixed t 3 ≈ 20 • (left) and t 3 (t) ≈ 22 • (right), respectively. The plots on the top depicts scales and those in the middle and the bottom depicts modes. We can remark the predominance of the mode p 4 (no heater is activated) by increasing the temperature of third room.

Conclusion

In this paper we have introduced CoSyMA, a tool that automatically synthesizes controllers for incrementally stable switched systems based on multi-scale discrete abstractions. We have illustrated by examples the synthesized controllers for safety and Figure 9: MLS controller for T 2 20,0.05 of three dimension and Q S ; Horizontal axis: t 1 (t) ; Vertical axis: t 2 (t) ; Top: t 3 (t) ≈ 20 ; Bottom: t 3 (t) ≈ 22; Left: actions of 20s are enabled (black); actions of 10s are enabled (gray); Middle: mode p 1 (black); mode p 2 (light gray); p 1 and p 2 (gray); Right: mode p 4 (black); mode p 3 (light gray); p 3 and p 4 (gray). time-bounded reachability problems. The benchmarks provide evidence that the use of multi-scale abstractions leads to a substantial reduction of synthesis time and size of the obtained controller while maintaining coverage of the state space.

A The tool Grammar

The grammar of CoSyMA configuration file is defined as follows.

  2}, be metric transition systems where L 1 = L 2 and O 1 = O 2 equipped with the metric d, and a precision ε

Figure 2 :

 2 Figure 2: Computation of the discrete abstraction:q ′ = succ 2 (q, (p, τ 2 )) = arg min m∈[R n ] 1 η ( x( τ2, q, p)m ) and q ′′ = succ 2 (q, (p, τ)) = arg min m∈[R n ] 0 η ( x(τ, q, p)m )

Definition 4 .

 4 We define the maximal lazy timebounded reachability controller based on Definition 3 as follows. The maximal lazy time-bounded reachability (MLBR) controller R m

, τ/ 2 )Figure 4 :

 24 Figure 4: The transition system C(T Q T )/V where V is a complete sub-controller of the MLBR controller R m shown in Figure 3 (right)Example 2.2. Let us consider the transition system C(T Q T )/R m of Fig.3(right), we can define the controller V such that C(T Q T )/V is the one shown in Fig.4. It is clear that V is a sub-controller for R m . The sub-controller V is complete because all the initial states {i ∈ I | ∃(i, 0) ∈ supp(R m )} (states {1, 5, 7, 9}) remain controllable in C(T Q T )/V. The states 1 and 7 are controllable with clocks greater than 0.

2 Initialization:

 2 is the set of transitions. begin 1 Invariants: dom(clock ∪ cutoff ) ⊆ dom(tmp_clock); dom(cutoff ) ∩ dom(clock) = / 0; initials contains tuples in {-} × I;

  14(k, p) := top(path) ; 15 call COPYTRANSITION((p, l, q)); the remaining labels to explore for q 19

3 4 →

 4 clock(p) := tmp_clock(p) ; R := → R \{(x, l, y) | x = p} ; //remove the outgoing transitions from p 5

  6). Now, we consider the time-bounded reachability specification (Q S , Q T , B) where Q S = [0.65, 1.65] × [4.95, 5.95], Q T = [1.1, 1.6] × [5.4, 5.9] and B = 20s. We synthe-Abstractions T N τ,η N = 0, τ = 0.5s, N = 2, τ = 2s, η = 10 -3 √ 2/4 η = 10

Figure 7 :

 7 Figure 7: The MLS controller for T 2 τ 2 ,η 2 and Q S . Top: mode 1 is activated (light gray); mode 2 is activated (black); modes 1 and 2 (gray); Bottom (Left): actions of 2s are enabled (light gray); Bottom (right): actions of 1s are enabled (light gray), actions of 0.5s are enabled (black).

Figure 8 :

 8 Figure 8: A static reachability controller for T 2 τ 2 ,η 2 and (Q S , Q T , B): mode p 1 (light gray); mode p 2 (black); both modes (medium gray)

n = 3 ,

 3 N = 2, n = 4, N = 2, n = 5, N = 1 η = 50 × 10 -3 η = 50 × 10 -3 η = 0.1 τ = 20s τ = 20s τ = 10s

  configuration ::= Switched id { Dimension = natural; Modes = enum_fileds; Constants : constants Variables : variables Trajectories : trajectories maybe_coefs VectorFields : vector_fields SamplingParameters : sampling_parameters spaces synthesis_specification plot_simulation} where id ::= string id_list ::= (id,) * id string_list ::= (string,) * string int_list ::= (int,) * int real_list ::= (real,) * real enum_fileds ::= (enum_filed,) * enum_filed enum_filed ::= 'string constants ::= (constants_by_value) * constants_by_value constants_by_value ::= id_list = xp; variables ::= (variables_by_type) * variables_by_type variables_by_type ::= id_list : variable_type; trajectories ::= (trajectory) * trajectory trajectory ::= id = [id_list]; maybe_coefs ::= Coefficients : coefs_types; coefs_values | <empty> coefs_types ::= (coefs_by_type) * coefs_by_type coefs_by_type ::= id_list : coef_type; coefs_values ::= (coefs_by_value) * coefs_by_value coefs_by_value ::= id_list : coef_value coef_value ::= matrix_of_xp | vector_of_xp | xp matrix_of_xp ::= xp array array vector_of_xp ::= xp array vector_fields ::= (trajectory_vfs) * trajectory_vfs trajectory_vfs ::= for id : (trajectory_vf ) * trajectory_vf trajectory_vf ::= simplified_vf | expanded_vf simplified_vf ::= simplified when enum_filed :: id concat "_dot" = xp; expanded_vf ::= expanded when enum_filed :: (var_diffeq) * var_diffeq var_diffeq ::= id concat "_dot" = xp; sampling_parameters ::= time = xp; space = xp; finer = natural; spaces ::= (space) * space space ::= id = (interval#) * interval; interval ::= [real .. real] synthesis_specification ::= safety_synthesis | reachability_synthesis safety_synthesis ::= SafetySynthesis(id,id,id,ds,natural,natural); reachability_synthesis ::= ReachabilitySynthesis(id,id,id,id,natural,ds,natural,natural); ds ::= Enum | Bdd plot_simulation ::= plot_2d | plot_2d3d | plot_2d4d | simulation plot_2d ::= Plot2D(dora,mors,[string_list],[string_list],real,real,[int_list]) plot_2d3d ::= Plot2Dof3D(dora,mors,[string_list],[string_list],real,real,[int_list],id, id,real) plot_2d4d ::= Plot2Dof4D(dora,mors,[string_list],[string_list],real,real,[int_list],id, id,real,real) Simulate ::= Simulate([real_list],real,real,int,real,real)InriaThe non-terminals xp (expressions), variable_type (types of variables), and coef_type (type of coefficients), are defined as follows.

  xp ::= @t | id | natural | int | real | (xp) | + xp |xp | xp + xp | xpxp | xp * xp | xp / xp | xp % xp | xp ^xp | unary_fun (xp) | log (xp,xp) unary_fun ::= sin | cos | tg | ctg | sec | csc | asin | acos | atg | sinh | cosh | tanh | sqrt | ln | exp variable_type ::= real | interval | real -> real | real -> interval coef_type ::= real -> real | real -> interval | Matrix (natural) | Vector (natural) RR n°8108 RESEARCH CENTRE GRENOBLE -RHÔNE-ALPES Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Table 3 :

 3 Experimental results for the MLS controller synthesis for the boost DC-DC converter

Table 4 :

 4 Experimental results for the static reachability controller for the boost DC-DC converter Now, we consider the time-bounded reachability specification (Q S , Q T , B) where Q S = [0.65, 1.65] × [4.95, 5.95], Q T = [1.1, 1.6] × [5.4, 5

RR n°8108

The ratio of controllable initial states for a controller S : Q →

L and a system T = (Q, L, →, O, H, I) is computed as |{q ∈ I|S(q) = / 0}|/|I|. Inria

The algorithm is a depth-first traversal of paths σ ∈ Φ(T N τ,η (Σ)) starting from initial states I until reaching Q T to keep trace of the clocks of the states reached by σ. From an operational point of view, the complete sub-controller V ⊆ R m , according to which the static reachability controller R ls V is defined, is computed according to the order of exploration of initial states I. By keeping the problem at its original dimension, the synthesis complexity is significantly reduced.

Configuration file

A configuration must be encapsulated in a preamble declared as follows Switched id { ... }


Adam Scheinherr

Lucie Bailly

Olivier Boiron

Thierry Legou

Antoine Giovanni

Georges Caillibotte

Marine Pichelin

Glottal motion and its impact on the respiratory flow

Keywords: Glottis, Breathing, Laryngofiberscopy, Image processing, CFD simulations

come

Introduction

Advantages of inhaled therapies as a priori targeted supply of drugs make them particularly convenient for the treatment of lung diseases. Nevertheless, several physical and anatomical factors can largely influence treatment efficiency. In particular, upper airways (UA) anatomic arrangement acts as an unwanted filter, which limits the amount of drug delivered to the lung [START_REF] Conway | Controlled, Parametric, Individualized, 2-D and 3-D Imaging Measurements of Aerosol Deposition in the Respiratory Tract of Healthy Human Subjects for Model Validation[END_REF]. More specifically, the glottis, defined by vocal folds aperture within the larynx, causes airways to narrow in a minimal transition cross-section. This anatomical singularity yields to a complex jet-like tracheal flow [START_REF] Brouns | Influence of glottal aperture on the tracheal flow[END_REF][START_REF] Katz | Threedimensional computational study of inspiratory aerosol flow through the larynx: the effect of glottal aperture modulation[END_REF][START_REF] Renotte | Numerical analysis of oscillatory flow in the time-varying laryngeal channel[END_REF], which can be determinant on particles deposition by inertial impaction. However, current studies are limited by two main issues: -Glottal dynamics during human breathing have been barely investigated so far, and despite a few reference in vivo studies [START_REF] Baier | Relationships among glottis opening, respiratory flow and upper airway resistance in humans[END_REF][START_REF] Brancatisano | Respiratory movements of the vocal cords[END_REF], the relationship between glottal area and inhaled airflow is still poorly understood.

-The aerodynamic influence of the glottal geometrical changes during breathing is often discarded in numerical works. Instead, a static glottis is often considered together with steady flow conditions [START_REF] Brouns | Influence of glottal aperture on the tracheal flow[END_REF][START_REF] Katz | Threedimensional computational study of inspiratory aerosol flow through the larynx: the effect of glottal aperture modulation[END_REF]. Therefore, the aim of the present study is (i) to characterize the glottal dynamics during human breathing in vivo using laryngofiberscopy and synchronized airflow recordings; (ii) to quantify the effects of a mobile glottis and unsteady flow conditions on laryngeal jet-flow dynamics using CFD modeling.

Methods

In vivo study. In vivo experiments were conducted in the ENT department of La Timone Hospital. One healthy female volunteer (LB, age 29) and one healthy male volunteer (OB, age 48) were recorded while performing two 30s breathing tasks: normal breathing (eupnea) and forced breathing (tachypnea). Laryngofiberscopic investigations were made using a flexible nasofiberscope (Storz 202220 20 tricam camera) with a continuous cold light source and a color CCD camera. Laryngeal images were captured with a camera frame rate of 25frames/s and an image resolution of 768x288 pixels. The oral airflow signal was simultaneously registered by means of a pneumotachograph placed at the mouth, EVA2 [START_REF] Ghio | Evaluation of the acoustic and aerodynamic constraints of a pneumotachograph for speech and voice studies[END_REF]. In silico study. As a first approximation of the glottal geometry was built a 2D rectangular moving constriction with a triangular dynamic mesh of 24000 el. The mesh density ensures grid-independent results. CFD simulations were conducted using Fluent 6.3.26 under laminar air-flow conditions, assuming an incompressible Newtonian gas of viscosity 1 equal to 1.789210 -5 kg/ms. A pressure outlet boundary condition was set to 0 Pa. Unsteady boundary conditions comprising the glottal width d g (t) and the velocity inlet were parametrically varied in agreement with the in vivo study. A no-slip shear boundary condition was applied at solid walls. Initially, zero velocities and pressures were assumed at all points. Equations were solved using a first-order time and spatial discretization schemes, and a time step set to 0.12s.

Results and discussion

Airflow rate. In vivo measurements yielded to about 30 respiratory cycles during eupnea and 40 during tachypnea. Every respiratory cycle was detected using a zero-tracking method developed in Matlab R2011b. Each airflow signal Q was normalized with respect to the maximum value achieved within the cycle, Q max . Time t was normalized by corresponding respiratory period, T = 2π/ω. The maximal period registered within each 30s sequence is noted T max . Signals were finally averaged into one mean flow-rate, herein noted <Q/Q max >. Figure 1a illustrates the typical flow-rate <Q/Q max > as a function of 3t, produced by subject OB during eupnea (maximal Q max = 1.9 dm 3 /s, T max = 0.62s) and tachypnea (maximal Q max = 4.3 dm 3 /s, T max = 4.8s). In comparison, a sinusoidal evolution is plotted. Conventionally, positive (resp. negative) flowrate values correspond to expiration (resp. inspiration) phase. During eupnea, the inspiration and expiration curves are similar (up to sign) and their durations are approximately equal, as commonly found in the literature [START_REF] Fenn | Handbook of Physiology, Section 3: Respiration[END_REF]. During tachypnea, the mean flow-rate curve deviates from the harmonic signal although inspiration and expiration durations remain roughly equal. A phase difference of about 22° in flow-rate maximal occurrences has been measured between sinusoid and eupnea curve, which corresponds to 6% of the breathing period. The phase difference increases during tachypnea up to 58°, namely 16% of the cycle duration. LB breathing shows similar characteristics. *Corresponding author. Email: Scheinherr@irphe.univ-mrs.fr b) Glottal motion. Glottal motion was extracted from the laryngoscopic images using Matlab Image Processing Toolbox and different phases: (i) Focus on a region of interest using a cross-correlation technique, (ii) Smoothing using a specific filter function [START_REF] Kroon | Hessian based Frangi Vesselness filter -Matlab code[END_REF], (iii) Detection of the glottal contours using a segmentation method [START_REF] Bernard | Variational B-Spline level-set: A linear filtering approach for fast deformable model evolution[END_REF] (iv) Measurement of the glottal anteroposterior diameter AP g , glottal area A g and glottal width d g (see Fig. 1b), (v) Normalization by comparing AP g diameter with the initial value AP g 0 , assumed as a geometrical invariant [START_REF] Higenbottam | Narrowing of glottis opening in humans associated with experimentally induced bronchoconstriction[END_REF], (vi) Conversion from pixels to millimeters, assuming AP g = 22mm [START_REF] Conway | Controlled, Parametric, Individualized, 2-D and 3-D Imaging Measurements of Aerosol Deposition in the Respiratory Tract of Healthy Human Subjects for Model Validation[END_REF]. Fig. 1c illustrates the evolution of the glottal dynamics measured during a typical eupnea cycle. Glottal area A g and glottal width d g are displayed as function of time, together with synchronized airflow-rate Q. It is shown that the glottis progressively widens during inspiration and narrows during expiration. Area A g varies in the range 90 -240mm 2 (mean value 165mm 2 ) during the cycle, while d g varies in the range 4 -13mm (mean value 7.8mm), thus achieving a peak value during inspiration nearly three times greater than that measured during expiration. This ratio is in line with previous studies [START_REF] Baier | Relationships among glottis opening, respiratory flow and upper airway resistance in humans[END_REF][START_REF] Brancatisano | Respiratory movements of the vocal cords[END_REF]. Above measurements allowed to assess a Reynolds number (Re={d g Q}/{1A g }) ranging from 500 to 2500. CFD simulations. Unsteady flow simulations were conducted considering Q and dg measured during the typical eupnea cycle (Fig. 1c). Fig. 2 shows the development of the glottal jet at 4 shot-instants at first during inspiration and consequently during expiration phase of breathing (crosses in Fig. 1c). Inertial effects associated to flow-rate variations yield to the jet instability and fluctuations of the reattachment area during the breathing cycle but the flow despite Re≅2500 stays laminar as already noticed in [START_REF] Boiron | Experimental and numerical studies on the starting effect on the secondary flow in a bend[END_REF].

Conclusions

The in vivo study showed that the glottis can be extremely variable during breathing and hence influence airflow characteristics. A glottal area widening was quantified during inspiration, with a typical ratio of 3:1 as compared to expiration. Airflow rate variations differ from harmonic signal during eupnea as well as tachypnea. The correlation between flow-rate and glottal area will be discussed and compared to previous clinical investigations. Preliminary 2D CFD simulations of the glottal jet were performed, based on the measured flow-rate and glottal changes during eupnea. Impact of unsteady flow conditions on the jet development is demonstrated.

Figure 1 :

 1 Figure 1: (a) Mean flowrate <Q/Qmax> measured during cycles of eupnea and tachypnea (subject OB) and comparison with sinusoid. (b) Illustration of glottal image post-processing.(c) Detected glottal area Ag and glottic width dg and flowrate Q.

Figure 2 :

 2 Figure 2: 2D unsteady simulations of velocity field through the moving glottis during a breathing cycle.