Nicolas Bousquet

APPROXIMATION OF LIMIT STATE SURFACES IN MONOTONIC MONTE CARLO SETTINGS

Keywords: limit state surface, computer model, monotonicity, convexity, structural reliability, rare event probability, random sets, support vector machines 2010 Mathematics Subject Classification: Primary 65C05 Secondary 62L12, 60K10

This article investigates the theoretical convergence properties of the estimators produced by a numerical exploration of a monotonic function with multivariate random inputs in a structural reliability framework. The quantity to be estimated is a probability typically associated to an undesirable (unsafe) event and the function is usually implemented as a computer model. The estimators produced by a Monte Carlo numerical design are two subsets of inputs leading to safe and unsafe situations, the measures of which can be traduced as deterministic bounds for the probability. Several situations are considered, when the design is independent, identically distributed or not, or sequential. As a major consequence, a consistent estimator of the (limit state) surface separating the subsets under isotonicity and regularity arguments can be built, and its convergence speed can be exhibited. This estimator is built by aggregating semi-supervized binary classifiers chosen as constrained Support Vector Machines. Numerical experiments conducted on toy examples highlight that they work faster than recently developed monotonic neural networks with comparable predictable power. They are therefore more adapted when the computational time is a key issue.

Introduction

Due to the increasingly powerful computational tools, so-called computer models g are developed to mimic complex processes (e.g., technical, biological, socio-economical, etc. [START_REF] Casti | Would-be worlds: How simulation is changing the frontiers of science[END_REF][START_REF] De Marchi | Computational and Mathematical Modeling in the Social Sciences[END_REF]). The exploration of multivariate deterministic black-box functions x → g(x) by numerical designs of experiments has become, over the last years, a major theme of research in the interacting fields of engineering, numerical analysis, probability and statistics [START_REF] Smith | Uncertainty Quantification: Theory, Implementation, and Applications[END_REF]. Numerical investigations based on Monte Carlo variance reduction techniques [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF] are often needed since, on the one hand, the complexity of g restricts the use of intrusive explorations, and in the other hand each run of the computer model can be very costly. A number of studies have for common aim to delimit the set of input situations x dominated by a limit state surface Γ defined by the equation g(x) = 0, where the dominance rule is defined by a partial order. In multi-objective optimization frameworks, where g is assimilated to a decision rule, Γ can be viewed as a Pareto frontier delimiting a performance space [START_REF] Figueira | Multiple criteria decision analysis -State of the art -Survey[END_REF]. In structural reliability it is often wanted to assess the measure of the set U -= {x; g(x) ≤ 0}, which can be interpreted as the probability p of an undesirable event when the input vector x is randomized [START_REF] Lemaire | Structural reliability[END_REF].

The Pareto dominance between the inputs x is the natural partial order needed to formalize an assumption of monotonicity placed on g. This hypothesis corresponds to a technical reality in numerous situations encountered in engineering [START_REF] De Rocquigny | Structural reliability under monotony: A review of properties of FORM and associated simulation methods and a new class of monotonous reliability methods (MRM)[END_REF] or decisionmaking [START_REF] Blaszczynski | Inductive discovery of laws using monotonic rules[END_REF], or a conservative approximation to reality in reliability studies. Therefore the exploration of monotonic time-consuming computer models by numerical means has been investigated by a growing number of researchers [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF][START_REF] Da Veiga | Gaussian process modelling with inequality constraints[END_REF][START_REF] Golchi | Monotone Function Estimation for Computer Experiments[END_REF]. The most immediate benefit of monotonicity is surrounding the isotonic limit state surface Γ by two Paretodominated sets delimited by the elements of the numerical design, the measures of which defining deterministic bounds for the probability p [START_REF] De Rocquigny | Structural reliability under monotony: A review of properties of FORM and associated simulation methods and a new class of monotonous reliability methods (MRM)[END_REF]. Recent works focused on the estimation of these bounds and the computation of statistical estimators of p based on Monte Carlo numerical designs [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF].

While random designs appear as powerful tools for surrounding Γ and building estimators of p, the stochastic behavior of the associated random sets and their measures has not been investigated yet. The aim of this paper is to fill in this gap by establishing first convergence results for these objects from random set theory. A corollary of these results, under a convexity argument, is the consistency of an estimator of Γ built from the hull of the numerical design, accompanied with its convergence rate. Such an estimator is proposed by a combination of Support Vector Machines (SVM), which appears as a faster alternative in large dimensions to neural networks specifically developed for isotonic situations [START_REF] Velikova | Solving Partially Monotone Problems with Neural Network[END_REF], while both techniques are usual in structural reliability frameworks [START_REF] Hurtado | An examination of methods for approximating implicit limit state functions from the viewpoint of statistical learning theory[END_REF]. This article is organized as follows. Section 2 details the main concepts and notations. The convergence of a dominated set is examined in Section 3. First we study conditions leading to convergence. Then, under some regularity assumption on the limit state, a rate of convergence for an estimator of Γ is provided in term of Hausdorff distance. Lastly, more precise results are given for two particular cases. Taking into account monotonicity properties, Section 4 focuses on a sequential sampling strategy.

A non naive acceptance-rejection method to simulate in the so-called non-dominated set is provided . Then the convergence rates of bounds is compared numerically with a standard Monte Carlo and a sequential sampling. These theoretical results are derived in Section 5 to build a consistent semi-adaptive SVM-based classifier of the limit state surface Γ respecting isotonic and convexity constraints. Numerical experiments conducted in Section 5.2 complete this analysis by comparing the properties of the classifier with the recently developed monotonic neural networks. The proofs of the new technical results are postponed to Section 7.

Framework

Let consider as in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF] the space U = [0, 1] d . Let X be a random vector uniformly distributed on U and g be a measurable application from U to IR. The function g discriminates two sets of input situations, leading to undesirable and safe events, respectively. These classes are denoted U -= {x ∈ [0, 1] d : g(x) ≤ 0} and U + = {x ∈ [0, 1] d : g(x) > 0}, with [0, 1] d = U -∪ U + . Hence the probability p that an undesirable event may occur is p = P(g(X) ≤ 0) = P(X ∈ U -).

Estimating p by a minimum number of calls to g is a prominent subject of interest in engineering. This number can be drastically diminished when g is assumed to be monotonous. Characterizing the monotonicity of g requires to use the Pareto dominance between two elements of IR d , recalled in next definition.

Definition 1. Let u = (u 1 , . . . , u d ) and v = (v 1 , . . . , v d ) ∈ IR d . We say that u is (resp. strictly) dominated by v, denoting u v (resp. u ≺ v), if for all i = 1, . . . , d, u i ≤ v i (resp. u i < v i ). Assumption 1. Let g : U ⊂ IR d → IR. It is assumed that g is globally increasing, i.e.
for all u, v ∈ U such that u v, then g(u) ≤ g(v).

Remark 1. Obviously as pointed in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF], any monotonic function can be reparametrized to be increasing with respect to all its inputs. This property has for consequence to simplify the exploration of U -and U + and provide deterministic bounds on p. Let x ∈ U -(resp. U + ) and y ∈ U. Since g is increasing, if y x (resp. y x) then y ∈ U -(resp. U + ). More generally, consider a set A of elements on [0, 1] d (for instance a design of numerical experiments) and define

U -(A) = x∈A∩U - {u ∈ [0, 1] d : u x}, U + (A) = x∈A∩U + {u ∈ [0, 1] d : u x}.
Then the monotonicity property implies the following result, since

U -(A) ⊂ U -⊂ [0, 1] d \U + (A).
On the convergence of random sets and measures in monotonic Monte Carlo frameworks 5

x8 x4 x2 x5 x7 x1 x3 x6 {1} d {0} d Γ U + U - Figure 1: Illustration for d = 2 of the designs U -(x 1 , x 2 , x 5 , x 6 ) and U + (x 3 , x 4 , x 7 , x 8 ). Proposition 1. Denote µ the Lebesgue measure on U. For any countable set A, µ(U -(A)) ≤ p ≤ 1 -µ(U + (A)).
Remark 2. In general, the set A will be {X 1 , . . . , X n }, the sets U -(A) and U + (A) would then be random sets and the inequality becomes an almost sure inequality.

The non-dominated subset

U(A) = [0, 1] d \(U -(A) ∪ U + (A)) contains necessarily Γ and
is the only subset of [0, 1] d for which a deeper numerical exploration is needed in view of estimating Γ and p. Formally

Γ = U -\ Ů-∩ U + \ Ů+
where E and E are respectively the closure and the interior of a set E. Next (mild) assumption is required to view the problem of estimating Γ as a binary classification problem with perfectly separable classes, and is needed to study the convergence of the dominated sets to U -and U + when the number of elements of the design increases. A two-dimensional illustration is provided on Figure 1.

Assumption 2. Assume that µ(Γ) = 0 and Γ is simply connex.

In the remainder of the article, a design of numerical experiments {X i , g(X i )} i=1,...,n is considered. To alleviate the notations denote for all n ≥ 1,

U - n = U -{X 1 , . . . , X n } , U + n = U + {X 1 , . . . , X n } ,
and the non-dominated set

U n = [0, 1] d \(U - n ∪ U + n ).
The Lebesgue measures (or hypervolumes) of sets (U - n , U + n ) are then denoted p - n and 1p + n , such that, from Proposition 1, we have with probability one

p - n ≤ p ≤ p + n . (1)
These two bounds can be computed using a sweepline algorithm described in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF] at an exponential cost. A faster approximation method is provided in Section 4.

Two situations are further investigated from the point of view of the convergence of dominated sets and bounds. First, the design is chosen static in [0, 1] d (usual Monte Carlo sampling). Then a sequential Monte Carlo approach is explored, assuming U 0 = [0, 1] d and X i being sampled within the current non-dominated set U i-1 for all i > 0. The convergence of the sets U - n and U + n will be studying via the Hausdorff distance defined as follows.

Definition 2. Let • q be the L q norm on IR d , that is for 0 < q < +∞, x q = ( d i=1 |x i | q ) 1/q , and for q = +∞, x ∞ = max i=1,...,d i=1,...,d

x i . Let (A, B) be two nonempty subsets of the normed vector space ([0, 1] d , • q ). The Hausdorff distance d H,q is defined by

d H,q (A, B) = max(sup y∈A inf x∈B x -y q ; sup x∈B inf y∈A x -y q ).
It is always finite since (A, B) ⊂ [0, 1] d 2 are bounded. The usual form of the Hausdorff distance is defined for q = 2. To alleviate notations, it is now denoted d H (., .) = d H,2 (., .) and . = . 2 .

Convergence results

Almost sure convergence for general sample strategy

The sequences (p - n ) n and (p + n ) n are respectively increasing bounded from above and decreasing bounded from below. Hence, there exists two random variables p - ∞ and p + ∞ such that almost surely p

- n -→ p - ∞ , p + n -→p + ∞ and p - ∞ ≤ p ≤ p + ∞ .
In this section we provide general conditions on the design of experiments that ensures that

p - ∞ = p + ∞ = p. Proposition 2. [Independant Sampling] Let (X k ) k≥1 be a sequence of independent random vectors on [0, 1] d . If there exists ε 1 > 0 such that for all x ∈ Γ ε1 = {u ∈ [0, 1] d , d(u, Γ) < ε 1 } there exists ε 2 > 0 such that n≥1 P(X n ∈ B(x, ε 2 )) = +∞, then (p - n , p + n ) a.
T n+1 = min{j > T n , Y j ∈ U Tn }.
Let (X n ) n≥1 be the sequence of random variables defined for all n ≥ 1 by X n = Y Tn .

Then, conditionally to X 1 , . . . , X n-1 , X n ∼ U(U Tn-1 ), and (p - n , p + n )

a.s.

→ (p, p), where

(p - n , p + n ) are obtained from the sequence (X n ) n≥1 .
It is now natural to quantify the rate of convergence. While it seems difficult to provide an explicit rate in the general case. We will then provide some explicit rate in particular cases.

