Jean-Antoine Désidéri

Application of MGDA to domain partitioning

Keywords: multiobjective optimization, descent direction, convex hull, Gram-Schmidt orthogonalization process optimisation multiobjectif, direction de descente, enveloppe convexe, processus d'orthogonalisation de Gram-Schmidt

. In this new report, the efficiency of the algorithm is tested in the context of a simulation by domain partitioning, as a technique to match the different interface components concurrently. For this, the very simple testcase of the finite-difference discretization of the Dirichlet problem over a square is considered. The study aims at assessing the performance of MGDA in a discretized functional setting. One of the main teachings is the necessiy, here found imperative, to normalize the gradients appropriately.

Introduction

Classically, in multi-objective optimization, several fundamental concepts are introduced: dominance in efficiency between design points, Pareto set, made of non-dominated solutions in design space, and Pareto front, its image in function space [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. The Pareto front provides the designer with the system maximum attainable performance. For complex systems, in particular those governed by partial-differential equations, a computational challenge is to devise algorithms permitting to identify numerically the Pareto set, or the most useful portions of it. In this respect, certain evolutionary strategies have been adapted to achieve this goal, and appear to provide the most robust algorithms. NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] is certainly one of the most widely-used methods for this purpose.

In the context of differentiable optimization, one would expect adequate strategies based on gradient evaluations to also be capable of capturing Pareto fronts, with less generality or robustness, but often far greater efficiency. However classical techniques, such as minimizing an agglomerated criterion, or one criterion at a time under the constraints of the others, are limited by hypotheses on the pattern of the Pareto front w.r.t. convexity and, or continuity. The Multiple-Gradient Descent Algorithm (MGDA), originally introduced in [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF], and again formalized in [START_REF] Désidéri | Multiple-gradient descent algorithm (mgda) for multiobjective optimization[END_REF], is based on a very general principle to define at each step a descent direction common to all criteria. This direction is the support of the minimum-norm element in the convex hull of the local gradients. Its efficacy to identify the Pareto front has been demonstrated in [START_REF] Zerbinati | Comparison between MGDA and PAES for Multi-Objective Optimization[END_REF] [START_REF] Désidéri | Multiple-gradient descent algorithm (mgda) for multiobjective optimization[END_REF] in a test-case in which the Pareto front was non-convex. The method was compared in efficiency with an evolutionary strategy, and was found to offer very promising performance.

More recently, a variant, MGDA II , has been proposed in which the descent direction is calculated by a direct procedure [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF], which provides a valuable simplification of implementation.

In this new report, the MGDA is tested in the context of a simulation by domain partitioning, as a technique to match the different interface components concurrently. For this, the very simple test-case of the finite-difference discretization of the Dirichlet problem over a square is considered.

The study aims at assessing the performance of MGDA in a discretized functional setting. One of the main teachings is the necessity, here found imperative, to normalize the gradients appropriately.

Dirichlet problem, domain partitioning and matching defects

We consider the problem of solving Laplace's equation,

-∆u = f (Ω) (1)
over the square

Ω = [-1, 1] × [-1, 1] (2)
subject to homogeneous boundary conditions:

u = 0 (Γ = ∂Ω) (3)
For this, the domain Ω is partitioned in four disjoint sub-domains:

         Ω 1 = [0, 1] × [0, 1] Ω 2 = [-1, 0] × [0, 1] Ω 3 = [-1, 0] × [-1, 0] Ω 4 = [0, 1] × [-1, 0] (4)
RR n°7968

with the following interfaces supported by the coordinate axes:

         γ 1 = { 0 ≤ x ≤ 1 ; y = 0 } γ 2 = { x = 0 ; 0 ≤ y ≤ 1} γ 3 = { -1 ≤ x ≤ 0 ; y = 0 } γ 4 = { x = 0 ; -1 ≤ y ≤ 0 } (5)
along which the following Dirichlet controls are applied:

         γ 1 : u = v 1 (x) γ 2 : u = v 2 (y)
γ 3 : u = v 3 (x) γ 4 : u = v 4 (y) (6)
(see Fig. 1).

-1 0 1 1 -1

Γ = ∂Ω : u = 0 γ 1 : u = v 1 (x) γ 3 : u = v 3 (x) γ 2 : u = v 2 (y) γ 4 : u = v 4 (y) Ω 1 Ω 2 Ω 4 Ω 3 x y
u = (u 1 , u 2 , u 3 , u 4 ) t (7)
(where the superscript t stands for transposition) is made based on a certain setting of the interface controls:

v = (v 1 , v 2 , v 3 , v 4 ) t (8)
by solving, possibly in parallel, the partial problems for i = 1, ..., 4:

         -∆u i = f u i = 0 u i = v i u i = v i+1 (Ω i ) (Γ ∩ ∂Ω i ) (γ i ) (γ i+1 ) (9)
In the above, and all throughout thereafter, by periodicity, the index i is understood modulo 4.

Since the interface controls are of Dirichlet type, the resulting compound solution u is continuous, and its derivative along each interface is also continuous. However, in general, unless the specified controls v i 's are equal to the restrictions of the global solution, the normal derivatives exhibit jump discontinuities, s i 's. Here, each interface is supported by a coordinate axis, and we adopt the following sign convention: on the interface γ i which is supported by the x (resp. y) axis for i = 1 and 3 (resp. 2 and 4), the jump, s i (x) (resp. s i (y)), is defined as the increment of the partial derivative ∂u/∂y (resp. ∂u/∂x) as y (resp. x) goes from 0 -to 0 + . Thus:

• over γ 1 (0 ≤ x ≤ 1 ; y = 0): s 1 (x) = ∂u ∂y (x, 0 + ) - ∂u ∂y (x, 0 -) = ∂u 1 ∂y - ∂u 4 ∂y (x, 0); • over γ 2 (x = 0 ; 0 ≤ y ≤ 1): s 2 (y) = ∂u ∂x (0 + , y) - ∂u ∂x (0 -, y) = ∂u 1 ∂x - ∂u 2 ∂x (0, y); • over γ 3 (-1 ≤ x ≤ 0 ; y = 0): s 3 (x) = ∂u ∂y (x, 0 + ) - ∂u ∂y (x, 0 -) = ∂u 2 ∂y - ∂u 3 ∂y (x, 0); • over γ 4 (x = 0 ; -1 ≤ y ≤ 0): s 4 (y) = ∂u ∂x (0 + , y) - ∂u ∂x (0 -, y) = ∂u 4 ∂x - ∂u 3 ∂x (0, y).
The above local measures of the defect in matching conditions can be associated with global measures defined as the functionals:

J i = γi 1 2 s 2 i w dγ i (10)
that is, explicitly:

J 1 = 1 0 1 2 s 1 (x) 2 w(x) dx J 2 = 1 0 1 2 s 2 (y) 2 w(y) dy (11)
J 3 = 0 -1 1 2 s 3 (x) 2 w(x) dx J 4 = 0 -1 1 2 s 4 (y) 2 w(y) dy (12)
Here, w(t) (t ∈ [0, 1]) is an optional weighting function, and w(-t) = w(t).

The jump s i depends on the partial solutions u i-1 and u i , which themselves, depend on (v i-1 , v i ) and (v i , v i+1 ) respectively. Hence, the integral J i depends on all four sub-controls except v i+2 . Nevertheless, these four integrals are thereafter considered as functionals of v.

The coordination problem is to devise a convergent iteration on the control v to satisfy in the limit the matching conditions

J 1 = J 2 = J 3 = J 4 = 0 (13)
To achieve this, the functional gradients are firstly established using the classical adjoint-equation approach, and several strategies are proposed and tested numerically.

Adjoint problems and functional gradients

A first calculation is made based on the four-component control v = (v 1 , v 2 , v 3 , v 4 ) t , resulting in the compound solution u = (u 1 , u 2 , u 3 , u 4 ) t , and the multi-component criterion J = (J 1 , J 2 , J 3 , J 4 ) t .

Then, one perturbs the control v of

v ′ = (v ′ 1 , v ′ 2 , v ′ 3 , v ′ 4 ) t = δv 1 (x), δv 2 (y), δv 3 (x), δv 4 (y) t (14)
Consequently, the compound solution u is perturbed of

u ′ = (u ′ 1 , u ′ 2 , u ′ 3 , u ′ 4 ) t = δu 1 (x, y), δu 2 (x, y), δu 3 (x, y), δu 4 (x, y) t ( 15
)
in which u ′ i is the solution of the following linearized system posed on sub-domain

Ω i :          ∆u ′ i = 0 u ′ i = 0 u ′ i = v ′ i u ′ i = v ′ i+1 (Ω i ) (Γ ∩ ∂Ω i ) (γ i ) (γ i+1 ) (16)
These state perturbations induce the following functional perturbations:

