Guillaume Deffuant

Isabelle Alvarez

AND Sophie Martin

Patrick Saint-Pierre

APPROXIMATING BACKSTOP VIABILITY KERNELS

Keywords: viability kernel, viability domain, resistar

This paper denes a specic viability domain called backstop viability kernel and proposes a new algorithm that approximates such sets in a compact constraint set of d dimensions. The algorithm uses a set approximation technique based on a sample chosen in a regular grid of n d vertices covering the considered constraint set. We show that, in some conditions, the result of the algorithm is the direct approximation of the backstop viability kernel by this approximation technique. This theoretical result is illustrated on examples for which the viability kernel is equal to the backstop viability kernel and can be derived analytically, using two set approximation methods: the nearest neighbour, ensuring a convergence rate in O(n -1 ), and a recently developed technique, the recursive simplex stars (resistars), ensuring a convergence rate in O(n -2 ).

1. Introduction. Viability theory [START_REF] Aubin | Viability theory[END_REF][START_REF] Aubin | Viability Theory: New Directions[END_REF] addresses the problem of maintaining a controlled dynamical system inside a given set of states, generally called the constraint set. This framework is particularly relevant for modelling sustainability problems in which the constraint set is interpreted as an acceptable or desirable property of the system that should be sustained [START_REF] Delara | Sustainable Management of Natural Resources[END_REF][START_REF] Chapel | Dening yield policies in a viability approach[END_REF][START_REF] Sabatier | Action versus result-oriented schemes in a grassland agroecosystem: a dynamic modelling approach[END_REF][START_REF] Mathias | Using the viability theory for assessing exibility of forest managers under ecological intensication[END_REF][START_REF] Oubraham | A survey of applications of viability theory to the sustainable exploitation of renewable resources[END_REF]. It appeared also relevant in a variety of engineering problems [START_REF] Tomlin | Computational techniques for the verication and control of hybrid systems[END_REF][START_REF] Kalisiak | Approximate safety enforcement using computed viability envelopes[END_REF][START_REF] Mesmoudi | Coupling geometric analysis and viability theory for system exploration: Application to a living food system[END_REF] and in nance and economics [START_REF] Doyen | Ecological-economic modelling for the sustainable management of biodiversity[END_REF][START_REF] Schulbauer | Economic viability and small-scale sheries: A review[END_REF]. In [START_REF] Martin | The cost of restoration as a way of dening resilience: a viability approach applied to a model of lake eutrophication[END_REF][START_REF] Chapel | Viability and resilience of languages in competition[END_REF][START_REF]Viability and Resilience of Complex Systems: Concepts, Methods and Case Studies from Ecology and Society[END_REF][START_REF] Rougé | Extending the viability theory framework of resilience to uncertain dynamics, and application to lake eutrophication[END_REF] viability theory is at the core of a mathematical denition of resilience, viewed as the capacity of the system to viably restore the property if it has been lost.

Recently, [START_REF] Heitzig | Topology of sustainable management of dynamical systems with desirable states: from dening planetary boundaries to safe operating spaces in the Earth system[END_REF] extended this view to a general theory of sustainable management.

One of the main concepts of viability theory is the viability kernel, the set of states from which the system can remain indenitely in the constraint set. From a state located outside the viability kernel, it is certain that the dynamical system will cross the limits of the constraint set after a nite time, whatever the chosen controls over time. The viability kernel is also important because it is easy to derive from it a variety of control policies keeping the system indenitely inside the constraint set (and actually also inside the viability kernel itself). Generally, it is not possible to determine a viability kernel analytically and several methods provide numerical approximations [START_REF] Bokanowski | An anti-diusive scheme for viability problems[END_REF][START_REF] Bonneuil | Computing the viability kernel in large state dimension[END_REF][START_REF] Coquelin | A dynamic programming approach to viability problems[END_REF][START_REF] Deffuant | Approximating viability kernel with support vector machines[END_REF][START_REF] Chapel | Dening yield policies in a viability approach[END_REF][START_REF] Delara | Sustainable Management of Natural Resources[END_REF][START_REF] Krawczyck | Viability Kernel Approximation, Analysis and Simulation Application -VIKAASAViability Kernel Approximation[END_REF][START_REF] Maidens | Lagrangian methods for approximating the viability kernel in high-dimensional systems[END_REF][START_REF] Alvarez | Viabilititree: a kd-tree framework for viabilitybased decision[END_REF]. In the frequent case of dynamics dened with ordinary dierential equations in a continuous state space, most algorithms start from a regularly distributed sample of points (vertices of a regular grid) covering the constraint set and another one covering the control space. The rst algorithm proposed by Saint-Pierre [START_REF] Saint-Pierre | Approximation of viability kernel[END_REF] uses a discrete approximation of the dynamical system on the vertices of a regular mesh. The algorithm computes a sequence of mesh subsets, until reaching a xed point. The nal set provides a discrete approximation of the viability kernel, which converges to the viability kernel when the time step tends to 0 and n, the number of points by grid axis, tends to innity.

In the particular case of viability kernels dened as epigraphs of a function, [START_REF] Cardaliaguet | Numerical schemes for discontinuous value functions of optimal control[END_REF] showed that the approximation error of the Saint-Pierre algorithm is in O(n -1 ), when choosing a time step in O(n -1 ). In the general case, [START_REF] Rieger | Shadowing and the viability kernel algorithm[END_REF] showed that the approximation error (dened with the Hausdor distance) of the Saint-Pierre algorithm is linear in n -1 and in the time step, if the problem satises some specic conditions (the shadowing property in particular). As far as we know, there is no other algorithm for which the convergence rate has been established, even in restrictive conditions. [START_REF] Alvarez | Viabilititree: a kd-tree framework for viabilitybased decision[END_REF] This manuscript is for review purposes only.

In this paper, we focus on approximating the backstop viability kernel, a subset of the viability kernel, with the aim to get a better convergence rate.

The rst contribution of this paper is theoretical. It denes the backstop viability kernel as the union of all viability domains satisfying a specic property. It introduces the extended discrete time in which the control can change at each time step (as in usual discrete time dynamics) and also when the system is about to leave the constraint set K, and extends the denition of the backstop viability kernel to the extended discrete time. It establishes that, in some conditions, the backstop viability kernels in extended discrete time and in continuous time are equal.

The second contribution is an algorithm approximating backstop viability kernels.

Approximating backstop viability kernels is easier than approximating viability kernels in general, because it boils down to testing long trajectories that change control at most once, starting from points located in a limited subset of the state space. The algorithm uses a set approximation technique which is assumed based on a sample derived from a regular grid of n d points covering K. In some conditions, the nal result is equal to the direct approximation of the exact backstop viability kernel by the set approximation technique. Therefore, the Hausdor distance between the backstop viability kernel and its approximation by our algorithm when n grows is the same as if the set approximation technique was directly applied to the backstop viability kernel itself.

The third contribution is a report on tests of the algorithm on viability problems for which the backstop viability kernel is equal to the viability kernel and can be derived analytically. The tests use two approximation techniques: the nearest neighbour and a recently developed method, the resistar surfaces [START_REF] Deffuant | Recursive simplex stars[END_REF]. These methods guarantee (if the backstop viability kernel holds some smoothness properties) a Hausdor approximation error respectively in O(n -1 ) and O(n -2 ). The results of the tests are in line with the theory.

The remaining of the paper is organised as follows: Section 2 denes the backstop viability kernel in continuous time and extended discrete time control and determines some conditions in which these sets are equal. Section 3 presents the approximation algorithm and the theoretical study of its convergence. Section 4 reports tests of the algorithm convergence on examples. The nal section discusses the contributions of the paper.

2. Backstop viability kernels.

Continuous time.

2.1.1. Viability problem and viability kernel. We consider a controlled dynamical system dened by its state x(t) ∈ R d which can be inuenced by a control

u(t), chosen in a compact set U ⊂ R m , dened by an ordinary dierential equation where ϕ is a continuous function from R d × R m to R d : (2.1) ẋ(t) = ϕ(x(t), u(t)) u(t) ∈ U.
A viability problem from point x 0 is to determine a measurable function u(t) : R + → U , such that the trajectory from x 0 when applying u(.) remains for all t in compact set K ⊂ R n .