Rate of convergence under some regularity assumptions

In this section convergence in Haussdorf distance sense of (U - n ) n towards U -is proven, and under regularity assumption for the set U -a rate of convergence is provided. The regularity considered here is the (α, γ)-regularity introduced in [START_REF] Audibert | Fast learning rates for plug-in classifiers[END_REF]. Definition 3. A set K is said (α, γ)-regular if there exist α, γ > 0 such that for all 0 < ε ≤ γ and for all x ∈ K one has

µ(B(x, ε) ∩ K) ≥ αµ(B(x, ε))
where B(x, r) is the open ball with center x and radius r.

This notion was introduced in [START_REF] Audibert | Fast learning rates for plug-in classifiers[END_REF] to prove minimax rates in classification problems. It has also been considered in [START_REF] Cuevas | A plug-in approach to support estimation[END_REF][START_REF] Rodríguez Casal | Set estimation under convexity type assumptions[END_REF] in order to state convergence rate when estimating sets. In particular Proposition 4 below can be seen as an application of Theorem 2 in [START_REF] Rodríguez Casal | Set estimation under convexity type assumptions[END_REF].

Roughly speaking, this condition excludes pathological sets U -, as those having infinitely many sharper and sharper peaks or presenting a jagged border Γ. For instance, it is indirectly proved in [START_REF] Dümbgen | Rates of convergence for random approximations of convex sets[END_REF] that if U -is convex and has non-empty interior then it is regular. In a structural reliability framework this condition is conservative since it traduces the hypothesis that any combination of inputs located between two situations leading to an undesirable event leads also to such an event. Described in [34], so-called r-convex sets that generalize the notion of convexity are also regular [START_REF] Rodríguez Casal | Set estimation under convexity type assumptions[END_REF].

It is traduced in U -by the assumption (or "rolling condition", cf. [34]) that for each boundary point x ∈ Γ, there exists y ∈ U -such that x ∈ B(y, r). In practice the assumption of r-convexity appears mild, as small values of r can produce a large variety of limit state surfaces presenting irregularities [START_REF] Rodríguez Casal | Set estimation under convexity type assumptions[END_REF]. Proposition 4. Let (X k ) k≥1 be a sequence of iid random variables uniformly distributed on [0, 1] d and ( Xk ) k≥1 be a sequence of iid random variables uniformly dis-

tributed on U -. Denote Ũ- n = U -( X1 , . . . , Xn ). Let (F n ) n≥1 be a sequence of measur- able subsets of [0, 1] d such that for all n ≥ 1, U - n ⊂ F n ⊂ [0, 1] d \U + n . Then (1) d H (F n , U -) a.s.
-→ n→+∞ 0 and µ(F n ) a.s.

-→ n→+∞ p.

(2) If U -is regular, then almost surely

d H ( Ũ- n , U -) = O (log n/n) 1/d .
(3) Furthermore, if U + is also regular, and if g is continuous, then almost surely

d H (F n , U -) = O (log n/n) 1/d . (2)
Remark 3. Note, that if in (2), we replace ( Ũn ) by (U - n ), the a.s. bound becomes O (log N 1 /N 1 )

1/d , where N 1 is the number of points X i such that g (X i )) < 0 and follows the binomial distribution with parameter n and p. Moreover the continuity assumption in (3) can be weakened, it is enough (see the proof) that one can cut [0, 1] d following g -1 ({0}) and glue it.

Convergence results in dimension 1

In dimension 1, we provide a complete description of the rate of convergence, in particular we prove that n(pp - n ) and n(p + np) are asymptotically exponentially distributed.

Proposition 5. Let (X n ) n≥1 be an iid sequence of random variables uniformly dis-

tributed on [0, 1] and p ∈ [0, 1]. Define p - n = max i=1,...,n (X i • 1 {Xi≤p} ) and p + n = 1 {maxi=1,...,n(Xi)≤p} + min i=1,...,n (X i • 1 {Xi>p} ). Then n(p -p - n ) L -----→ n→+∞ Exp(1), n(p + n -p) L -----→ n→+∞ Exp(1), E[p + n -p - n ] = 2 n + 1 - 1 n + 1 p n+1 + (1 -p) n+1 , (3)
where Exp(λ) is the exponential distribution with density

f λ (x) = λ exp(-λx)1 {x≥0} .
The cost of adopting a naive sequential strategy to decrease the volume of the nondominated set, by sampling a new design element within [0, 1] d , can be appreciated by coming back to the unidimensional case (d = 1). In this framework, denote

W n = min{j ≥ 1 : X n+j ∈]p - n , p + n [} where (X n+j ) j≥1
is an iid sequence of random variables uniformly sampled on [0, 1]. For all r ≥ 1, one has

P(W n = r) = P(X n+1 / ∈]p - n , p + n [, . . . , X n+r-1 / ∈]p - n , p + n [, X n+r ∈]p - n , p + n [), = E[P(X n+1 / ∈]p - n , p + n [, . . . , X n+r-1 / ∈]p - n , p + n [, X n+r ∈]p - n , p + n [|X 1 , . . . , X n )], = E[(1 -(p + n -p - n )) r-1 (p + n -p - n )].
From linearity of the expectation and Jensen inequality, it comes

E[W n ] = E[(p + n -p - n ) -1 ], ≥ E[p + n -p - n ] -1 , ≥ n + 1 2 from (3) (4)
which is (as expected) a prohibitive cost for large values of n.

More accurate results can be only obtained in some particular case. The following proposition provides a result of the expected Lebesgue measure of the non-dominated set for d = 1.

Proposition 6. Let p - 0 = 0, p + 0 = 1, X 1 ∼ U([p - 0 , p + 0 ]), ξ 1 = 1 {X1≤p} and F 1 = σ{X 1 }. For all n ≥ 1 define conditionally to F n = σ{X k , 1 ≤ k ≤ n}: X n+1 ∼ U([p - n , p + n ]), ξ n+1 = 1 {Xn+1≤p} p - n+1 = p - n + (X n+1 -p - n )ξ n+1 , p + n+1 = p + n -(p + n -X n+1 )(1 -ξ n+1 ).
Then, for n ≥ 1,

1 2 n ≤ E[p + n -p - n ] ≤ 3 4 n . Corollary 2. If p ∈ {0, 1}, E[p + n -p - n ] = 2 -n .
Remark 4. The shrinking convergence rate, which is inversely linear in n in a static Monte Carlo strategy, is significantly improved by becoming exponential when opting for a sequential strategy.

Asymptotic results when Γ = {1} d

Equivalent results are difficult to obtain for greater dimension and for a general limit state Γ. The followings proposition provides asymptotic results in the particulare case

Γ = {1} d . Proposition 7. Assume Γ = {1} d . Let (X k ) k≥1 be a sequence of iid random variables uniformly distributed on [0, 1] d . For 0 < q < +∞ denote A(1, q) = 1 and for d ≥ 2, A d,q = 1 dq d-1 d-1 i=1 B(i/q, 1/q) with B(a, b) = 1 0 t a-1 (1 -t) b-1 dt. For all n ≥ 1, let U - n = U -(X 1 , . . . , X n ).
(1) If 0 < q < +∞ then

(A d,q n) 1/d d H,q U - n , [0, 1] d L -----→ n→+∞ W(1, d).
(

) If q = +∞ then n 1/d d H,∞ U - n , [0, 1] d L -----→ n→+∞ W(1, d), 2
where W(1, d) is the Weibull distribution with scale parameter 1 and shape parameter d having cumulative density function F (t) = 1e -t d for all t ≥ 0.

Proposition 8. Under the assumptions and notations of Proposition 7 we have

E[µ([0, 1] d \U - n )] ∼ n→+∞ log(n) d-1 n(d -1)! . Remark 5. Denoting C ∞ (K) the convex hull of a set K, the result of Proposition 8
is tantamount to the following result provided in [START_REF] Bárány | Random polytopes, convex bodies, and approximation[END_REF]:

E[µ([0, 1] d \C ∞ ({X 1 , . . . , X n }))] ∼ n→+∞ log(n) d-1 n .
When d = 1, the convergence order obtained in Proposition 5 is retreived.

Sequential sampling

Using a naive methods to simulate in a constrained space can be time consuming.

Indeed, Equation 4gives the expected number of simulation needed to simulate in the non-dominated set in dimension 1. A method to simulate in the non-dominated set is provided in the first part of this section. In a second part, the measure of the non-dominated set is compared in function of the strategies of simulation adopted:

standard Monte Carlo or sequential sampling.

Fast uniform sampling within the non-dominated set

Simulating uniformly on [0, 1] d can be time consuming when the aim is to simulate in the non-dominated set. The strategy provided in this section involves to build a subset of [0, 1] d containing the non-dominated set then to simulate uniformly within this subset. After n calls to g let

x in U + n . From construction, the non-dominated set is in [0, 1] d \U + n (x)
, it remains to simulate uniformly on this subset of [0, 1] d . The choice of x is described further in this section. The proposed method use the following lemma coming from a simple acceptance-rejection method.

Lemma 1. Let A be a measurable subset of [0, 1] d with µ(A) > 0 and such that A = A 1 ∪ • • • ∪ A m with A i ∩ A j = ∅ if i = j. Let (X k ) k≥1 be an iid sequence of random vectors satisfying for all i, P(X k ∈ A i ) = µ(A i )\µ(A). Let C be a measurable subset of A, and T = inf{k ≥ 1, X k ∈ C} then X T is uniformly distributed on C.
The construction involves to split up [0, 1] d in hyperrectangles (see Figure 2) which are defined for all i = 0, . . . , 2 d -1 by

Q i (x) = I 1 i × I 2 i × • • • × I d i ,
with

I j i =      [x j , 1] if b i j = 0, [0, x j ] if b i j = 1, where b i = (b i 1 , . . . , b i d ) is equal to i coded in base 2. Since 0 is coded as 00 • • • 0 in base 2 with d numbers, one deduce that Q 0 (x) = [x 1 , 1] × [x 2 , 1] × • • • × [x d , 1] = U + (x). Let (X k ) k≥1 be an iid sequence of random vectors satisfying for all i = 1, . . . , 2 d -1 P(X ∈ Q i (x)) = µ(Q i (x)) 1 -µ(Q 0 (x)) . ( 5
) Define T = inf{k ≥ 1, X k ∈ U n }.Applying Lemma 1 with A = Q 1 (x) ∪ • • • ∪ Q 2 d -1 (x)
and C the non-dominated set U n , one deduce that X T is uniformly distributed on the non-dominated set. The point x is chosen in order to maximise the quantity in Equation ( 5). It must be noticed that more p is close to 0 more the points of Ξ + n are close to {0} d then the probability to have a simulation in the non-dominated set increases.

Numerical study

To complete the previous theoretical studies, a brief numerical comparison of convergence rates obtained for some static and sequential strategies is conducted in this paragraph. The impact of design strategies on the mean measure of the non-dominated set U n is examined through two repeated experiments conducted over the following model, defined such that U -= [0, 1] d , p = 1 and Γ = {1} d for various dimensions 1 shows that the gain provided by a sequential sampling strategy decrease as the dimension increase. For d = 10, the remaining volume obtain with the two strategies of simulations are close.

{1} d {0} d x * Q 0 (x * ) Q 1 (x * ) Q 3 (x * ) Q 2 (x * ) Figure 2: Splitting up [0, 1] d from x * in dimension 2. The hyperrectangle Q 0 (x * ) is in U + . d ∈ {1, 2,

Application : estimating Γ using Support Vector Machines (SVM)

The previous sections provide convergence results for different strategies of simulations without the estimation of the limit state Γ. This section provide an estimator of the sign of g under some convexity assumption on U -. The proposed estimator is in fact a classifier based on linear SVMs. It will then be compared on a toy example, with the monotonic neural networks recently developed in [START_REF] Velikova | Solving Partially Monotone Problems with Neural Network[END_REF].

Theoretical study

Dominated sets (U - n , U + n ) surrounding Γ can be improved by sampling within the non-dominated set U n . In view of improving a naive Monte Carlo approach to estimate p, as in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF], usual techniques like importance sampling or subset simulation should aim at targeting input situations close to Γ [START_REF] Bourinet | Assessing small failure probabilities by combined subset simulation and Support Vector Machines[END_REF]. Such approaches can be guided by a consistent estimation of Γ, under the form of a supervised binary classification rule calibrated from (U - n , U + n ). This classifier has to agree with the following (isotonic) ordinal property of the limit state surface Γ. Proposition 9. Under Assumption 2 For all u, v ∈ Γ such that u = v, u is not strictly dominated by v.