J ′ i = δJ i = γi s i s ′ i w dγ i (17) in which s ′ i = δs i . But:                            s 1 s ′ 1 = ∂u 1 ∂y - ∂u 4 ∂y ∂u ′ 1 ∂y - ∂u ′ 4 ∂y (x, 0) s 2 s ′ 2 = ∂u 1 ∂x - ∂u 2 ∂x ∂u ′ 1 ∂x - ∂u ′ 2 ∂x (0, y) s 3 s ′ 3 = ∂u 2 ∂y - ∂u 3 ∂y ∂u ′ 2 ∂y - ∂u ′ 3 ∂y (x, 0) s 4 s ′ 4 = ∂u 4 ∂x - ∂u 3 ∂x ∂u ′ 4 ∂x - ∂u ′ 3 ∂x (0, y) (18)
We now recall Green's formula for two functions φ and ψ ∈ H 2 (ω), for a simply-connected planar domain ω with smooth enough boundary. Since

ω φ ∆ψ = ω φ ∇.(∇ψ) = ∂ω φ ψ n - ω ∇φ.∇ψ ( 19
)
where n is outward unit vector normal to the boundary ∂ω, and ψ n = ∂ψ/∂n the normal derivative, and symmetrically,

ω ψ ∆φ = ω ψ ∇.(∇φ) = ∂ω ψ φ n - ω ∇ψ.∇φ (20)
where φ n = ∂φ/∂n is the normal derivative of φ, one has:

ω (φ ∆ψ -ψ ∆φ) = ∂ω (φ ψ n -ψ φ n ) (21)
Consider the following eight adjoint systems (two per sub-domain):

     ∆p i = 0 p i = 0 p i = s i w (Ω i ) (∂Ω i \γ i ) (γ i )      ∆q i = 0 q i = 0 q i = s i+1 w (Ω i ) (∂Ω i \γ i+1 ) (γ i+1 ) (22)

Inria

Then apply Green's formula, (21), to the eight cases corresponding to

ω = Ω i (i = 1, 2, 3, 4), φ = p i or q i , ψ = u ′ i ( 23
)
so that ∆φ = ∆ψ = 0 in ω giving:

∂Ωi φ ψ n = ∂Ωi ψ φ n (24)
On the boundary ∂ω = ∂Ω i :

• φ = 0 except for φ = p i = s i w along γ i , and φ = q i = s i+1 w along γ i+1 ;

• ψ = u ′ i =      v ′ i along γ i v ′ i+1 along γ i+1 0 along Γ ∩ Ω i
Hence (24) reduces to:

γi s i u ′ in w = γi p in v ′ i + γi+1 p in v ′ i+1 (25)
for φ = p i , and to:

γi+1 s i+1 u ′ in w = γi q in v ′ i + γi+1 q in v ′ i+1 (26)
for φ = q i . These two equations are now going to be particularized to sub-domains Ω i , for i = 1, ..., 4. For each sub-domain, attention must be paid to the orientation of the outward normal n along the two interfaces γ i and γ i+1 . For this, let ı and  be the unit vectors along the x and y axes.

Sub-domain Ω 1 : on γ 1 : n = n 14 = - ; on γ 2 , n = n 12 = -ı. Thus (25)-(26) write:

         1 0 s 1 (x) - ∂u ′ 1 ∂y (x, 0) w(x) dx = 1 0 - ∂p 1 ∂y (x, 0) v ′ 1 (x) dx + 1 0 - ∂p 1 ∂x (0, y) v ′ 2 (y) dy 1 0 s 2 (y) - ∂u ′ 1 ∂x (0, y) w(y) dy = 1 0 - ∂q 1 ∂y (x, 0) v ′ 1 (x) dx + 1 0 - ∂q 1 ∂x (0, y) v ′ 2 (y) dy (27) Sub-domain Ω 2 : on γ 2 : n = n 21 = + ı; on γ 3 , n = n 23 = - . Thus (25)-(26) write:          1 0 s 2 (y) ∂u ′ 2 ∂x (0, y) w(y) dy = 1 0 ∂p 2 ∂x (0, y) v ′ 2 (y) dy + 0 -1 - ∂p 2 ∂y (x, 0) v ′ 3 (x) dx 0 -1 s 3 (x) - ∂u ′ 2 ∂y (x, 0) w(x) dx = 1 0 ∂q 2 ∂x (0, y) v ′ 2 (y) dy + 0 -1 - ∂q 2 ∂y (x, 0) v ′ 3 (x) dx (28)
Sub-domain Ω 3 : on γ 3 : n = n 32 = +  ; on γ 4 , n = n 34 = + ı. Thus (25)-(26) write:

         0 -1 s 3 (x) ∂u ′ 3 ∂y (x, 0) w(x) dx = 0 -1 ∂p 3 ∂y (x, 0) v ′ 3 (x) dx + 0 -1 ∂p 3 ∂x (0, y) v ′ 4 (y) dy 0 -1 s 4 (y) ∂u ′ 3 ∂x (0, y) w(y) dy = 0 -1 ∂q 3 ∂y (x, 0) v ′ 3 (x) dx + 0 -1 ∂q 3 ∂x (0, y) v ′ 4 (y) dy (29)