Integrating equation 2.1, from initial state x 0 and for a chosen control function u(.), determines the successor of x 0 at time t as follows: 2 This manuscript is for review purposes only.

(2.2)

x x0,u(.) (t) := x 0 + t 0 ϕ(x x0,u(.) (t ), u(t ))dt .
We also denote all the successors (or the trajectory) of x 0 during time interval [0, t] when applying control function u(.) as follows:

(2.3)

X x0,u(.) (t) := t ∈[0,t]
x x0,u(.) (t ).

Definition 2.1. The viability kernel V iab(K) of set K under the dynamics dened by function ϕ and control set U is the set of states x 0 , for which there exists a control function u(.) such that all successors of x 0 when applying u(.) are in K:

(2.4)

V iab(K) = x 0 ∈ K, ∃u(.) : R + → U, ∀t ∈ R + , x x0,u(.) (t) ∈ K ,
where R + denotes the set of positive real numbers.

Determining the viability kernel is important because a variety of control functions u(.) keeping the trajectory indenitely within K can easily be derived from it.

However, the analytical determination of this set is generally impossible and it should therefore be approximated numerically. In this paper, we dene the backstop viability kernel, a subset of the viability kernel which is easier to approximate.

2.1.2. Backstop viability domain and backstop viability kernel. We need a few preliminary notations and denitions:

• For (x, y) ∈ K 2 , D(x, y) denotes the distance between x and y in K, which is the length of the shortest continuous path connecting x to y in K (K is supposed to be path-connected). Of course, if K is convex, D(x, y) = x -y , the usual Euclidean distance.

• For a set A ⊂ K, the distance from x ∈ K to A in K is: D(x, A) = inf y∈A D(x, y). By convention: D(x, ∅) = ∞.
• For any two sets (A, B) included in K, the distance from A to B in K is:

D(A, B) = inf x∈A D(x, B).
• R * + denotes the set of strictly positive real numbers.

• For a couple of sets A and B, such that B ⊂ A, A -B denotes the comple-

mentary set of B in A. • If u(t ) = u for all t ∈ [0, t], then x x0,u(.) (t) and X x0,u(.) (t) are respectively denoted x x0,u (t) and X x0,u (t).
Definition 2.2. The exit time from a set A for an evolution starting from x 0 when applying the control function u(.) is:

τ A (x 0 , u(.)) = inf{t ∈ R + , x x0,u(.) (t) / ∈ A}.
(2.5)

The exit time from set A and point x 0 obtained when applying a constant control u is denoted: τ A (x 0 , u).

Definition 2.3. For any couple (u, v) ∈ U 2 , the constant control from x 0 ∈ K with one possible bounce on the boundary of K, denoted (u, v) x0,K (.), is the function 3

This manuscript is for review purposes only.

from R + to U dened as follows:

(u, v) x0,K (t) = u for 0 ≤ t ≤ τ K (x 0 , u) v for t > τ K (x 0 , u).
(2.6)

The point reached at time t when applying this control function from x 0 is denoted

x x0,(u,v) K (t)
and the corresponding trajectory is denoted X x0,(u,v) K (t).

We can now dene a backstop viability domain.

Definition 2.4. The set V ⊂ K is a backstop viability domain of K under dynamics 2.1 for time threshold θ ∈ R * + i there exists a set B ⊂ V such that:

• D(B, K -V ) > 0;
• For all x 0 ∈ B, there exist

(u, v) ∈ U 2 such that X x0,(u,v) (θ) ⊂ V ;
• For all x 0 ∈ V -B, there exists (u, v) ∈ U 2 , such that:

either for all t ∈ R + , x x0,(u,v) K (t) ∈ V , or there exists t ∈ R * + , t ≥ θ | X x0,(u,v) K (t) ⊂ V and x x0,(u,v) K (t) ∈ B.
B is called a backstop in V for time threshold θ, V -B is called the catch of B in V , D(B, K -V ) the width of the catch of B in V . By convention, the empty set is a backstop viability domain for any catch width and any time threshold.

Proposition 2.5. Let V be a backstop viability domain of K and B a backstop

in V of catch width δ ∈ R * + for time threshold θ ∈ R * + . Let B ⊂ V be such that B ⊂ B and D(B , K -V ) = δ with δ ≤ δ. For any θ ∈ R * + such that 0 < θ ≤ θ, B
is a backstop in V for time threshold θ of catch width δ .

Proof.

Consider x 0 ∈ V -B . Because B ⊂ B , x 0 ∈ V -B.
Therefore, there exist (u, v) ∈ U 2 such that: Proof. Let D Bstop (δ, θ, K) be the set of all backstop viability domains in K including a backstop of catch width at least δ for time threshold θ. Let:

• either ∀t ∈ R+, x x0,(u,v) K (t) ∈ V , • or ∃t ∈ R * + , t ≥ θ such that X x0,(u,v) K (t) ⊂ V and x x0,(u,v) K (t) ∈ B. There- fore t ≥ θ and since B ⊂ B , x x0,(u,v) K (t) ∈ B . Now consider x 0 ∈ B . If x 0 ∈ B, there exists (u, v) ∈ U 2 such that X x0,u (θ) ⊂ V thus X x0,u (θ ) ⊂ V . If x 0 ∈ B -B, then x 0 ∈ V -B and
(2.7)

Consider x 0 ∈ V -B.
There exists V ∈ D Bstop (δ, θ, K) such that x 0 ∈ V , and for all V ∈ D Bstop (δ, θ, K), x 0 / ∈ B V , which implies x 0 ∈ V -B V . Therefore, there exist (u, v) ∈ U 2 such that:

• either ∀t ∈ R + , x x0,(u,v) K (t) ∈ V , then ∀t ∈ R + , x x0,(u,v) K (t) ∈ V , • or ∃t ∈ R * + , t ≥ θ such that X x0,(u,v) K (t) ⊂ V then X x0,(u,v) K (t) ⊂ V because V ⊂ V and x x0,(u,v) K (t) ∈ B V , then x x0,(u,v) K (t) ∈ B, because B V ⊂ B. Consider x 0 ∈ B. There exists V ∈ D Bstop (δ, θ, K) such that x 0 ∈ B V . Therefore, there exists (u, v) ∈ U 2 such that X x0,(u,v) K (θ) ⊂ V , hence X x0,(u,v) K (θ) ⊂ V .
Definition 2.7. If there exists (δ 0 , θ 0 ) ∈ (R * + ) 2 such that:

(2.8)

BSV iab(δ 0 , θ 0 , K) = (δ,θ)∈(R * + ) 2
BSV iab(δ, θ, K),

then
∈ K is denoted E θ (x 0 , K, U ) and
is the set of piece-wise constant functions u(.) : R + → U , dened by the sequence (u j , v j ) j∈N of pairs of elements of U such that, for j ∈ N:

(2.9)