Respecting this constraint a monotonic neural network classifier was recently proposed in [START_REF] Velikova | Solving Partially Monotone Problems with Neural Network[END_REF] and applied in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF] to structural reliability frameworks. While consistent, its computational cost remains high or even prohibitive when the size of the design

X 1 , . . . , X n defining (U - n , U + n ) increases.
Benefiting from a clear geometric interpretation, Support Vector Machines (SVM) offer an alternative to neural networks by their robustness to the curse of dimensionality [START_REF] Hurtado | An examination of methods for approximating implicit limit state functions from the viewpoint of statistical learning theory[END_REF]. A semi adaptive solution can be build from a combination of SVM when U -is convex. Conversely, it can be easily adapted when U + is a convex set.

Assuming U -is convex, any points x of U + can be separated from U -by a hyperplane h x (u) = α + β T

x u (see [START_REF] Rockafellar | Convex analysis[END_REF] Theorem 11.5) that maximises the minimal distance of h x to x and U -. It is also possible to construct h satisfying the ordinal property on Assumption 1 (see Appendix B for more details).

Given the numerical experiment

D n = (X i , y i ) 1≤i≤n ∈ U × {-1, 1}} where y i = 1 if g(X i ) > 0 and -1 otherwise. Let Ξ + n = {X 1 , . . . , X n } ∩ U + and Ξ - n = {X 1 , .
. . , X n } ∩ U -and for any x ∈ Ξ + n define h x as the hyperplane separating x from Ξ - n . The proposed classifier f n is defined by

f n : [0, 1] d → {-1, +1} y →      -1 if for all X ∈ Ξ + n , h X (y) ≤ 0 +1 otherwise.
Denote

F n = {x ∈ [0, 1] d : f n (x) = -1}
the set of all inputs classified as leading to undesirable events by f n .

Theorem 1. Assume U -is convex, then (1)
f n is globally increasing.

(2) For all X ∈ {X 1 , . . . , X n }, sign(g(X)) = f n (X).

(3) The set F n is a convex polyhedron.

(4) Furthermore if (X k ) k≥1 is a sequence of independent random vectors uniformly distributed on [0, 1] d , then

d H (F n , U -) a.s. -→ n→+∞ 0,
and almost surely,

d H (F n , U -) = O (log n/n) 1/d .
Updating the classifier given a new design element X n+1 found in U + can be done by a partially (semi) adaptive principle, illustrated on Figure 3 in dimension 2. To get

f n+1 it is enough to build the hyperplane h Xn+1 which separates X n+1 from Ξ - n+1 = Ξ - n (unfortunately, if X n+1 ∈ U -all
hyperplanes must be rebuild, but this occurs rarely with low probability pp - n when the design is sampled uniformly on [0, 1] d ). If X n+1 x for all x ∈ Ξ - n , the support vectors and the hyperplanes will not be modified. ). The dotted lines represents {x : h X (x) = 0} for some X in U + . Left: there is one point in U + then the dashed line is both {x : h X (x) = 0} and the frontier of the classifier. Right: Dotted lines represent the two sets {x : h X (x) = 0} for X represented by the two white points. Dashed line represent the frontier of the classifier.

Numerical experiments

This section examines first the gain of including the constraint of monotonicity in the calibration of SVM (detailed in Appendix B), with respect to usual SVM and the constrained neural networks proposed in [START_REF] Velikova | Solving Partially Monotone Problems with Neural Network[END_REF] when the associated code is globally monotonic.

First, the case of a linear limit state is studied. Next, the comparison will be made on a function which verifies the monotonicity and convexity properties. The rate of good classification obtained with the monotonic classifier and the monotonic neural networks are compared at each step of the algorithm. 5.2.1. Linear case. Investigating the efficiency of the new classifier, let start the numerical experiments with the linear case. The function g is a hyperplane defined by

h(x) = β 0 + β T x,
with positive parameters which are, in dimension d, uniformly distributed on the unit sphere

S + d = {x ∈ IR d + : x = 1}.
Various sample size N are used to examine the rate of good classification when the monotonicity is taking account or not. For each ordered pair (d, N ), a hyperplane is build with parameters uniformly distributed on S + d . Then, 200 data are simulated on the unit sphere to be predicted. The rate of good classification is compared for linear SVM and linear monotonic SVM. That step is repeated 100 times for each hyperplane, and also repeated with 100 differents hyperplane.

The results are summarised by Table 2 where there is respectively a proportion 0.1, 0.2, 0.3, 0.4 and 0.5 among N which are in {x : g(x) ≤ 0}. In general, taking account of the monotonicity provides better results. As expected, more N is great and d is low more the rate of good classification increase for the two methods. Results are comparable for d = 2 but for a fixed N the difference grows when d increase. Finally, the constrained SVM provides significantly better results for small N . 5.2.2. Convex case. In this section section, the classifier given in Theorem 1 is tested on a toy example. Let U = (U 1 . . . , U d ) be uniformly distributed on [0, 1] d , and denote Let q d,p be the p-quantile of Z d and define the function g(U) = Z dq d,p . This quantile can be deduced from Equation ( 13) in Section 7. Indeed, let t ∈ [0, 1], then

Z d = (U 1 U 2 • • • U d ) 2 .
P(Z d ≤ t) = P(U 1 U 2 • • • U d ≤ √ t) = √ t 0 (-log(x)) d-1 (d -1)! dx.
The function g is globally increasing and the set {u ∈ [0, 1] d , g(u) > 0} is convex. It is shown in this example that the semi-adaptive SVM-based classifier is less timeconsuming than the non-adaptive constrained neural networks for different dimensions.

Nonetheless, when the dimension increases the gain of time decrease.

Discussion

This articles initiates a theoretical counterpart to the increasing developments linked to the exploration of computer models that make use of their geometrical properties.

The original framework is related to structural reliability problems, but the obtained results can be adapted to multi-objective optimisation contexts, characterized by strong constraints of monotonicity. The main results presented here are the consistency and law convergence of the main ingredients (random sets, deterministic bounds, classification tools) used to solve such problems in the common case where the numerical design is generated by a Monte Carlo approach. Therefore such results should be understood as benchmark objectives for more elaborated approaches. For example, use the monotonic estimation of the limit state surface proposed in this article which presents good theoretical and practical properties.

These more elaborated approaches are, obviously, based on sequential sampling within the current non-dominated set. Importance sampling (or more generally variance reduction ) techniques could directly benefit from the current estimation of the limit state surface by the SVM-based tool : the importance distribution must be calibrated to optimize the sampling according to some criterion of volumic reduction.

A first step was done in this direction in [START_REF] Moutoussamy | Comparing conservative estimations of failure probabilities using sequential designs of experiments in monotone frameworks[END_REF]. Such approaches appear more difficult to study theoretically, especially in multidimensional settings. However, a straightforward argument consists in selecting the sequence of importance samplings as a (finite or infinite) subsequence of an infinite nested uniform sampling, such that the consistency properties be kept. The extension of the results presented here to sequential importance sampling, which was already suggested in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF], is the subject of a future work.

Developing specific techniques for fast sampling in a tortured set U n should also be part of this work.

Numerous other avenues for future research can be highlighted. Especially, monotonicity of computer models can be exploited to get deterministic bounds of output quantiles, instead of probabilities (associated to undesirable events too). Quantiles provided for a limit tail value are typical tools in risk analysis. A preliminary result is given beneath in dimension 1, in corollary of Proposition 5, to introduce this theme of research, which will be investigated in another future work.

Corollary 3. Assume that g is continuous, strictly increasing and differentiable at p. Denote q p = inf q∈I R (P(g(X) ≤ q) = p) the p-order quantile of g(X) and g ′ the derivative of g. Then

n(q p -g(p - n )) L -----→ n→+∞ Exp(1/g ′ (p)), n(g(p + n ) -q p ) L -----→ n→+∞ Exp(1/g ′ (p)). 7. Proofs Proof of Proposition 2. Set U - ∞ = k≥1 U - k , U + ∞ = k≥1 U + k and the non- dominated set U ∞ = [0, 1] d \ (U - ∞ ∪ U + ∞ ). We define p - ∞ = µ(U - ∞ ) and p + ∞ = µ([0, 1] d \U + ∞ )
. By inclusion and closure, the bounded sequences (p - n ) n≥1 and (p + n ) n≥1 are respectively increasing and decreasing, then p - n a.s.

→ p - ∞ and p + n a.s. → p + ∞ . Assume U ∞ is a non-empty open set such that µ(U ∞ ) = p + ∞ -p - ∞ > 0. There exist ε 1 > 0 and x 0 ∈ [0, 1] d such that B(x 0 , ε 1 ) ⊂ U ∞ . (6)
Hence no element of the sequence (X k ) k≥1 belongs to B(x 0 , ε 1 ). Now, we introduce the events A n = {X n ∈ B(x 0 , ε 1 )}, the are independent by construction and

n≥1 P(A N ) = n≥1 π d/2 ε d Γ(d/2 + 1) = +∞.
We can then apply to Borel-Cantelli's lemma that ensures that P(lim

sup n A n ) = 1.
Therefore it exists almost surely at least one X k ∈ B(x 0 , ε 1 ), which is contradictory with [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF]. Hence µ(U ∞ ) = 0 and necessarily p - ∞ = p + ∞ = p almost surely.

Proof of Corollary 1. From construction, one has

U - n ⊂ U -then from Proposition 2 d H (U - n , U -) = sup x∈U - n inf y∈U - x -y ≤ µ(U -) -µ(U - n ) = p -p - n , -→ n→+∞ 0.
Proof of Proposition 3. It is an alternative proof to the consistency result given in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF]. For any measurable set

A ⊂ [0, 1] d . If (Y n ) n is an i.i.d. sequence of variable
uniformly distributed on [0, 1] d and T = inf{n, X n ∈ A}. Then, it is a well known fact that Y T is uniformly distributed on A. Hence conditionnaly to Conditionally to X 1 , . . . , X n-1 one has

X n ∼ U(U Tn-1 ).
Denote (q - n ) n≥1 and (q + n ) n≥1 the sequences of bounds obtained from (Y n ) n≥1 . By construction, the sequences (p - n ) n≥1 and (p + n ) n≥1 are subsequences of (q - n ) n≥1 and (q + n ) n≥1 . Then

p + n -p - n ≤ q + n -q - n a.s. → n→+∞ 0.
Proof of Proposition 4. Proof of (1). By triangle inequality, one has

d H (F n , U -) ≤ d H (F n , U - n ) + d H (U - n , U -), ≤ d H ([0, 1] d \U + n , U - n ) + d H (U - n , U -), since U - n ⊂ F n ⊂ [0, 1] d \U + n ≤ d H ([0, 1] d \U + n , U -) + 2d H (U - n , U -)
, by a second triangle inequality.

We know from Corollary 1, that

d H ([0, 1] d \U + n , U - n ) a.s. -→ n→+∞ 0, d H (U - n , U -) a.s. -→ n→+∞ 0, hence d H (F n , U -) a.s. -→ n→+∞ 0. Besides µ(U - n ) ≤ µ(F n ) ≤ µ([0, 1] d \U + n ) by construction, then p - n ≤ µ(F n ) ≤ p +
n , and from Corollary 1, it can be deduced that

µ(F n ) a.s. -→ n→+∞ p,
that concludes the proof of (1). The proof of (2) and (3) are based on the following Theorem due to [START_REF] Rodríguez Casal | Set estimation under convexity type assumptions[END_REF].

