Vladimir Sinyakov

Antoine Girard

Controller Synthesis for Nonlinear Systems with Reachability Specifications Using Monotonicity

In this paper, we consider a control synthesis problem for a continuous-time nonlinear system. The problem under consideration consists in driving the state of the system to some target interval at a given time instant. We propose a solution based on candidate under-approximations of the backward reachable sets using multi-dimensional intervals. We show that a suitable controller can be designed by enforcing a monotonicity property of the closed-loop system on these intervals. For this purpose, we utilize the monotonicity conditions for nonlinear systems with inputs in the infinitesimal form. From these differential inequalities on the control strategy, we design some particular controllers which are time-varying, linear with respect to the state. The approach is illustrated by two examples.

I. INTRODUCTION

This paper is concerned with the control problem of steering a nonlinear system to a specified target set at a given time instant. The system under consideration is of general nonlinear time-varying type with time-varying uncertainty. The proposed solution strategy is constructed using candidate multi-dimensional interval under-approximations of the backward reachable (i.e. co-reachable) sets. A suitable controller, which keeps the trajectories inside these intervals, can be designed by enforcing a monotonicity property for the closed-loop system.

There exist several different approaches to control synthesis problem in the aforementioned setting. Our approach belongs to the class of methods which do not require state space discretization and do not directly approximate the righthand side [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF], [START_REF] Scott | Bounds on the reachable sets of nonlinear control systems[END_REF], [START_REF] Kostousova | On control synthesis for uncertain differential systems using a polyhedral technique[END_REF], [START_REF] Sinyakov | Method for computing exterior and interior approximations to the reachability sets of bilinear differential systems[END_REF]. Those methods usually consider simple-shaped regions like ellipsoids or parallelotops as basic under-approximations of the backward reachable sets. From these under-approximations, one can usually construct the corresponding controller (see, e.g. [START_REF] Subbotin | Generalized Solutions of First Order PDEs: The Dynamical Optimization Perspective[END_REF]). Such methods can yield families of approximations which may be combined into one covering a larger part of the actual co-reachable sets. Our paper follows this paradigm as well. Multi-dimensional intervals are obviously a particular case of parallelotops. However, we believe that the class of systems, to which our approach is applicable, and the proposed solution based on monotonicity, are different from what is present in the literature.

V. Sinyakov and A. Girard are with the Laboratoire des Signaux et Systèmes (L2S), CNRS, CentraleSupélec, Université Paris-Sud, Université, Paris-Saclay. 3, rue Joliot-Curie, 91192 Gif-sur-Yvette, cedex, France. email: vladimir.sinyakov@l2s.centralesupelec.fr

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 725144). This work is motivated by symbolic control, a control design methodology where discrete abstractions (also called symbolic models) are used to enforce complex behaviors of continuous dynamical systems [START_REF] Tabuada | Verification and Control of Hybrid Systems: a symbolic approach[END_REF], [START_REF] Belta | Formal Methods for Discrete-Time Dynamical Systems[END_REF]. In symbolic control, one of the methods for constructing a symbolic model from a continuous-time control system is to associate symbolic states with elements of a partition of the state space and to define transitions from one state to another if and only if there is a controller steering every trajectory initiating from the region associated with the first state to the one associated with the second state (see e.g. [START_REF] Caines | Hierarchical hybrid control systems: A lattice theoretic formulation[END_REF], [START_REF] Gol | Language-guided controller synthesis for linear systems[END_REF]). The problem of polytope-to-polytope control for nonlinear systems in relation with symbolic control has also been studied using algebraic properties of multi-affine and polynomial vector fields in [START_REF] Belta | Controlling a class of nonlinear systems on rectangles[END_REF], [START_REF] Ben Sassi | Control of polynomial dynamical systems on rectangles[END_REF], [START_REF] Sloth | Control to facet for polynomial systems[END_REF], or a combination of linear approximations with robust control in [START_REF] Girard | Synthesis for constrained nonlinear systems using hybridization and robust controllers on simplices[END_REF].