∀t ∈ [jθ, (j + 1)θ[, u(t) = (u j , v j ) xj ,K (t -jθ),
where x j = x x0,u(.) (jθ) and (u j , v j ) xj ,K (t -jθ) refers to denition 2.3.

Definition 2.9. The viability kernel of K in EDT of time step θ under dynamics 2.1, denoted V iab EDT (θ, K), is the set of states x 0 ∈ K for which there exists a control function u(.) in E θ (x 0 , U, K) such that the trajectory from x 0 applying control function u(.) remains indenitely in K:

(2.10) V iab EDT (θ, K) = {x 0 ∈ K, ∃u(.) ∈ E θ (x 0 , U, K) | ∀t ∈ R + , x x0,u(.) (t) ∈ K}.
We now dene a viability domain sewed in K with time step θ under dynamics 2.1. The image of sewing is suggested by considering the trajectory X x0,(u,v) K (θ) as a thread that should be included in K while the stitch x x0,(u,v) K (θ) should belong to the viability domain.

Definition 2.10. A set E ⊂ K is a viability domain sewed in K with time step θ under dynamics 2.1, i:

(2.11)

∀x 0 ∈ E, ∃(u, v) ∈ U 2 | X x0,(u,v) K (θ) ⊂ K and x x0,(u,v) K (θ) ∈ E.
Proposition 2.11. The viability kernel of K in EDT of time step θ is the largest viability domain sewed in K with time step θ.

Proof. For all x 0 ∈ V iab EDT (θ, K) there exists u(.) ∈ E θ (x 0 , U, K) such that for all t ∈ R + , X x0,u(.) (t) ⊂ K. Let u(.) be dened by the innite series {(u 0 , v 0 ), ... , (u j , v j ), ...} with, for j ∈ N, (u j , v j ) ∈ U 2 . We have thus X x0,(u0,v0) K (θ) ⊂ K. Let

x 1 = x x0,(u0,v0) K (θ). The function u (.) dened by {(u 1 , v 1 ), ..., (u j , v j ), ...} is such that for all t ∈ R + , X x1,u (.) (t) ⊂ K. Therefore, x 1 ∈ V iab EDT (θ, K). Therefore V iab EDT (θ, K) is a viability domain sewed in K with time step θ.
Let E ⊂ K be a viability domain sewed in K with time step θ. Consider

x 0 ∈ E. By denition, there exist

(u 0 , v 0 ) ∈ U 2 such that X x0,(u0,v0) K (θ) ⊂ K and x x0,(u0,v0) K (θ) ∈ E.
Consider now, for j ∈ N, {(u 0 , v 0 ), ..., (u j-1 , v j-1 )} ∈ (U 2 ) j , such that, denoting

x p+1 = x xp,(up,vp) K (θ), for p ∈ {0, ..., j -1}, we have: X xp,(up,vp) K (θ) ⊂ K and x xp,(up,vp) K (θ) ∈ E.
In particular,

x j = x xj-1,(uj-1,vj-1) K (θ) ∈ E. Therefore there exist (u j , v j ) ∈ U 2 such that x j+1 = x xj ,(uj ,vj ) K (θ) ∈ E and X xj ,(uj ,vj ) K (θ) ⊂ K because E is a viability domain sewed in K.
Therefore, there exists u(.

) ∈ E(x 0 , U, K) such that for all t ∈ R + , X x0,u(.) (t) ⊂ K, hence x 0 ∈ V iab EDT (θ, K). Therefore E ⊂ V iab EDT (θ, K).

Backstop viability kernel sewed in K.

Definition 2.12. V ⊂ K is a backstop viability domain sewed in K with time step θ under dynamics 2.1 i there exists a set B ⊂ V such that:

• D(B, K -V ) > 0;
• for all x 0 ∈ B, there exists

(u, v) ∈ U 2 such that X x0,(u,v) K (θ) ∈ K and x x0,(u,v) K (θ) ∈ V ;
• For all x 0 ∈ V -B, there exist (u, v) ∈ U 2 , such that:

either ∀j ∈ N, X x0,(u,v) K (jθ) ∈ K and x x0,(u,v) K (jθ) ∈ V , or ∃j ∈ N * , X x0,(u,v) K (jθ) ⊂ K and x x0,(u,v) K (jθ) ∈ B and ∀p < j, x x0,(u,v) K (pθ) ∈ V .
B is called a backstop of V in EDT time step θ and D(B, K -V ) is called the catch width of the backstop V .

It can immediately be seen that a backstop viability domain sewed in K with time step θ is a viability domain sewed in K with time step θ. Proof. The proof is similar to the one of proposition 2.6.

Connection between backstop viability kernels in continuous and

in extended discrete time.

Proposition 2.14. Let V θ be a backstop viability domain sewed in K with time step θ under dynamics 2.1, holding set B as a backstop of catch width δ ∈ R * + . There exists a backstop viability domain V of K holding B as a backstop for time threshold θ such that:

(2.12)

V θ ⊂ V. 6
This manuscript is for review purposes only.

Proof. Let B be a backstop of V θ of catch width δ and let C(B) be the set of points x 0 ∈ K for which there exists (u, v) ∈ U 2 such that:

• either for all t ∈ R + , X x0,(u,v) K (t) ⊂ K;

• or there exists

t ∈ R * + , t ≥ θ such that X x0,(u,v) K (t) ⊂ K and x x0,(u,v) K (t) ∈ B; Let V = B ∪ C(B). For all x 0 ∈ V θ -B, x 0 ∈ C(B), therefore, V θ ⊂ V . Therefore inf x0∈B D(x 0 , K -V ) ≥ δ.
V is thus a backstop viability domain holding B as a backstop of catch width greater or equal to δ, for time threshold θ.

Corollary 2.15. For any (δ, θ) ∈ (R * + ) 2 , we have:

BSV iab EDT (δ, θ, K) ⊂ BSV iab(δ, θ, K).
Proof. The proof comes directly from proposition 2.14 and from the denitions of BSV iab(δ, θ, K) and BSV iab EDT (δ, θ, K).

Proposition 2.16. Let V be a backstop viability domain of K under dynamics 2.1, holding a backstop B of catch width δ for time threshold θ and let:

(2.13)

M = max{ ϕ(x, u) , x ∈ K, u ∈ U }.
For all θ ≤ min(θ, δ 2M ), V is a backstop viability domain sewed in K with time step θ holding a backstop of catch width δ 2 .

Proof.

Consider θ ∈ R * + such that θ ≤ min(θ, δ 2M ). Let B = {x ∈ V, D(x, K - V ) ≥ δ 2 }. For all x ∈ B, D(x, K -B ) ≥ δ 2 . For x 0 ∈ B , because θ ≤ θ, there exists (u, v) ∈ U 2 such that X x0,(u,v) K (θ ) ⊂ K and x x0,(u,v) K (θ ) ∈ V . Consider x 0 ∈ V -B . Because B ⊂ B , x 0 ∈ V -B. Therefore, because B is a backstop in V ,
there exists (u, v) ∈ U 2 such that:

• Either for all t ∈ R + , x x0,(u,v) K (t) ∈ V , then ∀j ∈ N, x x0,(u,v) K (jθ ) ∈ V , and

X x0,(u,v) K (jθ ) ⊂ V , thus X x0,(u,v) K (jθ ) ⊂ K,
• Or there Corollary 2.17. Let M be dened by equation 2.13. For any (δ, θ)

exists t > 0, t ≥ θ | X x0,(u,v) K (t) ⊂ V and y = x x0,(u,v) K (t) ∈ B.
∈ (R * + ) 2 , for any (δ , θ ) ∈ (R * + ) 2 such that δ ≤ δ 2 and θ ≤ min(θ, δ 2M
), we have:

BSV iab(δ, θ, K) ⊂ BSV iab EDT (δ , θ , K).
Proof. The proof comes directly from proposition 2.16 and from the denitions of BSV iab(δ, θ, K) and BSV iab EDT (δ , θ , K).

Proposition 2.18. Let M be dened by equation 2.13. If BSV iab(K), the backstop viability kernel of K exists and BSV iab(K) = BSV iab(δ 0 , θ 0 , K), for (δ 0 , θ 0 ) ∈ (R * + ) 2 , then:

∀δ ≤ δ 0 2 , ∀θ ≤ min(θ 0 , δ 0 2M ), BSV iab EDT (δ, θ, K) = BSV iab(K).
(2.14)

Proof. The proposition is a direct consequence of corollaries 2.15 and 2.17.

Proposition 2.18 provides a connection between the backstop viability kernel and backstop viability kernels sewed in K (in EDT) which is important for the approximation algorithm described in the next section.

3. Approximating backstop viability kernels. Ĥi . The last iteration is reached when no point of the sample is eliminated. The algorithm also stops if

H 1 ∩ Γ(n) = K ∩ Γ(n).
The algorithm uses the following denitions:

• For any set H ⊂ K and any δ ∈ R + , the sets (H) δ and (H) -δ are respectively the erosion and the dilatation of size δ in K:

(H) δ := {x ∈ H | D(x, K -H) ≥ δ}, (H) -δ := {x ∈ K | D(x, H) ≤ δ}.
• The sets B i in and B i out approximate respectively ( Ĥi ) δ and ( Ĥi ) -δ :

B i in := A(( Ĥi ) δ ∩ Γ(n)), B i out := A(( Ĥi ) -δ ∩ Γ(n)).
• The sets H i are recursively dened by:

H 1 ←{x 0 ∈ K | ∃(u, v) ∈ U 2 , X x0,(u,v) K (θ) ⊂ K}, H i+1 ←{x 0 ∈ B i in | testInBstop(x 0 , B i in , B i out )}∪ {x 0 ∈ H i | x 0 / ∈ B i in and testInCatch(x 0 , B i in , B i out )}. 8
This manuscript is for review purposes only.

H 1 ∩ Γ(n) ← {ξ ∈ Γ(n), ∃(u, v) ∈ U 2 , X ξ,(u,v) K (θ) ⊂ K}; if H 1 ∩ Γ(n) = K ∩ Γ(n) then return K; end Ĥ1 ← A(H 1 ∩ Γ(n)); i ← 0; repeat i ← i + 1; B i in ← A(( Ĥi ) δ ∩ Γ(n)); B i out ← A(( Ĥi ) -δ ∩ Γ(n)); H i+1 ∩ Γ(n) ←{ξ ∈ B i in ∩ Γ(n) | testInBstop(ξ, B i in , B i out )}∪ {ξ ∈ H i ∩ Γ(n) | ξ / ∈ B i in and testInCatch(ξ, B i in , B i out )} Ĥi+1 ← A(H i+1 ∩ Γ(n)); until H i ∩ Γ(n) = H i+1 ∩ Γ(n); return ( Ĥi ); Algorithm 3.2 testInCatch(x 0 , B i in , B i out ). Test selecting points in B i out -B i in . Input: x 0 ∈ B i out -B i in , B i in , B i out . b ← ∃(u, v) ∈ U 2 | either ∀j ≤ j A , X x0,(u,v) K (jθ) ⊂ K and x x0,(u,v) K (jθ) ∈ B i out or ∃j < j A X x0,(u,v) K (jθ) ⊂ K and x x0,(u,v) K (jθ) ∈ B i in and ∀p < j, x x0,(u,v) K (pθ) ∈ B i out .
return (b);

Algorithm 3.3 testInBstop(x 0 , B i in , B i out ). Test selecting points in B i in . Input: x 0 ∈ B i in , B i in , B i out . b ← ∃(u, v) ∈ U 2 | X x0,(u,v) K (θ) ⊂ K and y ← x x0,(u,v) K (θ) | y ∈ B i in or y ∈ B i out -B i in and testInCatch(y, B i in , B i out )
return (b);

The basic properties of algorithm 3.1 require the following assumption.

Assumption 1. For any ∈ R * + , there exists n 1 ( ) ∈ N * , such that for all n > n 1 ( ), at any iteration i of the algorithm, for any set H ∈ {H i , ( Ĥi ) δ , ( Ĥi ) -δ }, we have:

(H) ⊂ Ĥ ⊂ (H) -and ( Ĥ) ⊂ H ⊂ ( Ĥ) -,
where Ĥ := A(H ∩ Γ(n)).

Proposition 3.1. Under assumption 1, for all n > n 1 ( δ 4 ), at each iteration of algorithm 3.1, we have:

(H i ) 3δ 2 ⊂ B i in ⊂ (H i ) δ 2 and (H i ) -δ 2 ⊂ B i out ⊂ (H i ) -3δ 2 . 9 Therefore testInCatch(x 0 , B i in , B i out ) = true. If x 0 / ∈ B i in , then x 0 ∈ H i+1 . If x 0 ∈ B i in , this implies testInBstop(x 0 , B i in , B i out ) = true, thus x 0 ∈ H i+1 . • x 0 ∈ (V ) 3δ 2 , then x 0 ∈ B i in (because of 3.1) and there exists (u, v) ∈ U 2 such that X x0,(u,v) K (θ) ⊂ K and y = x x0,(u,v) K (θ) ∈ V . Therefore y ∈ B i out . If y ∈ V -(V ) 3δ
2 , and as shown in the previous case, testInCatch(y,

B i in , B i out ) = true. Else y ∈ (V ) 3δ 2 hence y ∈ B i in therefore testInBstop(x 0 , B i in , B i out ) = true.
Finally, in all cases, x 0 ∈ H i+1 , therefore V ⊂ H i+1 . Assumption 2. For any ∈ R * + , there exists n 2 ( ) ∈ N, such that for all n > n 2 ( ), at any iteration i of the algorithm, for any set H i , we have:

Inclusion of the algorithm output in viability

∀x ∈ K -H i , ∃ξ ∈ (K -H i ) ∩ Γ(n), D(x, ξ) < ,
Definition 3.4. For j ∈ N, let H j be the set of points of K for which the exit time when applying the EDT control of time step θ is greater than jθ:

H j := {x 0 ∈ K, ∃u(.) ∈ E θ (x 0 , U, K), τ K (x 0 , u(.)) ≥ jθ}. (3.5)
We start by showing that the nal set H q dened by the algorithm is included in H j A . Then we distinguish the case where there exists j A ∈ N such that H j A = V iab EDT (θ, K) from the case where there is none.

Proposition 3.5. Under assumptions 1 and 2, for n > max(n

1 ( δ 4 ), n 2 ( δ 2
)), at the last iteration q of algorithm 3.1, we have:

H q ⊂ H j A .
Proof. At iteration i of algorithm 3.1, let

F i = B i in ∩ (K -H j A ). Assume F i = ∅
and, for all x 0 ∈ F i , testInBstop(x 0 , B i in , B i out ) = true. For each point x 0 ∈ F i , we dene U (x 0 ) ⊂ U 4 and j(x 0 ) ∈ N as follows:

• If there exists (u, v) ∈ U 2 such that X x0,(u,v) K (θ) ⊂ K and x x0,(u,v) K (θ) ∈ B i
in . Then U (x 0 ) := ∅ and j(x 0 ) := 0;

• Else there exists (u, v) ∈ K 2 | X x0,(u,v) K (θ) ⊂ K and y := x x0,(u,v) K (θ) ∈ B i out -B i
in , and, because x 0 / ∈ H j A , there exists (u , v ) ∈ U 2 and j ∈ N,

j < j A such that X y,(u ,v ) K (jθ) ⊂ K, x y,(u ,v ) K (jθ) ∈ B i in and for p ∈ N, 0 < p < j, x y,(u ,v ) K (pθ) ∈ B i out . Then U (x 0 ) := (u, v, u , v ) and j(x 0 ) := j.
We dene set E as follows:

E :=F i ∪ {x x x 0 ,(u,v) K (θ),(u ,v ) K (pθ), x 0 ∈ F i , (u, v, u , v ) = U (x 0 ) = ∅, 0 ≤ p ≤ j(x 0 )}.
It can easily be veried that E is a viability domain sewed in K with time step θ.

Therefore E ⊂ H j A . This is impossible because F i ⊂ E and F i ⊂ (K -H j A ), by denition.

Therefore, there exists x 0 ∈ F i , such that testInBstop(x 0 , B i in , B i out ) =false,

thus x 0 ∈ K -H i+1 . Because of assumption 2, there exists ξ ∈ (K -H i+1 ) ∩ Γ(n)
such that D(x 0 , ξ) < δ 2 therefore, ξ ∈ H i because x 0 ∈ B i in and assumption 1 ensures

B i in ⊂ (H i ) δ 2 , hence D(x 0 , K -H i ) ≥ δ 2 .
11

This manuscript is for review purposes only.

Consider now the case when F i = ∅, hence B i in ⊂ H j A , and consider a point

x 0 ∈ H i such that x 0 / ∈ H j A . For (u, v) ∈ U 2 , let y = x x0,(u,v) K (θ). y / ∈ H j A (otherwise we would have x 0 ∈ H j A ), therefore y / ∈ B i in . Assume y ∈ B i out -B i in .
Then, for any

(u , v ) ∈ U 2 , X y,(u ,v ) K (j A θ) ⊂ K is impossible because x 0 / ∈ H j A .
Suppose that there exist (u , v ) ∈ U 2 and j ≤ j A , such that X y,(u ,v ) K (jθ) ⊂ K and

x y,(u ,v ) K (jθ) ∈ B i in . B i in ⊂ H j A implies x y,(u ,v ) K (jθ) ∈ H j A which is impossible because x 0 / ∈ H j A . Therefore, testInBstop(x 0 , B i in , B i out )= false and, for all x 0 ∈ H i such that x 0 / ∈ H j A , x 0 / ∈ H i+1 . Therefore H i+1 ⊂ H j A .
To summarize, while

B i in is not included in H j A , H i+1 ∩ Γ(n) = H i ∩ Γ(n) hence
the algorithm does not stop and at the rst iteration i such that B i in ⊂ H j A , we have

H i+1 ⊂ H j A .
Assumption 1 ensures that the algorithm stops after a nite number of iterations q and that H q ⊂ H i+1 ⊂ H j A .

Definition 3.6. The supremum of steps that the system's trajectory remains in K when applying the EDT control of time step θ from points of K which are not viable for the EDT control of time step θ, is:

j EDT (θ, K) := sup x0∈K-V iab EDT (θ,K) {j ∈ N, ∃u(.) ∈ E θ (x 0 , U, K), X x0,u(.) (jθ) ⊂ K}.
For n ∈ N, the same supremum restricted to

(K -V iab EDT (θ, K)) ∩ Γ(n) is: j EDT (n, θ, K) := sup ξ∈(K-V iab EDT (θ,K))∩Γ(n) {j ∈ N, ∃u(.) ∈ E θ (ξ, U, K), X ξ,u(.) (jθ) ⊂ K}.
Proposition 3.7. Under assumptions 1 and 2, for n > max(n

1 ( δ 4 ), n 2 ( δ 2
)), and if j EDT (θ, K) < ∞ and j A > j EDT (θ, K), at the last step q of algorithm 3.1 we have:

H q ⊂ V iab EDT (θ, K).
(3.6)

Proof. The proof comes directly from proposition 3.5 and from the fact that

H j A = V iab EDT (θ, K) because j A > j EDT (θ, K).
Proposition 3.8. Assume j EDT (n, θ, K) < ∞. Under assumptions 1 and 2,

for n > max(n 1 ( δ 4 ), n 2 ( δ 2
)), and if j A > j EDT (n, θ, K), the last iteration q ∈ N of algorithm 3.1 denes set H q such that:

(3.7)

(H q ∩ Γ(n)) ⊂ (V iab EDT (θ, K) ∩ Γ(n)).
Proof. Proposition 3.5 ensures that for all Assumption 3. There exists (θ 0 , δ 0 ) ∈ (R * + ) 2 such that BSV iab(δ 0 , θ 0 , K) =

ξ ∈ Γ(n), if ξ ∈ H q then ξ ∈ H j A . Because j A > j EDT (n, θ, K), ξ ∈ H j A implies ξ ∈ V iab EDT (θ, K).
BSV iab(K) is the backstop viability kernel of K.

Definition 3.9. Under assumption 3, the maximum number of time steps in K for constant control with one possible bounce starting from points in K -BSV iab(K), the size of a population, which grows or diminishes with the evolution rate x 2 (t) that can be modied by the control. The populations should be kept within some bounds.

denoted j C (θ, K) is: j C (θ, K) = sup x0∈K-BSV iab(K) {j ∈ N, ∃(u, v) ∈ U 2 , X x0,(u,v) K (jθ) ⊂ K}.
The system is written as follows, m 1 , M 1 , m 2 , M 2 , m u , M u being positive parameters:

(4.1)

x 1 (t) = x 1 (t)x 2 (t) x 2 (t) = u(t), with -m u ≤ u(t) ≤ +M u , (x 1 (t), x 2 (t)) ∈ K = [m 1 , M 1 ] × [-m 2 , M 2 ].
In the extension to d dimensions, the model includes d -1 ≥ 2 populations

x 1 , x 2 , .., x d-1 and it uses variable x(t):

x(t) = d-1 i=2 (M 1 -m 1 ) 2 4 -x i (t) - M 1 + m 1 2 2 . (4.2)
The extended viability problem is:

(4.3)      x 1 (t) = y(t) (x 1 (t) + αx(t)) x i (t) = 0, i ∈ {2, .., d -1} x d (t) = u(t), with      -m u ≤ u(t) ≤ M u m 1 ≤ x i (t) ≤ M 1 , i ∈ {1, .., d -1}, -m 2 ≤ x d (t) ≤ M 2 .
Where α is a parameter. In this system, the dynamics of population x 1 nonlinearly depends on the abundance of the other populations which are constant.

The analytical denition of the viability kernel of the d-dimensional problem can be directly derived from the one of the 2D case (provided in [START_REF] Aubin | Elements of viability theory for the analysis of dynamic economics[END_REF]):

V iab(K) = {x ∈ K, ∂V -(x) ≤ x d ≤ ∂V + (x)}, with: (4.4)        ∂V + (x) = 2m u log Mx+αx x1+αx , ∂V -(x) = -2M u log x1+αx mx+αx .
(4.5)

The backstop viability kernel of this problem exists and is equal to the viability kernel.

Consumption problem. In the original 2D problem [3],

x 1 (t) represents the consumption of a primary good and x 2 (t) a critical level of consumption above which the prices can decrease and accelerate consumption and below which, on the contrary the prices increase and decrease the consumption. The critical level x 2 can be modied by a control within some bounds in order to maintain the consumption within some bounds. The system is written as follows (m u , M u , m 1 , M 1 and M 2 are positive):

(4.6)

x 1 (t) = x 1 (t) -x 2 (t) x 2 (t) = u(t), with -m u ≤ u(t) ≤ M u , (x 1 (t), x 2 (t)) ∈ K = [-m 1 , M 1 ] × [0, M 2 ].
Like the population model, we extend the consumption model to d dimensions by considering d -1 ≥ 2 consumption variables x 1 , x 2 , .., x d-1 . We also dene variable

x(t) with equation 4.2. The extended system is:

(4.7)      x 1 (t) = (x 1 (t) + αx(t)) -x d (t) x i (t) = 0, i = 2, .., d -1 x d (t) = u(t), with      -M u ≤ u(t) ≤ M u , -m 1 ≤ x i (t) ≤ M 1 , i ∈ {1, .., d -1}, 0 ≤ x d (t) ≤ M 2 .
The analytical denition of the viability kernel can easily be derived from its expression in the 2D problem (provided in [START_REF] Aubin | Elements of viability theory for the analysis of dynamic economics[END_REF]):

V iab(K) = {x ∈ K, ∂V -(x) ≤ x 1 ≤ ∂V + (x)} with: (4.8)    ∂V -(x) = x d -m u + m u exp -x d mu -αx, ∂V + (x) = x d + M u -M u exp x d -M1 Mu -αx. (4.9)
For this problem also, it can be veried that the backstop viability kernel exists and equals the viability kernel. 4.1.3. Variants of the models with oblique trajectories. In both population and consumption problems, the fact that x i (t) = 0 for i ∈ {2, .., d -1} can be seen as an easy particular case for the resistar approximation because in each 2dimensional grid dened by axes x 1 and x d , the problem to solve is the same as in 2 dimensions. Increasing the dimensionality requires only the interpolation by the resistars between these 2-dimensional classications.

In order to test how the approach performs on a more dicult problem, instead of keeping all trajectories in the 2 D spaces parallel to (b 1 , b d ), (b i , i ∈ {1, .., d} being the canonical basis vectors), for point x = (x 1 , ..., x d ) the trajectory is set in the 2D space generated by vectors (z(x), b d ), with:

z(x) = b 1 + d-1 i=2 z i (x)b i with, setting m(x 1 ) = 1 -β 2 + βx 1 , (4.10)
for i ∈ {2, ..., d -1}, z i (x) = β xi m(x1) , if x i ≤ m(x 1 ), β 1-xi 1-m(x1) , otherwise, (4.11)
where β is a parameter (0 ≤ β < 1) and the problem is rescaled so that

K = [0, 1] d . When x is located on the right line of direction b 1 + β d-1 i=2 b i which includes point ( 1-β 2 , ..., 1-β
2 ), it makes the maximum angle with b 1 . The components z i (x) equal 0

for x i = 0 or x i = 1.
This leads to population or consumption problems in the spaces This manuscript is for review purposes only.

(z(x), b d ) with a constraint set equal to [0, z(x) ] × [0, 1] instead of [0, 1] × [0,
(4.12)

         r = x 1 (t) 2 + x 2 (t) 2 , x 1 (t) = -x 2 (t) + ω(r -r 0 )x 1 (t), x 2 (t) = x 1 (t) + ω(r -r 0 )x 2 (t), x i (t) = -σx i (t), i ∈ {3, .., d}, with -1 ≤ x i (t) ≤ 1, i ∈ {1, .., d}.
Qualitatively, in the plane P generated by canonical basis vectors b 1 = (1, 0, ..., 0)

and b 2 = (0, 1, 0, ..., 0) and which includes point c = (0, .., 0), the dynamics turns around c in P and if the distance to c is higher than r 0 , then the trajectory is a spiral which increases its distance to c while the spiral goes towards c when the distance to c is smaller than r 0 . When the point is exactly at the distance r 0 from c, the radius is kept constant. The parameter ω rules the increase or decrease of the distance to c. When x does not belong to P , the dynamics is the combination of the spiral in the plane generated by (b 1 , b 2 ) which includes x and a translation towards P with a speed proportional to the distance from x to P (ruled by parameter σ). For sake of simplicity, there is no control in this problem (in our framework, the control set includes a single value).

In the spirals problem, j EDT (θ, K), the supremum of number of steps in K for the non-viable points of K, is innite, while it is nite for the two rst problems.

Indeed, consider x ∈ K such that x -c = r 0 + ( > 0); when tends to 0, the trajectory starting from x makes an indenitely increasing number of turns around c before exiting from the constraint set. In the 2D example, we set ω = 0.01 and the system makes a large number of rounds before exiting K even when it is moderately close to the circle of centre c and radius r 0 .

It can be veried that the backstop viability kernel exists and is equal to the viability kernel. The viability kernel can easily be dened analytically: Parameters of population and consumption models.

V iab(K) = {x ∈ K | x 2 1 + x 2 2 ≤ r 2 0 }. ( 4
Model Parameters of spirals model. In both tables the values are given for problems with modied axis scales so that K = [0, 1] d . We checked experimentally that these parameters are such that is dened from G n , the set of vertices of the regular grid covering K, as follows: of the cube. The complete description is available in [START_REF] Deffuant | Recursive simplex stars[END_REF]. that are likely to be the furthest to ∂V iab(K) in order to limit the computation time.

m 1 M 1 m 2 M 2 m u M u α β δ θ j A Popul.2D 0.
V iab(K) = BSV iab(K) = BSV iab EDT (θ,
x ∈ N n (H) ⇐⇒ ∃ξ ∈ G n ∩ H | x -ξ = min ξ ∈Gn x -ξ . ( 4
Theorem 4.2 ([11]). If ∂ K H is a (d -1)-dimensional manifold in K of reach 1 r such that r > √ 2dn -1 , if for all j-dimensional faces F of K, setting H F = H ∩ F , ∂ F H F is a (j -1)-dimensional manifold of reach r F > √ 2jn -1 ,
4.3.1. Resistar approximation. Among the vertices of the simplices dening

∂ K
Ĥq , the barycentres of all boundary points of a grid cube are likely to be the furthest from ∂ K V iab(K). The procedure estimates the distance from such vertices x to the projection of x on ∂ K V iab(K) parallel to the estimated normal vector to the resistar surface at x and it returns the maximum of these distances. It is indeed assumed that the normal vector to ∂ K Ĥq at x is close to the normal vector to ∂V iab(K) at the nearest point to x of ∂V iab(K). The procedure estimates the intersection y of ∂V iab(K) with the segment [x, x + νD M ] or with the segment [x, x -νD M ] (D M being a parameter) by performing successive dichotomies. The case when x + νD M or x -νD M is located outside K requires a specic management. The estimation of the normal vector to ∂ K Ĥq at point x, is based on the estimation (by successive dichotomies) of d -1 anely independent points of ∂ K Ĥq at a given distance from

x, from which the normal vector can be derived. The method includes a specic tratment for the cases when x is on the border of K or very close to it. using the normal to H estimated at the barycentre of the boundary points of the cube. Indeed, the centre of the cube is always on the boundary of Ĥq and is likely to be the point from which ∂V iab(K) is the furthest. The direction normal to ∂ K H at the barycentre of the boundary points of the cube is a reasonable approximation of the direction from the centre of the cube to its nearest point in ∂V iab(K).

Results

. Figure 1 shows the sets Ĥi for all the iterations of algorithm 3.1 applied on the 2D problems. The nal result can be visually compared with the theoretical viability kernel (quantied evaluations of the Hausdor distance are shown on gure 3). Note that the algorithm stops after 4 iterations for the population and consumption problems and after 6 iterations for the spirals problem (even with ω = 0.01). The number of iterations is similar in higher dimensionality. Figure 2 shows examples of nal results on the 3D problems. The smooth non-linearity along the are the same The axes are in a logarithmic scale.

Table 3 shows the estimation of the slopes of the logarithm of the Hausdor distance as a function of the logarithm of n, for the problems in 2 and 3 dimensions.

These results are in good agreement with the theoretical prediction of an Hausdor distance decreasing like n -1 for the approximation with the nearest vertex and like Figure 1. Approximations of the viability kernel in 2D for the three problems (rst row: population, second row: consumption, third row: spirals). Left column: nearest vertex approximation with n = 33 (grid size: 33 2 ), right column: resistar approximation with n = 9 (grid size: 9 2 ). The black curves are the boundaries of the theoretical viability kernel. The approximations Ĥi are represented in darker and darker grey as i increases. The darkest set is the output of algorithm 3.1. 20 This manuscript is for review purposes only.

2.2 x 1 -1 1 x 2 -1 1 x 3 -1 1 x 1 -1 1 x 2 -1 1 x 3 -1 1 Figure 2.
Approximations of the viability kernel in 3D for the three problems (rst row: population, second row: consumption, third row: spirals). Left column: nearest vertex approximation with n = 33 (grid size: 33 3 ), right column: resistar approximation with n = 9 (grid size: Slopes of the linear regression of the logarithm of the estimated Haussdorf distance between the viability kernel and its approximation as a function of the logarithm of n. The R 2 values are all superior to 0.98. n -2 for the approximation with the resistars. 3).

On gure 3, for a given value of n, the error does not change much when increasing the dimensionality. This observation is conrmed on Figure 4, panel (d) showing that the error of the resistar approximation for a grid dened by n = 5 does not vary signicantly when the dimensionality d varies from 3 to 9.

5. Discussion -conclusion. In general conditions, there exists a range of parameter values for which BSViabApp(n, δ, θ, j A ), the output of the approximation 22 This manuscript is for review purposes only. ), for respectively the population, consumption and spirals problems. Each panel represents the intersection of the resistar approximation with 3 hyperplanes (x 3 = 0.02, x 4 = 0.05 and x 5 = 0.08 in K = [0, 1] 6 ). Panel (d): Hausdor distance between viability kernel and its resistar approximation (y axis) for the population, consumption and spirals problems in a grid of size n = 5 and dimensionality varying from 3 to 9 (x axis). algorithm, satises:

BSV iab( 3δ 2 , θ, K) ∩ Γ(n) ⊂ BSViabApp(n, δ, θ, j A ) ∩ Γ(n) ⊂ V iab EDT (θ, K) ∩ Γ(n).
The second inclusion implies that all the points of BSViabApp(n, δ, θ, j A ) ∩ Γ(n) are viable. This is generally not guaranteed with the Saint-Pierre algorithm which provides an approximation of a set which contains the viability kernel. This dierence is important when the objective is to guarantee the viability of a system. When the backstop viability kernel BSV iab(K) exists, under general conditions, there exists a range of parameter values such that: BSViabApp(n, δ, θ, j A ) ∩ Γ(n) = BSV iab(K) ∩ Γ(n).

Overall, these results lead to a convergence rate of the algorithm to the backstop viability kernel which is the same as the convergence rate of the chosen set approximation technique.

When the viability kernel is equal to the backstop viability kernel, the conditions to get this convergence rate are more general. Even when using the nearest vertex approximation, we expect our algorithm to then outperform the current techniques approximating viability kernels for three reasons. Firstly, the convergence to the viability kernel is ensured without decreasing the time step to 0, which is a major dierence. Secondly, for a given time step, our algorithm requires a lower number of iterations, especially when the supremum of time steps in K for the non-viable points is innite (as illustrated on the spiral problem). Thirdly, our algorithm avoids cumulating the error of the successive set approximations taking place in standard algorithms.

When using resistars as set approximation technique, if the best conditions are satised, the convergence rate of our algorithm is like O(n -2 ) which signicantly increases the advantage over the standard methods, converging at best like O(n -1 ).

Indeed, in order to be as accurate as a resistar approximation using a grid of n d points, the standard methods need a grid of at least n 2d points. For instance, we have shown that it is possible to run resistars approximations in 5 dimensions using a grid of 25 5 (about 8 10 6 ) points. In order to reach the same accuracy as the one of these approximations, the standard methods would require a grid of at least (25 2 ) 5 (about 9.5 10 13 ) points, which is not manageable by current standard computers.

Proposition 2 .

 2 [START_REF]Viability and Resilience of Complex Systems: Concepts, Methods and Case Studies from Ecology and Society[END_REF]. For (δ, θ) ∈ (R * + ) 2 , the union of all backstop viability domains sewed in K with time step θ including a backstop of catch width greater or equal to δ is itself a backstop viability domain sewed in K with time step θ including a backstop of catch width greater or equal to δ. It is called the backstop viability kernel sewed in K with time step θ of catch width δ and denoted BSV iab EDT (δ, θ, K).

  kernel in EDT. The next propositions require assumption 2 and a new denition.

3. 3 .

 3 Convergence to the backstop viability kernel. The following proposition requires a new assumption and a new denition.

(3.8) 12 4. 1 .

 121 Tested viability problems. The tests are performed on three viability problems in d dimensions. 4.1.1. Population problem. In the original 2D version [3], x 1 (t) represents

  1]. Their viability kernel in the plane (z(x), b d ) can be derived directly from the viability kernels of the 2D problems. 4.1.4. Spirals problem. Equation4.12 denes the spirals problem in d dimensions, with 0 < r 0 < 1 and ω > 0..

16

 16

(d - 1 )

 1 -dimensional manifolds. The resistars are designed for approximating ∂ K H by hypersurfaces made of (d-1)-dimensional simplices. The rst step for deriving resistar surfaces is determining the boundary points B H (G n ) which are approximations of the intersections between ∂ K H and the edges of the grid. These points are computed on the edges [v, v ] of the grid such that one of the vertices is inside H and the other is outside. The estimation of a boundary point is done by successive dichotomies. The simplices dened in a cube share the barycentre of the boundary points located in the cube as a vertex, and their other vertices are dened similarly in the facets and faces

4. 3 .

 3 Evaluating the Hausdor distance between the viability kernel and its approximations. The procedure evaluates tho distance from Ĥq (result of the approximation algorithm) to V iab(K) (known analytically). It assumes that it is close to the distance from V iab(K) to Ĥq . It focuses on the points of the boundary ∂ K Ĥq

  x 2 axis ruled by parameters α and β appears in the viability kernel approximations of population and consumption problems. Panels (a), (b) and (c) of Figure 4 show the intersection with three chosen hyperplanes of viability kernel resistar approximations of algorithm 3.1 in 6 dimensions with n = 5 (grid size: 5 6 ). On Figure 3, panels (a), (b) and (c) show the estimated Hausdor distance between the viability kernel and its approximation (y axis) by the nearest vertex and by a resistar surface, in dimensionality 2, 3, 4 and 5, and for dierent values of the grid size n (x axis

Figure 3 .

 3 Figure 3. Estimation of the Hausdor distance between the viability kernel and its approximation (y axis) as a function of n dening the grid size as n d (x axis) for the population, consumption and spirals problems in d ∈ {2, 3, 4, 5} dimensions. The dashed lines are the linear regressions on the 2D values (slopes given in table3).

Figure 4 .

 4 Figure 4. Panels (a), (b) and (c): Resistar backstop viability kernel approximations in dimensionality d = 6 and n = 5 (grid size: 5 6), for respectively the population, consumption and spirals problems. Each panel represents the intersection of the resistar approximation with 3 hyperplanes (x 3 = 0.02, x 4 = 0.05 and x 5 = 0.08 in K = [0, 1] 6 ). Panel (d): Hausdor distance between viability kernel and its resistar approximation (y axis) for the population, consumption and spirals problems in a grid of size n = 5 and dimensionality varying from 3 to 9 (x axis).

  this set is called the backstop viability kernel of K under dynamics 2.1 and denoted

	BSV iab(K).
	Section 4.1 describes examples of viability problems for which the backstop via-
	bility kernel exists.
	2.2. Extended discrete time.
	2.2.1. Viability kernel in extended discrete time. We now dene the con-
	trol functions in extended discrete time (abbreviated as EDT), namely the control
	function for which it is possible to change the control at each clock tick (as usual
	in discrete time) and also just before the system leaves the constraint set K. This
	extended discrete time will be used in the algorithm approximating backstop viability
	kernels.

Definition 2.8. The EDT control function set of time step θ ∈ R * + , for constraint set K, dynamics 2.1 and starting from x 0

  3.1. Algorithm and its basic properties. Algorithm 3.1 takes as input: n ∈ N the number of points by axis of the grid, δ ∈ R * + a catch width value, θ ∈ R * + a time step value and j A ∈ N a maximum number of tested steps. It uses algorithm A approximating sets (or classication functions, this is equivalent) from a nite sample Γ(n) of points in K. This sample is based on a regular grid of n d points (n points by axis) covering K for the nearest vertex approximation and a more elaborated sample drawn on edges of the cubes dened by the grid for resistars (which are used in our tests of the algorithm reported in section 4). Other set approximation algorithms such as decision trees or support vector machines could be chosen. Considering a set

E ⊂ K, A(Γ(n) ∩ E) denotes the approximation of E derived from sample Γ(n) by algorithm A. Algorithm 3.1 builds sets H i ∩Γ(n) from which are derived Ĥi = A(H i ∩Γ(n)), the iterative approximations of the backstop viability kernel, with sets B i in as the approximation of their backstops. At each iteration, some points of the sample are eliminated from the denition of the next approximation if they do not pass testInBstop (for points in B i in ) or testInCatch (for points in H i -B i in ). Under some assumptions (see further), the sets B i out contain the sets Ĥi , thus the sets B i out -B i in are wide approximations of the catches. Actually, a point passes testInCatch if there exists a constant control with one possible bounce that denes a trajectory sewed in K with stitches in B i out during j A time steps, or reaching B in in less than j A time steps, without getting out from B i out . As will be shown further, in some conditions, at the last iteration of the algorithm, this test determines with full accuracy if a point belongs to the backstop viability kernel, independently from the precision of the approximation

  Algorithm 3.1 BSVApp(n, δ, θ, j A ). Backstop viability kernel approximation. Input: n ∈ N number of points by grid axis, δ ∈ R * + width of backstop catch, θ ∈ R * + time step, j A ∈ N number of tested steps.

Table 1

 1 .13) 4.1.5. Parameter values.Table 1 breaks down the parameter values used in the tests of the population and consumption models and Table 2 provides the values used in the tests of the spirals model.

  δ, K). Note that m u = 2 in the Population model in d dimensions in order to get a smoother boundary ∂V + (x), ensuring that the conditions of the convergence for resistars are satised in the considered range of values of n.4.2. Set approximation algorithms and their convergence rates. In this subsection, we assume K = [0, 1] d . The tests use two set approximation algorithms. 4.2.1. Nearest vertex. The nearest vertex approximation N n (H) of set H ⊂ K

  .14) Proposition 4.1.For any H ⊂ K such that, for all x ∈ H there exists ξ ∈G n ∩ H such that D(x, ξ) ≤ Proof. By denition, for all point x ∈ K such that x ∈ N n (H) there exists ξ ∈ (G n ∩ H) such that ξ is the point of G n which is the nearest to x,thus D(x, ξ) ≤ ∂H, the boundary of H and ∂ K H = ∂H -(∂H ∩ ∂K) are

		√ n , we have: d
	(4.15)	D h (H, N n (H)) = O(n -1 ).
	√ 2n . d 2n . Moreover, by hypothesis, for all x ∈ H there exists ξ ∈ G n √ d √ Therefore, D(x, H) ≤ such that D(x, ξ) < d n . Therefore, D(x, N n (H)) ≤ √ d n . We get: D h (H, N n (H)) ≤ √ d n .
	4.2.2. Recursive simplex stars (resistars). We assume that set H ⊂ K is a

d-dimensional manifold,

  and if all the boundary points are determined with at least log 2 (n) dichotomies, then the Hausdor distance between H and its resistar approximation decreases like O(dn -2 ).

  4.3.2. Nearest vertex approximation. The set Ĥq is now the nearest vertex approximation. The estimation of the Hausdor distance from Ĥq to V iab(K) is derived from the method dened for the resistar set approximation. It uses indeed the resistar surface denoted H , dened from the boundary points b = (v+v )/2 where

v and v dene a grid edge

[v, v ]

such that v ∈ Ĥq and v / ∈ Ĥq . It estimates the distance from the centres of the cubes containing boundary points of H to ∂V iab(K)

  ). For resistar approximations, the values of n are 7, 9, 13, 17, 25, 33,

	49, 65, 97, 129, 193, 257 for the 2D problems, 7, 9, 13, 17, 25, 33, 49, 65, 97 for the
	3D problems, 7, 9, 13, 17, 25, 33, 49 for the 4D problems and 7, 9, 13, 17, 25 for the
	5D problems. For nearest vertex approximations, the values 7 and 9 are not tested
	because they are too small for dening properly the sets B i in . The other values of n

  9 3 ).

		Nearest v.	Resistars
		2D	3D	2D	3D
	Population	-1.00 -0.97 -2.00 -1.99
	Consumption -1.01 -1.03 -2.01 -1.98
	Spiral	-1.02 -0.94 -2.02 -2.01