Theorem 2. ( [26], Theorem 2). Let K be a compact set on IR d and standard with respect to a measure ν. Let (X k ) k≥1 be a sequence of i.i.d. random variables uniformly distributed on K, and K n be a set such that for all n large enough one has almost surely (X 1 , . . . , X n ) ⊂ K n ⊂ K. Then almost surely

d H (K n , K) = O log n n 1/d .
The notion of standard set generalizes the notion of (α, γ)-regularity that is used here.

Hence any (α, γ)-regular set is also standard in the sense of Theorem 2.

Proof of (2). It is a direct consequence of Theorem 2.

Proof of (3). Application of Theorem 2. In order to do so, we will cut [0, 1] d following g -1 ({0}) and glue it on [0, 1] d-1 in the following way. Set Then for all x ∈ [0, 1],

G : [0, 1] d-1 -→ [0, 1] x →    inf{t ∈ [0, 1], g(x, t) = 0} if {t ∈ [0, 1], g(x, t) = 0} = ∅, 0 if {t ∈ [0, 1], g(x, t) = 0} = ∅,
(X i 1 {Xi≤p} ) ≤ x) = P(X 1 1 {X1≤p} ≤ x) n = 1 -P(X 1 1 {X1≤p} > x) n = (1 -P(p ≥ X 1 > x)) n = (1 -p + x) n . (7)
P(p - n ≤ x) =            0 if x < 0, (1 -p + x) n if 0 ≤ x ≤ p, 1 if x > p. (8)
Let x > 0 and n be large enough such that x ∈ [0, np], then

P(n(p -p - n ) < x) = 1 -(1 -x/n) n -----→ n→+∞ (1 -e -x ),
by recursion E[p + n ] = 2 -n . For p = 1, using X n+1 ∼ U([p - n , 1]) proves the result.

The proof of Proposition 7 relies on the following Lemma that gives the value of the probability density function of a sum of independent beta distributed random variables.

Lemma 2. Denote f d the probability density function (pdf ) of

d i=1 B i , supported over [0, d],
where the B i are iid random variables following the Beta(1/q, 1) distribution.

Then for all x ∈ [0, 1]

f d (x) = C d,q x d/q-1 (10)
with C 1,q = 1 and for d ≥ 2,

C d,q = 1 q d d-1 i=1 B(i/q, 1/q).
Proof of Lemma 2. We proceed by induction. For d = 1, the density given by [START_REF] De Marchi | Computational and Mathematical Modeling in the Social Sciences[END_REF] is the density of Beta(1/q, 1) distribution, that is

f 1 (x) = 1 q x 1/q-1 1 {0≤x≤1} .
Assume there exists k ∈ IN * such that for all 1 ≤ j ≤ k, f j is the density of j i=1 B i (only for x ∈ [0, 1]). Denote f ⋆ g the convolution of functions f and g. Then, for

x ∈ [0, 1],

f k+1 (x) = (f k ⋆ f 1 )(x) = x 0 f k (t)f 1 (x -t)dt, = x 0 1 q k k-1 i=1 B(i/q, 1/q) t k/q-1 1 q (x -t) 1/q-1 dt, = 1 q k+1 k-1 i=1 B(i/q, 1/q) x 0
x 1/q-1 t k/q-1 (1t/x) 1/q-1 dt, = C k+1,q x (k+1)/q-1 which proves the validity of [START_REF] De Marchi | Computational and Mathematical Modeling in the Social Sciences[END_REF].

Proof of Proposition 7. Set X i = (X 1 i , . . . , X j i ). Since U - n ⊂ U -, one has

d H,q (U - n , U -) = max( sup y∈U - n inf x∈U -x -y q ; sup x∈U - inf y∈U - n x -y q ), = sup x∈U - inf y∈U - n x -y q = inf y∈U - n {1} d -y q , = inf y∈{X1,...,Xn} {1} d -y q .
(1). Assume 0 < q < ∞. For t ∈ [0, d 1/q ], using the fact that if X is uniformly distributed on [0, 1] so is 1 -X we have

P(d H,q (U - n , U -) ≤ t) = 1 -   1 -P   d j=1 (X j 1 ) q ≤ t q     n .
Let (α n ) n≥1 be a sequence of real numbers such that α n → +∞ as n → +∞. Hence,

P(α n d H,q (U - n , U -) ≤ t) = 1 -1 -P d i=1 B i ≤ t q α q n n (11)
where each B i = X q i L ∼ Beta(1/q, 1). From Lemma 2, for u ∈ [0, 1],

P d i=1 B i ≤ u = u 0 C d,q x d/q-1 dx = q d C d,q u d/q , where C d,q = 1 q d d-1 i=1
B(i/q, 1/q). Therefore, for n large enough

P(α n d H,q (U - n , U -) ≤ t) = 1 -1 - q d C d,q t d α d n n . Denoting A 1,q = 1 and for d ≥ 2, A d,q = 1 dq d-1 d-1 i=1
B(i/q, 1/q), and choosing α n = (nA d,q ) 1/d it comes