Our approach is based on a monotonicity property. Intuitively, monotone systems are systems that preserve some partial ordering on their trajectories. The notion of monotone systems with inputs was introduced in [START_REF] Angeli | Monotone control systems[END_REF] where the definition depends on a particular partial order under consideration. In this paper, however, we focus on monotonicity with respect to a specific partial order on vector spaces given by less or equal relation between the components. In Section III, we establish sufficient conditions for monotonicity for a nonlinear time-varying system with inputs. Similar differential conditions are known for time-invariant systems with inputs [START_REF] Angeli | Monotone control systems[END_REF] and time-varying systems without inputs (e.g. [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear monotone systems[END_REF]). However, the authors are not aware of a previous result with both features.

In Section IV, we use the conditions for monotonicity to derive sufficient conditions on the controller to solve our reachability problem (Proposition 2). In Subsection IV.A, we then consider a class of monotone nonlinear control systems for which we verify using the obtained results that a particular control strategy solves the problem and is maximal in some sense. For that class of systems, the approach can be seen as an extension of [START_REF] Meyer | Robust controlled invariance for monotone systems: application to ventilation regulation in buildings[END_REF] where a robust controlled invariance problem was considered, which is essentially a problem of controlling a multi-dimensional interval into itself. Next, in Subsection IV.B we consider a more general class of (not necessarily monotone) nonlinear systems and propose a solution based on time-varying linear feedback control strategies. Finally, in Subsection IV.C we discuss how to construct a hybrid controller applicable in a larger region from several control strategies obtained for different multi-dimensional intervals.

The paper ends with two examples illustrating the approach. In the first example we construct one of the maximal by inclusion multi-dimensional intervals on which the reachability problem has a solution. The second example shows the application of our approach to an autonomous vehicle control problem. This example illustrates the intervalto-interval control synthesis using a number of elementary interval approximations.

Notations: For x ∈ R n , x ∞ = max i |x i | is the infinity norm. Let diag(x) denote a diagonal matrix with elements of x ∈ R n on the diagonal. Given vectors x, x ∈ R n , x x stands for x i ≤ x i for all i = 1, . . . , n. Using this partial order, we define multi-dimensional interval sets as follows:

for x, x ∈ R n , [x, x] = {x | x x, x x}. Similarly, given matrices M , M of the same dimensions, M
M stands for M i,j ≤ M i,j for all i, j. 2 X denotes the set of all subsets of X. For a set-valued map W :

[0, T ] -→ 2 R nw
, the space of all Lebesgue measurable functions w(•) on [0, T ] such that w(t) ∈ W (t) a.e. is denoted by

L ∞ ([0, T ], W (•)). For a set-valued map Z: [0, T ] -→ 2 R nx
, we denote the set of all pairs (t, x) ∈ [0, T ] × R nx such that x ∈ Z(t) by graph(Z). Given a function h: A -→ R, we denote by Arg min

II. PROBLEM FORMULATION

Consider a nonlinear system of the form:

ẋ = f (t, x, u, w), u ∈ R nu , w ∈ W (t) ⊆ R nw . (1)
Here u is the control and w is the disturbance. The problem under consideration is to synthesize a control strategy u : [0, T ] × R nx → R nu , which takes values in U ⊆ R nu and steers the system to the target set X 1 ⊆ R nx at time T > 0, for any disturbance w.

We assume that the vector field f is continuous in (t, x, u, w), continuously differentiable in (x, u, w), Lipschitz in (x, u). In particular, suppose we have

f (t, x, u, w) -f (t, y, u, w) ∞ ≤ L x -y ∞ . (2) A control strategy u : [0, T ] × R nx → R nu is admissible if u is continuous in (t, x), continuously differentiable and Lipschitz in x ∈ R nx . We assume that the set of all possible realizations of the disturbance w is L ∞ ([0, T ], W (•)) where the map W (•) is continuous.
Now let u be an admissible control strategy. Consider the closed-loop system

ẋ = f (t, x, u(t, x), w) = F (t, x, w), w ∈ W (t). (3)
Due to our assumptions, for all initial positions (t 0 , x 0 ) ∈ [0, T ) × R nx and all admissible disturbances w ∈ L ∞ ([0, T ], W (•)), the solution of this system exists on the whole interval [0, T ] and is unique. Let us denote this solution as x(•; t 0 , x, w).

Next, we give the definition of strong invariance.

Definition 1: A set-valued map Z: [0, T ] -→ 2 R nx is strongly invariant with respect to F if for all 0 ≤ t 1 < t 2 ≤
T , for all x ∈ Z(t 1 ) and all w ∈ L ∞ ([0, T ], W ) the solution x(•; t 1 , x, w) exists on [t 1 , t 2 ] and the following inclusion holds: x(t 2 ; t 1 , x, w) ∈ Z(t 2 ).

Formally the main problem we consider in this paper is formulated as follows.

Problem 1: Design an admissible control strategy u :

[0, T ] × R nx → R nu and a set-valued map X -: [0, T ] -→ 2 R nx
, strongly invariant with respect to F , such that

u(t, x) ∈ U ⊆ R nu for all (t, x) ∈ graph(X -),
and for every position (t 0 , x 0 ) ∈ graph(X -), the trajectory of the closed-loop system (3) which starts from this position reaches the target set at time T :

x(T ; t 0 , x 0 , w) ∈ X 1 for all w ∈ L ∞ ([0, T ], W (•)).
We call the set X -(t) the controllable set at time t. We would like to emphasize that if X -is strongly invariant with respect to F and such that X -(T ) = X 1 , then by Definition 1, it follows that x(T ; t 0 , x 0 , w) ∈ X 1 , for all

(t 0 , x 0 ) ∈ graph(X -), for all w ∈ L ∞ ([0, T ], W (•)).
For the following discussion let us introduce some additional related notion.

Definition 2:

The set of all states x for which there exist an initial state x 0 ∈ X 0 and an admissible input w(•) such that x(t; 0, x 0 , w) = x is called the reachable set of system (3) and is denoted X r (t; 0, X 0 ).

III. SUFFICIENT CONDITION FOR MONOTONICITY

Consider a nonlinear system of the form

ẋ = F (t, x, w), w ∈ W (t). (4)
Here W (t) = [w(t), w(t)] where w(•) and w(•) are continuous. Let F be continuous in (t, x, w), continuously differentiable in x and in w and Lipschitz in x. Let Z be a set-valued map, strongly invariant with respect to F .

Definition 3: The system (4) is monotone in graph(Z) if for all 0 ≤ t 1 < t 2 ≤ T , for all x, x ∈ Z(t 1 ) for all w, w ∈ L ∞ ([t 1 , t 2 ], W (•)) x x , and ∀t ∈ [t 1 , t 2 ], w(t) w (t) ⇒ ∀t ∈ [t 1 , t 2 ], x(t; t 1 , x, w) x(t; t 1 , x , w ). (5)
If Z(t) = R nx for t ∈ [0, T ] the system is simply called monotone.

Now consider Z(t) = [x(t), x(t)] where x(•) and x(•) solve the following system of equations:

ẋ(t) = F (t, x(t), w(t)), x(0) = x 0 , (6) ẋ(t) = F (t, x(t), w(t)), x(0) = x 0 . ( 7
)
Under our assumptions, the solution to these equations exists on the whole time interval [0, T ]. The main result of this section is the following.

Proposition 1: Consider system (4), initial value problem (6), [START_REF] Belta | Formal Methods for Discrete-Time Dynamical Systems[END_REF] and let the following conditions hold for all (t, x) in some neighborhood of graph(Z) and all w ∈ W (t):

∂F i ∂x j (t, x, w) ≥ 0, ∀i = j, (8)
∂F i ∂w k (t, x, w) ≥ 0, ∀i, k. (9)
Then Z(t) = ∅, for all t ∈ [0, T ] and Z is strongly invariant with respect to F . Remark 1: In the following proof we, in fact, only use the continuity and monotonicity of F i with respect to the corresponding components of the state space vector and the inputs. The monotonicity condition in the differential form above may be easier to use in the synthesis problem.

Proof: First, let us establish the claim of the proposition in the case when conditions ( 8), ( 9) hold for all x ∈ R nx .

To prove it we utilize the dynamic optimization approach. Let us define the value function:

V (t, x) = min w∈L ∞ ([0,t],W (•)) σ(x(0; t, x, w)) = min w∈L ∞ ([0,t],W (•)) max i max{x i (0; t, x, w) -x 0 i ,
x 0 i -x i (0; t, x, w), 0}. On one hand, the reachable set at time t satisfies the relation [START_REF] Kurzhanski | Dynamic optimization for reachability problems[END_REF]:

X r (t; 0, X 0 ) = {x | V (t, x) ≤ 0}, X 0 = [x 0 , x 0 ]. (10)
On the other hand, V is the viscosity solution of the Hamilton-Jacobi-Bellman equation (for instance, see [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF])

V t + max w∈W (t) V x , F (t, x, w) = 0, V (0, x) = σ(x). (11)
Consider now a function of the form

V -(t, x) = e -Lt max i max{x i -x i (t), x i (t) -x i , 0}. (12) For inequality V -(t, x) ≤ V (t, x) (13)
to hold it is sufficient for V -to be a subsolution of equation ( 11):

q + max w∈W (t) p, F (t, x, w) ≤ 0, (q, p) ∈ D + V -(t, x), (14)
V -(0, x) = σ(x), x ∈ R nx . ( 15
) Since V -is convex in x, its superdifferential D + V -(t, x)
is either empty or singleton for all (t, x). Therefore, it is sufficient for ( 14) to hold at every point of differentiability of V -. Let us consider the case when V -(t, x)e Lt = x i -x i (t):

-LV -(t, x)e Ltẋi (t) + max

w∈W (t) F i (t, x, w) ≤ 0 Since F i is monotone in each w k , we have ẋi (t) ≥ F i (t, x, w(t)) -LV -(t, x)e Lt .
Let ξ j (t, x) = x j , x j ≤ x j (t), x j (t), x j > x j .

, j = 1, . . . , n.

Then we estimate

F i (t, x, w(t)) -LV -(t, x)e Lt ≤ F i (t, ξ(t, x), w(t)) + L x -ξ(t, x) ∞ -LV -(t, x)e Lt .
Here we utilized [START_REF] Scott | Bounds on the reachable sets of nonlinear control systems[END_REF]. One may observe that the expression in the square brackets is equal to zero. Note that in the considered case ξ i (t, x) = x i (t). Using monotonicity of F i with respect to x j for i = j we obtain the relation:

ẋi (t) ≥ F i (t, x(t), w(t)).
Similar reasoning in the case when V -(t, x)e Lt = x i (t)x i gives us the differential inequality ẋi (t) ≤ F i (t, x(t), w(t)).

Thus, for ( 14) and [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear monotone systems[END_REF] to hold for all (t, x)

∈ [0, T ]×R nx it is sufficient that ẋi (t) = F i (t, x(t), w(t)), x(0) = x 0 , ẋi (t) = F i (t, x(t), w(t)), x(0) = x 0
which indeed hold. Therefore, a function V -given by ( 12) is well-defined on [0, T ] × R nx and satisfies inequality [START_REF] Girard | Synthesis for constrained nonlinear systems using hybridization and robust controllers on simplices[END_REF]. Observe from the definition of the value function that

V (t, x(t)) = V (t, x(t)) = 0. Therefore, V -(t, x(t)) = V -(t, x(t)) = 0.
But this is only possible if x(t) x(t). Thus, Z(t) = ∅. Now from [START_REF] Belta | Controlling a class of nonlinear systems on rectangles[END_REF] and [START_REF] Girard | Synthesis for constrained nonlinear systems using hybridization and robust controllers on simplices[END_REF] we have the over-approximation:

X r (t; 0, X 0 ) ⊆ Z(t) ≡ [x(t), x(t)].
From this relation it follows that for every x 0 ∈ X 0 and every w ∈ L ∞ ([0, T ], W ) the solution x(•; 0, x 0 , w) exists on [0, T ]. Moreover, repeating the whole argument for (t 1 , Z(t 1 )) instead of (0, X 0 ) we obtain the relation:

X r (t; t 1 , Z(t 1 )) ⊆ Z(t) ≡ [x(t), x(t)].
Therefore, from the definition of the reachable set we get the strong invariance of Z.

Observe now that the property Z(t) = ∅ depend on values of F in an arbitrary small neighborhood of {(t, x(t), w(t)), t ∈ [0, T ]} and {(t, x(t), w(t)), t ∈ [0, T ]}. In turn, the property of Z being strongly invariant depend on values of F in an arbitrary small neighborhood of graph(Z). This observation completes the proof of the proposition.

Corollary 1: The system (4) is monotone in graph(Z). If in Proposition 1 conditions (8), ( 9) hold for all x ∈ R nx then the system is monotone.

Proof: Consider x, x ∈ Z(t 1 ), x x and some

w, w ∈ L ∞ ([t 1 , t 2 ], W ), w(t) w (t), t ∈ [t 1 , t 2 ].
Since Z is strongly invariant and [x, x ] ⊆ Z(t 1 ), all trajectories starting from [x, x ] at time t 1 lie inside graph(Z). Then the assumptions of Proposition 1 hold for [x(t; t 1 , x, w), x(t, t 1 , x , w )]. Thus, this interval is nonempty for all t ∈ [t 1 , t 2 ] or, in other words, x(t; t 1 , x, w) x(t, t 1 , x , w ).

IV. CONTROL SYNTHESIS

In this section, we derive the conditions for control strategies u : [0, T ]×R nx → R nu to solve the problem of steering system (1) into the target set X 1 = [x 1 , x 1 ]. We assume that W (t) = [w(t), w(t)] as previously.

Consider two control inputs u 1 (•), u 2 (•) ∈ L ∞ ([0, T ]; U ). Let x(•) and x(•) be the trajectories of system (1) with disturbances w(•) and w(•) correspondingly:

ẋ(t) = f (t, x(t), u 1 (t), w(t)), x(T ) = x 1 , (16) ẋ(t) = f (t, x(t), u 2 (t), w(t)), x(T ) = x 1 . ( 17
)
Those equations are to be integrated backward in time on the interval [0, T ]. We assume that x(t) x(t), for all t ∈ [0, T ].

In our construction we use the interval [x(t), x(t)] as the controllable set X -(t). Proposition 2: Consider a control strategy u(t, x) which is continuous in (t, x), continuously differentiable and Lipschitz in x and such that u(t, x(t)) = u 1 (t) and u(t, x(t)) = u 2 (t) where the functions x(•) and x(•) are defined by equations ( 16) and ( 17) respectively with x(t)

x(t), t ∈ [0, T ]. Suppose that u(t, x) ∈ U for all t ∈ [0, T ] and x ∈ [x(t), x(t)]. Let conditions ∂f (t, x, u(t, x), w) ∂x + ∂f (t, x, u(t, x), w) ∂u ∂u(t, x) ∂x γI, ( 18
) ∂f (t, x, u(t, x), w) ∂w 0 ( 19
)
be satisfied for all t ∈ [0, T ], for all x in some neighborhood of [x(t), x(t)], all w ∈ W (t) and for some γ ∈ R. Here I is the identity matrix. Then u and X -given by X -(t) = [x(t), x(t)] solve Problem 1.

Proof: For this strategy to solve control synthesis problem of reaching the target set X 1 , it is sufficient that the closed-loop system satisfies the conditions of Proposition 1 since X -(T ) = X 1 . This is indeed the case since the conditions ( 18) and ( 19) coincide with ( 8) and ( 9) respectively written in a matrix form. It then follows that X -is strongly invariant and u enforces the monotonicity of the closed-loop system in graph(X -).

A. Special case

Let us now consider a nonlinear system with the following special structure (n u ≤ n x ):

ẋi = f i (t, x, u i , w), 1 ≤ i ≤ n u , ẋi = f i (t, x, w), n u < i ≤ n x . ( 20
)
Here u ∈ U = [u, u] and w ∈ W (t). Let us introduce some notations:

x(t) = 1 2 (x(t) + x(t)), x(t) = 1 2 (x(t) -x(t)), û = 1 2 (u + u), ũ = 1 2 (u -u). ( 21
)
Proposition 3: Let the monotonicity conditions

∂f i ∂x j ≥ 0, ∂f i ∂u i ≥ 0, ∂f i ∂w k ≥ 0 hold for all 1 ≤ i, j ≤ n x , 1 ≤ k ≤ n w , i = j, for all t ∈ [0, T ],
for all u ∈ U and all w ∈ W (t).

Consider an interval X -(t) = [x(t), x(t)] defined by the equations

ẋi (t) = f i (t, x(t), u i , w(t)), 1 ≤ i ≤ n u , ẋi (t) = f i (t, x(t), w(t)), n u < i ≤ n x , ẋi (t) = f i (t, x(t), u i , w(t)), 1 ≤ i ≤ n u , ẋi (t) = f i (t, x(t), w(t)), n u < i ≤ n x .
with the terminal conditions

x(T ) = x 1 , x(T ) = x 1 .
This is a maximal by inclusion interval which solves Problem 1 for some strategy u(t, x): every multi-dimensional interval Z ⊃ X -(t) contains uncontrollable states.

Proof: First, observe that the strategy u * defined by

u * i (t, x) = ûi - x i -xi (t) xi (t) ũi (22)
solves Problem 1 for the controllable set

X -(t) = [x(t), x (t)
]. Indeed, conditions of Proposition 2 are clearly satisfied in this case. Now let u(•) be an arbitrary control input from L ∞ ([0, T ], U ). Consider the following system

ẋi (t) = f i (t, x(t), u i (t), w(t)), i = 1, . . . , n u , ẋi (t) = f i (t, x(t), w(t)), n u < i ≤ n x
with the initial condition

x(t 0 ) = x 0 .
Assume that x 0 x(t 0 ). Then due to monotonicity of [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF] and uniqueness of its solution, x(T ) x 1 . Therefore, the system is not controllable from (t 0 , x 0 ).

Similar reasoning for the case x 0 ≺ x 1 completes the proof.

B. Linear feedback case

We build up on the previous case considering now more general linear feedback control strategies for system [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF] with

U = [u, u]. Let Γ(t) = (γ ki (t)) be an n u × n x matrix such that ni i=1 γ ki (t) ≤ 1, γ ki (t) ≥ 0, ∀k.
Then we may define a linear-feedback control strategy u by the formulas:

u(t, x) = û -diag(ũ)Γ(t)[diag(x(t))] -1 (x -x(t)), ẋ(t) = f (t, x(t), û -diag(ũ)Γ(t)e, w(t)), x(T ) = x 1 , ẋ(t) = f (t, x(t), û + diag(ũ)Γ(t)e, w(t)), x(T ) = x 1
where x(t), x(t), û and ũ are the same as in ( 21) and e = (1, . . . , 1) T . If x i (t) < x i (t) for all 1 ≤ i ≤ n x and all t ∈ [0, T ] then this strategy is well-defined in the region graph(X -). Indeed, one may directly check that u(t, x) ∈ [u, u] in this region. Controllers of this form (for parallelotopic regions) were considered in [START_REF] Kostousova | On control synthesis for uncertain differential systems using a polyhedral technique[END_REF].

The monotonicity condition above gives us a relation for

Γ(t): ∂f (t, x, u(t, x), w) ∂x diag(x(t))- ∂f (t, x, u(t, x), w) ∂u diag(ũ)Γ(t) γI, ∀t ∈ [0, T ], ∀x ∈ R nx , ∀w ∈ W (t)
for some γ ∈ R.

C. Discontinuous control strategies

Normally, a different choice of inputs u 1 and u 2 results in different interval sets [x(t), x(t)]. In general, among those sets there are maximal ones: the intervals which are not included into any other such intervals at some t. Computing different such intervals we may then construct a discontinuous controller for a larger controllable set as follows.

Let u s (t, x), s = 1, . . . , N be basic controllers for the controllable intervals X - s (t) = [x s (t), x s (t)] and let V - s (t, x) be the corresponding functions of type [START_REF] Sloth | Control to facet for polynomial systems[END_REF]. Then u(t, x) can be formally defined by

s * (t, x) = min Arg min 1≤s≤N V - s (t, x) , xs (t) = (x s (t) + x s (t))/2, xs (t) = (x s (t) -x s (t))/2 u(t, x) = u s * (t,x) (t, x)
with the corresponding controllable set

X -(t) = 1≤s≤N X - s (t).
Effectively, we chose the number s such that x(t 0 ) ∈ X - s (t 0 ) once and then use the corresponding control strategy u s (t, x).

The controller u(t, x) above solves a version of Problem 1 where the admissible strategy is just Lebesgue measurable in (t, x).

V. EXAMPLES

Example 1. Consider a two-dimensional nonlinear system:

ẋ1 = ux 1 (1 + x 2 ) + w 1 , ẋ2 = -1 1+x 2 1 - 2x 3 2 1+x 2 2 + w 2 .
Here u ∈ [0, 1] is the control, w ∈ [-1, 0] × [-1, 1] is the disturbance. This system is of the type (20) above. Let us compute the strategy and the corresponding controllable set of Proposition 3. Figure 1 depicts the target set in red, X -(0) in green, x(•) and x(•) in black and a simulated trajectory x(•) of the closed-loop system in blue.

Example 2. Examples of monotone systems naturally arise in many applications. One such area is the control of autonomous vehicles. In [START_REF] Saoud | Contract based design of symbolic controllers for interconnected multiperiodic sampled-data systems[END_REF], for instance, an efficient solution to the vehicle platooning problem, which is based on monotonicity of the system, was proposed.

Here we would like to consider the following takeover maneuvre problem. Taking the kinematic car model as a starting point we then perform an input transformation which gives the following nonlinear control system ( [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF], p. 371):

ẋ = u 1 , ẏ = u 1 tan θ, θ = u 1 tan φ, φ = u 2 .
Here u 1 ∈ [0, 1], u 2 ∈ [-1, 1] are the controls. Variables x and y describe the position of the car, θ is the angle between the car body and x axis, φ is the steering angle with respect to the car body. This model is applicable for small absolute values of θ and φ. Therefore, we pose the following statespace constraints |θ| ≤ θ 1 , |φ| ≤ φ 1 < π/2.

In the following we only consider the last three equations of this system. Let z = [y, θ, φ] T . The goal of control here is to steer the system from any state in

X 0 = [0.01, 0.03] × [-0.01, 0.01] × [-0.001, 0.001] to some state in X 1 = [-0.01, 0.01] × [-0.01, 0.01] × [-0.001, 0.001].
Physically, if the road is aligned with the x axis, this may correspond to a transition from one lane to another.

As we consider potential applications of the approach to the symbolic control framework, the regions X 0 and X 1 may, for instance, represent two states of a symbolic model. If using our approach we obtain a controllable set which includes X 0 then the symbolic model would have the corresponding transition between X 0 and X 1 . Therefore, we obtain the relations: Let u 1 (•) and u 2 (•) be arbitrary control inputs with u 1 1 (t) = u 2 1 (t), let z(•) = (y(•), θ(•), φ(•)) T and z(•) = (y(•), θ(•), φ(•)) T form the solution of the corresponding system ( 16), [START_REF] Kurzhanski | Dynamic optimization for reachability problems[END_REF] and define a control strategy as follows: u(t, z) = u 1 (t) .

∂u 1 (t,
Here φ(t) = 1 2 (φ(t)+φ(t)) and φ(t) = 1 2 (φ(t)-φ(t)). Then u(t, z) is an admissible control strategy which satisfies the monotonicity condition and solves Problem 1 in X -(t).

Constructing a discontinuous control strategy as shown above we are able to solve the control problem from any state in X 0 . In this construction we compute the trajectories of ( 16), [START_REF] Kurzhanski | Dynamic optimization for reachability problems[END_REF] checking that the condition z(t) z(t) holds for all t ∈ [0, T ].

Figure 2 depicts the controllable sets X -(0) for different basic control strategies (in black) which correspond to different control inputs u 1 (•) and u 2 (•). The initial state X 0 is given in red.

VI. CONCLUSION

The considered approach gives a new solution to the control synthesis problem for a certain class of nonlinear control systems with uncertain parameters. The condition which must be satisfied by a control strategy ultimately is a system of linear partial differential inequalities. Although, finding a solution to it may be very difficult in general, we presented certain classes of systems and examples where it is possible. The last example illustrates interval-to-interval control synthesis which could be useful from the symbolic control perspective. We are planning to expand in this direction in our future work. We also plan to extend the results to the continuous-time analogue of so called mixedmonotone systems.

  a∈A h(a) the set of all minimizers {a * ∈ A | h(a * ) = min a∈A h(a)}.

2 Fig. 1 .

 21 Fig. 1. Simulated trajectory of the closed-loop system.

Fig. 2 .

 2 Fig. 2. Cover of the initial region X 0 by the controllable regions of different basic control strategies.


K Makki

B Borotikar

M Garetier

S Brochard

D Ben Salem

F Rousseau

HIGH-RESOLUTION TEMPORAL RECONSTRUCTION OF ANKLE JOINT FROM DYNAMIC MRI

Keywords: High-resolution reconstruction, dynamic MRI, intensity-based registration, locally rigid transformation

Cerebral palsy is the leading cause of motor disabilities affecting children. The ankle is the most common equine musculoskeletal strain in children with cerebral palsy. Despite multiple medical and surgical therapies, postoperative recurrence rate is still very high (48%). A major reason for therapy failure is the lack of knowledge of the ankle joint biomechanics. Dynamic MRI can be used to acquire high resolution static data and low resolution temporal images. However, spatial and temporal data should be combined to provide the most comprehensive point of view to study joint motion. In this paper, we first present an intensity-based registration method to accurately estimate the rigid motion of the ankle bones. Second, we investigate the use of the log-euclidean framework to reconstruct a four-dimensional (3D+time) high-resolution dynamic MRI sequence from a low-resolution dynamic sequence and one high resolution static MR image. The proposed approach has been applied and evaluated on in vivo MRI data acquired for a pilot study on child motor disability. Results demonstrate the robustness of the proposed pipeline and very promising high resolution visualization of the ankle joint.

INTRODUCTION

Cerebral palsy is the leading cause of child motor disabilities [1]. This is a term for permanent developmental disorders of movement and posture, responsible of activity limitations caused by nonprogressive damage occurring during early brain development. Approximately 90% of the deformities caused by cerebral palsy occur in the ankle and foot regions. The muscle contractions associated with the defect bone growth impact directly on the posture. These disorders are a frequent causative factor of moral distress for the patient and his/her family. Logistical and financial costs are thus added to the psychological cost.

Management of spastic equinus deformity is still a clinical challenge. There is a general agreement that surgical intervention is the most viable option for fixed equinus deformity with the goal to stabilize the lower extremity, and improve patient walk abilities. However, post-surgical statistics and the long term follow-up studies show up to 48% of recurrence rate and shows a general lack of understanding of spastic equinus deformity in children.

In this context, dynamic MRI is a non-invasive method that can provide spatial and temporal in vivo information of the joint [2]. Acquired MRI data usually consist in a high-resolution static scan containing very accurate anatomical information and a series of lowresolution dynamic scans capturing slow joint motion. However, the combined use of high resolution static data and temporal sequences is limited. Joint analysis of these in vivo data could provide major insights of biomechanical studies of human locomotion.

Joint analysis of high resolution static data and temporal sequences has been used by Gilles et al. in [3] to select of the best dynamic imaging protocol and to track bone motion from real-time dynamic MRI. A 2D/3D registration was used to estimate the hip joint kinematics of healthy volunteers. Clarke et al. have described in [START_REF] Ec Clarke | A non-invasive, 3d, dynamic MRI method for measuring muscle moment arms in vivo : demonstration in the human ankle joint and achilles tendon[END_REF] a non-invasive 3D dynamic MRI method for measuring muscle moment arms in vivo. In this work, which is close to our study, the authors have proposed a 3D mesh-based registration method using an Iterative Closest Point (ICP) algorithm. Bone motions were estimated by the registration of high-resolution joint geometry from MRI scans of the stationary joint with low-resolution geometries from ultrafast MRI scans of the slowly moving joint. The main limitation is the need of manual segmentation of each low-resolution time-frame in addition to the segmentation of the high-resolution static scan.

In this article, we focus on the estimation of ankle joint dynamics of children from MRI. Contrary to the previous work of Clarke et al., we propose to use an automatic intensity-based registration technique to avoid the need of dynamic MRI segmentation. Moreover, we also investigate in this work the use of log-euclidean polyrigid framework developed by Arsigny et al. in [START_REF] Arsigny | Polyrigid and polyaffine transformations : a novel geometrical tool to deal with non-rigid deformations-application to the registration of histological slices[END_REF] for ankle joint analysis. In [START_REF] Papademetris | Articulated rigid registration for serial lowerlimb mouse imaging[END_REF], the authors prove that this type parameterization is well-adapted to the registration of articulated structures such as the hip and knee joint for serial x-ray CT mouse images. In [START_REF] Á Martín-Fernández | A log-euclidean polyaffine registration for articulated structures in medical images[END_REF], an extension of the Log Euclidean polyrigid registration was presented to deal with human articulated structures from 2D radiographs. In our study, we propose to use the log-euclidean framework to estimate high-resolution temporal MRI data. Using a set of bone-dependent rigid transforms, a dense temporal deformation field can be estimated by fusing multiple rigid transforms through the dynamic sequence.

CLINICAL DATA

MRI data have been acquired using a 3T MR scanner (Achieva dStream, Philips Medical Systems, Best, NL). An MRI safe orthotic fixture has been specifically designed in such way that it can be used for performing either passive plantar-dorsiflexion or voluntary active ankle joint movements. These two cases will be called "passive motion" and "active motion" in the next sections. Data were acquired as part of a pilot study for children between 7 and 14 years old at CHRU of Brest, France, under a protocol approved by the local ethical committee. The scanning protocol included one highresolution static 3D scan of the stationary ankle joint with a resolution of 0.26 × 0.26 × 0.8mm, and three low-resolution dynamic scans (two passive sequences for two knee angle positions and one active sequence per child). Each dynamic sequence is composed of 15 time frames with a resolution of 0.57 × 0.57 × 8mm.

METHOD

The proposed approach for dynamic MRI analysis consists in two steps : 1) rigid motion estimation of each bone from dynamic sequence, 2) computation of dense deformation fields by fusing multiple rigid transforms.

Motion estimation

The first step of the proposed method is the motion estimation of each bone. In the context of ankle MRI, we focus on three bones : calcaneus, talus and tibia. In this work, bone motion is modeled with a 3D rigid transform (i.e. rotations and translations). The objective of this step is two-fold : 1) to estimate a mapping for each bone between high resolution static MRI data and dynamic sequences, 2) to estimate the relative motion of each bone using the dynamic sequence.