Sub-domain Ω 4 : on γ 4 : n = n 43 = -ı; on γ 1 , n = n 41 = +  . Thus (25)-(26) write: 27)-( 30) are injected in ( 18) and ( 17) to get:

         0 -1 s 4 (y) - ∂u ′ 4 ∂x (0, y) w(y) dy = 0 -1 - ∂p 4 ∂x (0, y) v ′ 4 (y) dy + 1 0 ∂p 4 ∂y (x, 0) v ′ 1 (x) dx 1 0 s 1 (x) ∂u ′ 4 ∂y (x, 0) w(x) dx = 0 -1 - ∂q 4 ∂x (0, y) v ′ 4 (y) dy + 1 0 ∂q 4 ∂y (x, 0) v ′ 1 (x) dx (30) Lastly, (
                                                                   J ′ 1 = 1 0 s 1 (x) s ′ 1 (x) w(x) dx = 1 0 s 1 (x) ∂u ′ 1 ∂y - ∂u ′ 4 ∂y (x, 0) w(x) dx = 1 0 ∂(p 1 -q 4 ) ∂y (x, 0) v ′ 1 (x) dx + 1 0 ∂p 1 ∂x (0, y) v ′ 2 (y) dy + 0 -1 ∂q 4 ∂x (0, y) v ′ 4 (y) dy J ′ 2 = 1 0 s 2 (y) s ′ 2 (y) w(y) dy = 1 0 s 2 (y) ∂u ′ 1 ∂x - ∂u ′ 2 ∂x (0, y) w(y) dy = 1 0 ∂q 1 ∂y (x, 0) v ′ 1 (x) dx + 1 0 ∂(q 1 -p 2 ) ∂x (0, y) v ′ 2 (y) dy + 0 -1 ∂p 2 ∂y (x, 0) v ′ 3 (x) dx J ′ 3 = 0 -1 s 3 (x) s ′ 3 (x) w(x) dx = 0 -1 s 3 (x) ∂u ′ 2 ∂y - ∂u ′ 3 ∂y (x, 0) w(x) dx = - 1 0 ∂q 2 ∂x (0, y) v ′ 2 (y) dy + 0 -1 ∂(q 2 -p 3 ) ∂y (x, 0) v ′ 3 (x) dx - 0 -1 ∂p 3 ∂x (0, y) v ′ 4 (y) dy J ′ 4 = 0 -1 s 4 (y) s ′ 4 (y) w(y) dy = 0 -1 s 4 (y) ∂u ′ 4 ∂x - ∂u ′ 3 ∂x (0, y) w(y) dy = - 1 0 ∂p 4 ∂y (x, 0) v ′ 1 (x) dx - 0 -1 ∂q 3 ∂y (x, 0) v ′ 3 (x) dx + 0 -1 ∂(p 4 -q 3 ) ∂x (0, y) v ′ 4 (y) dy (31)
These formulas are of the form:

J ′ i = 4 j=1 γj G i,j v ′ j dγ j (i = 1, ..., 4) (32)
in which the kernels, {G i,j }, are partial gradients given in terms of the partial derivatives of the eight adjoint states {p i , q i } (i=1,...,4) . These equations can also be collected in the following symbolic matrix form:

[J ′ ] = (J ′ 1 , J ′ 2 , J ′ 3 , J ′ 4 ) t =                 1 0 ∂(p 1 -q 4 ) ∂y 1 0 ∂p 1 ∂x 0 0 -1 ∂q 4 ∂x 1 0 ∂q 1 ∂y 1 0 ∂(q 1 -p 2 ) ∂x 0 -1 ∂p 2 ∂y 0 0 - 1 0 ∂q 2 ∂x 0 -1 ∂(q 2 -p 3 ) ∂y - 0 -1 ∂p 3 ∂x - 1 0 ∂p 4 ∂y 0 - 0 -1 ∂q 3 ∂y 0 -1 ∂(p 4 -q 3 ) ∂x                        v ′ 1 (x) dx v ′ 2 (y) dy v ′ 3 (x) dx v ′ 4 (y) dy       
(33) in which the partial derivatives w.r.t. x (resp. y) are evaluated at (0, y) (resp. (x, 0)).

Remark : to calculate the co-state variables p i and q i over sub-domain Ω i , firstly, one constructs two continuous functions pi (x, y) and qi (x, y) that satisfy the same boundary conditions respectively. These functions are listed below along with their Laplacians:

                             p1 (x, y) = (1 -y)s 1 (x); φ 1 := ∆p 1 = (1 -y)s ′′ 1 (x) q1 (x, y) = (1 -x)s 2 (y); ψ 1 := ∆q 1 = (1 -x)s ′′ 2 (y) p2 (x, y) = (1 + x)s 2 (y); φ 2 := ∆p 2 = (1 + x)s ′′ 2 (y) q2 (x, y) = (1 -y)s 3 (x); ψ 2 := ∆q 2 = (1 -y)s ′′ 3 (x) p3 (x, y) = (1 + y)s 3 (x); φ 3 := ∆p 3 = (1 + y)s ′′ 3 (x) q3 (x, y) = (1 + x)s 4 (y); ψ 3 := ∆q 3 = (1 + x)s ′′ 4 (y) p4 (x, y) = (1 -x)s 4 (y); φ 4 := ∆p 4 = (1 -x)s ′′ 4 (y) q4 (x, y) = (1 + y)s 1 (x); ψ 4 := ∆q 4 = (1 + y)s ′′ 1 (x) (34)
Then one lets:

p i = pi + λ i q i = qi + µ i (i = 1, ..., 4) (35)
so that, the new unknown variables, λ i (x, y) and µ i (x, y) (i = 1, ..., 4), are the solutions of the problems:

λ i , µ i ∈ H 1 0 (Ω i ) : -∆λ i = φ i -∆µ i = ψ i (i = 1, ..., 4) (36)
Thus, the unknowns λ i and µ i are the solutions of problems of the same form as the Poisson problem for u, except that the functions φ i and ψ i respectively are assigned to the right-hand side in place of f .

Discretization

For purpose of numerical treatment, we assume that each sub-problem is discretized by standard centered finite-differences over a uniform (sub-)mesh of dimension N X × N Y rectangular cells.

For example, for sub-domain Ω 1 , the partial solution u 1 is approximated by the solution of the system:

- u j-1,k -2u j,k + u j+1,k h 2 X - u j,k-1 -2u j,k + u j,k+1 h 2 Y = f j,k (j = 1, ..., N X -1, k = 1, ..., N Y -1) (37) where h X = 1/N X , h Y = 1/N Y , f j,k = f (jh,
kh), and only the nodal indices (j, k) are indicated for notational simplicity. The homogeneous boundary condition u 1 = 0 on the portion Γ ∩ ∂Ω 1 of the outer boundary is enforced by letting:

u j,NY = 0 (j = 0, ..., N X ) u NX ,k = 0 (k = 0, ..., N Y ) (38)
The controls on the interfaces γ 1 and γ 2 are enforced by letting:

γ 1 : u j,0 = v 1,j (j = 1, ..., N X -1) γ 2 : u 0,k = v 2,k (k = 1, ..., N Y -1) (39)
By passing the controls to the right-hand side, the system for the interior unknowns writes:

A h u h = f h (40)
in which the interior nodal values are chosen to be ordered by columns:

u h = u 1,1 , u 1,2 , ..., u 1,NY -1 1st column , u 2,1 , u 2,2 , ..., u 2,NY -1 2nd column , ..., u NX -1,1 , u NX -1,2 , ..., u NX -1,NY -1 NX -1st column , t (41)

RR n°7968

The vector f h contains the corresponding nodal values of the function f appropriately altered by the application of the boundary conditions, and the matrix A h is a Krönecker sum:

A h = T X ⊕ T Y = T X ⊗ I Y + I X ⊗ T Y (42)
in which I X and I Y are the identity matrices of order N X -1 and N Y -1 respectively, and T X and T Y are the tridiagonal matrices of order N X -1 and N Y -1 respectively associated with the (one-dimensional) centered second-difference operator when homogeneous boundary conditions are applied at both limits:

T X = 1 h 2 X        2 -1 -1 2 -1 . . . . . . . . . -1 2 -1 -1 2        (NX -1) T Y = 1 h 2 Y        2 -1 -1 2 -1 . . . . . . . . . -1 2 -1 -1 2        (NY -1)
(43) where the orders of the matrices are indicated in parentheses at their lower right for convenience.

The diagonalization of the above matrices is well-known (see for example [START_REF] Désidéri | Modèles discrets et schémas itératifs. Application aux algorithmes multigrilles et multidomaines[END_REF]):

T X = Ω t X Λ X Ω X T Y = Ω t Y Λ Y Ω Y (44)
in which the orthogonal matrix Ω X (resp. Ω Y ) is associated with the discrete sine-transform:

       (Ω X ) j,mX = 2h X sin jθ (mX ) X = 2h X sin jm X π N X (j, m X = 1, ..., N X -1) (Ω Y ) k,mY = 2h Y sin kθ (mY ) Y = 2h Y sin km Y π N Y (k, m Y = 1, ..., N Y -1) (45)
Here, θ

(mX ) X = m X π/N X (m X = 1, ..., N X -1) and θ (mY ) Y = m Y π/N Y (m Y = 1, ..., N Y -1
) are frequency parameters. Hence, in this case where Dirichlet conditions are applied, the orthogonal matrices Ω X and Ω Y are also symmetric:

(Ω X ) t = Ω X , Ω 2 X = I X , (Ω Y ) t = Ω Y , Ω 2 Y = I Y (46)
The corresponding eigenvalues are given by:

             λ X ,mX = (Λ X ) mX ,mX = 2 -2 cos θ (mX ) X h 2 X = 2 -2 cos (m X π/N X ) h 2 X (m X = 1, ..., N X -1) λ Y ,mY = (Λ Y ) mY ,mY = 2 -2 cos θ (mY ) Y h 2 Y = 2 -2 cos (m Y π/N Y ) h 2 Y (m Y = 1, ..., N Y -1)
(47) Consequently, the diagonalization of the matrix A h is also known:

A h = (Ω X ⊗ Ω Y ) (Λ X ⊕ Λ Y ) (Ω X ⊗ Ω Y ) (48)
This writing corresponds to the discrete form of the separation of variables. It permits the direct inversion of (40) by repeated application of one-dimensional operators:

u h = (Ω X ⊗ Ω Y ) (Λ X ⊕ Λ Y ) -1 (Ω X ⊗ Ω Y ) f h (49)
For more general boundary-value problems, and in particular when the mesh is not tensorial, this direct inversion is not possible and sub-problems such as (40) are then solved by more general techniques (relaxation, multigrid, etc).

The partial solutions u 2 , u 3 and u 4 associated with the sub-domains Ω 2 , Ω 3 and Ω 4 are obtained by similar discretizations and solution methods, except that the indices run within other bounds. The definition of the discrete solution u h = {u j,k } is then extended to include all four partial solutions by letting the index j (resp. k) run between ±N X (resp. ±N Y ). Lastly, the value of u h at the center of the square, origin of the coordinate system (x = y = 0 ), is calculated in terms of the values at the 4 surrounding points to satisfy the discrete Laplace equation:

- u -1,0 -2u 0,0 + u 1,0 h 2 X - u 0,-1 -2u 0,0 + u 0,1 h 2 Y = f (0, 0) := f 0,0 (50)
giving:

u 0,0 = u 1,0 + u -1,0 h 2 X + u 0,1 + u 0,-1 h 2 Y + f 0,0 2 h 2 X + 2 h 2 Y = v 1,1 + v 3,N X -1 h 2 X + v 2,1 + v 4,N Y -1 h 2 Y + f 0,0 2 h 2 X + 2 h 2 Y ( 51
)
This condition imposes discretely the satisfaction of the elliptic PDE at the origin. Hence it is viewed as a regularity condition, and the partial derivatives at this point are considered smooth (no jumps).

Once the compound solution u h calculated, the jumps at the interfaces are evaluated. Here, for each interface, we need approximations of the partial derivative in the direction normal to the interface on both sides of the interface. Second-order one-sided finite-differences are used for this purpose. For example, for the interface γ 1 , recalling the definition of s 1 (x), this gives at x = x j = jh X (j = 1, ..., N X -1):

s 1,j := s 1 (x j ) = ∂u 1 ∂y - ∂u 4 ∂y (x j , 0) . = -3u j,0 + 4u j,1 -u j,2 2h Y - 3u j,0 -4u j,-1 + u j,-2 2h Y = -6u j,0 + 4 (u j,1 + u j,-1 ) -(u j,2 + u j,-2 ) 2h Y (52)
at endpoints, as mentioned above, the jumps are equal to 0:

s 1,0 = s 1,N X = 0 (53)
For the interface γ 3 , we apply a shift of -N X on the index j:

s 3,j := s 3 (x j-NX ) . = -6u j-NX ,0 + 4 (u j-NX ,1 + u j-NX ,-1 ) -(u j-NX ,2 + u j-NX ,-2 ) 2h Y (54)
At endpoints:

s 3,0 = s 3,N X = 0 (55)
Similarly, for k = 1, ..., N Y -1:

s 2,k := s 2 (y k ) . = -6u 0,k + 4 (u 1,k + u -1,k ) -(u 2,k + u -2,k ) 2h X (56)
and, at endpoints:

s 2,0 = s 2,N Y = 0 (57)
and:

s 4,k := s 2 (y k-NY ) . = -6u 0,k-NY + 4 (u 1,k-NY + u -1,k-NY ) -(u 2,k-NY + u -2,k-NY ) 2h X (58)
and, at endpoints

s 4,0 = s 4,N Y = 0 (59)
Consistently with this discretization, the functionals {J i } (i = 1, ..., 4) are approximated by numerical quadrature; specifically, the trapezoidal rule is used:

J 1 . = h Z 2 NX -1 j=1 s 1 2 ,j w j , J 2 . = h Y 2 NY -1 k=1 s 2 2 ,k w k , J 3 . = h X 2 NX -1 j=1 s 3 2 ,j w j , J 4 . = h Y 2 NY -1 k=1 s 4 2 ,k w k . (60)
Once all criteria {J i } (i=1,...,4 ) evaluated, to prepare the calculation of their gradients, two adjoint problems per sub-domain Ω i (i = 1, ..., 4) need be solved to get the co-state functions p i and q i . Using again second-order central differencing, the unknown variables λ i and µ i are calculated by the same direct solver for i = 1, ..., 4, and the co-state variables p i and q i are then computed straightforwardly.

Then, the partial derivatives of p i and q i along the two controlled interfaces are approximated by second-order one-sided finite-differences. Consequently, nodal values along these interfaces are known for the gradients G i,j ,k (i, j = 1, ..., 4; k = 1, .., N Z -1 where Z stands for X or Y ) in (32). We also set G i,j ,k = 0 at endpoints (k = 0 and N Z ). The integral in (32) is approximated by the trapezoidal rule:

J ′ i . = 4 j=1 NZ k=1 G i,j .k v j ′ ,k h Z (i = 1, ..., 4) (61)
Hence, the discrete gradient of the criterion J i w.r.t. the nodal values of the control v j is given by:

∂J i ∂v j ,k . = G i,j .k h Z (62)
For each criterion J i , four such discrete gradients are calculated (one per control v j ), except that one of them is equal to 0. These four vectors are assembled in one, thereafter denoted ∇J i , of dimension 2(N X + N Y -2).

Now, knowing (second-order approximations of) the criteria {J i } (i=1,...4) and their gradients {∇J i } (i=1,...,4) w.r.t. the 2(N X + N Y -2) nodal controls, we need to set up a strategy to iterate on these controls to satisfy the matching conditions at convergence.

Gradient-based coordination iterations

Our main objective is to compare the standard steepest-descent method with the Multiple-Gradient Descent Algorithm (MGDA) as potential iterative methods to satisfy the matching conditions by driving the defect functionals to 0.

Conventional steepest-descent method

In the conventional approach, one considers a single matching defect measure, treating all interfaces as one:

J = 4 i=1 J i (63)
The discrete gradient is then simply equal to the sum of the individual contributions of the interfaces:

∇J = 4 i=1 ∇J i (64)
The above global criterion can then be driven to 0 by the classical steepest-descent method [START_REF]WikipediA: The Free Encyclopedia[END_REF] [7]: at iteration ℓ, the control v is updated proportionally to (the opposite of) the discrete gradient:

v (ℓ+1) = v (ℓ) -ρ ℓ ∇J (ℓ) (65)
Inria for some appropriate positive step-size ρ ℓ (see below), and a new compound solution u (ℓ+1) is calculated, the defect-functional and its gradient reevaluated, and so on until a satisfactory convergence is achieved.

Strictly speaking, in the standard steepest-descent method, once the direction of search is identified, by the calculation of the gradient ∇J (ℓ) , the step-size ρ ℓ is defined via a one-dimensional minimization:

ρ ℓ = Argmin ρ j(ρ) ; j(ρ) := J v (ℓ) -ρ∇J (ℓ) (66)
This minimization is usually carried out by a numerical procedure. However here, we know of an additional information: the targeted value of J is known: J = 0. An estimation of the variation of J is given by the differential:

δJ = ∇J (ℓ) .δv (ℓ) = -ρ ℓ ∇J (ℓ) 2 (67)
Hence, the step-size expected to diminish J (ℓ) of the amount δJ = -εJ (ℓ) is estimated to be:

ρ ℓ = εJ (ℓ) ∇J (ℓ) 2 (68)
In particular for ε = 1, we get the quasi-Newton method since the employed discrete gradient is only approximately equal to the gradient of the discrete J.

Multiple-gradient descent algorithm (MGDA )

In this subsection, we propose an alternative coordination algorithm in which the matching of the sub-solutions is treated as a multi-objective optimization problem, considering that all defectfunctionals J i 's should be driven to 0 concurrently.

In the Multiple Gradient Descent Algorithm MGDA (see [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF] for a detailed definition and convergence proof), once the individual discrete gradients,

u i = ∇J i (i = 1, ..., 4) (u i ∈ R N ) (69)
are known, one considers the convex hull U of these vectors, and identifies its minimum-norm element ω:

       ω = Argmin u∈U u 2 U = u ∈ R N / u = 4 i=1 α i u i ; α i ≥ 0 (∀i) ; 4 i=1 α i = 1 (70)
In our problem, the dimension N is the number of nodal controls:

N = 2(N X + N Y -2) (71)
In Appendix A, a special parameterization of the convex hull is proposed to facilitate the determination of the element ω by a numerical optimization procedure.

Then, once the element ω is determined, if ω = 01 , the current iterate is, or is treated as Pareto stationary. But here, the Pareto front is made of only one point corresponding to J i = 0 for all i. This situation corresponds to full convergence of the coordination algorithm. Otherwise (ω = 0), -ω is a descent direction for all criteria simultaneously. Thus, (65) is replaced by:

v (ℓ+1) = v (ℓ) -ρ ℓ ω (ℓ) (72)
Here again, we propose to adjust the step-size ρ ℓ according to (68). However, it is not clear here that the proper scaling corresponds to ε ∼ 1.

6 Numerical experimentation

Test-case

Let a and b be two adjustable constants, and:

                   ψ = (x + a) 2 + (y + b) 2 (-1 ≤ x ≤ 1) (-1 ≤ y ≤ 1) τ = 1 ln(a 2 + b 2 ) φ = τ ln ψ f (x) = 1 -x 2 f (y) = 1 -y 2 u e = f (x) f (y) φ (73)
As a result, φ is a harmonic function:

∆φ = 0 (74)
and this permits us to simplify somewhat the expression of the Laplacian of u e :

∆u e = ∆ f (x) f (y) φ + 2∇ f (x) f (y) .∇φ + f (x) f (y) ∆φ = -2 f (x) f (y) φ + 2∇ f (x) f (y) .∇φ (75) = -f e (76)
where

f e = 2 f (x) f (y) φ + 8τ ψ x(x + a)f (y) + y(y + b)f (x) (77)
Hence, for f = f e , the exact solution of the continuous problem is u = u e .

The constants a and b have been introduced to destroy the symmetry in the solution. More specifically, the following settings were made: a = 5 4 and b = 3 4 . The corresponding problem has been discretized and solved using either one domain to establish a reference, or four to experiment multi-criterion optimization algorithms.

The single-domain discrete solution u h is depicted in Fig. 2 as a surface in 3D, and the corresponding contour map is given more precisely in Fig. 3.

Quasi-Newton steepest descent

In a first series of numerical experiments, the steepest-descent method was applied to drive to 0 the global criterion J. After a few trials, it appeared that best convergence was achieved by setting ε to 1 in (68), which corresponds to the quasi-Newton method.

Two experiments are reported presently. They differ in the setting of the initial interface conditions. In the first, the v i 's are initially assigned the restriction to the corresponding interface of the exact solution u e , which differs from the discrete solution by truncation errors. In this case, iterative errors are initially very small, which permits the asymptotic convergence to be assessed. In the second experiment, the controls are initially set to 0 order to assess the global convergence.

Asymptotic convergence. The convergence history of the global criterion J as well as its individual parts, {J i } (i=1,...,4) is represented in Fig. 4. The criterion J, in 20 iterations, goes from 2.9 × 10 -2 to a level below 10 -4 . Note that J 4 is somewhat smaller in magnitude and more subject to oscillations.

Global convergence. The convergence history of the global criterion J as well as its individual parts, {J i } (i=1,...,4) is represented in Fig. 5. The criterion J, in 200 iterations, is reduced by 8 orders of magnitude. The different criteria, apart from small oscillations, converge at essentially the same rate. In a linear convergence process, this rate is imposed by the most persistent mode, present in all criteria when the initial condition is arbitrary.

Convergence of the gradients. The evolution of the four gradients {∂J/∂v i } (i=1,...,4) over 200 iterations is given on Fig. 6-Fig. 9. They appear as high-frequency modes. Each one vanishes asymptotically over the corresponding interface.

Discrete solution. The four-domain discrete solution is found perfectly smooth, in fact even smoother than the single-domain discrete solution. This is due to a higher degree of iterative convergence. In the experiments of this section, at each iteration, the 3 parameters c 1 , c 2 and c 3 of (101) have been discretized uniformly by step of 0.01, and ω was set equal to the vector of minimum norm among the essentially 10 6 corresponding candidates. The resulting vector was therefore a very coarse approximation of the actual vector ω.

Asymptotic convergence For this experiment, the discretized continuous solution is again imposed initially at the interfaces. The convergence history of the above basic algorithm is indicated in Fig. 12 for 20 iterations. After some initial adjustment, the trend is towards decaying, but at a deceiving slow rate. In an attempt to explain this poor convergence, the following observation was made: suppose the gradients of the individual criteria, {∂J i /∂v} (i=1,...,4) , are very different in magnitude. Remember that in a linear iterative process, unless initial conditions are very special, all quantities converge at the same rate, say Cρ (iter) , where ρ is the spectral radius, and C a constant which depends on the quantity considered. For example, in the previous experiment, J 4 itself was observed to be somewhat smaller than the other criteria, and so was its gradient. Then, the convex-hull minimum-norm element ω is paradoxically dominated by the gradient of smallest magnitude, since in the convex combination of the gradients, putting the largest weight on the smallest has the effect of reducing the norm of the combination. But this is not efficient, since this gradient corresponds to the already small criterion for which minimization is the least necessary. This observation has led us to calculate the direction ω as the minimum-norm element in the convex hull of normalized gradients. Here, the normalization was made by scaling each gradient to the corresponding value of the individual criterion. In other words, logarithmic gradients were considered. In fact, that is exactly what Newton's method does with the global criterion J.

The above experiment was then repeated, using logarithmic gradients to determine the vector ω. The corresponding convergence history is indicated in Fig. 13. This new result is now found very similar to the analogous result previously achieved by the quasi-Newton method (Fig. 4). The importance of scaling is therefore confirmed. Global convergence. All interface controls v i 's are initially set to 0. The resulting convergence history over 200 iterations is indicated in Fig. 14 for the basic algorithm, and in Fig. 15 for the scaled algorithm based on logarithmic gradients. The first algorithm seems to convergence, but at a very slow rate. The second seems to subject to accidents and to experience a difficulty to enter the asymptotic convergence phase. The criteria stagnate. This may be caused by many factors on which future investigation will focus:

• the insufficiently accurate determination of ω;

• the non-optimality of the scaling of gradients;

• the non-optimality of the step-size, the parameter ε in (68) being maintained equal to 1 throughout;

• the large dimension of the design space, here 76 (4 interfaces associated with 19 d.o.f.'s).

S ′ i -1 J ′ i , ω = λ 2 (∀i) (88)
that is, the same positive constant [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF].

As a result of this direct and fast construction, the vector ω is usually different from the former definition, except in special cases, as for example, when n = 2 and the angle between the two gradients is obtuse. Nevertheless, the new ω also provides a descent direction common to all criteria, scaled essentially as initially prescribed.

Automatic rescaling: MGDA II b. The examination of the case n = 2 has led us to propose a slightly different handling of the scales. Here, one lets

A i = S i - i-1 k=1 c i,k (89)
only this number is strictly-positive. Otherwise, the scale S i is redefined ("automatic rescaling") as follows:

S i = i-1 k=1 c i,k (90)
and one sets:

A i = ε i S i .
The result in (88) is still valid, and it now provides an information on gradients that have been weighted as prescribed whenever S i > i-1 k=1 c i,k , and otherwise by the procedure itself. This rescaling procedure is certainly perfectible.

Convergence experiments and discussion. MGDA II has been tested on the partitioning problem in the four possible options corresponding to "no scaling prescribed" vs "prescribed scaling", and "automatic rescale off" vs "on". In MGDA II b, when k<i c i,k was found greater or equal the prescribed S i (=1 or J i ), ε i was set to 0.01 (and maintained to 0 otherwise).

A first observation was made: the new procedure for determining ω is much faster, and MGDA II seems to be less sensitive to round-off errors.

In Fig. 16 and Fig. 17, the automatic rescale is off, and the effect of scaling alone is evaluated. Over the first 200 iterations, the result is about the same. However the scaled version indicates a trend to convergence acceleration to be confirmed.

In Fig. 18 and Fig. 19, the automatic rescale is on, and the option of prescribed scaling is off/on. Again a better convergence is achieved when scales are prescribed.

In order to confirm these results, the best option "prescribed scales and automatic rescale" is compared with the basic method in Fig. 20 and Fig. 21 over 500 iterations. The trends indicate a linear convergence for the first method, and a seemingly-quadratic convergence for the second. Compared to the quasi-Newton method of Fig. 5, MGDA II b is grossly-speaking twice slower, but it indicates a more definite trend to asymptotic convergence acceleration.

One last remark: in these experiments, we observe that scaling has the effect of making the convergence curves associated with the different criteria closer to one another. In this report, various versions of the Multiple-Gradient Descent Algorithm (MGDA ) have been tested numerically over a domain-partitioning problem treated as a multi-objective problem in which matching defect integrals at the different interfaces are to be minimized concurrently, and in fact, all driven to 0.

The major objective of this experimentation was to assess the potential of the MGDA to handle multi-objective problems in which the finite-dimensional setting was the result of discretization, thus approaching a more general functional setting. In this respect, the demonstration was made. Indeed convergence was achieved by MGDA . However the quasi-Newton method applied to the agglomerated criterion was globally found more efficient, but in the most sophisticated version (MGDA II b), the algorithm seems to demonstrate a promising asymptotically-quadratic convergence.

Thus, if the convergence was not always found satisfactory, several observations should temper this conclusion, and many promising directions of improvement can be envisaged:

• in the problem under study, the Pareto set was reduced to the single point corresponding to all criteria equal to 0 associated with the unique solution of the discretized Poisson problem; this situation is really atypical of standard multi-objective problems; additionally, the criteria to be minimized were not really antagonistic, since they all converged at almost the same rate with the quasi-Newton method, leaving little possibility of improvement from the start; for these two reasons MGDA has been tested in a very straining situation for which it was not devised originally;

• the large dimension of the design space, here 76 (4 interfaces associated with 19 d.o.f.'s), was probably a handicap;

• a robust procedure to define the step-size should be devised; in our experiments, the parameter ε was maintained equal to 1 throughout;

• the determination of ω in the basic method should be made more accurately;

• the scaling of gradients was found important; alternatives to the logarithmic gradients should be analyzed and rationalized; more generally, preconditioning remains an open question;

• the MGDA II variant was found faster and more robust;

• at present, our most sophisticated algorithm, MGDA II b, also involves an automatic rescaling procedure; it indicates a definite trend to asymptotic convergence acceleration (quadratic convergence). Usually, but not necessarily : N ≥ n.

A Generalities about convex hulls and parameterization

Parameterization -The convex hull may parameterized by identifying the set of allowable coefficients {α i } (i=1,...,n) . To satisfy the positivity condition automatically, one lets:

α i = σ 2 i (i = 1, ..., n) (92)
Then, the equality constraint

n i=1 α i = n i=1 σ 2 i = 1 (93)
states that

σ = (σ 1 , σ 2 , ..., σ n ) ∈ S n (94)
where S n is the unit sphere of R n , and precisely not R N . This sphere is easily parameterized using trigonometric functions of n -1 independent arcs φ 1 , φ 2 , ..., φ n-1 :

                 σ 1 =
that is:

σ i = sin φ i-1 . n-1 j=i cos φ j (96)
with φ 0 = π 2 . It is sufficient to consider the portion of the sphere corresponding to φ i ∈ [0, π 2 ] for all i ≥ 1 since the sign of the σ i 's makes no difference.

The usage of trigonometric functions is not really necessary, since one can let:

c i = cos 2 φ i (i = 1, ..., n) (97)
and get:

                           α 1 = c 1 .

Figure 1 :

 1 Figure 1: Partition of a square in sub-domains {Ω i } (i=1,...,4) to solve the Dirichlet problem A first calculation of the compound solution

Figure 2 :Figure 3 :

 23 Figure 2: Single-domain discrete solution u h

Figure 4 :Figure 5 :Figure 6 : 1 -Figure 7 : 2 InriaFigure 8 : 3 -Figure 9 :Figure 11 :

 45617283911 Figure 4: Quasi-Newton steepest descent -convergence history of criteria (discretized continuous solution imposed initially at interfaces)

Figure 12 :

 12 Figure 12: Basic MGDA -asymptotic convergence history of criteria

Figure 13 :

 13 Figure 13: MGDA based on logarithmic gradients -asymptotic convergence history of criteria

Figure 15 :

 15 Figure 14: Basic MGDA : global convergence

Figure 16 :Figure 17 :Figure 18 :Figure 19 :Figure 20 :Figure 21 :

 161718192021 Figure 16: MGDA II , no scaling prescribed

α

  i u i ; {u i } (i=1,...,n) : given family of vectors in R N constraints: ∀i : α i ≥ 0 ;

c 2 . c 3 1 . ( 1 -

 2311 . ... . c n-1 α 2 = (1 -c 1 ) . c 2 . c 3 . ... . c n-1 α 3 = 1 . (1 -c 2 ) . c 3 . c n-1 )

  cos φ 1 . cos φ 2 . cos φ 3 . ... . cos φ n-1 σ 2 = sin φ 1 . cos φ 2 . cos φ 3 .

								...	. cos φ n-1
	σ 3 = . . .	1	. sin φ 2 . cos φ 3 .	... . . .	. cos φ n-1
	σ n-1 =	1	.	1	.	...	. sin φ n-2 . cos φ n-1
	σ n =	1	.	1	.	...	.	1	. sin φ n-1

RR n°7968

In the numerical implementation, this condition is relaxed to be: ||ω|| < T OL, for a given tolerance T OL.RR n°7968

Inria

In this section, J ′ i denotes the gradient of J i , and no longer the first variation.

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

MGDA II

Recently, a variant, MGDA II , has been proposed in which the descent direction is calculated by a direct procedure, which provides a valuable simplification of implementation [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF]. This new algorithm is now presented again along with a new variant (MGDA II b), and tested on the DDM problem.

Basic definition, scaling. In MGDA II , the possibility to prescribe scales for the gradients, {S i } (i=1,...,n) (S i > 0 (∀i)) is offered. In the following experiments, at a given iteration, these scales are either set all equal to 1 ("no scaling prescribed"), or equal to the current values of the criteria ("prescribed scaling"):

the latter implies that the descent direction is based on logarithmic gradients.

Assuming that the gradients form a linearly-independent family, an assumption never contradicted in the numerical experiments, a family of orthogonal, but usually not orthonormal vectors are formed {u i } (i=1,...,n) according to the following:

where A 1 = S 1 , and, for i = 2, 3, ..., n:

where:

and

for some arbitrary, but small ε i (0

The minimum-norm element in the convex of the family {u i } (i=1,...,n) is given by:

and one finds:

which confirms that ω does belong to the interior of the convex hull, so that:

and:

Convening that ε i = 0 in the regular case (S i = k<i c i,k ), and otherwise by modifying slightly the definition of the scaling factor according to

that is:

with c 0 = 0, and c i ∈ [0, 1] for all i ≥ 1.

In this way, the constraints on the coefficients {α i } have been replaced by the bounds 0 and 1 on the new parameters {c i }, independently of one another. However, the criterion to be minimized,

is now a polynomial of possibly large degree, namely 2(n -1), of the new parameters {c i }.

In the particular case of the coordination of 4 sub-domains by MGDA , and independently of the degree of refinement of the spatial discretization controlled by the integers N X and N Y , n = 4, and once the 10 scalar products (u i , u j ) (i, j = 1, ..., 4) calculated, the determination of the minimum-norm element ω is equivalent to minimizing a 6th-degree polynomial of (c 1 , c 2 , c 3 ) in the R 3 unit cube (limits included):


Fatima Nasser

Zhongyang Li

Nadine Martin

Michelle Vieira

Philippe Gueguen

Seismic Response Analysis of Different Buildings using Time-Invariant and Time-Variant Damping Coefficients

Seismic signals are characterized by strong excitations, short durations, non-linearity and non-stationarity having both the amplitude and frequency content vary as a function of the recorded time. Various classical detection and estimation techniques, like the time-frequency representations and the Fourier-based techniques have been used to analyze such signals, but these techniques have variety of limitations and they fail to correctly estimate the concerned signals. The Damped-Amplitude and Polynomial-Frequency Model has been introduced to help in adapting to the seismic signals where the amplitude is damped. This model is based on approximating the frequency by loworder polynomials and the amplitude by damped exponentials. Its amplitude in turn is characterized by a damping coefficient; which was firstly assumed to be time-invariant. However the results of the studied signals showed rapid amplitude fluctuations and frequency content variations of each component that could be justified by the fact that the dynamic response of the structure is highly sensitive to the dynamic characteristics of the ground motion. Accordingly and to be more adapted with the physical model of the building motion that is characterized by damped exponential functions, the damping coefficient was then assumed to be time-variant leading to the foundation of a new model that keeps the same approximation for the frequency like the aforementioned model, and changes that of the amplitude by approximating its damping-coefficient by low-order polynomials. This model was then named Polynomial Damping Function Model. Results on different seismic signals show that the time-variant assumption is more efficient than the time-invariant one.

Introduction

Seismic signals are non-stationary; having both the amplitude and frequency content vary considerably as a function of recording time. Consequently, for a more reliable representation of such signals, both the amplitude and frequency variations of the recorded time histories should be accounted for. Recently, a large number of methods have been proposed in this context. For example, Reine et al. (2009) [START_REF] Reine | The robustness of seismic attenuation measurements using fixed-and variable-window time-frequency transforms[END_REF] compared timefrequency and time-scale methods in seismic data processing and interpretation, the Stransform method was also applied in seismic data analysis by Pinnegar and Mansinha (2003) [START_REF] Pinnegar | The s-transform with windows of arbitrary and varying shape[END_REF] , and the reassignment method was proposed for the small frequency variations in civil engineering structures under weak and strong motions by Michel and Gueguen (2009) [START_REF] Michel | Time-Frequency Analysis of Small Frequency Variations in Civil Engineering Structures Under Weak and Strong Motions Using a Reassignment Method[END_REF] . In the context of the signals under study in this paper, these techniques present certain drawbacks. The analysis of (1, 2, and 3) is based on the time-frequency representation which is limited in time-frequency resolution by the choice of a window length. Particular drawback in [START_REF] Michel | Time-Frequency Analysis of Small Frequency Variations in Civil Engineering Structures Under Weak and Strong Motions Using a Reassignment Method[END_REF] is that the reassignment of the energy yields wrong amplitude estimations. Accordingly, and in order to approach the dynamic evolution of the modulations, a high time resolution is required, and to have better performance in multi-component case, a high frequency resolution is needed.

Generally, to get away from the constraints and drawbacks of such non-parametric classical methods, and to correctly reconstruct the signals under study, a model should be set. This model should be as general as possible to be applicable in different domains.

In previous papers, we already proposed a polynomial frequency model with polynomial amplitude in [START_REF] Jabloun | Estimation of the amplitude and the frequency of nonstationary short-time signals[END_REF][START_REF] Jabloun | Local orthonormal polynomial decomposition for both instantaneous amplitude and frequency of highly nonstationary discrete signals[END_REF][START_REF] Jabloun | A AM/FM Single Component Signal Reconstruction using a Nonsequential Time Segmentation and Polynomial Modeling[END_REF] , or with damped amplitude in [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] . In this paper, a new model named polynomial damping function model is proposed for the amplitude modulation. This model is of great interest since it is still able to deal with multicomponent signals of very short durations, and with the rapid non-linear amplitude and frequency variations of each component of a signal. Moreover, compared to the damped amplitude model in [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] , the damping coefficient of its amplitude is assumed to be time variant, so that it is less constrained and better convenient for damping factor estimation of seismic signals.

The model proposed in this paper needs a parameter estimation method, and a procedure to solve the non-linearity problem of the likelihood function similar to those of [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] , which are the Maximum Likelihood Estimation (MLE) and the Adaptive Simulated Annealing (ASA) techniques respectively. The objective of this paper is then to apply the method proposed herein to estimate the amplitude and frequency modulations (AM/FM) of seismic vibration signals of buildings, recorded by different accelerometers placed on the top of the buildings.

The remainder of the paper is structured as follows. In section 2, the physical model for the seismic vibration signals of interest is formulated. The signal model for the AM/FM estimations and the parameter estimation method are introduced in section 3. The obtained results over the analysis of simulated signals computed from the physical model and real-world seismic data are discussed in section 4. Finally, the paper ends with final conclusions in section 5.

Physical model of the building motion

In [START_REF] Michel | Seismic risk scenario in Grenoble (France) using experimental dynamic properties of buildings[END_REF] the dynamic behavior of a building is deduced from the matching of the building with a continuous Timoshenko beam. In this model, the mass of each storey of the building is considered to be mostly concentrated at its floor, and a lumped mass modeling is assumed for this structure. Therefore, the Duhamel integral gives the elastic motion of the building at each floor ( ) fl U t by only knowing the mass of each floor, the modal parameters (namely the mode shapes , frequencies and damping ratios), and the ground motion ( ) P t . The mode shape of the th fl floor [ ] k fl φ depends on the number of the floor fl and it participates as the weighting factor to generate the final response of the structure, as will be shown later in this section. The impulse response of a single degree of freedom system for a k mode is

( ) [ ] mod 1 ( ) sin , 1, , -ξ ω - = ω ∀ ∈ ω k k t k Dk e Dk h t e t k N (1) with 2 2 1 ω = π -ξ Dk k k f
the pseudo-pulsation of mode k . This system being excited by the ground solicitations, the time-varying response of the mode k is obtained from the convolution with the ground solicitation ( )

P
- = - = * ∫ (2)
with k p the participation factor that is defined as In fact, a building is more complex than a single degree of freedom system, so that the motion of the building is then modelized as a linear superposition of the motion of all the mod e N modes. Therefore, the dynamic response for a given floor with

[ ] [ ] 1
[ ] mod 1 ( ) ( ) e N fl k k k U t y t fl φ = = ∑ . (3)
From this point, it is worthwhile mentioning that the non-stationary dynamic response of the structure is highly sensitive to the dynamic characteristics of the ground motions.

In the next two sections, two cases will be studied, the ambient ground excitations and the seismic ground excitations.

Accordingly, as an initial step to clarify this physical model, and to introduce the next two sections, let's consider one example with the continuous physical model of a commercial building, located in Sherman-Oaks, California, assuming that the mass of each floor is constant, and the number of floors being 13 (see Fig. 1 left). From the spectrogram and the Fourier spectrum of the data recorded at the top of this building when excited by the real-world earthquake of Northridge January 17, 1994 ( 6.

4

L M = )
(see Fig. 1 middle and right respectively), the number of modes mod e N could be clearly seen as 3 located at 0.4, 1.2 and 2.2 Hz respectively. Thus, knowledge of the first two modal frequencies 1 f =0.4 Hz and 2 f =1.2 Hz, necessary in the Timoshenko beam model [START_REF] Michel | In situ experiment and modelling of RC-structure using ambient vibration and Timoshenko beam[END_REF] , allows evaluating the third modal frequency at 2.2 Hz which is perfectly matched with what is presented by the spectrogram and the Fourier spectrum of this signal. For each mode k , the damping ratio k ξ is assumed to be fixed at 0.05, accordingly the pseudo pulsations Dk ω of this building will be

Under Ambient ground excitations

Ambient vibrations in buildings are produced by the wind (low frequencies < 1Hz), internal sources (machinery, lift at high frequencies) and seismic noise (broadband).

The seismic noise being prevailing in case of ambient ground excitations, the input in model ( 2) is then assumed to be a stationary white noise with zero mean and constant variance.

For the example of Sherman-Oaks building of Fig. 1, when excited by such a ground excitations, the result of Eq. ( 2) for simulating modal responses ( ) k y t for the 3 modes, and the result of Eq. ( 3) for simulating the full dynamic response will be as shown in Fig ( 3).

Under Seismic ground excitations

Seismic ground excitations are more difficult to be modeled, and their spectral shape is no longer white. In this paper, these seismic excitations will be a real-world ones recorded at the ground level of the Sherman-Oaks building of Fig. 1 during the earthquake of Northridge January 17, 1994 ( 6.4

L M =

). For the most energetic part of this signal between 2 and 12 s, as shown by the spectrogram in Fig. 4 left, the spectral shape is no longer white. The log-log presentation of the spectrum of this energetic part, (see Fig. 4 right), corresponds actually to the Brune model [START_REF] Brune | Tectonic Stress and the Spectra of Seismic Shear Waves[END_REF] with a flat frequency band between 0.6 c f = Hz and max 5 f = Hz.

Such an input in Eq. ( 2) provides the modal response ( ) k y t of each mode and the full dynamic response of the Sherman-Oaks building to the earthquake considered (see Fig. 5). Generally, different types of excitations will affect the building response differently. Whereas for the ambient ground excitation, all the modal responses in Fig. 3 (1 st , 2 nd and 3 rd lines) are affected in the same way by ground motions during the whole recorded time series, where for seismic ground excitations case (see Fig. 5), the response increases in presence of the seism between 2-12s, then return to the amplitude variations due to the ambient ground excitations.

Comparing the two spectrograms, (Fig. 3 (5 th line) and Fig. 5 (5 th line)) show that the mode at 0.4 Hz is more energetic than the other two at 1.2 and 2.2 Hz. These last two modes have also been modified by the earthquakes. Consider the energy of the mode k k E calculating in the most energetic part between 2-8s as follows ( ) ( )

8 2 2 s k k t s E y t = = ∑ [ ] mod 1, e k N ∀ ∈ (4)
where ( ) k y t is the modal response of the mode k . Table 1 shows an energy calculation of all the modes for both ambient and seismic excitations during the seismic period between 2-8s.

[ ] [ ] [ ] y n s n e n = +

where [ ] [ ] [ ]

1 i K j n i i s n A n e Φ = = ∑ . ( 5
)
The time reference is set to the center of the window as 2 2

N N n -≤ ≤ , K is the number of components, j is the complex number verifying [START_REF] Pinnegar | The s-transform with windows of arbitrary and varying shape[END_REF] 1

j = -, [ ] i
A n is the time-varying amplitude, and

[ ] i n Φ
is the phase of the th i component, obtained by numerical integration of the Instantaneous Frequency IF, and centered in the middle of the observation window in order to minimize the estimation error [START_REF] Jabloun | Estimation of the amplitude and the frequency of nonstationary short-time signals[END_REF] , thus making [ ] ,0

0 i i ϕ = Φ , [ ] [ ] [ ] ( ) 0 ,0 / 2 / 2 2 n i i i i k N k N n F k F k ϕ π =- =- Φ = + - ∑ ∑ (6)
To assure the uniqueness of the model, the Instantaneous Amplitude (IA) should be strictly positive [ ]

0 i A n > ,
[ ] 0 / 2 i s F n F < <
where s F is the sampling frequency [START_REF] Jabloun | Local orthonormal polynomial decomposition for both instantaneous amplitude and frequency of highly nonstationary discrete signals[END_REF][START_REF] Jabloun | A AM/FM Single Component Signal Reconstruction using a Nonsequential Time Segmentation and Polynomial Modeling[END_REF] .

Damped amplitude and polynomial frequency model

For the frequency modulation model, like in [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] , the IF is approximated by low order polynomials, in order to track the strong local variations, In this paper, the amplitude modulations are approximated by two types of damped functions, one with time-invariant damping coefficient like in [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] and another, which is the proposition of this paper, with time-variant polynomial damping coefficient. Table 2 presents these two studied signal models and their characteristics.

[ ] [ ] , 0 f M i m i m m F n f g n = = ∑ , (7)

Time-invariant damping coefficient

In this section, the damping coefficient of the structure is assumed to be time-invariant, so that

[ ] i n i i A n e α β - = , (8)
where i β and i α are the initial amplitude and the damping coefficient respectively, characterizing the amplitude of the th i component of this model. With regard to the realworld data, i β is always constrained to be positive [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] .

In summary, the P parameters to be estimated are ,0 ,0 , , , , , ,

f i i i i i i M f f θ β α ϕ   =   K with 4 f P M = + .

Time-variant polynomial damping coefficient

The real seismic signals are characterized by rapid amplitude fluctuations and frequency variations of the energy contents of each of their components, from this point the damping coefficient was assumed to be time-variant and thus a new model known as "Polynomial Damping Function Model" is proposed to approximate the amplitude fluctuations. This model is more adapted to the civil engineering buildings that are characterized by the exponentially damping functions, and thus it is in well accordance with the physical model presented in section 2.

For this model, the amplitude is defined as follows Accordingly, the P parameters to be estimated are ,0 , ,

[ ] [ ] [ ] ( ) 0 / 2 / 2 n i i k N k N k k i i A n e α α β =- =- - - ∑ ∑ = where [ ] [ ] , 0 M i m i m m n g n α α ρ = = ∑ , ( 9
i i i M i i i M f f α θ β ρ ρ ϕ   =   K K with 4 f P M M α = + + . 0 ,0 , , , , , , , , f i
[ ] i n i i A n e α β - = [ ] [ ] [ ] ( ) 0 / 2 / 2 n i i k N k N k k i i A n e α α β =- =- - - ∑ ∑ = (IF) [ ] [ ] , 0 f M i m i m m F n f g n = = ∑ [ ] [ ] , 0 f M i m i m m F n f g n = = ∑ ( ) i α i α constant [ ] [ ] , 0 M i m i m m n g n α α ρ = = ∑ ( ) P 4 f M + 4 f M M α + + ( ) i θ ,0 ,0 ,
, , , , ,

f i i i i i i M f f θ β α ϕ   =   K ,0 , , 0 ,0 , , , , , , , ,
f i i i i M i i i M f f α θ β ρ ρ ϕ   =   K K

Parameter estimation method

A maximum-likelihood approach as in [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] is used to estimate the model parameter vector i θ . As the noise is assumed to be a white Gaussian process; the MLE is equivalent to minimization of a least Square approach, ( )

ˆarg min K P LS R l θ θ θ × ∈ = , with ( ) [ ] [ ] 2 / 2 / 2 N LS n N l y n s n θ =- = - ∑ , (10)
with [ ] y n the noisy observation, and [ ] ŝ n the signal model computed by substituting

[ ] i A n , [ ] i F n and [ ] i n Φ in [ ] s n .
Analytical solutions are not applicable on (10) due to its non-linearity, the same stochastic optimization technique as in [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] is then applied. This technique speeds up the convergence process to the global optima and induces a gain in the computing time. The detailed steps of this technique can be found in (13) .

Applications on seismic signals

Firstly in this section, the aforementioned two signal models of section 3 will be applied over the real-world seismic data recorded at the top of three different buildings:

-1 st building: The 7-storey hotel building located in Van Nuys, California, Fig. 6 (left), which was severely damaged by the Northridge earthquake of January17, 1994, ( 6.4

L M =

). The signal is of 60 seconds (3000 samples) length in total, sampled at 50Hz and decimated at 5Hz. -2 nd building: The 13-storey governmental office San José building located in the Bay area of California, San Francisco, Fig. 6 (right) which was severely damaged as well by the earthquake of Loma Prieta of October 17, 1989 ( 6.9

L M =

). The signal is of 120 seconds (6000 samples) length in total, sampled at 50Hz and decimated at 10Hz.

-3 rd building: The 13-storey commercial building, located in Sherman-Oaks, California, already presented in Sec. 2 and Fig. 1. This building has been damaged by the earthquake of Northridge January 17, 1994 ( 6.

4

L M = ),
The signal is of 60 seconds (3000 samples) length in total, sampled at 50Hz and decimated at 10Hz. The two signal models considered will be evaluated over these real-world seismic data and compared all together in the aim of highlighting the interest of the model proposed in this paper when dealing with such data. Then at the end of this section, the proposed model will be applied over the simulated seismic signals of section 2.

Analysis of real-world seismic data

For each of the three real signals presented in the previous section, an appropriate segment of 6 s length is studied. For the case of the time-invariant damping coefficient, the frequency is estimated at the 2 nd order. While for the case of time-variant damping coefficient, the polynomial damping function is estimated at the 3 rd order.

Time-invariant and time-variant damping coefficient

In [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] , the initial amplitude i β and the damping coefficient i α are constrained to be strictly positive (Sec. 4.1). Fig. 7 presents the results of the damped amplitude and polynomial frequency model, with i α constrained to be strictly positive, when studied over the three real-world seismic data of the buildings presented in section 4. This model has shown its efficiency when applied over ambient vibration signals [START_REF] Li | Damped-amplitude & polynomial frequency model for short-time ambient vibrations of a building[END_REF][START_REF] Li | Parameter estimation of shorttime multi-component signals using damped-amplitude & polynomial-frequency model[END_REF] . Our aim in this paper is mainly to test this model over real-world seismic data, which are characterized by more intense vibrations. The results reported in Fig. 7 proved that this model is no longer reliable when dealing with seismic data. The amplitude modulations of the three chosen studied signals are not correlated with the spectrogram, which indicates that the local variations of the concerned signals weren't tracked correctly using this model. From this point, and in the aim of improving the performance of this model, i α was then assumed to be without constraints. The results of this hypothesis, studied but not presented here, didn't improve the performance at all, neither the normalized errors were minimized, nor the amplitude variations were in agreement with the spectrogram. Fig. 8 presents the results of the polynomial damping function model, when studied over the three real-world seismic data of the buildings presented in section 4. The improved performance could be clearly noticed, the amplitude variations are in agreement with the spectrogram, the frequency variations are acceptable as well, and the normalized error has been decreased noticeably.

It has been shown that the amplitude and frequency content of the seismic data are well extracted when applying the polynomial damping function model; this model has proven to be a good choice when dealing with seismic data.

Analysis using simulated seismic signals

In this section, the polynomial damping function model is applied over the range 3-9 s of the simulated seismic signal presented in section 2, where the seismic ground excitations correspond to those in Fig. 4. , which justifies why this mode is less influenced by the seism. Accordingly, and in order to highlight the effect of the seism on the simulated signals as well as the real ones, we applied a high-pass filter with 0.6 c f Hz = to eliminate the component at 0.4Hz . From Fig. 9, we can conclude that the two components and their variations are correctly identified for both cases when applying the polynomial damping function model. Meanwhile, this model performs better with the simulated seismic signal than with the real-world seismic one, indeed the normalized root mean square error of the simulated signal is about two times smaller than the seismic data.

The results of the real-world seismic data could be justified because for such data the modal parameters (such as natural frequency, damping ratio, and mode shapes) are always time-variant, however, for the simulated seismic signal the modal parameters are constant and are a priori known. Nevertheless, the polynomial damping function model has proven its efficiency even when dealing with the physical modeling of the building motion.

Conclusions

Unlike input-output techniques for which building response is evaluated by knowing both the input signal (for example, seismic ground motion) and the output signal corresponding to building motion (for example, motions recorded on the roof), we propose in this paper an Output-Only Modal Analysis method, considering the input as unknown. The model proposed in this paper is the polynomial damping function model based on the time-variant damping coefficient and applied in the context of seismic vibrations, which are characterized by non-linearity, non-stationarity and shortdurations.

The proposed model has led to significantly improved estimations of both the amplitude and the frequency modulations as compared to the damped amplitude & polynomial frequency model, the variations of the multi-component seismic signals were tracked correctly over the very short duration of 6 seconds, and the normalized error for each of the studied segments for the concerned signals was noticeably low. Furthermore, this model is more adapted to the civil engineering buildings that are characterized by exponential damping functions, and it is in well accordance with the physical model of the building motion.

In future we intend to study the applicability of this Output-Only signal model in extracting the dynamic properties of the buildings, such as the natural frequency and the damping ratios.

Figure 1 :

 1 Figure 1: The 13-storey commercial building located in Sherman-Oaks, California (left), the spectrogram calculated with Hann window of 12.8 s (64 points) (middle), and the Fourier spectrum (right), of the data recorded at the top of the Sherman-Oaks building when excited by the realworld earthquake of Northridge January 17, 1994 ( L M = 6.4 ) Figs. 2 (left and right) present respectively the mode shapes of the 13 th floors and the impulse response of the 3 modes of the Sherman-Oaks building.

Figure 2 :

 2 Figure 2: The mode shapes of the 13 th floor (left) and the impulse response of the 3 modes (right) of the Sherman-Oaks building

Figure 3 :

 3 Figure 3: Dynamic response model at the top of Sherman-Oaks building when excited by white Gaussian ground excitations: modal responses k y (t) of the 3 modes, mode k=1, 2 & 3 (1 st line top, 2 nd line, 3 rd line respectively), the dynamic response of the 13 th floor U fl (t) (4 th line), and its spectrogram (Hann window of 51.2 s-256 points) (5 th line)

Figure 4 :Figure 5 :

 45 Figure 4: The spectrogram calculated with Hanning window of 3.2 s (16 points) (left), and the loglog presentation of the Fourier spectrum around the most energetic part of the seismic data [2-12s] corresponds to Brune Model (right) of the seismic data recorded at the ground level of the Sherman-Oaks building

  ) with i β being the initial amplitude, M α is the approximation order of the damping coefficient, [ ] m g n is an orthonormal polynomial of order m, and , m i ρ are the amplitude parameters.

Figure 6 : 7 -

 67 Figure 6: 7-storey hotel building located in Van-Nuys, California (top left), 13-storey San-José building located in the Bay area, California (top right), and the spectrograms of the real-world seismic data recorded at the top of the Van-Nuys and San José buildings calculated with Hann window of 3.2 s (32 points) and 6.4 s (32 points) (bottom left and bottom right respectively)

Figure 7 :

 7 Figure 7: Analysis of the real-world seismic data recorded at the top of Van-Nuys hotel (1 st column left), San-José building (2 nd column middle) and Sherman-Oaks building (3 rd column right) when using the damped amplitude and polynomial frequency model ( i α constant): The spectrograms calculated with Hann window of 15 points (1 st line top), the estimated amplitude modulations (2 nd line middle), and the estimated frequency modulations (3 rd line bottom)

  Based on these conclusions and contrary to the damped amplitude and polynomial frequency model that is concentrated on the time-invariant damping-coefficient, in this paper, a new model namely the polynomial damping function model that accepts i α to be time-variant is proposed (c.f. section 3.1.2).

Figure 8 :

 8 Figure 8: Analysis of the real-world seismic data recorded at the top of Van-Nuys hotel (1 st column left), San-José building (2 nd column middle) and Sherman-Oaks building (3 rd column right) when using the polynomial damping function model ( i α variable): The spectrograms calculated with Hann window of 15 points (1 st line top), the estimated amplitude modulations (2 nd line middle), and the estimated frequency modulations (3 rd line bottom)

Figure 9 :

 9 Figure 9: Analysis of the simulated seismic data recorded at the top of the Sherman-Oaks building when excited by seismic ground excitations (1 st column) and the real-world seismic data recorded at the top of Sherman-Oaks building during the Northridge earthquake of Jan. 17, 1994 (2 nd column): The spectrogram calculated with Hann window of 32 points (1 st line top), the time-frequency representation (2 nd line), the estimated frequency modulations (3 rd line), and the estimated amplitude modulations (4 th line)

  t and is given by the Duhamel integral

	( ) y t k	0 t p P h t ( ). ( k k τ	) τ τ d	( ) p P t k	ω	1 Dk	e	-	k k ξ ω	t	sin( ω	Dk	t	),

Table 1 : Energy calculation for each mode of the ambient and seismic excitations

 1

		E = Energy of the Ambient excitations	' E = Energy of the Seismic excitations	E ∆ =	E E	'
	Mode 1 at 0.4 Hz	5.24	0.1141	45.9
	Mode 2 at 1.2 Hz	0.0468	1.044* 4 10 -	448.2
	Mode 3 at 2.2 Hz	0.0278	9.351* 5 10 -	297.3262

3. Signal models in view of AM/FM estimations

In this paper, we intend to estimate the AM/FM of these seismic signals. A signal model [ ] y n is defined as a discrete time signal with time varying amplitude and frequency, having short time duration, and composed of a multi-component deterministic process [ ] s n embedded in an additive white Gaussian noise [ ] e n with zero mean and variance 2 σ .

Table 2 : The two studied signal models

 2

	Time-Invariant Damping Coefficient	Time-Variant Damping Coefficient
	-Damped Amp.& Polynomial Freq. Model -	-Polynomial Damping Function Model-
	(IA)

Acknowledgements

This work has been supported by French Research National Agency (ANR) through RISKNAT program (project URBASIS ANR-09-RISK-009).