P (A d,q n) 1/d d H,q (U - n , U -) ≤ t -→ n→+∞ 1 -e -t d .
(2). Assume now q = ∞. Let (β n ) n≥1 be a sequence of real numbers such that β n → +∞ as n → +∞. Then for all t ∈]0, 1[

P(β n d ∞ H (U - n , U -) > t) = (1 -P( X 1 ∞ ≤ t/β n )) n = 1 -P( max j=1,...,d X j 1 ≤ t/β n ) n , = 1 -t d /β d n n , choosing β n = n 1/d , it comes P(n 1/d d ∞ H (U - n , U -) > t) -→ n→+∞ e -t d , and
n 1/d d ∞ H (U - n , U -) L -----→ n→+∞ W(1, d).
Proof of Proposition 8. Let (X n ) n≥1 be a sequence of iid random variables uniformly distributed on [0, 1] d . Let U = (U 1 , . . . , U d ) be uniformly distributed on [0, 1] d and independent of the sample (X n ) n≥1 , then

E[µ([0, 1] d \U - n )] = E[E[1 U∈Un |U]], = E[E[1 U X1,...,U Xn |U]], = E[E[1 U X1 |U] • • • E[1 U Xn |U]], = E[E[1 U X1 |U] n ] = E[ 1 -(1 -U 1 ) • • • (1 -U d ) n ], = E[ 1 -U 1 • • • U d n ]. (12)
Using the fact thatlog(U i )

L ∼ Exp(1), it is easy to see that d i=1 U i L = exp(-G d ),
where G d ∼ Gamma(d, 1) with density function

f (x) = x d-1 e -x (d -1)! .
We then easily get that the density function of d i=1 U i is given by

f d (x) = (-log(x)) d-1 (d -1)! 1 {x∈[0,1]} . (13)
From ( 12) and ( 13), it comes

E[µ([0, 1] d \U n )] = E[ 1 -U 1 • • • U d n ] = 1 0 (1 -u) n (-log(u)) d-1 (d -1)! du.
One has

n(d -1)! log(n) d-1 E[µ([0, 1] d \U n )] = n log(n) d-1 1 0 (1 -u) n (-log(u)) d-1 du, and the substitution u = x/n gives n log(n) d-1 1 0 (1 -u) n (-log(u)) d-1 du = 1 log(n) d-1 +∞ 0 (1 -x/n) n (-log(x/n)) d-1 1 {x≤n} dx, = +∞ 0 (1 -x/n) n (1 -log(x)/ log(n)) d-1 1 {x≤n} dx, -→ n→+∞ 1.
The last equation holds by the dominated convergence theorem apply to

f n (x) = (1 -x/n) n (1 -log(x)/ log(n)) d-1 1 {0≤x≤n} ≤ exp (-x)1 {x≥0} . Proof of Proposition 9. Let u ∈ Γ and v ∈ [0, 1] d such that u ≺ v. Assume v ∈ Γ.
There exist ε > 0 such that B((u + v)/2, ε) ⊂ Γ. That implies µ(Γ) > 0, which is impossible by Assumption 2. Then v / ∈ Γ.

Proof of Theorem 1. (1). Obvious since h X (u) ≤ h X (v) for all X ∈ Ξ + n and for all u, v ∈ [0, 1] d such that u v.

(2), [START_REF] Bárány | Random polytopes, convex bodies, and approximation[END_REF].

By construction U - n ⊂ F n ⊂ [0, 1] d \U + n . ( 4
). Since U -is regular (by convexity) and U + is regular (by Lemma 3 below), the result is a consequence of Equation (2). Lemma 3. Let K be a compact, convex subset of [0, 1] d such that for all y ∈ K there is no x ∈ K c = [0, 1] d \K such that x y. Then K c is (α, γ)-regular.

Proof of Lemma 3. Since K is a compact convex set, there exists α, γ > 0 such that K is (α, γ)-regular [START_REF] Audibert | Fast learning rates for plug-in classifiers[END_REF]. Given x ∈ K c , since K is closed convex the projection of

x on K, denoted P K (x), is unique. Let 0 ≤ ε ≤ γ and x ∈ K c , then µ(B(x, ε) ∩ K c ) ≥ µ(B(P K (x), ε) ∩ K c ), ≥ µ(B(P K (x), ε) ∩ K), since K is convex ≥ αµ(B(P K (x), ε)), since K is (α, γ)-regular ≥ αµ(B(x, ε)), then K c is (α, γ)-regular.
Proof of Corollary 3. Since g is differentiable at p, the Delta method and Proposition 5 imply that

n(g(p) -g(p - n )) L -----→ n→+∞ Exp(1/g ′ (p)), n(g(p + n ) -g(q)) L -----→ n→+∞ Exp(1/g ′ (p)).
Since q p is the p-order quantile of g(X), it comes that P(g(X) ≤ q p ) = p. Since g is continuous and strictly increasing, then P(g(X) ≤ q p ) = P(X ≤ g -1 (q p )).

Moreover X ∼ U([0, 1]), then P(X ≤ g -1 (q p )) = g -1 (q p ) = p. It can be deduced that q p = g(p). Then, the two last equations becomes

n(q p -g(p - n )) L -----→ n→+∞ Exp(1/g ′ (p)), n(g(p + n ) -q p ) L -----→ n→+∞ Exp(1/g ′ (p)).
[33] Walther, G. and Lüts D. (1996). Rates of convergence for random approximations of convex sets.

[34] Walther, G. (1999). On a generalization of Blaschke's Rolling Theorem and the smoothing of surfaces, Mathematical methods in the applied sciences, 22, 301-316.

Appendix A. Computing hypervolumes (deterministic bounds)

The computation of bounds (p - n , p + n ) defining in (1) can be conducted exactly or using simulation, in function of the dimension.

A.1. An exact method in dimension d = 2

Consider a design x 1 , . . . , x n ∈ [0, 1] 2 . The first stage is to order the element according to their first (or second) coordinate, such that

x 1 1 < • • • < x 1 n .
Since no design element is dominated by another one, then x 2 i > x 2 i+1 for i ∈ {1, . . . , n -1}. The point x 1 delimit a first rectangle P 1 with the following vertices:

        0 0 x 1 1 0 0 x 2 1 x 1 1 x 2 1         , such that µ(P 1 ) = x 1 1 x 2 1 .
For all i ∈ {2, . . . , n} a new rectangle P i can be defined with the following vertices:

        x 1 i-1 0 x 1 i-1 x 2 i x 1 i 0 x 1 i x 2 i         , such that µ(P i ) = (x 1 i -x 1 i-1 )x 2 i .
The second stage consist to compute the volume of each rectangles:

p - n = x 1 1 x 2 1 + n i=2 (x 1 i -x 1 i-1 )x 2 i .
can be perfectly separated according to y by an hyperplane h defined by

h(x) = β 0 + β T x, (14)
with (β 0 , β) ∈ IR×IR d . In this linear framework, an infinity of hyperplanes can separate perfectly the data. Therefore it is needed to add some constraints to make a unique (optimal) choice of hyperplane. Define the distance of x from Γ : ∆(x, Γ) = inf for i = 1, . . . , n y i (w 0 + w T X i ) ≥ 1.

(

) 16
A number d + 1 of parameters associated to n linear constraints has to be estimated.

Since the problem ( 16) is quadratic with linear constraints, a unique solution can be found. The dual form of the problem is given by the following optimisation problem

                 max a n i=1 a i - 1 2 n i,j=1
a i a j y i y j x T i x j n i=1 a i y i = 0

a 0. ( 17
)
where a = (a 1 , . . . , a n ) is a Lagrange multiplier. Let a = (a 1 , . . . , a n ) be the solution of problem [START_REF] Golchi | Monotone Function Estimation for Computer Experiments[END_REF]. The computer model g can be estimated by

g n (x) = w 0 + n i=1 a i y i X T i x. (18)
such that the data can be linearly separated. Then g can be written as g(x) = β 0 + β T h(x), and using ( 18) and [START_REF] Hastie | The elements of statistical learning[END_REF], it can be estimated by: g n (x) = w 0 + n i=1 a i y i K(X i , x).

On the convergence of random sets and measures in monotonic Monte Carlo frameworks 41

The optimization problem [START_REF] Golchi | Monotone Function Estimation for Computer Experiments[END_REF] becomes

                 max a n i=1 a i - 1 2 n i,j=1
a i a j y i y j K(X i , x j ) and K i,j=1,...,n = y i y j K(X i , X j ), with I n the identity matrix.

B.2. Monotonic SVM

Let g be defined by [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] and assume now that g is globally increasing. Let x = (x 1 , . . . , x d ) ∈ [0, 1] d and x i = (x 1 , . . . , x i-1 , x i + η, x i+1 , . . . , x d ) with η > 0, then x x i and g(x) ≤ g(x i ). It comes that

g x i -g(x) ≥ 0 ⇒ β i ≥ 0,

Figure 3 :

 3 Figure3: Construction of the classifier based on SVM. The plain line represent Γ and black (resp. white) points are in U -(resp. U + ). The dotted lines represents {x : h X (x) = 0} for some X in U + . Left: there is one point in U + then the dashed line is both {x : h X (x) = 0} and the frontier of the classifier. Right: Dotted lines represent the two sets {x : h X (x) = 0} for X represented by the two white points. Dashed line represent the frontier of the classifier.

  The SVM-based classifier and the constrained neural networks are built from the points of a sequential design used to delimit the non-dominated set. At step n, the number of points which delimits this non-dominated set is denoted N n . The comparison is conducted as follows. At each step n ≥ 1, M = 500 random vectors are uniformly distributed on the non-dominated set U n-1 . The rate of good classification on these M random vectors and the time needed to build the two classifiers are stored. Then non-dominated set is updated with a random vector uniformly distributed on it. The comparison is conducted for different dimensions: 2, 3 , 4 and 5 and for p = 0.01 and have been averaged on 40 independent experiments. In Figure 4 the rate of good classification is compared in function of N n . The red values represent the number of times the situation where the non-dominated set is delimited by N n points. For some n ≥ 1, this number is greater than 40. Indeed, the sequence (N n ) n≥1 is not monotonic: at any step n, a new simulation employed to update the non-dominated set U n-1 can dominate one of the point of its frontier. Then the value of N n can be lower than N n-1 . At step n, having N n = n means no points of the sequential design is dominated by one another. The results show that the mean rate of good classification is equivalent for the two classifier and for the tested dimension. But, when the dimension increases the SVM-based classifier is more stable than the constrained neural networks. Denote (T SV M n ) n≥1 and (T N N n ) n≥1 respectively the cumulative time needed to build the SVM-based classifier and the constrained neural networks until step n. The ratio to construct these estimator are compared on Figure 5.

Figure 4 :Figure 5 :

 45 Figure 4: Boxplots of the rate of good classification for d ∈ {3, 4, 5, 6}, in function of n. In red are the sample size used for each boxplot. Left: Monotonic SVM. Right: Constrained neural networks. The results have been averaged on 40 independent experiments.

  and define Γ as the graph of the application x → G(x) -1. Now letΨ : [0, 1] d -→ IR d (x, y) →    (x, y) if g(x, y) ≤ 0, (x, y -1) if g(x, y) > 0. . We define now Ũ = Ψ([0, 1] d ) Γ (we have cut the space [0, 1] d following g -1 ({0}) and glue together [0, 1] d-1 × {0} and [0, 1] d-1 × {1}, see Figure 6). Now by assumption Ũ is a compact set and the random variables (Ψ(X i )) i are i.i.d uniformly distributed in Ũ. The result follows applying Theorem 2. Proof of Proposition 5. Let x ∈ [0, p]. Then P(p - n ≤ x) = P( max i=1,...,n

Figure 6 :

 6 Figure 6: Illustration of the cut given in proof of Proposition 4. Left: Illustration of [0, 1] d after some simulations. Right: Representation of Ũ delimited by the set in gray, Γ and Γ. Ψ(x i ) = x i for i ∈ {3, 4, 7, 8}.

  y∈Γxy = |β 0 + β T x|/ β , from[START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF]. The chosen hyperplane is the solution of the following problem: w = β m β and w 0 = β0 m β[START_REF] Hastie | The elements of statistical learning[END_REF], the problem (15) becomes:

1 2

 1 a T Kad T a Aa * ≥ c, where the symbol * ≥ means that the first constraint is an equality constraint and the other are inequalities constraints, and

for all i = 1 , 2 for i = 1 ,. ( 20 )a. 1 Table 4 :

 1212014 . . . , d. The problem[START_REF] Gadat | Classification with the nearest neighbor rule in general finite dimensional spaces[END_REF] becomes: . . . , n y i (w 0 + w T X i ) ≥ 1 for i = 1, . . . , d w i ≥ 0The associated Lagrangian M(w, w 0 , a, b) can be written asM(w, w 0 , a, b) = 1 2 w 2 -n i=1 a i [y i (w 0 + w T X i ) -1]b T w,where a and b are the Lagrangian multipliers. The dual problem becomes                         max i a j y i y j X T i X jb T n i=1 a i y i X i -Table4shows there is always less constraints and parameters to estimate in the primal problem than in the dual problem. In practice, it is more interesting to solve the primal problem. primal dual Number of parameters d + 1 d + n Number of constraints d + n d + n + Numbers of constraints and dimension of parameters to be estimated for the primal and dual problems in the monotonic SVM framework.The problem[START_REF] Hurtado | An examination of methods for approximating implicit limit state functions from the viewpoint of statistical learning theory[END_REF] can then be rewritten. On the convergence of random sets and measures in monotonic Monte Carlo frameworks 43 n×(d+1) (IR), A 2,1 = I d 0 d ∈ M d×(d+1) (IR), and I d the identity matrix in IR d , 1 T n = (1, . . . , 1) ∈ IR n and 0 T d = (0, . . . , 0) ∈ IR d .

  s.

	→ n→+∞	(p, p).
	Corollary 1. Under Assumption 2 and the conditions of Proposition 2, then
		d H (U -n , U -), d H (U + n , U + )	a.s. -→ n→+∞	(0, 0).
	Example 1. If a sequence (X n ) n∈I N * that is i.i.d. and uniformly distributed on [0, 1] d
	then it satisfies the assumption of Proposition 2. More generally any sequence of iid
	random variables with bounded from below densities on [0, 1] d satisfies the assumption
	of Proposition 2.
	Since the only useful experiments are those which fall into the non-dominated set U n ,
	a sequential Monte Carlo strategy can be carried out, as detailed in next proposition.
	Proposition 3. [Sequential sampling] Let (Y n ) n≥1 be a sequence of iid random vari-
	ables uniformly distributed on [0, 1] d . Let (T n ) n≥1 be a sequence of random variables
	on N defined by T 1 = 1 and for all n ≥ 1

Table 1

 1

	. As

3, 4, 5, 10}. That means, after n simulations, the volume of x∈{X1,...,Xn} {u ∈ [0, 1] d : u x}, is compared, where X 1 , . . . , X n are respectively obtained from a uniform and a sequential sampling strategy. The model is run n = 100 times per experiment and repeated 100 times. The results of these experiments are summarized on

Table 1 :

 1 Mean measures (volumes) of non-dominated set U n after n = 100 runs.

	Remark 6. Table

Table 2 :

 2 Rate

		/78.40	58.14/62.72	53.05/57.65	50.80/53.53	50.49/52.32	50.09/50.88
	N =	84.04/84.19	63.77/68.99	56.66/62.20	52.37/56.83	51.80/55.85	50.40/52.83
	N =	91.22/91.34	71.65/74.13	62.64/67.38	56.25/61.44	54.59/59.91	51.68/55.75
	N =	93.04/93.07	73.92/76.41	65.04/70.05	57.43/63.04	56.16/61.65	52.36/56.87
	N =	92.57/92.74	78.72/80.34	68.95/72.59	60.24/66.02	58.71/63.99	53.73/58.58
	N = 100	94.41/94.41	83.36/84.43	73.66/76.95	64.50/69.03	62.11/67.26	55.54/60.96
	N = 200	97.09/97.09	88.87/89.30	82.09/83.55	72.83/76.07	69.55/73.22	61.24/66.50
	N =	84.86/85.73	62.79/68.05	57.78/63.67	53.08/58.46	52.75/57.92	50.95/54.41
	N =	90.41/90.44	70.27/74.69	61.98/68.49	57.64/63.62	55.88/62.30	52.56/57.87
	N =	93.65/93.67	78.89/80.90	70.51/74.88	62.48/68.38	60.23/66.30	55.44/61.93
	N =	95.51/95.51	82.06/83.78	72.30/76.28	64.84/70.50	62.62/68.35	57.01/63.39
	N =	96.28/96.29	85.36/86.99	77.63/80.25	69.12/73.65	66.41/71.86	59.88/66.15
	N = 100	96.87/96.87	88.44/89.18	81.04/82.94	72.39/76.21	69.71/74.30	62.02/67.80
	N = 200	98.48/98.49	93.37/93.63	88.37/89.19	81.10/83.18	78.02/80.40	69.65/74.13
	N =	86.55/87.42	67.76/73.2	60.98/68.26	56.49/63.74	55.67/62.95	52.86/59.30
	N =	92.39/92.68	74.94/78.36	67.68/73.32	61.13/68.43	58.87/66.38	55.40/62.69
	N =	95.27/95.40	83.20/85	75.22/78.75	66.73/72.70	65.15/71.56	59.24/66.67
	N =	95.35/95.37	85.12/86.46	77.81/81.04	69.39/74.81	66.72/72.89	61.00/68.37
	N =	97.06/97.06	88.88/89.42	82.24/84.51	73.64/78.11	71.40/76.19	63.80/70.19
	N = 100	97.90/97.90	91.55/91.79	85.68/86.71	77.47/80.77	74.69/78.78	66.59/72.88
	N = 200	98.73/98.73	95.16/95.26	91.53/91.94	85.40/86.99	82.95/85.09	74.62/79.08
	N =	88.53/89.14	68.99/74.31	64.20/71.35	59.36/67.83	58.06/66.78	55.23/64.17
	N =	92.10/92.22	77.54/81.02	70.60/75.69	63.92/71.54	62.88/70.56	58.04/66.88
	N =	95.25/95.26	85.06/86.52	77.75/81.31	70.62/76.28	68.48/74.74	62.33/70.43
	N =	97.20/97.20	87.02/87.74	80.91/83.52	72.58/77.60	70.60/75.88	63.98/71.50
	N =	97.51/97.51	90.99/91.54	84.84/86.84	76.90/80.92	74.86/79.32	67.58/74.37
	N = 100	97.98/97.99	92.86/93.06	87.65/88.57	80.11/83.03	77.88/81.74	70.57/76.39
	N = 200	98.81/98.81	95.95/96.07	93.50/93.83	88.27/89.44	86.55/88.18	78.6/82.36
	N =	88.91/89.52	70.12/75.03	64.23/71.39	59.93/68.41	59.17/67.97	56.88/66.48
	N =	93.44/93.59	77.65/80.22	71.63/76.71	65.03/72.34	63.77/71.59	59.74/68.46
	N =	95.85/95.85	86.33/87.23	78.95/81.88	71.58/77.11	69.16/75.35	64.00/71.87
	N =	96.28/96.36	88.25/89.00	81.36/84.56	74.22/79.23	71.72/77.56	65.78/73.53
	N =	97.51/97.51	91.01/91.61	85.96/87.60	78.60/82.51	75.85/80.09	69.66/76.18
	N = 100	98.04/98.04	93.18/93.46	89.03/89.81	81.92/84.43	79.55/82.73	72.08/78.06
	N = 200	99.07/99.08	96.18/96.27	93.57/94.12	89.05/90.12	87.16/88.70	80.09/83.76

of good classification for usual SVM (left) and monotonic SVM (right) in function of d and N . From up to down, there is respectively a proportion of 0.1, 0.2, 0.3, 0.4, 0.5 in the set {x : g(x) ≤ 0}.

Table 3

 3 suggest to solve the primal problem if d ≤ n, and the dual problem in the

	other case.
		primal dual
	Number of parameters d + 1	n
	Number of constraints	n	n + 1

Table 3 :

 3 Numbers of constraints and dimension of parameters to be estimated for the primal and dual problems in the usual SVM framework.

B.1.2. Nonlinear situation. Assume now that D n cannot be separated by a hyperplane but by a nonlinear surface. Let K be a symmetrical positive definite kernel. From Mercer's theorem there exists a transformation h : [0, 1] d → H where H is an Hilbert space with the inner product h(x), h(y) H = K(x, y).

Acknowledgement

The authors express their grateful thanks to Fabrice Gamboa (Institut de Mathématiques de Toulouse) and Bertrand Iooss (EDF R&D and Institut de Mathématique de Toulouse) for their great comments, technical advices and support during the redaction of this article.

Exp [START_REF] Audibert | Fast learning rates for plug-in classifiers[END_REF]. From [START_REF] Cuevas | A plug-in approach to support estimation[END_REF], it comes:

Denote p - n (p) a random variable such that for all x ∈ [0, 1],

. Then, p + n have the same distribution as 1p - n (1p). For all x ∈ [0, 1], one has

, and

It can be deduced that for all n,

Proof of Proposition 6. F n will denote the natural filtration o One has

. By recursion, it can be deduced that

Since

, it comes by recursion that

, it can be deduced from ( 9) that

By recursion, it comes that

To compute p + n it is enough to transform each x i into 1x i , to compute p+ n with the previous approach. Then p + n = 1 -p+ n .

A.2. An accelerated Monte Carlo method in upper dimensions

A sweepline algorithm described in [START_REF] Bousquet | Accelerated Monte Carlo estimation of exceedance probabilities under monotonicity constraints[END_REF] allows to compute the bounds in any dimension, but at an exponential cost. An alternative approach is using a Monte Carlo method. Considering an iid uniform sample ŪN = U 1 , . . . ,

The computation can be strongly accelerated by adapting the order of evaluation to the monotonic context, using the following algorithm (easily adaptable to estimating where X i ∈ [0, 1] d and y i = sign(h(X i )), with sign(x) = 1 if x > 0 and -1 otherwise,


Zhenhua Wang

Mickael Rodrigues

Didier Theilliol
email: didier.theilliol@univ-lorraine.fr

Yi Shen

Fault estimation filter design for discrete-time descriptor systems

Keywords: Actuator faults, Descriptor systems, Discrete-time, Fault estimation, Filter design, LMI

published or not. The documents may come

Fault estimation filter design for discrete-time descriptor systems

Introduction

The model-based fault diagnosis techniques have received considerable attention during the past decades. Various methods have been proposed in the literature, e.g. observer-based fault detection method [START_REF] Seliger | Fault diagnosis by disturbance decoupled nonlinear observers[END_REF][START_REF] Hou | Fault detection and isolation observers[END_REF], parity space approach [START_REF] Gertler | Fault detection and isolation using parity relations[END_REF], H ∞ fault detection filter [START_REF] Zhong | An LMI approach to design robust fault detection filter for uncertain LTI systems[END_REF][START_REF] Zhao | H∞ fault detection for linear discrete time-varying systems with delayed state[END_REF], and the recently presented finite frequency fault detection method [START_REF] Wang | A finite frequency domain approach to fault detection for linear discrete-time systems[END_REF][START_REF] Yang | Fault Detection for T-S Fuzzy Discrete Systems in Finite-Frequency Domain[END_REF][START_REF] Long | Fault detection for a class of networked control systems with finite-frequency servo inputs and random packet dropouts[END_REF]. For more comprehensive knowledge on model-based fault diagnosis techniques, see the monographs [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF][START_REF] Ding | Model-based fault diagnosis techniques: Design schemes, algorithms, and tools[END_REF][START_REF] Isermann | Fault-Diagnosis Systems: An Introduction from Fault Detection to Fault Tolerance[END_REF], survey papers [START_REF] Frank | Fault diagnosis in dynamic systems using analytical and knowledge based redundancyA survey and some new results[END_REF][START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF][START_REF] Dai | From model, signal to knowledge: A data-driven perspective of fault detection and diagnosis[END_REF] and references therein.

On another research front, descriptor systems appear in many fields such as power systems, electrical networks, and mechanical systems [START_REF] Dai | Singular Control Systems[END_REF][START_REF] Duan | Analysis and design of descriptor linear systems[END_REF]. In the last few decades, many approaches have been proposed to design observers for descriptor systems [START_REF] Dai | Observers for discrete singular systems[END_REF][START_REF] Darouach | Design of observers for descriptor systems[END_REF][START_REF] Hou | Observer design for descriptor systems[END_REF][START_REF] Lu | Full-order and reduced-order observers for Lipschitz descriptor systems: The unified LMI approach[END_REF][START_REF] Wang | Observer design for discrete-time descriptor systems: An LMI approach[END_REF][START_REF] Hammouri | An observer design for a class of implicit systems[END_REF]. References [START_REF] Dai | Observers for discrete singular systems[END_REF][START_REF] Darouach | Design of observers for descriptor systems[END_REF][START_REF] Hou | Observer design for descriptor systems[END_REF] have studied the observer design problem in linear descriptor systems. Observer design methods for nonlinear descriptor systems have been reported in [START_REF] Lu | Full-order and reduced-order observers for Lipschitz descriptor systems: The unified LMI approach[END_REF][START_REF] Wang | Observer design for discrete-time descriptor systems: An LMI approach[END_REF][START_REF] Hammouri | An observer design for a class of implicit systems[END_REF].

Compared with the results on observer design for the fault-free descriptor systems, there are limited papers on fault diagnosis methods for descriptor systems [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF][START_REF] Duan | Robust fault detection in descriptor linear systems via generalized unknown input observers[END_REF][START_REF] Gao | Fault estimation and fault-tolerant control for descriptor systems via proportional, multiple-integral and derivative observer design[END_REF][START_REF] Wang | Actuator fault estimation for a class of nonlinear descriptor systems[END_REF][START_REF] Marx | Design of observers for Takagi-Sugeno descriptor systems with unknown inputs and application to fault diagnosis[END_REF][START_REF] Chadli | H-/H∞ fault detection filter design for discrete-time Takagi-Sugeno fuzzy system[END_REF][START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF]. [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF] has proposed a fault estimation method for linear descriptor systems by using on-line learning methodology. In [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF], the state variables are required to be fully measurable. However, in many situations, it is too expensive or even impossible to measure all the states. [START_REF] Duan | Robust fault detection in descriptor linear systems via generalized unknown input observers[END_REF] has considered the robust fault detection problem for linear descriptor systems via generalized unknown input observers. Nevertheless, as pointed out in [START_REF] Duan | Robust fault detection in descriptor linear systems via generalized unknown input observers[END_REF], these generalized observers are difficult to implement due to the fact that the derivatives of known signals may be involved in realizing the generalized observers. More recently, [START_REF] Gao | Fault estimation and fault-tolerant control for descriptor systems via proportional, multiple-integral and derivative observer design[END_REF] has studied fault estimation and fault-tolerant control for descriptor systems by using proportional, multipleintegral and derivative observer. [START_REF] Wang | Actuator fault estimation for a class of nonlinear descriptor systems[END_REF] has proposed an actuator fault estimation method for a class of nonlinear descriptor systems via neural adaptive observer design. It should be noticed that most of the aforementioned results focus on fault diagnosis in continuous-time descriptor systems, only [START_REF] Marx | Design of observers for Takagi-Sugeno descriptor systems with unknown inputs and application to fault diagnosis[END_REF][START_REF] Chadli | H-/H∞ fault detection filter design for discrete-time Takagi-Sugeno fuzzy system[END_REF][START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF] have considered fault diagnosis methods for discrete-time descriptor systems. In [START_REF] Marx | Design of observers for Takagi-Sugeno descriptor systems with unknown inputs and application to fault diagnosis[END_REF],

an LMI-based unknown input observer design method is proposed and applied to fault detection and isolation for Takagi-Sugeno descriptor systems. [START_REF] Chadli | H-/H∞ fault detection filter design for discrete-time Takagi-Sugeno fuzzy system[END_REF] has proposed an H -/H ∞ fault detection filter for a class of discrete-time descriptor systems described by Linear Parameter-Varying (LPV) form with a globally Lipschitz term. However, the methods in [START_REF] Marx | Design of observers for Takagi-Sugeno descriptor systems with unknown inputs and application to fault diagnosis[END_REF] and [START_REF] Chadli | H-/H∞ fault detection filter design for discrete-time Takagi-Sugeno fuzzy system[END_REF] only address the fault detection problem but the fault magnitudes cannot be estimated. As pointed out in [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF], accurate and timely fault estimation is an important basis for active control reconfiguration. For regular systems, several model-based techniques such as parity space approach, adaptive observer and sliding mode observer have been proposed to achieve fault estimation [START_REF] Zhong | Parity space-based fault estimation for linear discrete time-varying systems[END_REF][START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF][START_REF] Zhang | Fast fault estimation and accommodation for dynamical systems[END_REF][START_REF] Shen | Actuator fault-tolerant robust control for a linear system with adaptive fault identification[END_REF][START_REF] Alwi | Sliding mode estimation schemes for incipient sensor faults[END_REF][START_REF] Jiang | Fault estimation and accommodation for linear MIMO discrete-time systems[END_REF]. However, little work has been done to fault estimation for descriptor systems [START_REF] Gao | Fault estimation and fault-tolerant control for descriptor systems via proportional, multiple-integral and derivative observer design[END_REF][START_REF] Wang | Actuator fault estimation for a class of nonlinear descriptor systems[END_REF][START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF], especially for the discrete descriptor systems. In [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF], a fault estimation method based on Proportional-Integral (PI) observer is proposed to deal with fault estimation for discrete-time LPV descriptor systems. However, the method in [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF] requires the faults to be constant or slow varying, which is a restrictive condition.

In this paper, a novel fault estimation filter design method without the constant fault assumption is proposed for discrete-time linear descriptor systems. A restricted system equivalent model of the considered descriptor system is firstly established. We firstly construct a restricted system equivalent model of the considered descriptor system, and then design a fault estimation filter based on the restricted system equivalent representation. The filter design method is inspired by [START_REF] Keller | Fault isolation filter design for linear stochastic systems[END_REF] such that the ith component of the fault estimation vector asymptotically converges to the ith fault with a time delay equals to its fault detectability index. Moreover, this paper proposes a robust filter design method to deal with the effect of the unknown disturbance, which is not considered in [START_REF] Keller | Fault isolation filter design for linear stochastic systems[END_REF]. In this paper, design conditions for the proposed filter are formulated as a linear matrix inequality (LMI) feasibility problem, which can be easily solved by using standard LMI toolboxes. The main novelty of this paper lies on the following aspects. First, the proposed approach does not assume the faults to be constant. Therefore, it has a larger application scope than the method based on PI observers. Secondly, the form of the proposed filter is similar with the normal state observer, which is easy to implement. Moreover, this paper proposes a robust filter which is able to attenuate the effect of unknown disturbance and to provide a robust fault estimation. In addition, the proposed filter design method provides a dynamic detection threshold for the residuals in the fault-free situation. This can be used to achieve robust fault detection.

The rest of this paper is organized as follows. The problem formulation is presented in Section 2. The fault estimation filter is designed in Section 3. In Section 4, a numerical example is given to illustrate the effectiveness of the proposed approach. Conclusions are given in Section 5.

Problem statement and preliminaries

Consider the following discrete-time descriptor system with actuator faults

       Ex(k + 1) = Ax(k) + Bu(k) + F n(k) + Dw(k) y(k) = Cx(k) (1)
where x(k) ∈ R n is the state vector, u(k) ∈ R p is the input vector, w(k) ∈ R d denotes the unknown disturbance and y(k) ∈ R m is the output vector. E ∈ R n×n may be rank deficient, i.e. rank(E) = r ≤ n, A, B, C and D are known matrices with appropriate dimensions.

n(k) = [n 1 (k) n 2 (k) • • • n q (k)] T ∈ R q is the vector of fault magnitudes and F = [f 1 f 2 • • • f q ] ∈ R n×q
is the fault distribution matrix. In the following of the paper, we assume that rank(C) = m, rank(F ) = q, q ≤ m, and rank

    E C     = n (2)
In this paper, a filter is designed to detect and estimate the faults in descriptor system [START_REF] Seliger | Fault diagnosis by disturbance decoupled nonlinear observers[END_REF].

To this end, some usefull preliminaries are presented in this section. First, recall the definition of restricted system equivalent.

Definition 1 [START_REF] Dai | Singular Control Systems[END_REF]. The following descriptor systems

       Ex(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (3) and        Ẽ x(k + 1) = Ãx(k) + Bu(k) y(k) = C x(k) (4)
are called Restricted System Equivalent (RSE) if there exist two nonsingular matrices Q and P such that

x(k) = P x(k) (5)
and

QEP = Ẽ, QAP = Ã, QB = B, CP = C (6)
In this paper, it is required that T ∈ R n×n is a nonsingular matrix satisfing the following equiation

T E + N C = I n ( 7
)
where N ∈ R n×m is an arbitrary matrix.

In order to ensure the existence of such a matrix T , the following Lemma is introduced.

Lemma 1. Under the assumption (2), there exist a nonsingular matrix T ∈ R n×n and a matrix

N ∈ R n×m satisfying T E + N C = I n (8)
Proof. Since rank(E) = r, there exist two nonsingular matrices Q 0 and P 0 such that

    x 1 (k) x 2 (k)     = P -1 0 x(k), Q 0 EP 0 =     I r 0 0 0     , Q 0 AP 0 =     A 11 A 12 A 21 A 22     , Q 0 B =     B 1 B 2     , Q 0 F =     F 1 F 2     , Q 0 D =     D 1 D 2     , CP 0 = [ C 1 C 2 ] (9)
then the system (1) is restricted system equivalent to

                       Q 0 EP 0     x 1 (k + 1) x 2 (k + 1)     = Q 0 AP 0     x 1 (k) x 2 (k)     + Q 0 Bu(k) + Q 0 F n(k) + Q 0 Dw(k) y(k) = CP 0     x 1 (k) x 2 (k)     (10)
Under the assumption (2), it is easy to shown that rank

    Q 0 EP 0 CP 0     = rank         I r 0 0 0 C 1 C 2         = rank     I r 0 C 1 C 2     = n (11) Let H = [0 I m ] T , we have Q 0 EP 0 + HCP 0 =     I r 0 C 1 C 2     (12)
Since P 0 is nonsingular, the following matrix

G Q 0 E + HC =     I r 0 C 1 C 2     P -1 0 ( 13
)
is also a nonsingular matrix.

Let

T = G -1 Q 0 , N = G -1 H (14)
we get a nonsingular matrix T and a matrix N such that (8) holds.

Remark 1. Lemma 1 guarantees the existence of a nonsingular matrix T which also satisfies the constraint [START_REF] Long | Fault detection for a class of networked control systems with finite-frequency servo inputs and random packet dropouts[END_REF]. However, it is not very computationally efficient. Therefore, the following Lemma is introduced to provide a convenient alternative to find such a matrix T .

Lemma 2 [START_REF] Ben-Israel | Generalized inverses: Theory and applications[END_REF]. Given matrices

X ∈ R a×b , Y ∈ R b×c , and Z ∈ R a×c . If Y = c, the general solution of X Y = Z is given by X = ZY † + S[I b -YY † ] ( 15
)
where S ∈ R a×b is an arbitrary matrix. Herein and throughout this paper, the symbol † is use to denote the pseudo-inverse of a matrix.

Since rank

    E C     = n,
there exist matrices T and N satisfying [START_REF] Yang | Fault Detection for T-S Fuzzy Discrete Systems in Finite-Frequency Domain[END_REF]. Using Lemma 2, the general solution for T and N is given by

T = MR T , N = MR N ( 16
)
where

M =     E C     † + S      I n+m -     E C         E C     †      , R T =     I n 0     , R N =     0 I m     (17)
and SR n×(n+m) is an arbitrary matrix, which provides us the degrees of design freedom. By choosing matrix S, we can determine a nonsingular matrix T which satisfies equation ( 7).

Since T is a nonsingular matrix, descriptor system (1) is RSE to the following representation

       T Ex(k + 1) = T Ax(k) + T Bu(k) + T F n(k) + T Dw(k) y(k) = Cx(k) (18)
Motivated by [START_REF] Keller | Fault isolation filter design for linear stochastic systems[END_REF], this paper introduces the following definitions of fault detectability indexes and fault detectability matrix, which will be used in the sequel.

Definition 2. The fault detectability indexes of system [START_REF] Darouach | Design of observers for descriptor systems[END_REF] are defined as

ρ = {ρ 1 , ρ 2 , • • • , ρ q }
where

ρ i = min{v : C(T A) v-1 T f i ̸ = 0, v = 1, 2, • • •} (19)
Without loss of generality, it is assumed that the system (18) has finite fault detectability indexes and s min{ρ 1 , ρ 2 , • • • , ρ q } is used to represent the maximum value of fault detectability indexes.

Definition 3.

With the assumption that the system (18) has finite fault detectability indexes, the fault detectability matrix is defined as

Ψ = [(T A) ρ 1 -1 T f 1 (T A) ρ 2 -1 T f 2 • • • (T A) ρq-1 T f q } (20)
Besides, the discrete-time Bounded Real Lemma will be used in this paper.

Lemma 3 [START_REF] Hsiung | Lyapunov inequality and bounded real lemma for discrete-time descriptor systems[END_REF]. Consider the state-space system

       x(k + 1) = Ax(k) + Bw(k) y(k) = Cx(k) + Dw(k) (21)
For a given scalar γ > 0, the following statements are equivalent:

(i) A is a stable matrix and ∥C(zI -A) -1 B + D∥ < γ.

(ii) There exists a symmetric P satisfying P > 0 and the following matrix inequality    

A T P A -P + C T C A T P B + C T D B T P A + D T C -γ 2 I + B T P B + D T D     < 0 (22)

Fault estimation filter design

Based on the RSE representation [START_REF] Darouach | Design of observers for descriptor systems[END_REF], the following filter is constructed for descriptor system

(1)                z(k + 1) = T Ax(k) + T Bu(k) + L(y(k) -C x(k)) x(k) = z(k) + N y(k) r(k) = W (y(k) -C x(k)) (23)
where z(k) ∈ R n is an intermediate state vector, x(k) ∈ R n is the state estimation vector and r(k) ∈ R q denotes the residual vector. In [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF], matrices T and N can be determined by using [START_REF] Duan | Analysis and design of descriptor linear systems[END_REF], L ∈ R n×m and W ∈ R q×m are matrices to be designed.

Define the state estimation error and output estimation error as

e(k) = x(k) -x(k) (24)
and

ϵ(k) = y(k) -C x(k) (25)
respectively. Then, by using ( 7), ( 18) and ( 23), the error dynamic system is obtained as

       e(k + 1) = (T A -LC)e(k) + T F n(k) + T Dw(k) ϵ(k) = Ce(k) (26)
In order to analyse the estimation error, the error system ( 26) is split into two parts, i.e.

       e f (k + 1) = (T A -LC)e f (k) + T F n(k) ϵ f (k) = Ce f (k) (27)
and

       e w (k + 1) = (T A -LC)e w (k) + T Dw(k) ϵ w (k) = Ce w (k) (28)
where

e(k) = e f (k) + e w (k) ( 29
)
ϵ(k) = ϵ f (k) + ϵ w (k) ( 30
)
ϵ f (k) = Ce f (k), ϵ w (k) = Ce w (k) (31)
Remark 2. The idea of separating ( 26) as ( 27) and ( 28) is inspired by [START_REF] Zhang | Adaptive observer for multiple-input-multiple-output (MIMO) linear time varying systems[END_REF], where a very similar methodology was used to analyze and design an adaptive observer. With this treatment, the output estimation error ϵ(k) is split into two components, ϵ f (k) and ϵ w (k). The first component ϵ f (k) reflects the impact of fault on the output estimation error, while the other part ϵ w (k) contains the influence of the disturbance. Note that this separation is only used to facilitate the analysis of error dynamic [START_REF] Wang | Actuator fault estimation for a class of nonlinear descriptor systems[END_REF] and the design of filter [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF]. It should be noticed that only the out estimation error ϵ(k), which is available in practice, will be used in the proposed fault estimation method.

Therefore, this separation does not involve any practical realization problem.

Then, the following Theorem is proposed to constrain the structure of matrix L in the filter [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF]. This Theorem is of essential importance in filter design.

Theorem 1. If there exists a matrix L satisfying

(T A -LC)Ψ = 0 ( 32
)
then ϵ f (k) can be expressed as

ϵ f (k) = C(T A -LC) k e(0) + Υ [n 1 (k -ρ 1 ) n 2 (k -ρ 2 ) • • • n q (k -ρ q )] T , k > k f + s ( 33
)
where Υ = CΨ, e(0) is the initial estimation error and k f denotes the time instant when faults occur.

Proof. From ( 27), it is easy to derive that

ϵ f (k) = C(T A -LC)e f (k -1) + CT F n(k -1) = C(T A -LC) 2 e f (k -2) + C(T A -LC)T F n(k -2) + CT F n(k -1) . . . = C(T A -LC) k e f (0) + CG k-k f T F n(k f ) + • • • + CG s T F n(k -s) + • • • + CG 1 T F n(k -1) (34)
where k > s + k f , matrix G j is given by

G j = (T A -LC) j-1 , j = 1, • • • , k -k f (35)
From the definitions of T and n(k), we have

CG j T F n(k -j) = CG j T f 1 n 1 (k -j) + • • • + CG j T f i n i (k -j) + • • • + CG j T f q n q (k -j) (36)
Now, we should discuss about CG j T f i . By using the definition of fault detectability indexes, we have

CG j T f i = C(T A) j-1 T f i = 0, j < ρ i ( 37
)
CG j T f i = C(T A) ρ i -1 T f i , j = ρ i ( 38
)
CG j T f i = C(T A -LC) j-ρ i (T A) ρ i -1 T f i = 0, j > ρ i (39)
Using [START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF], we obtain

(T A -LC)(T A) ρ i -1 T f i = 0 (40)
Substituting ( 40) into (39) yields

CG j T f i = 0, j > ρ i ( 41
)
To sum up, we have

CG j T f i =        C(T A) ρ i -1 T f i j = ρ i 0 j ̸ = ρ i j = 1, • • • , k -k f (42)
It follows that

CG j T f i =        C(T A) ρ i -1 T f i j = ρ i 0 j ̸ = ρ i j = 1, • • • , k -k f (43)
Substituting ( 43) into [START_REF] Shen | Actuator fault-tolerant robust control for a linear system with adaptive fault identification[END_REF] gives

ϵ f (k) = C(T A -LC) k e(0) + C(T A) ρq-1 T f q n q (k -ρ q ) + • • • + C(T A) ρ 1 -1 T f 1 n 1 (k -ρ 1 ) = C(T A -LC) k e(0) + Υ [n 1 (k -ρ 1 ) n 2 (k -ρ 2 ) • • • n q (k -ρ q )] T ( 44
)
where k > s + k f . This completes the proof.

Theorem 1 implies that the vector ϵ f (k) will exhibit the form of (33) if the matrix L satisfies the equation [START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF]. Since the output estimation error ϵ(k) is composed of ϵ f (k) and ϵ w (k), the effect of w(k) should be considered in the design of matrix L.

Based on Theorem 1, the following Theorem is proposed to design the fault estimation filter [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF]. First, the gain matrix L in ( 23) is required to have a special structure which makes the equation (32) satisfied. Moreover, the matrix L is designed so that the residual r(k) is robust against the disturbance w(k).

Theorem 2. For a given scalar γ > 0, if the following conditions are satisfied:

(i) rank(Υ) = q (ii) there exist a symmetric positive definite matrix P > 0, P ∈ R n×n and a matrix

Y ∈ R n×m such that         -P + (Υ † C) T (Υ † C) 0 A T 1 P -C T 1 Y T * -γ 2 I d (T D) T P * * -P         < 0 ( 45
)
where

A 1 = T A -T AΨΥ † C, C 1 = (I m -ΥΥ † )C (46)
Then, by choosing

W = Υ † ( 47
)
and

L = T AΨΥ † -P -1 Y (I m -ΥΥ † ) (48)
the residual r(k) in the filter ( 23) satisfies the following robust performance

∥ξ∥ ≤ √ γ 2 ∥w∥ 2 + N (e(0), k) ( 49
)
where

ξ(k) = r(k) -n d (k) (50) n d (k) [n 1 (k -ρ 1 ) n 2 (k -ρ 2 ) • • • n q (k -ρ q )] T (51)
and N (e(0), k) is a function satisfying lim k→∞ N (e(0), k) = 0.

Proof. If rank(Υ) = q, then we have

Υ † Υ = I q (52)
By letting W = Ψ † and considering [START_REF] Zhang | Fast fault estimation and accommodation for dynamical systems[END_REF], the residual r(k) is given by

r(k) = W ϵ(k) = Υ † ϵ f (k) + Υ † ϵ w (k) = Υ † C(T A -LC) k e(0) + [n 1 (k -ρ 1 ) n 2 (k -ρ 2 ) • • • n q (k -ρ q )] T + Υ † ϵ w (k) = Υ † C(T A -LC) k e(0) + n d (k) + Υ † ϵ w (k) (53)
It follows that

ξ(k) = ξ 0 (k) -ξ w (k) ( 54
)
where

ξ 0 (k) = Υ † C(T A -LC) k e(0), ξ w (k) = Υ † ϵ w (k) ( 55
)
Using triangle inequality gives

∥ξ∥ = √ ∥ξ 0 ∥ 2 + ∥ξ w ∥ 2 ≤ √ (∥Υ † C∥∥(T A -LC)∥ k ∥e(0)∥) 2 + ∥ξ w ∥ 2 (56)
Let

N (e(0), k) = ( ∥Υ † C∥∥(T A -LC)∥ k ∥e(0)∥ ) 2 (57)
It is obvious that N (e(0), k) satisfies lim k→∞ N (e(0), k) = 0 if matrix T A -LC is stable, i.e. ∥(T A -LC)∥ < 1 Moreover, the inequality (49) is fulfilled if ∥ξ w ∥ < ∥w∥ holds. Now, it is ready to show that T A -LC is stable and ∥ξ w ∥ < ∥w∥ holds if (45) is satisfied.

Applying Lemma 3 to system [START_REF] Chadli | H-/H∞ fault detection filter design for discrete-time Takagi-Sugeno fuzzy system[END_REF], it is obtained that matrix T A-LC is stable and ∥ξ w ∥ < ∥w∥ holds if the following matrix inequality is satisfied

    (T A -LC) T P (T A -LC) -P + (Υ † C) T (Υ † C) (T A -LC) T P D * -γ 2 I + (T D) T P T D     < 0 ( 58
)
Using the Schur complement Lemma [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF], the inequality (58) is equivalent to

        -P + (Υ † C) T (Υ † C) 0 (T A -LC) T P * -γ 2 I d (T D) T P * * -P         < 0 (59)
Since rank(Υ) = q, we can apply Lemma 2 to (32), then the general solution for L is obtained as

L = T AΨΥ † + R(I m -ΥΥ † ) (60)
where R ∈ R n×m is an arbitrary matrix.

Substituting (60) into T A -LC gives

T A -LC = T A -T AΨΥ † -R(I m -ΥΥ † ) = A 1 -RC 1 (61)
It follows that inequality (59) is equivalent to

        -P + (Υ † C) T (Υ † C) 0 (A 1 -RC 1 ) T P * -γ 2 I d (T D) T P * * -P         < 0 (62)
It is obvious that inequality (62) becomes (45) by letting Y = P R. Now, it has been shown that T A -LC is stable and ∥ξ w ∥ < ∥w∥ holds if (45) is satisfied.

If there exist a symmetric positive definite matrix P and a matrix Y such that (45) holds. Then R is given by

R = P -1 Y (63)
Substituting (63) into (60) yields (48). This completes the proof.

Remark 3. It is known that fault detection is an essential task in fault diagnosis. To achieve the fault detection goal, a residual evaluation procedure is necessary. A commonly used strategy is the so-called norm based residual evaluation [START_REF] Ding | Model-based fault diagnosis techniques: Design schemes, algorithms, and tools[END_REF], i.e. the norm of the residual r(k) is compared with a threshold to detect the occurrence of a fault. In this paper, the robust performance criterion (49) can be used to compute the detection threshold. Considering the fault-free situation (n(k) = 0) and using (49), it comes

∥r∥ = ∥ξ∥ ≤ √ γ 2 ∥w∥ 2 + N (e(0), k) (64)
Without loss of generality, it is assumed that the bounds of ∥w∥ and ∥e(0)∥ are w and ē0 , respectively. Then, we can obtain the following fault detection scheme

       Fault -free if ∥r∥ ̸ = θ(k) Fauly if ∥r∥ > θ(k) (65)
where

θ(k) = √ γ 2 w2 + N (ē 0 , k) (66)
is a dynamic detection threshold.

Remark 4. After the fault is detected, fault isolation and estimation should be implemented in order to obtain further information about the fault. It is noted that the proposed method is able to achieve fault isolation and estimation simultaneously. Note that the performance criterion in (49) implies that the residual r(k) robustly tends to the delayed fault magnitudes vector n d (k).

First of all, the residual r(k) can be viewed as an estimate of fault n(k). Therefore, the fault estimation goal is achieved. From the point of view of fault isolation, the ith component r i (k) of the residual can be used as a fault indicator for the ith actuator because r i (k) robustly tends to the ith fault n i (k -ρ i ) without affecting by other faults. In other words, the ith residual component r i (k) can be used to determine whether the ith actuator is fault.

Finally, for the sake of readability, the overall procedure for the proposed fault estimation filter design method is summarized in Algorithm 1.

Simulations

A simulation example is given in this section to demonstrate the design procedure and the effectiveness of the proposed fault estimation filter.

ρ i ← min { v : C(T A) v-1 T f i ̸ = 0, v = 1, 2, • • • } Calculate the fault detectability matrix: Ψ ← [ (T A) ρ 1 -1 T f 1 • • • (T A) ρq-1 T f q ] Υ ← CΨ γ ← γ 0 repeat
Solve the LMI problem (45) ◃ To design the gain matrix L if the LMI (45) has a solution

R ← P -1 Y else γ ← γ 0 + ∆γ end if until the LMI (45) is solved L ← T AΨΥ † -R(I m -ΥΥ † ) W ← Υ † ◃ To obtain the matrix W
Suppose that the matrices of system (1) are given by

E =             1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0             , A =            
0.9 0.005 -0.095 0 0.005 0.995 0.0997 0 0.095 -0.0997 0.99 0

1 0 1 1             B = F =             0.1 0 1 1 -0.1 1 -1 0             , D =             1 1 1 0             , C =         0 1 0 0 0 0 1 0 0 0 0 1        
By choosing S in [START_REF] Duan | Analysis and design of descriptor linear systems[END_REF] as

S =            
1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1

           
It comes the following matrices

T =             1 0 0 0 0 0.5 0 0 0 0 0.5 0 0 0 0 1             , N =             0 0 0 0.5 0 0 0 0.5 0 0 0 1            
Note that CT f 1 ̸ = 0 and CT f 2 ̸ = 0, the fault detectability indexes are ρ 1 = 1 and ρ 2 = 1.

Then Ψ = [T f 1 T f 2 ] =            
0.9 0 0.0025 0.5 0.0475 0.5

1 0             , Υ = CΨ =         0.0025 0.5 0.0475 0.5 1 0        
It is obvious that rank(Υ) = 2, i.e. the condition (i) in Theorem 2 is satisfied.

By the definitions of matrices A 1 and C 1 given in Theorem 2, it follows

A 1 =             0
            , C 1 =         0
L =             2
           
In the following simulation, the input signal is u(k) = [2sin(k) 2sin(k)] T , the initial value of state vector is x(0) = [0.5 1 0 -0.5] T , while the initial state estimation is x(0) = [0 1 0 -0.5] T . The disturbance w(k) in the simulation is set as a random noise drawn from the uniform distribution on the interval (0, 0.1). The reason to choose a random noise with uniform distribution is that this kind of noise exhibits both bias and stochastic properties. Herein, it is assumed that we can obtain the following priopri: w = 0.11, x0 = 0.55.

Consider the following actuator faults

n 1 (k) =        0 k < 80 1 k ≥ 80 , n 2 (k) =        0 k < 100 1.2sin(0.1k) k ≥ 100
In this situation, the fault detection result is shown in For comparative purpose, a PI observer proposed in [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF] has been designed and simulated.

Fig. 3 shows the actuator fault estimation results of the PI observer. As shown in Fig. 3, periodic vibration appears in the fault estimation of n 1 (k). This phenomenon illustrates that the fault estimation of n 1 (k) is disturbed by the dynamic of n 2 (k). In fact, the dynamics of fault estimations and state estimations vectors are coupled in the PI observer. Therefore, the fault estimation errors and state estimation errors are affected by each other and deteriorate the fault estimation performance. Moreover, it can be seen from Fig. 2 and Fig. 3 that our presented filter has faster convergence speed and better estimation accuracy than the PI observer from [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF].

Conclusion

This paper proposes an actuator fault estimation filter design method for the discrete-time linear descriptor systems. The proposed fault estimation filter is able to detect and estimate actuator faults. Moreover, the filter is designed such that the residual is robust against the unknown disturbance. The filter design conditions have been derived as an LMI feasibility problem, which can be solved efficiently.

Algorithm 1

 1 Fault estimation filter design for discrete-time descriptor systems Given the parameters of the descriptor system (1): E, A, B, F , D, C Choose a matrix S Solve the equation (7) and obtain a matrix M ◃ To design matrices T and N T ← MR T , N ← MR N Calculate the fault detectability indexes:

  By repeated applications of Theorem 2, we obtain a minimal γ = 1.5 and the corresponding L is

Fig. 1 .Figure 1 :

 11 Figure 1: The fault detection result of the proposed filter

Figure 2 :Figure 3 :

 23 Figure 2: The fault estimation result of the proposed filter

Acknowledgements

This work was financially supported by the National Natural Science Foundation of China (Grant No. 61273162). The authors would like to thank the associate editor and the anonymous reviewers for their helpful remarks, which have improved the quality of this paper.