In this work, it is assumed that the bones of interest have been segmented in the high resolution static image (see Figure 1). In order to reduce manual interventions as much as possible, an intensitybased registration is used to propagate the initial segmentation from the high resolution static image through the dynamic low resolution sequence. The first step is to estimate the transformation from the highresolution static image to the dynamic sequence. This leads to the estimation of a rigid transformation for each bone between images that undergo non-rigid deformations (due to the movement of the joint). In order to make this step robust, the registration is carried out in two steps. First, a rigid registration between the static image and all the images in the dynamic sequence is performed. Then, initialized from the previous step, a rigid transformation is estimated between the static image and each image of the dynamic sequence for each bone. Finally, to maximize the robustness of the algorithm, a image from the dynamic sequence whose the content is the closest to the static image, is selected. To this end, the image of the sequence for which the overlap between regions of interest before and after bone-dependent registration is maximum is selected as the reference image for dynamic motion estimation.

The second step focuses on the estimation of the motion of each bone during the dynamic sequence. To do this, we propose to take advantage of the temporal regularity of the data by propagating the masks of each bone in the dynamic sequence. For each bone, a rigid transform is estimated between successive images. All the registrations are performed using FSL-FLIRT (FMRIB's Linear Image Registration Tool [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]) with correlation ratio as similarity measure. A voxel-wise weighting based on the segmentation masks is used when estimating bone-dependent motion.

The overall algorithm for bone motion estimation is described in Algorithm 1, using the following notations : S is the high-resolution static image, {D k } k=1,...,K is the set of low resolution dynamic images, {M i }i=1,...,N is the set of mask of bones of interest, TA→B is the rigid transform from image A to image B, T i A→B is the rigid transform from image A to image B for the bone i, Dice(A,B) is the DICE coefficient which measures the overlap between segmented regions (Dice(A, B) = 2|A∩B| |A|+|B| ).

Algorithm 1 Bone motion estimation

Fusion of locally rigid transformations

Given a set of rigid transforms for each bone at every time frame of the dynamic MRI sequence, a high resolution dynamic MRI sequence can be estimated by computing a set of dense deformation fields corresponding to each time frame. To this end, each dense deformation field is obtained by fusing rigid transforms corresponding to the bones of interest. For each bone, a weighting function is computed based on a distance to the mask [START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF]. Each weighting function reflects the local influence of the rigid transform of each bone of interest. In this study, following previous work of Commowick et al. [START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF], for each voxel x, the weighting functions are computed as follows :

w i D k (x) = 1 1 + α dist(x, M i D k ) (1)
where dist(x, M i D k ) is the distance between the point x and the mask of the bone i propagated onto the dynamic image D k , α is set to 0.5 according to [START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF] to ensure smooth interpolation.

Given a point x in the high-resolution static image, the target location of this point in the dynamic image D k can be computed within the log-euclidean polyrigid framework proposed by Arsigny et al. in [START_REF] Arsigny | Polyrigid and polyaffine transformations : a novel geometrical tool to deal with non-rigid deformations-application to the registration of histological slices[END_REF], with the following equation :

TS→D k (x) = exp N i=1 wi D k (x) log(T i S→D k ) .x (2)
where TS→D k is the dense deformation field from the static image S to the dynamic image D k , wi D k is a normalized weighting function (i.e.

N i=1 wi D k (x) = 1, ∀x ∈ D k ).

RESULTS

Each high resolution static image of the dataset has been carefully manually segmented to extract accurate region of interest of the three bones considered in this work (calcaneus, talus and tibia, see Figure 1). In order to take into account the partial volume effect due to anisotropy of the image resolution, a Gaussian blur with a standard deviation set to 2 voxels has been applied onto each bone mask.

To evaluate the proposed registration approach for bone motion estimation, the propagated mask of each bone through the dynamic sequence is compared with a manually delineated bone on the first and the last dynamic image (i.e. D1 and DK ). The accuracy of the estimations of the successive rigid transforms is assessed by computing the following DICE coefficient Dice(B i k , B i m k ) where B i k is the propagated mask of the bone i on dynamic image D k and B i m k is a manually delineated mask of the same bone on the dynamic mask D k . A DICE value close to one indicates that the bone masks have been well propagated through the entire dynamic sequence.

Results on bone motion estimation are reported in Table 1 for six subjects of the pilot study. Results for passive and active motion are separated. In 33 of 36 cases, the DICE coefficient is greater than 0.8 indicating accurate propagations of bone masks over the dynamic sequence, for both passive and active motions. Although based on a successive estimation of transformations, the proposed approach appears to be robust enough for almost all subjects. It can be noticed that the proposed approach can be used for both passive and active motion of the ankle joint.

From the set of rigid transforms of each bone over the dynamic sequence, it is then possible to compute a temporal dense deformation field to synthesize a high resolution temporal image of the ankle joint. Figure 3 illustrates, for two subjects, such high resolution reconstructed data. The use of a diffeomorphic approach for polyrigid transform fusion leads to the estimation of realistic high resolution temporal data to study in vivo joint motion through dynamic MRI.

CONCLUSION

Our work focuses on dynamic MRI analysis for child ankle joint study. An intensity-based registration pipeline have been proposed to estimate bone motion through the dynamic sequence. The robustness of this approach for both passive and active motion has been evaluated on in vivo data acquired for a pilot study on child motor disability due to cerebral palsy. We have also investigated the use of log-euclidean framework to estimate dense temporal deformation fields from multiple rigid transforms. The polyrigid approach provides a way to combine high resolution spatial information with temporal dynamics of joints. Visual promising results have been presented showing the potential of this framework for detailed in vivo biomechanical analysis of joint motion from dynamic MRI.

In future works, the proposed approach will be applied to a larger dataset to compare ankle joint dynamics of children with cerebral palsy with age-matched healthy children. We will further explore this scheme to extract fine biomechanical parameters of tendon and ligaments. Dynamic MRI and dedicated image processing techniques can open a new way to study in vivo human joints.

[1] A Colver, C Fairhurst, and P-OD Pharoah, "Cerebral palsy,"

The Lancet, vol. 383, no. 9924, pp. 1240-1249, 2014.

[2] Won C Bae, T Ruangchaijatuporn, and CB Chung, "New techniques in mr imaging of the ankle and foot," Magnetic Resonance Imaging Clinics, vol. 25, no. 1, pp. 211-225, 2017.

[3] B Gilles, R Perrin, N Magnenat-Thalmann, and J-P Vallee, "Bone motion analysis from dynamic MRI : acquisition and tracking," Academic Radiology, vol. 12, no.

= 1 k = K k = 1 k = K k = 1 k = K k = 1 k = K k = 1 k = K k = 1 k = K Calcaneus 0.
= 1 k = K k = 1 k = K k = 1 k = K k = 1 k = K k = 1 k = K k = 1 k = K Calcaneus 0.

Fig. 1 :

 1 Fig. 1: Bones of interest : calcaneus (yellow), talus (green) and tibia (red). (a) : Mid-sagittal image from the high-resolution static scan ; (b) : Three dimensional rendering of segmented bones.

1:

  Input : segmentation of bones of interest in S. 2: Mapping estimation from S to {D k } : -Estimate rigid transforms {TS→D k } k=1,...,K -Estimate rigid transform {T i S→D k } k=1,...,K using bone mask M i , initialized with {TS→D k } k=1,...,K -Select D k in {D k } k=1,...,N such that k = argmax k ( N i=1 Dice(TS→D k (M i ), T i S→D k (M i ))) 3: Motion estimation : -Estimate forward successive rigid transforms {T i D k →D k+1 } k = k ,...,K-1; i=1,...,N . -Estimate backward successive rigid transforms {T i D k →D k-1 } k = k ,...,1; i=1,...,N .

Fig. 2 :

 2 Fig. 2: Normalized weighting functions wi D k (x). (a) : corresponding high resolution reconstructed image, weighting function for the calcaneus (b), talus (c), tibia (d) with the following color coding : blue if greater than 0.5, red if lesser than 0.5, white if equal to 0.5. The use of polyrigid framework leads to smooth diffeomorphic dense deformation fields.

  k = 1 corresponds to the first time frame of the dynamic sequence while k = K corresponds to the last time frame.

Fig. 3 :

 3 Fig. 3: High resolution temporal reconstruction for one subject (subject3). First row : original dynamic images (time frames : 1, 3, 6, 9 and 13) , second row : corresponding high resolution reconstructed image. Each column corresponds to one time frame. The first time frame contours shows the reconstruction accuracy (column 1) and the joint motion (columns 2 to 5)

Table 1 :

 1 3D DICE coefficients between manual and automatic bone segmentation maps for the six subjects of the pilot study.

		88	0.9	0.95	0.92	0.97	0.91	0.88	0.9	0.87	0.9	0.92	0.94
	Talus	0.8	0.88	0.9	0.9	0.94	0.92	0.93	0.95	0.52	0.45	0.4	0.99
	Tibia	0.95	0.98	0.96	0.9	0.97	0.97	0.99	0.95	0.95	0.96	0.98	0.97

* . The research leading to these results has received funding from