21

This manuscript is for review purposes only.

Table 3

 3

This manuscript is for review purposes only.

The reach of ∂ K H is the supremum of ρ such that for any point x of K for which D(x, ∂ K H) = ρ, there is only one point y ∈ ∂ K H such that x -y = ρ[START_REF] Federer | Curvature measures[END_REF].

Proof. For ∈ R * + , < δ 2 , for n > n 1 ( ), at any iteration i, we have:

This implies:

((H i ) ) δ ⊂ ( Ĥi ) δ hence: (H i ) δ+ ⊂ ( Ĥi ) δ . (3.1) and assumption 1 ensures also:

which implies:

Therefore:

Applying assumption 1 to set ( Ĥi ) δ , and setting B in = A(( Ĥi ) δ ∩ Γ(n)) we get:

in ⊂ (( Ĥi ) δ ) -hence (eq: 3.3):

Applying the same reasoning to B i out (mutatis mutandis) and taking = δ Proof. Proposition 3.1 guarantees that, for n > n 1 ( δ 4 ), at each iteration i of algorithm 3.1, we have:

is nite, the procedure reaches

after a nite number of iterations.

Proposition 3.3. Under assumption 1, for n > n 1 ( δ 4 ), at any step i of algorithm 3.1 we have:

Assume now V ⊂ H i and consider x 0 ∈ V . Two cases occur:

