Didier Rullière
email: didier.rulliere@univ-lyon1.fr

Stéphane Loisel
email: stephane.loisel@univ-lyon1.fr

The win-first probability under interest force

Keywords: Ruin probability, hazard rate, upper absorbing barrier, constant interest force, risk-return indicator, win-first probability. JEL Classification codes: G22, C60

In a classical risk model under constant interest force, we study the probability that the surplus of an insurance company reaches an upper barrier before a lower barrier. We define this probability as win-first probability. Borrowing ideas from life-insurance theory, hazard rates of the maximum of the surplus before ruin, regarded as a remaining future lifetime random variable, are studied, and provide an original derivation of the win-first probability. We propose an algorithm to efficiently compute this risk-return indicator and its derivatives in the general case, as well as bounds of these quantities. The efficiency of the proposed algorithm is compared with adaptations of other existing methods, and its interest is illustrated by the computation of the expected amount of dividends paid until ruin in a risk model with a dividend barrier strategy.

Introduction

In this paper, we propose a way to compute the probability that a risk process reaches an upper barrier (representing a goal or a threshold for a dividend policy) before crossing a lower barrier (representing the ruin of the company, or a threshold for insolvency penalties). We define this probability as win-first probability.

We consider the compound Poisson risk model with a constant instantaneous interest force δ. The surplus of an insurance company at time t is modeled by the process R t , where R 0 = u and R t satisfies the stochastic differential equation:

dR t = cdt -dS t + δR t dt. Here, u is the initial surplus, c the premium income rate, and the cumulated claims process S t is a compound Poisson process given by the Poisson parameter λ and the distribution function F W of the individual claim amount W , with mean m. Assume that c > λm. Denote by T u and T v u the respective times to lower or upper barrier, with initial surplus u,

T u = inf {t, R t < 0} and T v u = inf {t, R t ≥ u + v} , with T u = +∞ if ∀t ≥ 0, R t ≥ 0 and T v u = +∞ if ∀t ≥ 0, R t < u + v.
The non-ruin probability within finite time t is ϕ δ (u, t) = P (T u > t) , and the eventual non-ruin probability and ruin probability are respectively ϕ δ (u) = P (T u = +∞) and ψ δ (u) = 1 -ϕ δ (u) .

As c > λm, ct -S t a.s.

→ +∞ as t → ∞. If δ = 0 (no interest force), for any (u, v) ∈ R 2 , T v u is an almost surely finite stopping time and one can determine whether or not T u > T v u . However, if δ > 0, P(R t → +∞ as t → +∞) = 1, because there exists a threshold y < 0 such that, if for some t > 0, R t < y, then surely ∀s > t, R s < 0. This corresponds to the definition of ruin under interest force of [START_REF] Gerber | An introduction to mathematical risk theory[END_REF]. This phenomenon causes many generalizations of the classical risk model to fail. Nevertheless, if for all t ≥ 0, R t ≥ 0, then R t a.s.

→ +∞ as t → ∞. This will be very important to compute the win-first and the lose-first probabilities with constant interest force, respectively defined as

WF (u, v) = P (T v u < T u ) , LF (u, v) = 1 -WF (u, v) .
These probabilities may provide risk and profit indicators with the same unit: subjectivity is reduced to the choice of the lower bound u, which represents the event "lose", and the upper bound v, which represents the event "win". Without upper barrier, one drawback of the probability of ruin is that its minimization often prescribes the cession of the whole activity by the insurer to the reinsurer. Besides, it does not give any information about the possible profit, even for very small ruin probabilities. It is interesting to combine it with a return indicator, and one of the simplest compromises is to consider the probability WF(u, v) to reach a level u + v from initial surplus u before being ruined. It has the advantage not to require constrained optimization techniques. Risk and return indicators can be built from the win-first probability, such as the initial surplus required to avoid a failure, u ε (v) = inf {u, 1 -WF (u, v) ≤ ε}, the objective level v and confidence level ε being given, or the maximal objective level that is reasonably achievable v ε (u) = sup {v, WF (u, v) ≥ 1 -ε}, u and ε being given. The two barriers thus help to define synthetic risk-return indicators having the same unit, like (u ε (v), v) and (u, v ε (u)), useful to compare reinsurance or investment strategies. Other quantities involving win-first probabilities can be considered, such as E ((T u -T v u ) + ), E ((T v u -T u ) + ) ... Double barrier problems have been studied in the compound Poisson model without interest force by [START_REF] Segerdahl | On some distributions in time connected with the collective theory of risk[END_REF], Dickson and Gray (1984a,b), [START_REF] Wang | Some characteristics of a surplus process in the presence of an upper barrier[END_REF]. We first give properties of win-first probabilities in subsection 2.1, including a differential equation and a direct adaptation of a result of [START_REF] Segerdahl | Über einige risikotheoretische Fragestellungen[END_REF]. We thus obtain the win-first probability as a quotient of two non-ruin probabilities. A first way to tackle the problem of numerically compute win-first probabilities would be to use existing methods [START_REF] Brekelmans | Approximating the finite-time ruin probability under interest force[END_REF], [START_REF] Sundt | Ruin estimates under interest force[END_REF], 1997[START_REF] De | Numerical finite-time ruin probabilities by the Picard-Lefèvre formula[END_REF]) of computing ruin probabilities for some particular claim amount distributions, or for small δ, and to take the quotient. For exponentially distributed claim amounts, the probability of ruin under constant interest force is well-known (see [START_REF] Segerdahl | Über einige risikotheoretische Fragestellungen[END_REF], or [START_REF] Sundt | Ruin estimates under interest force[END_REF]). For general claim size distribution, bounds and Lundberg coefficients have been derived by Sundt andTeugels (1995, 1997), and several others. [START_REF] Sundt | The adjustment function in ruin estimates under interest force[END_REF] obtain bounds for the adjustment function. Konstantinides et al. (2002) obtain an asymptotical two-sided bound for heavytailed claim size distribution from generalizing results of the classical case δ = 0 to the general case. It is possible to use these bounds to get a two-sided bound for the win-first probability with interest force with heavy-tailed claim size distribution. However, we do not need in our problem to compute ruin probabilities, and we shall introduce an original method which is adapted to the present framework and more suitable in the general case and for general interest force δ than the method consisting in computing the two corresponding ruin probabilities.

The formulation of the problem, and the quotient of survival probabilities suggest the possibility to study a ruin-related survival function of some defective random variable θ, inspired from life-insurance theory. We study in subsection 2.2 its hazard rate function and propose an algorithm to compute the win-first probabilities and its derivatives, and a bound of the numerical error. A particular property of the hazard rates of θ (see theorem 5) is the key argument which makes the method so efficient. The algorithm and reasons for someone to want to use it are detailed in section 3. In section 4, numerical examples are given to demonstrate the accuracy of the algorithm and applications are proposed. In particular, computing expressions like E[WF(u -W, W )] which involve win-first probabilities are of real interest in models with dividends. For example, [START_REF] Frostig | The expected time to ruin in a risk process with constant barrier via martingales[END_REF] and [START_REF] Gerber | On the time value of ruin[END_REF], considered risk models with a dividend barrier, and computed the expected amount of dividends until time t and until ruin, or optimal dividend strategies. These quantities are expressed in subsection 4.1 in terms of win-first probabilities, which correspond in this framework to the probability that the dividends are positive. We compare our method with the one using [START_REF] Sundt | Ruin estimates under interest force[END_REF] in subsection 4.3.

Win-first probability

In this section, we first adapt classical results of ruin theory to our framework. There is no essentially new idea in subsection 2.1. This is the reason why we only state the results we shall need later. The proofs are similar as in the case δ = 0. We introduce in subsection 2.2 the new method we propose to compute the win-first probabilities in the general case.

Adaptation of classical results and methods of ruin theory

Note that WF (u, v) is nondecreasing with respect to u, nonincreasing with respect to v, and that WF (u, v) = 0 for all u < 0, and

WF (u, v) = 1 for all u ≥ 0, v ≤ 0.
Remark 1 In the special case δ = 0, R t = u + ct -S t corresponds to the classical risk process, and R t -R 0 does not depend on R 0 = u. In this case, u is not necessarily the initial reserve, and WF(u, v) corresponds to the probability that the surplus process R t reaches R 0 + v before reaching the barrier R 0 -u, and does not depend on R 0 .

Theorem 1 For v ≥ 0, w ≥ 0,

WF (u, v + w) = WF (u, v) • WF (u + v, w) .
(1)

Proof : For u ≥ 0, v > 0, w > 0, from stationarity and Markov property of R t , earning v + w before losing u may be decomposed into: earning v before losing u and then earning w before losing u + v. If v = 0 or w = 0 equality is obvious. For u < 0, both terms are equal to 0.

Theorem 2 For u ≥ 0, v > 0,

∂ ∂u WF (u, v) - ∂ ∂v WF (u, v) = λ c + δu (WF (u, v) -E [WF (u -W, v + W ))] , (2)
∂ ∂u WF (u, v) - ∂ ∂v WF (u, v) = λ c + δu WF (u, v) • (1 -E [WF (u -W, W )]) . ( 3
)
Proof : From Poisson process properties, we get

WF (u, v) = (1 -λ∆t) • WF ue δ∆t + c δ (e δ∆t -1), v -u + c δ (e δ∆t -1) +λ∆t•E WF ue δ∆t + c δ (e δ∆t -1) -W, v -u + c δ (e δ∆t -1) + W +o (∆t) .
This heuristic argument shows that equation (2) may be derived with classical ruin theory tools. For u = 0, we take the convention that ∂ ∂u WF (u, v) is the right derivative of WF (u, v). Note that in this case, the last term of equation (2) disappears. Starting from (2), a direct application of (1) leads to