or there exists j ∈ N * | X x0,(u,v) K (jθ) ⊂ K and x x0,(u,v) K (jθ) ∈ (V ) 3δ Proposition 3.10. Let M be dened by equation 2.13. If:

• Assumptions 1, 2 and 3 are satised and n > max(n 1 ( δ 4 ), n 2 ( δ 4 )),

• j C (θ, K) < ∞ and j A > j C (θ, K), then for all δ ∈ R * + such that 3δ 2 ≤ δ0

2 , for all θ ≤ min(θ 0 , δ0 2M ), at the last step q of algorithm 3.1 we have:

Proof. The structure of the proof is similar to the one of proposition 3.5. Let

and j(x 0 ) ∈ N as follows:

in and for all p < j,

we can dene U (x 0 ) and j(x 0 ) like we did in the proof of proposition 3.5.

Let E be dened as follows:

By construction, E is a backstop viability domain sewed in K with time step θ holding

in as a backstop of catch width greater or equal to δ 4 . E is not included in V because

4 and by hypothesis,

This is impossible because the choice of θ ensures that V = BSV iab EDT ( δ 4 , θ, K) and thus includes all backstop viability domains sewed in K of catch width greater or equal to δ 4 with time step θ, by denition.

Therefore, there exists

= false, or there exists