WF (u -W, v + W ) = WF (u -W, W ) • WF (u, v) ,
which provides the second equation. Inequalities between win-first probabilities and some finite-time ruin probabilities may be derived.

Proposition 1 For any u ≥ 0, v ≥ 0, we have

ϕ δ (u) ≤ WF (u, v) ≤ ϕ δ (u, τ δ (u, v)) , (4)
where

τ δ (u, v) = 1 δ ln 1 + v u + c/δ if δ > 0 , and τ 0 (u, v) = v/c.
Proof : For u ≥ 0, v ≥ 0, if T u = +∞ then the insurer earns almost surely v before losing u, because R t a.s → +∞ as t → ∞. It follows WF (u, v) ≥ P (T u = +∞) = ϕ δ (u). Now, if the insurer earns v before losing u, time needed to earn v is necessarily greater than the solution τ δ (u, v) of equation in t:

ue δt + c δ e δt -1 = u + v,
and

T u > τ δ (u, v). So, WF (u, v) ≤ P[T u > τ δ (u, v)].
Finally, considering lim v→∞ WF (u, v), enables us to express WF (u, v) as a quotient of survival probabilities.

Theorem 3 For u ≥ 0, v ≥ 0, WF (u, v) = ϕ δ (u) ϕ δ (u + v) . (5)
In the special case δ = 0, this result has been recently developed by [START_REF] Wang | Some characteristics of a surplus process in the presence of an upper barrier[END_REF], and had also been treated previously by Dickson and Gray (1984b) and [START_REF] Segerdahl | On some distributions in time connected with the collective theory of risk[END_REF]. The idea is here exactly the same, and we omit the proof of the extension, which is rather direct.

From equation ( 5), it is possible to derive an exact formula for WF(u, v) in the case of exponentially distributed claim amounts (see [START_REF] Segerdahl | Über einige risikotheoretische Fragestellungen[END_REF], or [START_REF] Sundt | Ruin estimates under interest force[END_REF]), and asymptotical equivalents and bounds for general claim size distribution, as mentioned in the introduction.

Hazard rates of θ and applications

In this section we present an interesting interpretation of WF (u, v). Let us change our notation for an instant and write

v p u = WF (u, v) .
Property (1) can be written

v+w p u = v p u • w p u+v ,
and corresponds to a simple classical formula, expressed in International Actuarial Notation (see Actuarial Mathematics), stating that for a positive future lifetime θ,

P (θ ≥ u + v + w|θ ≥ u) = P (θ ≥ u + v|θ ≥ u) • P (θ ≥ u + v + w|θ ≥ u + v) .
This formula, based on elementary conditioning, illustrates the fact that someone aged u survives v + w years, if he first survives v years, and, being then aged u + v, survives w more years. So, it seems logical to look for a nonnegative random variable θ such that WF (u, v) = P (θ ≥ u + v|θ ≥ u). Let θ be the positive, defective random variable

θ = sup {R t , t ≤ T 0 | R 0 = 0} .
Define the survival function of θ by S (x) = P (θ ≥ x) , x ∈ R + , and its hazard rate by

µ x = - S (x) S(x) = - ∂ ∂x ln(S(x)).
Theorem 4 For u ≥ 0, v ≥ 0, the win-first probability can be written as

WF (u, v) = P (θ ≥ u + v|θ ≥ u) = S (u + v) /S (u) , (6) with S (x) = P (θ ≥ x) = ϕ δ (0) /ϕ δ (x) , x ≥ 0 . Proof : Let us first consider the case u = 0, v ≥ 0. If T 0 = +∞, R t a.s.
→ +∞ as t → ∞, and upper barrier v is reached after an almost surely finite time

T v u < T 0 . In this case, given that T 0 = +∞, WF (0, v) = 1 = P (θ ≥ v). If T 0 < +∞, upper barrier v is reached if and only if θ ≥ v, and WF (0, v) = P (θ ≥ v). In every case WF (0, v) = P (θ ≥ v), v ≥ 0. Consider now u ≥ 0, v ≥ 0. We have seen that T 0 = +∞ implies θ ≥ u. So P (θ ≥ u) ≥ ϕ δ (0) > 0. Starting from property (1), we have WF (u, v) = WF(0, u + v)/WF (0, u) = P (θ ≥ u + v) /P (θ ≥ u).
And the result is obvious since v ≥ 0.

For u ≥ 0, v ≥ 0, note that

µ u+v = - ∂ ∂v ln WF (u, v) .
This rate is finite and only depends on the sum u + v. In the case of integervalued claim amounts, we will see that µ u is continuous and derivable at each u ∈ R + N. For u ∈ N, µ u will be only right-continuous and rightdifferentiable, so that we will take the convention that each derivative of µ is its right derivative. We will take the same convention for derivatives in u of WF(u, v). Given that θ ≥ u, the conditional density of θ is

f θu (x) = ∂ ∂x P (θ < u + x|θ ≥ u) = WF (u, v) • µ u+v . Hence, for example, LF (u, v) = P (θ < u + v|θ ≥ u) = v 0 WF (u, s) µ u+s ds.
In the sequel, since we will use common actuarial tools, we will most often preferably write probabilities with standard actuarial notations, using t p x instead of WF(x, t), and will also write:

µ (i) u = ∂ i ∂u i µ u , S (i) u = ∂ i ∂u i S(u), t p (i) x = ∂ i ∂x i t p x , w p (i) u-w = ∂ i ∂u i w p u-w .
Note that, due to these definitions, we do not have an equality between S (i) u and w p

(i) u-w when w = u. Let us denote by C k n the binomial coefficient for integers k and n, 0 ≤ k ≤ n.
Proposition 2 For u, v ≥ 0, we have

WF (u, v) = exp - u+v u µ s ds, ( 7
) t p (1) x = t p x (µ x -µ x+t ), t p (k+1) x = k i=0 C i k t p (i) x (µ (k-i) x -µ (k-i) x+t ), k ≥ 0. ( 8
)
Proof : (7) holds directly from theorem 4. Differentiations are straightforward.

Proposition 3 A general link between unconditional survival function and hazard rate is given for x ≥ 0, k ≥ 0 by

S (1) (x) = -µ x S(x), S (k+1) (x) = - k i=0 C i k µ (i) x S (k-i) (x). ( 9
)
Theorem 5 The hazard rate of θ and its right derivatives are as follows:

µ u = α u (1 -E( W p u-W )), µ (k) u = α (k) u -k j=0 C j k α (k-j) u E( W p (j) u-W ), with α (k) u = k!λ(-δ) k (c + δu) -(k+1) , α u = α (0) u .
for u ≥ 0, k ≥ 0. ( 10)

Proof : direct from (3) and from ( 7).

Note that first equation in previous relations could also be written:

µ u = λ c + δu 1 -E 1 W ≤u exp - u u-W µ s ds , µ u = λ c + δu ((1 -F W (u)) + E [ 1 W ≤u • LF (u -W, W )]) .
In particular, suppose W is a continuous random variable. It is clear that

WF (u -W, W ) = 0 if W > u. It follows from (5) that µ 0 = λ c and that ∀u ≥ 0, 0 ≤ µ u ≤ µ 0 . Since WF (0, v) = exp -v 0 µ s ds > 0 for each v > 0, µ +∞ = lim u→+∞ µ u = 0.
Furthermore, differentiation of µ u follows immediately from ( 7) and ( 5). Hence, when W is a continuous random variable, the hazard rate µ u is a continuous, decreasing function of u, such that

µ 0 = λ c , lim u→+∞ µ u = 0, µ 0 = - λδ c 2 and lim u→+∞ µ u = 0 . ( 11
)
Remark 2 For δ = 0, differentiation of WF (u, v) makes sense, and computing µ u , u > 0 in terms of ϕ 0 (u) leads to

µ u = ϕ 0 (u) /ϕ 0 (u) .
We also check that, in the special case δ = 0, formula (5) is a version of the classical risk theory formula

ϕ 0 (u) = λ c ϕ 0 (u) - λ c E [ϕ 0 (u -W )] .

Algorithm

The recursive determination of hazard rate µ u and its derivatives, for successive values of u, gives a set of values of S(u) and its derivatives up to a given order. Despite the purpose is here to find values of win-first probabilities, this will eventually give results on ψ δ (u) when ψ δ (0) is known.

The proposed iterative algorithm allows to re-use previous computed quantities, and then reduce the complexity of the determination of the whole function S. It gives many derivatives of µ u and S(u), and all numerical errors will be bounded in a further section. No assumption is made neither on claim amounts nor on the interest force δ, making the context different from studies using small δ (see [START_REF] Sundt | Ruin estimates under interest force[END_REF] and section 4.3), and from the one using particular distributions for claim amounts (see Konstantinides et al. (2002) and [START_REF] Brekelmans | Approximating the finite-time ruin probability under interest force[END_REF])

Approximations

In the sequel, W is assumed to be a random variable taking values in the set N * of positive integers. Define

π i = P (W = i), i ∈ N.
This hypothesis is not so stringent: in practice, we may approach any continuous random variable by a discrete one, and the discretization step may be chosen as small as necessary. Instead of taking this step smaller than 1, we choose this step equal to 1 and change the monetary unit.

The restriction π 0 = 0 can be easily eliminated: if π 0 > 0, one may replace π 0 , π 1 , π 2 ... with 0, [START_REF] De | Numerical finite-time ruin probabilities by the Picard-Lefèvre formula[END_REF].

π 1 / (1 -π 0 ) , π 2 / (1 -π 0 ) ... and λ with λ (1 -π 0 ) (see
The main assumption we shall use for approximations is :

Assumption H ε r : µ is locally polynomial of order r on intervals [kε, kε + ε[, k ∈ N.
We shall see further that even a choice like r = 2 and ε = 0.5 gives numerically quite good results (see section 4.2), and the precision of the algorithm increases rapidly for a better choice of these two parameters. In section 3.3, we will derive bounds for each approximated quantity in the algorithm, in order to ensure the numerical validity of this assumption. Under H ε r , we get from the r first derivatives of µ:

∀s ∈ [0, ε[, µ(x + s) = r i=0 µ (i) (x) i! s i , For x ∈ εN, S(x + ε) = S(x) exp - ε 0 r i=0 µ (i) (x) i! s i ds , S(x + ε) = S(x) exp - r i=0 µ (i) (x) (i + 1)! ε i+1 . ( 12
)
We can also derive S(x + ε) from derivatives of S(x), but we choose to use the single hypothesis H ε r . In life insurance, the hypothesis of constant hazard rate is often considered for survival lifetimes, and corresponds here to the order r = 0. In practice, it is possible to get higher order derivatives of µ u since computation of µ (r+1) u may be replaced with an approximation like µ left u

(r) = µ (r) u -µ (r)
u-ε /ε. However, one should keep in mind that, if W takes values in N, each x ∈ N is a point of discontinuity for function µ (see figure 4). So, using this approximation will give good results, except for u ∈ N. Nevertheless, since numerical results are fine enough, and since the parameter r could be chosen, we did not use this approximation.

Proposition 4 (approximation algorithm) Under hypothesis H ε r , the following algorithm computes recursively the values of S(u), µ u , E[ W p u-W ] and all their derivatives up to a given order r. With S(0) = 1, and for

u ∈ εN, u ≤ u max , k ∈ N, k ≤ r, w p (k) u-w = 1 {k=0} S(u) S(u -w) + 1 {k≥1} k-1 i=0 C i k-1 w p (i) u-w (µ (k-1-i) u-w -µ (k-1-i) u ), w = 1..[u], µ (k) u = α (k) u - k j=0 C j k α (k-j) u E( W p (j) u-W ), S(u + ε) = S(u) exp - r i=0 µ (i) (u) (i + 1)! ε i+1 .
From the second equation, quantities µ

(0) 0 = α (0) 0 and µ (k)
0 are given by the recursion, as the E( W p (j) u-W ), j ≤ k are derived from previously computed quantities w p k+1) , and that for w > u, w p (j) u-w = 0. The previous algorithm gives derivatives of µ from order 0 to r. It also gives for each u ≤ u max , u ∈ εN, S(u) and eventually ψ δ (u) = 1 -ϕ δ (0)/S(u). To obtain derivatives of WF and S, we can use (8) for ∂ k ∂u k WF(u, v) and following relations:

(j) u-w . Note that α (k) u = k!λ(-δ) k (c + δu) -(
∂ k ∂v k WF(u, v) = S (k) (u + v)/S(u), S (k+1) (x) = - k i=0 C i k µ (i) x S (k-i) (x) , k ≥ 0. ( 13
)
Let i be a positive integer. We have seen that S (i) = ϕ δ (0) /ϕ δ (i). In the special case δ = 0, it is known that ϕ 0 (i) can be exactly computable by classical formulae (see [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] and De Vylder (1999)):

ϕ 0 (i) = 1 - λm c i j=0 h j (j -i) , with h j (τ ) = λτ cj j k=1 kπ k h j-k (τ ) and h j (0) = e -λτ c .
This formula has the advantage to give exact values if W is integer-valued. Let us compare the number of loops involved in the computation of S (i) , i = 1...x by algorithm (4) and the number of loops involved in the computation of ϕ 0 (i) , i = 1...x by the Picard and Lefèvre (1997) algorithm. Computing S (i) , i = 1...x implies r 2 /2 loops for i = 1..x/ε, j = 1...iε, so that complexity of algorithm 4 is quite proportional to r 2 x 2 /ε. Computing ϕ δ (i) , i = 1...x requires loops for i = 1...x, j = 0...i, k = 1...j, so that complexity of Picard-Lefèvre formula is quite proportional to x 3 . To approximate a continuous distribution W by W d taking values in dN, time needed by algorithm ( 4) is proportional to r 2 x 2 / (d 2 ε) against x 3 /d 3 for the [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] formula. Noting that hypothesis H ε 1 means that µ is linear on intervals of length dε, one can use ε = 1 and r = 1 if d is small. As both formulae lead to an approximation of values obtained for a continuous W , the algorithm may be of a practical interest even in the case δ = 0. Moreover, we will see in section 3.4 that the complexity of the algorithm can be reduced in this case.

Convergence for parameters r and n

The highest order of derivatives that are computed by the algorithm is r, and n = 1/ε is an integer that represents the number of sub-periods in one unit of time. The hypothesis in approximation algorithm is that, on each sub-period, µ u is locally polynomial of order r. The precision of the algorithm, at one step, is given by η, which represents the number of decimal digits that one aims at obtaining. More precisely, 10 -η represents the error in the approximation of µ u+ε by the Taylor expansion of order r. To improve the local precision of the algorithm, we can increase either n or r; this may have different effects on the complexity of the algorithm. We only give here informal considerations for the choice of the couple (n, r) to minimize the complexity of the algorithm. It would be possible to get more rigorous results for that choice of parameters, but they are omitted here in the interest of conciseness. Note first that the remaining part in the Taylor expansion behaves like µ (r+1) (r+1)! ε r+1 . To simplify further calculation, take u = 0, since µ (r+1) (0) is known, equal to α (r+1) (0). In absolute value, the error is then comparable to λ( δ cn ) r+1 . If this last quantity is set to be equal to 10 -η , then a link appears between r and n:

r = η ln(10) + ln(λ) ln(cn/δ) -1, n = δ c (λ10 η ) 1/(r+1) .
For a given u the local complexity of the algorithm is then proportional to

c(r) = nr 2 = δ c λ10 η 1 r+1 r 2 .
Trying to find r 0 that minimizes c(r), we find, in the case where ln(λ10 η ) ≥ 8

r 0 = 1 4 (ln(λ10 η ) -4 + ln(λ10 η )(ln(λ10 η ) -8) , n 0 = δ c λ10 η 1 r 0 +1 .
Since n 0 is here a real number, and should better be an integer greater than 1, and since r must also be an integer, we may choose the following parameters to ensure that required precision on µ is reached at the first point following u = 0:

n opt = max(2, [n 0 ] + (0 or 1)
) and r opt = η ln(10) + ln(λ) ln(cn opt /δ) + (0 or 1).

As an example, take δ = 100%, so that we do not suppose that δ is close to 0. For λ = 1, c = 1 and η = 12 decimal digits, we get n opt = 9 and r opt = 12. With η = 16 decimal digits, we get n opt = 9 and r opt = 16, so that the complexity is multiplied by something less than 1.8 to reach 4 more decimal digits.

Bounds for µ u , WF(u, v) and their derivatives

The algorithm makes only one approximation by replacing µ u+ε with its Taylor expansion. Nevertheless, this approximation is used recursively, so that even if the error is locally bounded, we cannot ensure that the global result will be precise enough. For this reason, we must give exact bounds for the values we approximate.

For a function f u of u, we will use the following notations:

f u [-1] and f u [+1] will be bounds of f u , such that f u ∈ f u [-1] , f u [+1]
. We define by this way µ (k)

u [σ] ,
S(u) [σ] and

β (k) u [σ] , with β (k) u = E W p (k)
u-W , for σ ∈ {-1, +1}. For two bounded quantities a and b, we will use following arithmetic, that might be simplified when the signs of a and b are known.

(a + b) [σ] = a [σ] + b [σ] , (a -b) [σ] = a [σ] -b [-σ] , (ab) [σ] = max σ 1 ,σ 2 ∈{-σ,σ} {a [σ 1 ] b [σ 2 ] }, σ ≥ 0 (ab) [-σ] = min σ 1 ,σ 2 ∈{-σ,σ} {a [σ 1 ] b [σ 2 ] }, σ ≥ 0.
Note first that, when u = 0, S(0) [+1] = S(0) [-1] = 1. From ( 8), we can bound w p (k) u-w from bounds of µ (j) u and w p (j) u-w , j < k, w ≤ u. From (10), we can also S(u) [-1] and S(u) [+1] give bounds for S(u), not for its approximation S(u) [0] , which can be outside the interval.

Proposition 5 (bounding algorithm) Bounds for µ u , E[ W p u-W ], S(u) and their derivatives up to order r are given by following algorithm, with initialization values S(0) [-1] = S(0) [+1] = 1. For u = 0..u max by step ε, for k = 0..k + 1, and for σ 0 = +1, +2, S(u) [+2] = S(u) [+1] , S(u) [-2] = 0

w p (0) u-w [σ] = S(u) [σ] S(u -w) [-σ] , σ = ±σ 0 w p (k) u-w [σ] = k-1 j=0 C j k-1 w p (k-1-j) u-w (µ (j) u-w -µ (j) u ) [σ] , u ≥ 1, k ≥ 1, σ = ±σ 0 µ (k) u [σ] = α (k) u [σ] -1 {u≥1} k j=0 C j k α (k-j) u E[ W p (j) u-W ] [-σ] , σ = ±σ 0 S(u + ε) [σ] = S(u) [σ] exp - r k=0 µ (k) u [-σ] (k + 1)! ε k+1 - µ (r+1) u [-2σ] (r + 2)! ε r+2 , σ = ±1.
This algorithm is quite similar to the first one we proposed. Some remarks can be done for its practical implementation. First, we had better use only integer arguments, so that for n = 1/ε,n ∈ N, we preferably replace u with an index i = 0..nu max , where i denotes nu.

Second, for each value of u, we do not use previous values w p (k)

u 0 -w and E[ W p (k)
u 0 -W ], u 0 < u. In the algorithm, these quantities do not need to depend on u, and that spares stocking memory. Third, many quantities, like E[ W p (k) u-W ] or like Taylor integrated approximation in the exponential, can be computed in previous sums giving respectively w p (k) u-w and µ (k) u . We may check at each step if the precision of the computer is high enough. If not, it is possible to change lower and upper bounds in order to include, at each step, the maximum numerical computer error. At last, bounds for derivatives of S and WF(u, v) with respect to u and v are given for σ ∈ {-2, -1, 1, 2} by:

S (k+1)[σ] (x) = - k i=0 C i k µ (i) x S (k-i) (x) [σ] . ∂ k ∂v k v p u [σ] = S (k) (u + v) [σ] /S(u) [-σ] , k ≥ 0, ∂ k ∂u k v p u [σ] = k-1 i=0 C i k-1 v p (k-1-i) u (µ (i) u -µ (i)
u+v ) [σ] , k ≥ 1.

Further results and improved algorithm

We have seen that, given the survival function S(x) for x ∈ [0, u], it is possible to deduce exactly as many derivatives of µ u and w p u-w as wanted, and to get then an approximation of S(u + ε). The previous algorithm was constructed on this idea. For k varying from 0 to a given derivative order r, let us recall here equations that are used in this exact differentiation step (w ≤ u):

w p (k) u-w = 1 {k=0} S(u) S(u -w) + 1 {k≥1} k-1 i=0 C i k-1 w p (i) u-w (µ (k-1-i) u-w -µ (k-1-i) u
),( 14)

µ (k) u = α (k) u - k j=0 C j k α (k-j) u E( W p (j)
u-W ), ( 15)

α (k) u = k! λ c + δu -δ c + δu k . ( 16
)
This step was of complexity proportional to ur 2 . We will see here that it is sometimes possible to reduce this complexity to something proportional to ln(u)r 2 . To do so, we shall denote by Ω ν a random variable distributed as

W * 2 ν = W 1 + ... + W 2 ν , with Ω 0 = W .
The law of Ω ν can be easily constructed for integer claim amount W , since for k ∈ N,

P[Ω 0 = k] = P[W = k], P[Ω ν+1 = k] = k i=0 P[Ω ν = i]P[Ω ν = k -i], k ∈ N, ν ≥ 0.
Remark also that if S is given on [0, u], we can easily deduce Ων p u-Ων from S.

We will see that since W ≥ 1 we will only need law of Ω ν when Ω ν ≤ u, i.e. for ν ≤ ln(u) ln(2) . We previously gave derivatives for almost all relations, except an important one:

Proposition 6 By derivation of actuarial property of win-first probabilities,

s+t p x = s p x • t p x+s , s+t p (k) x = k j=0 C j k s p (j) x t p (k-j) x+s . ( 17
)
Consider first the case δ = 0. In this case, α (k) u = λ c 1 {k=0} . As w p u-w = 0 when w > u, injecting (15) into ( 14) gives

W 1 p (k+1) u-W 1 = λ c k i=0 C i k W 1 p (i) u-W 1 E W 2 [ W 2 p (k-i) u-W 2 ]- λ c k i=0 E W 2 C i k W 2 p (k-i) u-W 1 -W 2 W 1 p (i) u-W 1 .
Using proposition 6, we get then the following theorem, reducing the complexity of differentiation step to something proportional to ln(u)r 2 :

Theorem 6 When δ = 0, and if S(x) is given on x ∈ [0, u], then all derivatives of Ων p u-Ων are given by the following recursion: for ν from [ln u/ ln 2] down to 0, for k from 0 to r,

β (k+1) u,ν = 1 {2 ν ≤u} - λ c β (k) u,ν+1 + λ c k i=0 C i k β (i) u,ν β (k-i) u,ν , k ≥ 0, δ = 0 with β (k) u,ν = E Ων p (k) u-Ων .
Consider now the case δ > 0. In this case, note that

C j k α (j) u α (k-j) u = α (0) u α (k) u , with α (0) u = λ c + δu . ( 18
)
Assume that u is given and define μ(k)

x = µ (k) x /α (k) u , t px (k) = t p (k) x /α (k) u and α(k) x = α (k)
x /α (k) u . Substituting ( 18) in ( 14) and ( 15) yields:

w p(k+1) x-w = γ k k i=0 w p(i) x-w μ(k-i) x-w -μ(k-i) x , with γ k = -λ (k + 1)δ , ( 19
) μ(k) x = α(k) x -α (0) u k j=0 α(k-j) x E( W p(j) x-W ), (20)
and from 17,

s+t p(k) x = λ c + δu k j=0 s p(j) x t p(k-j) x+s . ( 21
)
Equations ( 19) to ( 21) could be useful for computations and for further analysis, since they avoid to compute binomial terms in the recurrence or factorials in the Taylor's expansion:

r k=0 µ (k) u (k + 1)! ε k+1 = λε c + δu r k=0 μ(k) u k + 1 -δε c + δu k .
This improvement being quite simple, the resulting algorithms for approximation and bounds are omitted here. Other extensions may be found for δ > 0 by similar arguments as in the case δ = 0. Injecting ( 20) into (19), and using the actuarial property, we can get an expression depending on quantities α(k) u-w , w ≤ u. Since these quantities are bounded, with α(k)

u-w ∈ [1, (1 + δ u c ) k+1 ], we can derive recursive bounds for β(k) u,ν = E[ Ων p(k) u-Ων ] as a function of β(k) u,ν+1
, for k ≥ 0 and δ > 0. We can thus construct bounds for βu = β

(k) u,0 , μ(k)
u and µ (k) u . The complexity of the differentiation step is then proportional to ln(u)r 2 instead of ur 2 , but the obtained bounds are less precise than in previous bounding algorithm. Nevertheless, this approach might be useful when looking for analytic bounds of β(k) u,ν .

4 Applications and numerical results

An example of application : payment of dividends

Let us now modify our process R t with an horizontal dividend barrier strategy. Starting from u, if the surplus reaches the upper barrier u + v, all the premium income and the interests (at rate δ) are paid as dividends until the next claim, i.e. during an exponentially distributed time ξ, with parameter λ.

We shall show here that it is possible to determine the total amount of dividends that will be paid until the process reaches the lower barrier 0, and that this cumulative amount of dividends depends on win-first probabilities and on quantities which are computed in the previous algorithm. The total expected amount of dividends is given here as a simple example, and depending on the purpose of the study, one may introduce either a discounting factor or other parameters. We shall keep in mind that the total dividend amount might be here represented by a defective random variable. Denote by D i the cumulative amount of dividends that is paid during the i th period of payment, distributed as D = ξ 0 e δs ds, where is N the number of payment periods, and T is the total amount of dividends T = N i=0 D i . We also use N 0 and T 0 , the random variables distributed as N and T given that N > 0. For any random variable X, we will denote respectively by F X and f X its distribution and density function. From the memoryless property of the modified risk process, N 0 -1 is a geometric random variable with parameter

β u+v = E[WF(u + v -W, W )],
and it is easy to get the classical results :

P[N = 0] = 1 -WF(u, v) and P[N = k] = WF(u, v)β k-1 u+v (1 -β u+v ), k ≥ 1,
and if λ > δ and β u+v < 1, supposing that u and v are fixed, with the notation ω = WF(u, v) and β = β u+v , the distribution function and the mean of the cumulated amount of dividends are as follows :

P[T ≤ x] = (1 -ω) + ω(1 -β) ∞ k=1 β k-1 F D * k (x) , E[T ] = WF(u, v) (1 -β u+v )(λ -δ) .
The expected value only depends on λ, δ, WF(u, v) and β u+v , which we are able to compute with as much precision as necessary.

Proposition 7 If there exists R such that E[e RD ] = 1/β, then we can get, by application of Smith's theorem,

lim x→∞ P[T ≤ x] = (1 -ω) + ω lim x→∞ e -Rx ∞ 0 ā(y)dy ∞ 0 (1 -Ḡ(y))dy , with ā(y) = (1 -β)e Ry F D (y), d Ḡ(y) = βe Ry f D (y), F D (y) = 1 -(1 + δy) -λ δ , and f D (y) = λ(1 + δy) (-λ δ -1) .
This simple example shows how quantities β = E[WF(u + v -W, W )] and ω = WF(u, v) naturally appear in the computation of the expected dividends.

Numerical results

The results presented hereafter have been obtained for λ = 1 and c = 1.05. W is first exponentially distributed with parameter 1, and then discretized with F d (id) defined on each interval [id, id + d[, such that:

F d (id) = 1 d [id,id+d[ F W (x) dx.
We have taken d = 1. As explained in section 3.1, in order to cancel π 0 = P(W = 0), the Poisson parameter λ has been modified into λ(1 -π 0 ), and the π i have been changed too. This discretization procedure if fully described [START_REF] De | Numerical finite-time ruin probabilities by the Picard-Lefèvre formula[END_REF]. This explains values for x = 0 in figures 2 and 3. The observation of the evolution of µ x , x > 0 for integer-valued claim amounts confirm that it is nonincreasing, but not continuous. This is a classical fact in ruin theory, and it explains that we usually observe discontinuity points that really exist, even if exact computations are carried out. Let us explain, for example, that if δ = 0,

µ x = µ 0 for each x < 1. Note that θ = sup {R t , t ≤ T 0 | R 0 = 0}.
Starting from 0, the random variable θ keeps growing, as a survival lifetime, until the first claim. If the claim occurs before θ reaches the value 1, then ruin occurs since the claim amount is a positive integer. As long as θ < 1, for δ = 0, the probability that θ stops growing is directly linked with the hazard rate of the time of the first claim, which is constant and equal to the modified λ. If θ has reached x > 1, situation is more complex, since the probability that θ stops growing will also depend on the claim amount. The analysis of derivatives of hazard rates (see figures 4 and 5) may be important to understand approximations that are made in the proposed algorithm. Replacing µ (k+1) iε with the approximation µ left iε

(k) = µ (k) iε -µ (k)
(i-1)ε /ε gives good results for iε / ∈ N, but it must be done keeping in mind the discontinuity of µ x and of its derivatives on atoms of the distribution of the claim amount (see tables 1 and 2). Despite discontinuities of hazard rates of θ (see figures 2 and 3), survival function S (x) is continuous (see figure 6), and tends to ϕ δ (0) as x → +∞. This function is sufficient to obtain all values of WF (u, v) = S (u + v) /S (u) , u, v > 0. Of course, in the special case δ = 0, computation of probabilities of ruin and non-ruin are already well-known, and may be computed for example with classical formulae (see [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] or [START_REF] Rullière | Another look at the picard-lefèvre formula for finite-time ruin probabilities[END_REF]). We retrieve S (u) by computing the ratio ϕ 0 (0) /ϕ 0 (u) (see table 3). We shall remember that, for u > 0, although the computation of ϕ 0 (u) is exact, it does not use previous computations of ϕ 0 (x) , x = 1, ..u -1. This implies, especially if discretization of W is really accurate, a computation time that could be important. It is thus interesting to propose another way to determine WF (u, v) that would help to understand the structure of θ. In table 3, we see that approximation algorithm gives quite precise results for small values of convergence parameters n and r, and that precision increases rapidly when r becomes larger.

To give an idea of the convergence of the bounding algorithm, we have taken convergence parameters n = 2 and r = 100. Keeping c = 1.05, λ = 1, we obtain quantities in tables 4, 5, 6, 7 and 8, in both cases δ = 0.05 or δ = 1.2. Rather than proposing very near bounds for each quantity, we preferred showing only decimals that were in common in lower and upper bounds. The great number of correct digits shows that the algorithm gives very thin bounds when r becomes large. It may help measuring quality of analytical approximations, and also helps comparing precision of the algorithm with the existing one in the literature. This last point will be developed in section 4.3. In table 6, we gives bounds for the 10 first order derivatives of S and µ. When r is large enough, bounds remain very thin also for these derivatives, and are far much precise than the one that could be obtained by successive finite differences on Table 6 Exact decimal digits of right derivatives S(u) (k) and µ thin intervals of length ε. This result comes directly from the fact that (5) gives at point x exact values of these derivatives, when S is given on [0, x].

Only approximations on S, that are numerically very precise, have an impact on these derivatives. At last, we have bounded, for u = 5 and v = 4, terms appearing in differential equation (3):

∂ ∂u WF (u, v) - ∂ ∂v WF (u, v) = λ c + δu WF (u, v) • (1 -E [WF (u -W, W )]) .
Only decimals that are in common in lower and upper bounds are written in tables 7 and 8. As λ was modified to eliminate the mass P [W = 0], we easily verify this differential equation, in both cases δ = 0.05 and δ = 1.2. Convergence parameters n = 2 and r = 100 give in both cases, for the equality, a better precision than the 120 decimal digits we used for calculations. An interesting result of the algorithm is that it also gives all derivatives up to a given order, with respect to u or in v, of WF(u, v). To give a concise illustration of section 4.1, figure 7 draws the evolution of average cumulative dividends that may be paid each time the process reaches the upper barrier without having reached the lower one. This simple, natural example is based on quantities computed in approximation or bounding algorithm. It is given here in a simplified environment, and introduction of other economical parameters, such as a discounting factor, would require further analysis. [START_REF] Sundt | Ruin estimates under interest force[END_REF] proposed several methods to compute ψ δ (u). Each one is based on the value of ψ δ (0). If these methods are used to compute S(u) = WF(0, u) = 1-ψ δ (0) 1-ψ δ (u) , then the result obtained depends on the value of ψ δ (0). [START_REF] Sundt | Ruin estimates under interest force[END_REF] proposed for example a recursive algorithm that we rewrite with our notations:

Comparison with other methods

Sundt

ϕ δ (hk) [-1] = γ -   cϕ δ (0) + k j=1 ϕ δ (h(k -j)) [-1] f + j   , ϕ δ (hk) [+1] = γ +   cϕ δ (0) + k-1 j=1 ϕ δ (h(k -j)) [+1] f - j   .

Figure 1 .

 1 Figure 1. Sample path of R t , with θ and θ u = (θ -u) + (δ = 20%).

Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 2345 Figure 2. Aspect of µ for integer-valued W and δ = 0.

Figure 6 .

 6 Figure 6. Aspect of survival function S(u) = WF(0, u) for δ = 0 and δ = 0.05.

  δ = 0.05). quantity value for u = 5, v = 4, δ = 0.05 WF(u, v) .8230914532618470298053719011486982427784142006271074... ∂ ∂u WF(u, v) .0550920169557785632478932959551552772017674284516602... ∂ ∂v WF(u, v) -.01814985623171288470150842179406380159358909721724867... E[WF(u -W, W )] .816998441718484612223534760053005148559101965324086718334... λ/(c + δu) .486246583714275137234212484491183948118606822283255511917...

Figure 7 .

 7 Figure 7. Average cumulative dividends as a function of premium rate c (δ = 0.05).

Table 2

 2 Some values of derivatives of µ for δ = 0.05 and ε = 0.01.

	1
	0.9
	0.8
	0.7
	0.6
	0.5
	0.4
	0.3
												S with δ=5%
	0.2
	0.1											S with δ=0%
	0
	0	1	2	3	4	5	6	7	8	9	1 0	1 1	1 2

Table 5 Exact

 5

	k	S(u) (k)	(k) µ u
	0	.244757282698191919476064786081132301397790822287	.088983882528775672523830328545016100434528133542
	1	-.021779453291678248139853343964504022264773000051	-.027629483673356369324950664995212902914462559697
	2	.008700537659492416889348321750556619101482134957	.004781759780650901141178540920815767739921848509
	3	-.003148088249731182436118591683531570999288924430	.001184056174728555633066782660597992807695785290
	4	.001023929174136961425000957228656082809071769453	.000079184281921029225092433587328484379390687953
	5	-.000604885236309820279361901997938225188501618568	-.000402623545415252800623317732679845596127644524
	6	.000349960969975695293643389220473108020035197533	-.000386370000451382731645846489531581540120338400
	7	-.000098690210126223893537006123193270538997311567	-.000179180225484452079551481731074750159874505135
	8	.000162013429988748369733435193279071491692638059	.000076722640836791688533098696105266989982210445
	9	-.000075506148302808719759688754740222191140782554	.000264991262234335480487532290438883180373071878
	10 -.000008701456970729554583982507829971722968348809	.000285704192102520471952411764658858649021463062

decimal digits for S(u) and µ u by bounding algorithm (δ = 1.2).

Table 7

 7 Exact decimal digits of quantities in differential equation by bounding algorithm (δ = 0.05).

Table 8

 8 Exact decimal digits of quantities in differential equation by bounding algorithm (δ = 1.2).

bound µ (k) u from bounds of w p (j) u-w , j ≤ k. We get then w p (k) u-w [σ] and µ (k) u [σ] for σ ∈ {-1, +1}.

We will now use for a function f u of u the following notations: f u [-2] and f u [+2] will be bounds of f u , such that ∀s < ε, f u+s ∈ f u [-2] , f u [+2] .

Note that, when u = 0, S(0) [+2] = 1, and since S is decreasing, S(0) [-2] ≤ S(ε) [-1] . The sign of α (k) u is the same as the one of (-1) k . Since α (k) u is thus either increasing or decreasing in u, depending on k mod 2, we remark that

u+ε .

Using such α (k) u

[σ] , we can easily check that ( 8) and ( 10) can be adapted to get bounds w p

The knowledge of bounds of µ (k) u+s , s < ε will allow us to derive bounds of the derivative form of Taylor's remainder, and then bounds of S(u). For s ∈ [0, ε[, µ u+s is a continuous and r + 1 times differentiable function of s. We have

is bounded, we can bound R (r) u,s , and then S(u + ε):

S(u+ε) [σ] = S(u) [σ] exp -

The only difficulty to build the bounding algorithm is the following: since S(u) is decreasing in u, a good lower bound for S(x), x ∈ [u, u + ε[ is given by the lowest value of S(u + ε), so that we can propose S(u) [-2] = S(u + ε) [-1] . Nevertheless, the calculation of S(u + ε) [-1] from S(u) [-1] uses µ (r+1)

, that is then calculated from S(u) [-2] . Using such a bound gives then S(u) [-2] as a computable function of itself. We have built both a formal computation algorithm, in order to get the root value of S(u) [-2] , and also a fixed-point algorithm, starting from S(u) [-2] = 0. Nevertheless, since the last term of Taylor expansion becomes very small for large values of r, such precise bounds of S(u) [-2] could be replaced with S(u) [-2] = 0. The great acceleration resulting of this choice can be exploited to increase r or n = 1/ε, for example, and thus the precision of the algorithm. We will see with numerical figures that this approximation is sufficient to get very precise results. Indeed, it only changes bounds for the r + 1 th derivative order of µ u and has an impact comparable to µ (r+1) u [+2] ε r+2 /(r + 2)!. The problem does not hold for S(u) [+2] since the better bound we can propose is S(u) [+2] = S(u) [+1] . Note that, by construction,

, and, in the special case of integer-valued claim amounts,

Note that the corresponding formulae for this quantities in [START_REF] Sundt | Ruin estimates under interest force[END_REF] (top of page 12) have to be switched. Let us try to minimize

Hence, after some omitted computations, with u = hk,

As an example, in the case δ = 0, we can get from table 3 values for the right member of (22). With same numerical parameters, the minoration of ∆ ϕ h k changes from values 10 -3 to 10 -5 when hk ∈ [1, 10]. To get the same precision level 10 -η as in table 4, would require an h smaller than 10 -η+5 , and a much higher complexity in 1/h 2 than with our method. One must add to this problem the error possibly made in ϕ δ (0), which was supposed to be avoided, and the propagation error due to ∆ ϕ h k-j , j ≤ k. Other methods might not be more efficient, except in the case where ϕ δ (0) is precisely bounded, for as well small or large values of δ. Besides, the adjustment functions do not in general provide directly two-sided bounds for ruin probabilities. They are particularly adapted to the case of large initial reserve, which does not correspond to the assumption made here. This shows that the method consisting in taking the quotient of two non-ruin probabilities (computed with methods efficient for that problem) is not adapted to our framework, and that the algorithm proposed in section 3 is more convenient to this problem.


Nicolas Simond
email: nicolas.simond@inria.fr

Obstacle Detection from IPM and Super-Homography

We present in this article a simple method to estimate an IPM view from an embedded camera. The method is based on the tracking of the road markers assuming that the road is locally planar. Our aim is the development of a freespace estimator which can be implemented in an Autonomous Guided Vehicle to allow a safe path planning. Opposite to most of the obstacle detection methods which make assumptions on the shape or height of the obstacles, all the scene elements above the road plane (particularly kerbs and poles) have to be detected as obstacles. Combined with the IPM tranformation, this obstacle detection stage can be viewed as the first stage of a free-space estimator dedicated to AGV in the complex urban environments.

I. INTRODUCTION

A. Related works

For two years now, most of the automotive manufacturers propose on their new cars Advanced Driving Assistance Systems that allow active security like Automatic Cruise Control, Lane Crossing Detection, obstacle avoidance, etc. The obstacle detection systems can be divided into different groups according to the types of obstacle the system detects, the ranges, the refresh rate, the reliability.

The range-finders based on radars, Laser or Lidar provide high spatial resolutions data at high scanning speeds for long ranges but remain too much expensive for automotive applications. Furthermore, the interpretation of the beam impacts on the environment is difficult in dynamic and complex environments. Cheaper sensors like sonars are nevertheless used for difficult maneuver as parking but are useless at higher speed. The drawbacks of range-finder methods are improved as soon as they are coupled with a vision system to perform a segmentation of the scene.

The vision sensor appears particularly well adapted to segment complex scenes thanks to the rich information it captures in a single frame. Several methods have been implementing for years now: they essentially differ according to the applications. From highways to urban environment, many factors vary like the density of traffic, the field of view, the types of obstacles (vehicles, trucks, cycles, pedestrians, others), the numbers of road lanes, the presence or not of road markers. In [START_REF] Sun | On-road vehicle detection: a review[END_REF], where an exhaustive recent review of the vehicle detection is presented, all the methods are visionbased except those that are coupled with a range finder.

Methods using a mono camera are generally based on apriory knowledge of the scene and particularly on the type of obstacles. These methods are dedicated to detect either vehicles according their vertical edges, shadows, symmetry, lights [START_REF] Steux | FADE: a vehicle detection and tracking system featuring monocular vision and radar data fusion[END_REF] or pedestrians and crowd according to human shape and behaviour constraints. Some original methods track sets of features points ( [START_REF] Rabie | Active-Vision based traffic surveillance and control[END_REF]) or edges ([4]) on obstacles because their relative motion is different from the camera motion. In a general way, the segmentation of the elements of a scene can be performed with a dense estimation of the optical flow [START_REF] Stein | A robust method for computing vehicle ego-motion[END_REF] and the introduction of colour, texture constraints ( [START_REF] Heisele | Motion-based object detection and tracking in color images sequences[END_REF]) or a parametrization of the static world [START_REF] Braillon | Real-time moving obstacle detection using optical flow models[END_REF].

An other technique common to mono and stereo-vision consist on considering all the elements which are lying above the road as obstacles. In a first step, the projection of the road surface is segmented according to a-priory knowledges (road model, planarity assumption). The estimation of camera motion between two poses (relative motion of the static road) in a second step, allows, in a third step, the detection of obstacles according to their elevations which don't verify the road surface motion (optical flow) or they are distorted (homography).

The stereo-vision allows by triangulation the estimation of the depth of homologous pixels matched in a pair of images if the stereo-rig is calibrated. When the images are rectified, the computation of the V-disparity, an original method developed by Labayrade in [START_REF] Labayrade | Real time obstacle detection in stereo-vision on non flat road geometry through "v-disparity" representation[END_REF], allows with real-time constraints the segmentation of obstacles whatever the longitudinal profile of the road and the discrimination between vehicles and trucks. Many authors perform the method with introducing U-disparity [START_REF] Hu | A Complete U-V-Disparity Study for Stereovision Based 3D Driving Environment Analysis[END_REF], infrared sensors to detect pedestrians [START_REF] Bertozzi | Stereo Vision-based approaches for Pedestrian Detection[END_REF] or multi-baseline sensors for off-road environment [START_REF] Broggi | A Single Frame Stereo Vision System for Reliable Obstacle Detection during Darpa Grand Challenge 2005[END_REF].

B. Motivations

The method we present in this paper is dedicated to Autonomous Guided Vehicles where the estimation of the free-space in front of the mobile robot is more crucial than the identification of the obstacles. In other words, all the elements which are lying above the road surface represent an obstacle, whatever their shape or height. Hence, the kerbs and the green strips which generally highlight the roads boundaries have to be segmented as obstacles if they have a non-null elevation.

Our aim is also the detection of all the obstacles which reduce the free-space to allow an overtaking maneuver if it is necessary. This study represents the last improvement of a feature-based method we have been developing for some years ( [START_REF] Simond | Homography from a Vanishing Point in Urban Scenes[END_REF]) to locate a mobile robot when GPS data are corrupted. The method is inspired by the works of Okutomi [START_REF] Okutomi | Robust Estimation of Planar Regions for Visual Navigation Using sequential Stereo Images[END_REF], which proposed to identify the road surface with the homography induced in images captured by a calibrated stereo-rig. The computation is performed with a dense method whereas we use points and lines extracted with a Harris and a Canny detector.

The main restriction of the method remains on the assumption of flat road that seems not so restrictive in real conditions due to the limited field of view in case of urban traffic. We hence assume that the free space in front of the vehicle does not exceed some meters, that corresponds to the security distance with the preceding vehicle. Nevertheless, the proposed method has the significant interest not to depend to the calibration of the stereo-rig.

This paper is organized as follows: in Section II, we propose a method to compute easily a bird eye view of the road plane. The expression in the IPM image of the features lying above the road plane is presented in section III. Some results obtained on video sequences are shown in Section IV. The conclusion and a list of possible research directions are given in Section V.

II. INVERSE PERSPECTIVE MAPPING A. Interest

We present in this part a simple methodology to compute an Inverse Perspective Mapping (IPM) image [START_REF] Mallot | Inverse perspective mapping simplifies optical flow computation and obstacle detection[END_REF] or a Virtual Projection Plane, as it is called by some authors. The IPM transformation consist on modifying the angle of view under which a scene is acquired to remove the perspective effect. In a practical way, the IPM transformation consist to set the angles of rotation (mainly the tilt one) to values such as the normal to the plane n and the focal axis of the virtual camera z v are opposite: n.z v = -1, according to the Fig. 1. The virtual camera follows the vehicle path but has a constant ordinate z v along its z axis, whatever the slope of the road s. The method only requires the extraction of two coplanar parallel lines. Due to the perspective effect, the 3D lines which are parallel have projections in image which converge to a Vanishing Point. The IPM transformation also consist on determining the point of view of a virtual camera where the projections of the two lines appear parallel.

In case of structured road, the navigation lanes are usually highlighted with parallel painted markers. Hence, the extraction of two parallel lines along a straight path is relatively easy. Many authors propose several models of the road projection that mainly differ on the free field of view, the presence of road markers or not, the number of degrees of liberty the road model assumes. We personally use the simplest model: the road borders are linearized at the foreground and converge to the Dominant Vanishing Point (DVP) x whatever the road curvature is. The reader can retrieve in [START_REF] Simond | Vehicle Trajectory from an Uncalibrated Stereo-Rig with Super-Homography[END_REF] the tracking process we have developed to segment the projection of the road in images.

In the following, we detail the method to compute an IPM while in the second part, we discuss on the interest of the method.

B. Methodology

Lets now considering R i and R v the frameworks relative to the camera and the virtual camera. We assume in the following that the original and virtual images have the same dimensions [W, L] and 256 gray levels.

We also introduce two horizontal lines l t and l b which have respectively v t and v b v-ordinates. In the original image, the intersections of these two lines with the pencil of road markers (restricted to the first and the last median) and the image vertical boundaries generate two quadruplets of points. The first one p πik with k ∈ {1, 2, 3, 4} forms a trapezoid, as in Fig. 2, while the second one p eik which is the bottom part of the original image forms a rectangle.

If the road has a constant width in the scene, its projection in the virtual framework is also constant w πv0 , that means the four points p πvk form a parallelogram in R v whereas the image is represented with a truncated cone. The IPM transformation can hence be resumed as a homography H vi such as:

p πvk ≃ H vi p πik (1)
p evk ≃ H vi p eik (2)
where ≃ represents a projective equality. 1) Initialization of the IPM process: At this point, the coordinates of the two quadruplets in the virtual framework R v are unknown. We only suppose that the aperture angle of the camera lens α is supposed to be known and the elevation z v of the virtual camera relative to the plane is constant. That means if we fixed the location of the cone summit c, the projection of the original image in R v is invariant whatever the camera pose in R i . To maintain a coherence in the projective transformation, we consider that the top and the bottom lines are invariant:

l bv ≃ H -t vi l bi ≃ l b (3) l tv ≃ H -t vi l ti ≃ l t (4)
The quadruplet p evk also forms an isosceles trapezoid whose the little base is supported by l b and the vertical axis represents a symmetry axis.

In an other hand, the cross-ratio and the incidence are the lonely properties left invariant by a projective transformation. For the two quadruplets of points, the consistency of the cross-ratio is checked 1 :

CR(p ei1 , p ri1 , p ri2 , p ei2 ) = CR(p ev1 , p rv1 , p rv2 , p ev2 )(5) CR(p ei4 , p ri4 , p ri3 , p ei3 ) = CR(p ev4 , p rv4 , p rv3 , p ev3 )(6)
To entirely constraint the system form with the two quadruplets p πvk and p evk , we now have to determine the v-ordinates v c of the cone summit (the u-ordinate is assumed fixed u c = W/2) and v b of the line l b .

While there exist an infinity of solutions to the IPM problem, a trivial solution consist to fix the width of the plane projection w πv0 that determines the elevation z v of the virtual camera and induces the v-ordinate v c of the cone summit. At this point, the pan angle γ between the camera axis and the road direction is unknown. Therefore, we personally prefer to fix w πv = 200 pixels to compute the locations of c, p ev3 , p ev4 . The width of the little base w eb of the trapezoid is also computed thanks to the consistency of the cross-ratio (6): the points p ev4 , p ev3 have invariant locations and form an isosceles triangle with the cone summit c. The v-ordinate v c is such as:

v c = v b + w eb /2 tan(alpha/2) ( 7
)
where v b = L is the height of the virtual image. The width of the great base w et is computed thanks to the consistency of the cross-ratio (5), according to:

w πv = w πvb = w πvt = p πv3 -p πv4 = p πv2 -p πv1 (8)
The v-ordinate v t of the top line l t can also be computed as:

v t = w πv /2 tan(alpha/2) -v b (9)
and the road width w πv0 is linked to the camera orientation γ as:

w πv0 = w πv cos(γ) (10)
1 Whatever the quadruplet of aligned points {A, B, C, D},

CR(A, B, C, D) = B-A C-B / D-A C-D = cte with p πv1 -p πv4 = (v b -v t ) tan(γ).
The homography H vi is now totally defined thanks to the quadruplet p eik and their image p evk . We show in Fig. 3 an IPM image computed with this procedure. Although the camera is not calibrated, we compute up to a scale factor the orientation γ and the lateral position relatively to the road direction and width. We propose on our web site2 the result of this transformation along a video sequence where the pan and tilt rotation motions of the camera vary within large ranges. 2) Simplifications: The consistency of the cross-ratio is still checked if an affine transformation is applied to the quadruplet. For the lines l b and l t , the IPM transformation corresponds to a translation with a scale factor. To compute w eb then w tb , we had better to check directly the following ratio rather than use ( 5) and [START_REF] Heisele | Motion-based object detection and tracking in color images sequences[END_REF], which lead to quadratic equation where we have to determine the right solution:

w πi w ei = w πv w ev (11)
where w ei = W , the image width.

As soon as the location of the optical axis of the virtual camera c is fixed, the determination of the top v-ordinate v t highly depends on the location of the DVP x. The Thales theorem applied to the triangles {x, p πi4 , p πi3 } and {c, p ev1 , p ev2 }, crossed by the parallel lines l b and l t , allows the following equalities:

w πit w πib = v t -v x v b -v x (12)
w evb w evt = v c -v t v c -v b (13)
According to v c , v b , w evb are now fixed and w πib , v x are measured in the original image. The other variables w πit , w evt depend on the v t ordinate. Parallely, a new equality can be introduce from ( 8) and [START_REF] Broggi | A Single Frame Stereo Vision System for Reliable Obstacle Detection during Darpa Grand Challenge 2005[END_REF]:

w πit w πib = w πit w eit . w eib w πib = w πvt w evt . w evb w πvb = w evb w evt (14)
The introduction of ( 12) in the last equation implies:

v t -v x v b -v x = v c -v t v c -v b (15)
Hence, the v-ordinate of the l t line originally fixed with (9) becomes:

v t = v c + v x -v b . (16)

C. Discussion

The main drawback of the method is that the IPM image doesn't have squared pixels. The method is simple and only requires the aperture angle of the camera lens. But due to the perspective projection, the bottom line l b doesn't have a projection in R v on the same v-ordinate v t but rather some pixels above. The trouble is we can't estimate how many. Therefore a calibration stage is necessary to determine the two different scale factors which have to be applied on each axis to authorize metric measures.

As soon as two coplanar parallel lines can be detected, the method provides a solution to the IPM transformation. In the other cases, the IPM transformation can be computed thanks to the estimation of the camera motion from the last pose where the IPM was estimated. Lets consider H vi (k 0 ) the last IPM transformation available at the frame k = k 0 . At the frame k = k 0 +N , the camera motion can be computed from the camera pose at k = k 0 as the sum of the N elementary motions between each camera pose:

H k0+N,k0 ≃ H k0+N,k0+N-1 • • • H k0+2,k0+1 . H k0+1,k0 (17)
The IPM transformation at the frame k = k 0 + N is also computed as:

H vi (k 0 + N) ≃ H k0+N,k0 . H vi (k 0 ) . H -1 k0+N,k0 (18)
Thanks to the computation of the Super-Homography (cf. [START_REF] Simond | Vehicle Trajectory from an Uncalibrated Stereo-Rig with Super-Homography[END_REF]), we can afford to reconstruct the road plane with a mosaicing of warped images all expressed in the first one. The reference image was a bird-eye view of the first image of the video sequence of 312 frames long that corresponds to a path greater than 250m. We hence would like to stress that if the vehicle path reach an unstructured road or an intersection, the loss of the two parallel coplanar lines during a short period is not crucial.

III. OBSTACLE DETECTION THANKS TO THE SUPER-HOMOGRAPHY

A. Methodology

The accuracy of the homographies extracted from the computation of the Super-Homography leads us to detect all the elements which are above the considered plane. Thanks to the feature-based method we have developed, which identifies coplanar homologous lines and points matched in several views of the same plane, we had the choice to detect the obstacles with points and/or lines. We nevertheless use edges because the matching process for points correlation is time consuming and the obstacles have sometimes uniform textures.

We present in this section a two-stages method where the coplanar edges are first identified to be removed, then the remaining edges are projected with an IPM transformation to determine the regions free from obstacles. We also only focus on edges which have at least one extremity under the horizon line, assumed as the line which has the same vordinate as the DVP.

B. Identification of edges above the horizon line

The edges which are crossing the horizon line have projections in the IPM which are passing through the projection center c of the virtual camera. With the real data all the vertical edges don't converge to the same point but they generally have extremities out of the image boundaries as in Fig. 4. The identification of the edges above the horizon line is also straight forward thanks to the IPM transformation we detailed in the last section.

C. Identification of coplanar edges with homography

Lets considering two views of a same planar scene and H 21 the homography induced by the plane between the images I 1 and I 2 . The projection according H 21 of the features of one image into the other one allows the identification of the coplanar features because only the coplanar features are overlapping with their homologous. The main difficulty we also face is the edges detected with a Canny detector are not strictly overlapping. A 3D edge can be plotted with a long segment in one image and with zero, one, two or more in the other.

The matching of edges between two images thanks to the homography H 21 be nevertheless performed with a two steps process. The pre-selection stage consist on assuming as potential homologous the couples of edges whose the extremities, computed by retro-projection, are close to their homologous with a tolerance of 3 pixels. The pre-selection is checked for edges of I 1 in I 2 and inversely with H -1 21 . As soon as the couples have been identified, the selection stage consists on verifying the planarity assumption with checking if the couple of edges have common u-or v-ordinates when they are expressed in the same image.

The edges, shown in Fig. 5(a), issued from the left (red) and right (green) cameras, check the overlapping criterion. As expected, the majority of outliers (above the road plane) are horizontal edges whose the pan angle between the camera framework is not large enough to discriminate them with such a method whatever the camera motion. The best way is certainly do not take account the horizontal edges because the determination of their coplanarity is not reliable.

We actually obtain similar results if we only consider the camera motion between two consecutive poses. The trouble is the rate of ouliers increases when the camera motion is not large enough. With such a configuration, the camera frame-rate should be synchronized with the vehicle speed to maintain a critical change of view. The best method certainly consists on tracking the coplanar edges with two couples of views: stereo and consecutive. But as far as we are concerned, the results provided by the stereo-images are sufficiently right to allow a reliable discrimination even some coplanar edges can't be be segmented as if when they are occluded by an obstacle in one image.

D. Identification of coplanar edges with correlation

The last sequential operations allow the detection of the majority of the coplanar edges. The left coplanar edges don't fit precisely the contours or don't have homologous in the second view due to occlusion or illumination effects. The rejection of these edges requires one additional process to verify their planarity. Generally, Sum of Absolute Difference (SAD) image computed from the difference between the warping of one image into an other according to a transformation, highlight the scene regions where the assumption (planarity, same relative motion) is checked because their gray levels are the darkest.

The discrimination between edges lying on or above the road with SAD image is not easy due to the illumination conditions the areas close to the edges often appear clearest, whatever the quality of the estimation of homography. Therefore we prefer using normalized criteria on 3 pixels large region centered on each remaining edge in the original and warped images. We fix two new criteria based on the correlation coefficient to consider edges as coplanar: it has to be either greater than 0.75 or greater than 0 with a variance less than 400 and the meaning gray level have to be less than 24 for a 256 gray levels image.

Although some edges have no homologous in the other image, all the (blue) coplanar edges are rejected with the application of the correlation criterion. The (red) remaining edges of the Fig. 5(b) have to be considered as lying above the road plane.

IV. RESULTS

A. Camera rotation motions estimated with IPM

The decomposition of the homography H πvi requires the knowledge of the internal parameters of the camera K, an estimation of the normal to plane n = [0; -1; 0] t (close to the vertical axis) and an estimation of the camera elevation d in the absolute framework, linked to the road plane. A coarse calibration of the stereo-rig allows us to estimate the three rotation R and translation t motions required to compute the image and the evolution of the normal of the plane n along the vehicle path according to:

H πvi = K[R + t.n t d ]K -1 (19)
During the Antibes sequence, the test-vehicle follows a straight path, stops before a parked vehicle, overtakes it then continues its straight motion. We present the decompositions of the IPM homography in Fig. 6. The smoothness of the different time diagrams is relevant: the IPM homography is only computed with four couples of matched points. The numerical values of the rotation motions are not correct: assuming that the optical axis z is abroad 20 • under the horizon line, we expect a rotation around the x axis higher than -17 • . Such an error is mainly due to the scale factors we introduce with our coarse calibration. The decomposition of the translation motions seem reliable: to obtain an IPM view, the translation is mainly support by the z and in a second time by the y, as in Fig. 1.

B. Free-space segmentation

In the preceding section, we explain how to discriminate between edges lying on or above a plane. The best way to interpret this result is using an IPM transformation to produce a polar chart of the obstacle edges. We can hence segment in the 2D scene the areas which are free from obstacles in front of the vehicle thanks to a bird-eye view. The constant constraints (road width and aperture angle) on IPM images allow a correct segmentation of the obstacles in the road framework at two scale-factors due to the calibration of the stereo-rig is not available.

Opposite to most of the methods developed on the reconstruction of the 3D scene based on (V-)disparity, the proposed method does not make assumption on the obstacles elevation. The limitation of the method is the length of the edges which have to be greater than 10 pixels, that allows us to detect the majority of the white small bollards which are lying on left road boundary of the Antibes sequence like in Fig. 7. That is our goal for an implementation in an AGV.

Free-space at the frame 0908

V. CONCLUSIONS AND FUTURE WORKS

We present in this paper a methodology to compute reliable IPM transformation as soon as two parallel lines can be extracted from an image. The scene elements whose edges don't check the planarity criteria are segmented as obstacles. The expression of these remaining edges in the IPM view allow a coarse segmentation of the obstacles location in front of the vehicle.

The described methodology have to be improved to perform an estimation of the free-space. The stereo-vision system has to be calibrated to represent the edges in a polar chart and the region free from obstacles has to be unified and precisely segmented from the bottom image to the intersection of the obstacles with the road plane.

Fig. 1 .

 1 Fig.1. Principle of the IPM transformation: the virtual camera Rv has a constant elevation zv along the third zv axis which is parallel to the normal to the plane n whatever the camera pose R i and the road slope s while the aperture angle α is maintained.

Fig. 2 .

 2 Fig. 2. Principle of the IPM transformation: the road boundaries converge to the DVP x in R i and appear parallel in Rv. The intersections of the image boundaries and the road borders with the two parallel lines l b and l t respectively generate the quadruplets {p e1 , p e2 , p e3 , p e4 } and {p π1 , p π2 , p π3 , p π4 }.

Fig. 3 .

 3 Fig.3. The IPM transformation is computed thanks to the intersections of the first and the third blue medians of the road boundaries with the green l b and l t lines. The second median is shown for convenience.

Fig. 4 .

 4 Fig. 4. Identification with an IPM transformation of the edges whose extremities are above and under the horizon line.

  Identification of coplanar edges which are overlapping by retro-projection between the left and right images. Identification of coplanar edges (blue) which have similar neighborhood in original and warped images.

Fig. 5 .

 5 Fig. 5. Sequential search of the coplanar edges in the left (red) image thanks to the right (green) image and the stereo homography Hst.

Fig. 6 .

 6 Fig.6. Orientation of the camera (brown) and decomposition of the H πvi homography along the Antibes sequence. The numerical values are given up to a scale factor due to the pixels of the IPM image are not squared and the camera is not calibrated (d = 1m).

Fig. 7 .

 7 Fig. 7. The free-space in front of the camera is the lowest region of the IPM image under the obstacles edges. The edges which are crossing the horizon line are removed according to the methodology developed in III-B.

http://www.lara.prd.fr/Nicolas.Simond

VI. ACKNOWLEDGMENTS

The author would like to thank Akihito Seki from the Okotumi Lab. for the interesting exchanges of these last months. This work is supported by the E.C. CyberCars 2 project.


Anthony Giacoma

David Dureisseix

Anthony Gravouil

Michel Rochette

A two-scale LATIN-PGD for efficient frictional contact problem solving

come   L'archive ouverte pluridisciplinaire

Contact mechanics is a fundamental element governing an interaction boundary phenomenon of one or multiple bodies. Due to its non-linear and non-smooth behavior, the frictional contact constitutive law is difficult to handle. Within the linear elastic isotropic mechanical framework, we assume small displacements hypothesis and quasi-static evolution. We consider non-regularized laws for the contact behavior (Signorini's conditions and Coulomb's model). The finite element method is used to discretize the mechanical problem.

Industrial complex problems involving large contact zones may lead to prohibitive computational time. Unfortunately, classic incremental solvers (e.g. Newton-Raphson) for contact problems could defeat the deployment of model reduction techniques well-known for their efficiency to reduce computing time and memory usage [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF][START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF][START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF].

Even if reduced basis methods seem to be inappropriate to non-smooth phenomenon like frictional contact, solutions of such problems show a great scale separability (both in space and in time). Thus, these solutions can be compressed into a small reduced number of vectors depicting different scales of the problem [START_REF] Giacoma | A multiscale large time increment/FAS algorithm with time-space model reduction for frictional contact problems[END_REF][START_REF] Giacoma | Techniques de rduction de modles espace-temps a posteriori et a priori pour les problmes de contact frottant[END_REF].

We propose to tackle the contact problem with the non-linear and non-incremental twoscaled iterative LATIN (Large Time Increment) solver [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF]. Then, unknowns are sought into a separated form (a time-space rank-one approximation) using a PGD strategy. The two-scale feature is dealt with a FAS-multigrid like algorithm. Thus, an efficient strategy is proposed combining the accelerating PGD method still keeping the strong physical meaning of the LATIN method. Moreover, a size control is proposed allowing an optimality mastering at each enriching stage of the growing PGD basis.

On academic cases, we show the efficiency of the proposed method. This new strategy requiring definitely less computational charge and memory usage may overtake full-order LATIN solver for given levels of accuracy.

Ongoing investigations concerning multiparametric studies and multiscale properties will be alluded.

  11th. World Congress on Computational Mechanics (WCCM XI) 5th. European Conference on Computational Mechanics (ECCM V) 6th. European Conference on Computational Fluid Dynamics (ECFD VI) July 20 -25, 2014, Barcelona, Spain

	A TWO-SCALE LATIN-PGD FOR EFFICIENT
	FRICTIONAL CONTACT PROBLEM SOLVING
	A. Giacoma 1,2, * , D. Dureisseix 1 , A. Gravouil 1,3 and M. Rochette 2
	1 Université de Lyon, INSA-Lyon, LaMCoS, CNRS UMR 5259, F-69621, Villeurbanne (France)
	2 ANSYS France, Villeurbanne (France)
	3 Institut universitaire de France
	Key words: Non-smooth mechanics, Multigrid, Multiscale, Reduced Order Model.