In both cases, there exists a point

Suppose now

in and for all p < j, x x0,(u,v) K (pθ) ∈ B i out . Then the set V ∪ {x x0,(u,v) K (pθ), p ∈ {0, .., j}} is a backstop viability domain of time threshold θ admitting B i in as a backstop of catch width at least δ 4 and this backstop viability domain is not included in V . This is impossible. Therefore, testInCatch(x 0 , B i in , B i out ) = false. This implies:

Because of proposition 2.18, the choice of δ and θ ensures V = BSV iab EDT ( 3δ 2 , θ, K).

Proposition 3.3 can be applied because assumption 1 is satised thus V ⊂ H i+1 and

This manuscript is for review purposes only.

Overall, while

the algorithm does not stop. Once (H i ) δ 4 ⊂ V , the algorithm stops after at most one iteration and H i+1 = BSV iab(K). (3.9) Assumption 4. There exists a function (n) :

such that, for any set H satisfying some smoothness properties (see examples in section 4), there exists n 4 ∈ N such that for n > n 4 :

(3.10) Proposition 3.11. Assume the conditions of proposition 3.10. Under assumption 4 and if BSV iab(K) satises the required smoothness conditions and n is large enough, at the nal iteration q of algorithm 3.1, we have:

Proof. Let V = BSV iab(K). Proposition 3.10 ensures H q = V and because of assumption 4,

Proposition 3.12. Assume that BSV iab(K) = BSV iab(δ 0 , θ 0 ) is the backstop viability kernel and is equal to the viability kernel: BSV iab(K) = V iab(K). Assume that the conditions of proposition 3.8 are satised and j A > j EDT (n, θ, K). Under assumption 4 and if V iab(K) satises the required smoothness conditions and n is large enough, then for all δ ∈ R * + such that 3δ 2 ≤ δ0

2 , for all θ ≤ min(θ 0 , δ0 2M ), at the nal iteration q of algorithm 3.1, we have:

(3.12)

Proof. Let V = V iab(K). In the considered conditions,

Note that the case V iab(K) = BSV iab(K) is particularly interesting, because the convergence rate (n) can be ensured even for j EDT (θ, K) = ∞ and j C (θ, K) = ∞, whereas if V iab(K) = BSV iab(K), it requires j C (θ, K) < ∞. The next section reports tests performed in this case.

4. Tests of the convergence rate when viability kernel and backstop viability kernel are equal. In the tests, the trajectories X x0,(u,v) K (jθ) are approximated by the Runge and Kutta method with a time step smaller than θ ensuring that several trajectory points are computed in each grid cube. Hence the approximation error on the trajectory (of the order of n -4 ) is negligible with respect to the set approximation error (at best of the order of n -2 ).