Haiyang Ding
email: haiyang.ding@gipsa-lab.fr

Mazen Alamir
email: mazen.alamir@gipsa-lab.fr

Ahmad Hably
email: ahmad.hably@gipsa-lab.fr

A distributed cooperative control scheme with optimal priority assignment and stability assessment

In this paper, a distributed partially cooperative control framework is proposed for a network of linear interconnected subsystems. It is assumed that each subsystem in the network possesses its own objective and a corresponding nominal interaction-free state feedback law. The proposed framework enables each subsystem to compute an additional control term in order to help maintaining the integrity of the overall network. As this cooperation-like behavior involves relative priority assignment, a communication aware heuristic is proposed with an associated stability assessment that is based on the closed-loop network matrix's spectrum monitoring. Illustrative examples are used to assess the effectiveness of the proposed scheme including a distributed load frequency problem.

INTRODUCTION

The design of distributed control framework for interconnected systems has received a great attention in the last recent years due to the fact that centralized frameworks are often very expensive in terms of communication if not impossible to be settled. On the other hand, when the interactions between subsystems are strong, the fully decentralized control schemes show bad performance [START_REF] Ho | On centralized optimal control[END_REF], [START_REF] Cui | Performance limitations in decentralized control[END_REF]. To overcome these limits, distributed control schemes have been proposed. The main idea is that each subsystem of the network has to cooperate with its neighbors by sharing partial information in order to contribute to the global objective of the network (or to manage its own objective) while avoiding destabilizing the whole network. The literature related to the distributed/cooperative control of interconnected networks has witnessed an impressive growth these last years. This makes quite difficult an exhaustive relative positioning of the existing propositions. An efficient way towards a comprehensive positioning is to consider the control paradigm being addressed, the information shared between network's subsystems and the coupling between them, and at the end the stability enhancement and assessment.

• In [START_REF] Dunbar | Distributed receding horizon control for multivehicle formation stabilization[END_REF], Richards and How [2004a], [START_REF] Franco | Cooperative constrained control of distributed agent with nonlinear dynamics and delayed information exchange: A stabilizing receding horizon approach[END_REF], the paradigm lies in performing a distributed solution of a unified optimal control problem that involves all the subsystems. • The amount of information items that is communicated or shared among the subsystems is another important issue in distributed control design. In T. [START_REF] Borrelli | Predictive control approach to autonomous vehicle steering[END_REF], Brett T. [START_REF] James | Cooperative distributed model predictive control[END_REF], the dynamic models of all neighbors are completely shared while in [START_REF] Dunbar | Distributed receding horizon control of dynamically coupled nonlinear systems[END_REF], A.N. Venkat and Wright [2008], Richards and How [2004a], the whole states and control vector and/or trajectories are transmitted. • The coupling between subsystems in the networks can be limited to terms involving only control vectors [START_REF] Camponogara | Distributed optimization for model predictive control of lineardynamic networks[END_REF]), state vectors [START_REF] Jia | Distributed model predictive control[END_REF], E. [START_REF] Krogh | Distributed model predictive control[END_REF], [START_REF] Dunbar | Distributed receding horizon control of dynamically coupled nonlinear systems[END_REF]) or general coupling terms that do not induce instability [START_REF] Dunbar | Distributed receding horizon control of vehicle platoons: Stability and string stability[END_REF], Richards and How [2004b], T. [START_REF] Borrelli | Predictive control approach to autonomous vehicle steering[END_REF]). • The stability assessment results can be obtained by monitoring the spectrum of the closed-loop network's matrix either through LMI design or by checking the eigenvalues [START_REF] Michel | Stability analysis of discretetime interconnected systems via computer-generated lyapunov functions with applications to digital filters[END_REF], [START_REF] Jia | Distributed model predictive control[END_REF], E. [START_REF] Krogh | Distributed model predictive control[END_REF], T. Li and [START_REF] Li | Synchronization control for arrays of coupled discrete-time delayed cohengrossberg neural networks[END_REF]).

In the present work, we consider linear networks involving general potentially destabilizing interconnection terms that involve both states and control. It is assumed that each subsystem has its own objective, its own nominal feedback law and its own Lyapunov function. This paradigm is industrially relevant as practitioners often prefer rather modular design framework in which networkrelated issues and solutions are added as an additional layer that does not question the old existing widely assessed local controllers. The information exchanged between subsystems preserves the details of the local decisions and set-ups from being totally exchanged and only aggregated quantities are blindly transmitted as far as this transmission helps preserving the stability of the overall network. Such kind of parsimonious information exchange has been already adopted in [START_REF] Farina | Distributed predictive control: A non-cooperative algorithm with neighbor-toneighbor communication for linear systems[END_REF], [START_REF] Namerikawa | Distributed load frequency control of electrical power networks via iterative gradient methods[END_REF], [START_REF] Richards | Robust distributed model predictive control[END_REF]) where only the interaction-related terms are exchanged between neighbors. For the stability assessment, some priority coefficients optimized off-line have been used. More precisely, for a given set of these priority coefficients, the controls at the subsystems' level are completely defined.

Only once these coefficients have been successively optimized by the network's control designer who has the complete knowledge of the overall network then the distributed cooperative control scheme stability can be guaranteed.

The paper is organized as follows. In section 2, the control problem is stated and the network equations and related assumptions are given. In section 3, the distributed cooperative control scheme is progressively introduced highlighting the role of the priority matrix. More precisely, it is shown how each subsystem in the network constructs an optimal control problem involving both its own objective and the stability-related concerns of its neighbors. The latter are transmitted through dedicated exchanged items.

The way terms coming from neighbors impacts the local optimal control problem heavily depends on the priority of the neighbors as viewed by the local subsystems. The way the priority matrix is optimized using off-line dedicated nonlinear programming is detailed in section 4. Section 5 shows the effectiveness of the proposed scheme with stability assessment using two examples including a well known frequency distributed control problem in order to illustrate the different features of the proposed framework and to assess its efficiency in yielding communication-aware stabilizing cooperative control scheme. Finally section 6 concludes the paper and gives hints for further investigation.

The following notations are used throughout the paper:

• N : total number of subsystems in the network • N → j : number of systems that are affected by system j • n u i and n i : dimension of input and state vector of subsystem i • O: array with only zero elements of appropriate dimensions • → and ← indicate which system is affecting the other.

For example, a variable with (j → i) indicates that this variables is involved in the way subsystem j influences subsystem i • and indicate the direction of the information circulation. For example, a variable with (j i) indicates that it is sent from subsystem j to subsystem i • I → j : A set containing indices of subsystems that are affected by subsystem j • I ← j : A set containing indices of subsystems that affect system j

PROBLEM STATEMENT AT A SUBSYSTEM LEVEL

Let us consider an interconnected network in which a subsystem i is governed by the following Ordinary Differential Equation (ODE):

ẋi = A i x i + B i u i + j∈I ← i (A j→i x j + B j→i u j ) (1)
where x i , u i are respectively the state and the control vector of subsystem i while the summation term describes the effects of interaction with neighbors on subsystem i. Based on the above argument, it is assumed that subsystem i can be stabilized by a local controller: u i = -K i x i such that there are a Lyapunov function V i (x i ) =

x T i P i x i and some positive definite function S(x i ) such that the following inequality holds in the nominal case:

Vi

(x i )| ẋi=A cl i xi ≤ -S i (x i ) (2)
where P i is a symmetric positive definite matrix and A cl i = A i -B i K i is used to denote the nominal closedloop matrix of subsystem i. When interaction with other subsystems is present, a cooperative control term v i should be designed and added to the nominal term by which subsystem i can keep its own objective while helping its neighbors to maintain the global stability. Thus, 1) can be rewritten in the following more compact form:

u i = -K i x i + v i (3) Consequently, equation (
ẋi = A cl i x i + j∈I ← i (A j→i x j + B j→i u j ) (4)
in which i ∈ I ← i meaning simply that subsystem i is impacting itself. This convention is used in the reminder of the paper. Once the interaction is taken into account, the inequality (2) becomes:

Vi

(x i ) ≤ -S i (x i )+ j∈I ← i L j→i (x i )•(A j→i x j +B j→i u j ) (5)
where L j→i (x i ) = 2x i T P i . Based on the above expressions, one has the following straightforward lemma:

Lemma 1. [Stability constraint for subsystem i]
The subsystem i is asymptotically stable provided that the following inequality holds:

j∈I ← i L j→i (x i ) • (A j→i x j + B j→i u j ) < µS i (x i ) (6)
for some µ ∈ [0, 1[. △

Now this requirement has to be sent, in a rather condensed form to all the neighbors affecting i (namely, having indices j ∈ I ← i ). For this reason, the following assumption is needed regarding communication rules as well as mutual a priori knowledge available at each subsystem regarding the neighbors characteristics:

Assumption 1. [Knowledge/Communication rule & Network topology]
(1) The structure of interaction in the network is bidirectional, namely:

{i ∈ I ← j } ⇐⇒ {j ∈ I ← i } (2) Each subsystem j that affects i (j ∈ I ← i
) is aware of the interaction matrices A j→i and B j→i involved in (4).

(3) At each time step, subsystem i sends to its neighbors j ∈ I ← i (affecting it) the following quantity:

W ij := L j→i (x i ) (7)
It should be pointed out that each subsystem j (j = i) neither knows the dynamics of W ij nor the explicit value of the state x i for all i that is affected by it (i ∈ I → j ). The information being sent concerns the instantaneous value of L j→i (x i ). Note however that as soon as the subsystem j receives the value of W ij , it can consider the left hand side of (6) as a linear combination of x j and u j that has to be reduced if subsystem i needs to be stabilized. But since there are as many subsystems i as there are elements in the set of indices I → j , there is clearly need to define a relative priority weighting viewed by subsystem j. This is the aim of the following section.

The priority vector

In many networks, some subsystems are more critical than others. In such situations, it seems reasonable to consider the stability constraints sent by such critical subsystems with higher weight (or priority) that the stability constraints sent by other less critical subsystems. Since only neighbors belonging to I → j are considered by j, the relative priority weights are intimately linked to the connection topology. To better explain this feature, let us consider the following example where a network with 4 subsystems coupled in a chain graph as shown in the right drawing of Fig. 1. An associated priority matrix can be given in the 8), Π and π j denote respectively the priority matrix of the network and the priority vector of subsystem j. π j i determines how important subsystem i is viewed by subsystem j. Roughly speaking, the higher π j i is, the more important subsystem i is considered by subsystem j.

Π = [ π 1 π 2 π 3 π 4 ] =          π 1 1 π 1 2 0 0 π 2 1 π 2 2 π 2 3 0 0 π 3 2 π 3 3 π 3 4 0 0 π 4 3 π 4 4          ; π j i ≥ 0 (8) In (
On receiving the information W ij subsystem j defines a cost function to be minimized based on a given priority vector in order to cooperate with all subsystems having index i ∈ I → j . Namely min

vj { v j 2 Rj + i∈I → j π i j [W ij (A j→i x j + (9)
B j→i (-K j x j + v j ))]} Note that the above optimization problem is defined given a pre-assigned value of the priority vector π j .

SOLVING THE OPTIMIZATION PROBLEM (9)

In order to solve (9) and derive the additional control term v j , an extended dynamic system is constructed:

żj = A j z j + B j v j = A cl j G 1 j . . . G N → j j O Γ z j + B j v j (10)
where

z j :=       x j π 1 j W T 1j . . . π N → j j W T N → j j       , B j = B j O , Γ = diag(-ǫ) with ǫ is a positive constant, G i j = 1 2π i j • (A i→j -B i→j K i ) • (P i ) -1
while assuming a hypothetical exponentially decaying dynamic of W T ij 1 .

Rewriting of the cost function

By expressing π i j W ij and x j in terms of the extended state z j , and by using the following assumption: Assumption 2. [Shared Knowledge] For any subsystem i that are affected by j (i ∈ I → j ), subsystem j knows the matrix

(A i→j -B i→j K i ) • (P i ) -1 ∈ R nj ×ni △
The optimization problem ( 9) can be put in the following form: min

vj {v j T R j v j + z j T Q coop j z j + z j T N coop j v j } (11)
with

Q coop j =      0 A j→1 . . . A j→N → j -B j→1 K j 0 . . . 0 . . . 0 . . . 0 -B j→N → j K j 0 . . . 0     
and

N coop j =      O n * nu B j→1 . . . B j→N → j     
By now, the linear dynamic model ( 10) and the quadratic cost function (11) completely define a quadratic programming problem that can be solved in a discrete-time setting to compute the matrix gain K v j that is used by subsystem j in order to compute the additional control v j according to:

v j = -K v j z j (12)
The use of sampling control here enables the non necessarily negligible communication time to be taken into account. Note that the above optimization problem is defined given a pre-assigned value of the priority vector π j . In section 4, it is shown how the priority vectors can be computed such that the resulting overall interconnected cooperative network is asymptotically stable. However, let us first show how (9) can be solved at the subsystem j level.

STABILITY ENHANCEMENT THROUGH THE OPTIMIZATION OF THE PRIORITY VECTORS

The preceding sections propose a heuristic that results in the complete definition of the distributed feedback control once given a priority matrix Π [see ( 8)]. In other words, regardless how justified is the heuristic proposed above, it leads to a rigorously defined terms v j , j ∈ {1, . . . , N } leading to rigorously defined controls u j , j ∈ {1, . . . , N }. In this section, it is shown that based on these resulting control expressions (issued from heuristically defined procedure), one can rigorously compute the spectrum of the network's dynamic closed loop matrix as a function of the priority matrix Π. Therefore, Π can be optimized by the network control designer (in the off-line design step) in order to guarantee the stability of the network. This is explained in the following sections.

The complete state space model

A complete state space model of the distributed scheme is built in which the states of all the subsystems are involved. The rationale behind this is that although W i→j is known to subsystem j as an information sent by subsystem i, it is known from a designer's view as:

W ij = L j→i (x i ) = 2x i T P i (13) Since W ij
is involved in the extended state z j and (13) is linear in x i , z j can be expressed as linear combination of the state vectors {x i |i ∈ I → j }. Once all the extended states for each subsystem in the network are expressed in this manner, the complete state space model that is equivalent to the distributed scheme can be established. This results in the following compact autonomous form:

X + = E(Π)X (14) with X = x 1 T x 2 T • • • x N T T
in which Π is the priority matrix defined above.

Stability analysis and off-line optimization of the priority matrix Π

Let us use the following notation in order to refer to the spectrum of the network matrix for a given priority vector Π:

α(Π) = max i |λ i (E(Π))| (15)
where λ i is the i-th eigenvalue of E(Π). Therefore, the stability of the distributed control based on the heuristic explained in the preceding sections can be guaranteed as soon as a priority matrix Π can be found such that α(Π) < 1. Note also that E(Π) can be written in the following form:

E(Π) = π1 E 1 + • • • + πNp E Np ( 16
) in which N p is the number of nonzero elements in the connection matrix C while πi 's are the corresponding elements of the priority matrix. The problem of minimizing the radius of the spectrum of ( 16) is a well known problem in the matrix-pencil literature [START_REF] Burke | Optimizing matrix stability[END_REF]). Several methods can be used to minimize the radius α(Π):

min Π∈P {α(Π)} (17)
where P ∈ R Np is a compact subset of admissible values to defined according to the context. In the sequel, the MATLAB subroutine fmincon is used to solve the above optimization problem leading to the optimal priority matrix Π.

Remark 2. When considering a large scale network, the problem of defining the form of the priority matrix is important since the matrix pencil optimization with high dimension is hard to perform. One potential solution is to reduce the number of degrees of freedom in the definition of the priority matrix by using reduced dimensional parametrization of the form Π = Π(ρ). The dimension of the optimization problem is therefore defined by the dimension of the parameter vector ρ. ♦

ILLUSTRATIVE EXAMPLES

In this section, some illustrative examples are given to show the efficiency of the proposed framework. In all the examples, evolutions labeled with coop refer to the distributed cooperative method while those labeled with noncoop correspond to scenarios where the distributed cooperative control is not activated (i.e. u i = -K i x i ; v i = 0).

Example 1: Stability enhancement of a numerical case study

In this first example, a network is built randomly using 3 subsystems with 2 states and 1 control each. The subsystems are interconnected with each other in a cyclic graph as in Fig. 1 through states and control inputs. The continuous-time equation of the subsystems are given as follows: ẋi

= A i x i + B i u i + β j∈I ← i (A j→i x j + B j→i u j ) (18)
in which A i is the open-loop state matrix governing subsystem i. The nominal local controllers are designed with LQR method:

A cl i = A i -B i K i .
The cooperative control is computed as presented in the previous sections. Note that β expresses the strength of the coupling between subsystems. The values of the matrices can be found in the appendix. At each β, three different control schemes are applied to the resulted network: (a) only nominal decentralized controllers are activated, (b) the cooperative control with equal priority and (c) the cooperative control with optimal priority. In the cooperative control with equal priority, all the priority vectors are set to π i = [ 1 1 1 ] T for each subsystem. The optimal priority is obtained using the method described in section 4. The radiuses of the spectrum α(Π) of the three complete state space models when β grows from 0 are shown in Fig. 2. Fig. 3 is a zoom-in of Fig. 2. The zoom-in figure shows that nominal decentralized control can stabilize the network when β < 1.62. Beyond this value, subsystems with cooperation (using equal priorities) preserves the network stability as far as β < 1.82. The figure on the right shows that the cooperative control with optimal priority is able to stabilize the subsystems until β reaches 15. The fluctuations of curve is due to the stopping criterion of the algorithm.

Example 2: Load Frequency control problem

In this example, the distributed control scheme is applied to the load frequency control problem described in [START_REF] Namerikawa | Distributed load frequency control of electrical power networks via iterative gradient methods[END_REF]). The network consists of 4). The parameters of the network are taken according to [START_REF] Namerikawa | Distributed load frequency control of electrical power networks via iterative gradient methods[END_REF]) and the states of each area are influenced by a white noise with covariance matrix 10 -3 × I. An additional coefficient β is added to tune the strength of coupling between areas. Four case studies are presented to illustrate the advantages of the proposed scheme. The equations of the network are given in appendix C.

Case study 1 -Cooperative control with equal priority:

In the first simulation, the coupling parameter β = 1 is used (leading to the same model as [START_REF] Namerikawa | Distributed load frequency control of electrical power networks via iterative gradient methods[END_REF])). When the cooperative control is not activated, the network spectrum radius is 0.99 which means that the network is stable under the decentralized nominal control. Nevertheless, when looking at the behaviors of the Lyapunov function of the cooperative control V coop and that of the decentralized nominal control V noncoop , it can be seen on Fig. 5 that the firing of the cooperation algorithm leads to far better overall performance. Note that for the purpose of better illustrating the evolution of the Lyapunov functions (which we choose to evaluate the performances of the areas), the initial conditions for all states of the ares are set to 1 (since all of them are stable). V coop -V noncoop indicating that cooperation with equal priority enhance the overall performance of the network even when decentralized control is sufficient to achieve the stability (weak coupling).

Case study 2 -Cooperative control with higher priority on certain subsystem: In this simulation the coupling parameter is also chosen as β = 1, and we apply the priority matrix defined by:

Π =    10 10 0 0 1 1 1 0 0 1 1 1 0 0 1 1   
which means that Area 1 is given a higher priority viewed by both itself and its neighbor Area 2. The behavior of the Lyapunov function V coopArea1 is compared to that of the cooperative control with equal priority V coop . One can see form Fig. 6 that the behavior of Area 1 is improved while for other areas, since they are defined to be less important, their performances are a little sacrificed. V coopArea1 -V coop indicating that when higher priority is applied to Area1, its performance is better than in the equal priority situation at a cost of the reduction of its neighbors' performances (weak coupling)

Case study 3 -Cooperative control with optimal priority: With the coupling parameter remaining as β = 1, we explore the performance of the network when the optimal priority matrix is applied. At this step, each element in the priority matrix is chosen between π j i ∈ [0.1, 10] in the matrix pencil optimization. The resulted optimal priority matrix is given by: Π =    9.25 4.67 0 0 5.26 9.56 3.70 0 0 9.42 9.64 2.48 0 0 8.51 9.52

  
Similarly, the performance of the network is verified by looking at the behaviors of the Lyapunov functions of each area V coopOptimal . V coopOptimal is compared to V coop obtained in section 5.2.1. Fig. 7 depicts the improvement of the behaviors of all four areas. V coopOptimal -V coop indicating that, in this case when optimal priority assignment is applied to all areas, their performances are all better than in the equal priority situation (weak coupling)

Case study 4 -Cooperative control under destabilizing interaction: In this simulation, the coupling parameter is set as β = 2 and this changes the interaction between the areas to a destabilizing one. To stabilize the network despite of the destabilizing interconnection, we design the optimal priority matrix as presented in section 4 and each element is chosen between [0.1, 1000]. The optimal priority matrix is:

Π =   
23.6 0.1 0 0 1000 0.1 1000 0 0 1000 0.1 1000 0 0 1000 67.4

  

The initial conditions for all areas are zero and the evolutions of frequency deviation of each area under decentralized nominal control and cooperative control are shown in Fig. 8 and Fig. 9 respectively. The two figures show that the behaviors of the areas are unstable and unacceptable under decentralized nominal control while the cooperative control with the presence of optimal priority selection is able to stabilize all areas and the resulting frequency deviation is well controlled. The fluctuation is caused by a white noise to all states with covariance matrix 10 -3 × I.

CONCLUSION

In this paper, distributed control framework with limited communication has been proposed for interconnected subsystems with potentially destabilizing coupling. While each subsystem in the network has its own objective through an associated Lyapunov function, subsystems may cooperate when the integrity of the network is threatened. The choice of priority vector represents a standard matrix pencil optimization problem which can be solved offline using non convex optimization algorithms that may becomes interesting when high dimensional problems are tackled. The proposed framework is particularly attractive as a first step in the design of distributed control for nonlinear networks. Indeed, unlike many works where the degrees of freedom are represented by the linear feedback gain matrices which loose their relevance as soon as nonlinear dynamics are involved, the concept of relative priority coefficients can be used in both linear and nonlinear design [see (M. Alamir and Ding [2011])]. Moreover, the optimal relative priority set that is computed using the linearized model is likely to remain relevant even for the nonlinear network (at least as far as deviations are not very important). They can therefore be used as relevant design parameters for the nonlinear distributed control design.

Investigating such situations and comparisons between such algorithms are currently under progress and represent obvious perspectives for the present paper. The network consists of four identical subsystems (areas) that are interconnected in a chain graph as shown in Fig. 4 . Each subsystem has 7 states: ∆P tiei , ∆f i , ∆P gi , ∆x vgi , ∆P Ei , ∆P Hi and U ARi . ∆P tiei is the tie-line power flow of subsystem i. ∆f i is the frequency deviation. ∆P gi is the output of gas turbine generator while ∆x vgi is its input. ∆P Ei , ∆P Hi are the outputs of the battery electric storage system and the electric water heater respectively. The regional demand is AR i = ∆P tiei -B i ∆f i and we define U ARi = AR i dt. The interconnected term ∆P tiei of subsystem i is expressed as ∆P tiei = T ij (∆f j -f i ) where j stands for the interconnecting subsystem with i. a g , a E and a H are chosen to satisfy the ratio of the capacity of the corresponding gas turbine, battery electric storage and electric water heater. Moreover, they also satisfy the following equality a g +a E +a H = 1. B i , T ij and R g are the frequency bias, synchronizing coefficient and regulation constant respectively. The state space models of each subsystem are given in the following form:

ẋi = A i x i + B i u i + ω + β j∈I ← i A j→i x j (C.1)
ω is a white noise added to all states (with covariance 10 -3 in this paper). The states and the system matrices can be written as:

x i = [ ∆Ptie i ∆fi ∆Pg i ∆xvg i ∆P E i ∆P H i U AR i ] T (C.2) B i =              
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 T gi 0 0 0 0 0 0 0 1 T Ei 0 0 0 0 0 0 0 1 T Hi 0 0 0 0 0 0 0 0

              (C.3) A i =                    0 - j∈I ← i Tij 0 0 0 0 0 1 M i - D M i 1 M i 0 1 M i - 1 M i 0 0 0 - 1 T di 1 T di 0 0 0 0 - 1 T gi R gi 0 - 1 T gi 0 0 a g K i T g 0 - 1 T Ei R gi 0 0 - 1 T Ei 0 a g K i T E 0 - 1 T Hi R gi 0 0 0 - 1 T Hi a g K i T H 1 -B i 0 0 0 0 0                    (C.4) A j→i =        
0 T ij 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

        , j ∈ I ← i (C.5) A j→i = O ∈ R 7×7 , j / ∈ I ← i (C.6)
The numerical values of the parameters are taken from [START_REF] Namerikawa | Distributed load frequency control of electrical power networks via iterative gradient methods[END_REF]) and are shown on

Fig. 1 .

 1 Fig. 1. Example of a coupled networks following form:

Fig. 3 .

 3 Fig. 2. Example 1: The evolution of the spectrum of the network under growing interconnection strength.

Fig. 4 .

 4 Fig. 4. The structure of the power network N = 4 subsystems (areas) and each of them has n = 7 states which are coupled in a chain-like graph (see Fig.4). The parameters of the network are taken according to[START_REF] Namerikawa | Distributed load frequency control of electrical power networks via iterative gradient methods[END_REF]) and the states of each area are influenced by a white noise with covariance matrix 10 -3 × I. An additional coefficient β is added to tune the strength of coupling between areas. Four case studies are presented to illustrate the advantages of the proposed scheme. The equations of the network are given in appendix C.

Fig. 5 .

 5 Fig. 5. Example 2 -case study 1: Evolution of ∆V r =V coop -V noncoop indicating that cooperation with equal priority enhance the overall performance of the network even when decentralized control is sufficient to achieve the stability (weak coupling).

Fig. 6 .

 6 Fig. 6. Example 2 -case study 2: Evolution of ∆V r =V coopArea1 -V coop indicating that when higher priority is applied to Area1, its performance is better than in the equal priority situation at a cost of the reduction of its neighbors' performances (weak coupling)

Fig. 7 .

 7 Fig. 7. Example 2 -case study 3: Evolution of ∆V r =V coopOptimal -V coop indicating that, in this case when optimal priority assignment is applied to all areas, their performances are all better than in the equal priority situation (weak coupling)

Fig. 8 .Fig. 9 .

 89 Fig. 8. Example 2 -case study 4: Frequency deviation without cooperative control in the case β = 2 is used (destabilizing interaction)

  Table. C. Table C.1. Values of the parameters

	Parameters	Symbol	Value	Unit
	inertia constant	M	0.20	puM W • s/Hz
	damping constant	D	0.26	puM W/Hz
	governor time constant	Tg	0.20	s
	gas turbine constant	T d	5.0	s
	BESS time constant	T E	0.20	s
	HP time constant	T H	4.5	s
	Regulation constant	Rg	2.5	Hz/puM W
	Synchronizing coefficient	T ij	0.50	puM W

Appendix A. FINDING THE MATRICES As for the matrices G i j , i = 1, . . . , N → j , they are derived as follows (resulting in the equation (A.5) below) : Let us first recall the equation of subsystem j:

In order to build the dynamics that enables the additional control v j to be computed, subsystem j assumes that the cooperative actions v i are zero for all subsystems is that are affected by j (i ∈ I ← j ). The bi-directional interconnection relationship required in Assumption 1 is used now, and one can compare (A.2) and the first n j rows (the dynamics corresponding to ẋj ) of ( 10) to find that:

and since W ij = 2x i T P i , one has:

Therefore the matrices G i j involved in the extended model (10) that is used by subsystem j is given by:


J-C Beugnot

S Lebrun

G Pauliat

H Maillotte

V Laude

T Sylvestre

Observation of surface acoustic wave Brillouin scattering in optical microfibers

We report the observation of surface-acoustic-wave Brillouin scattering (SAWBS) in a subwavelength-diameter fiber and show that this effect relies on the generation by the electrostrictive force of phonons confined at the surface of the microwire.

Brillouin scattering between sound and light waves has been extensively investigated in optical fibers. It is wellknown that silica optical fibers support many transverse and longitudinal elastic waves over a broad frequency range from 50 MHz to 11 GHz, that lead to guided-acoustic-wave Brillouin scattering (GAWBS) and stimulated Brillouin scattering (SBS), respectively [START_REF] Sylvestre | Nonlinear Fiber Optics[END_REF].

In this work, we demonstrate, for the first time to our knowledge, the generation of a new class of surface acoustic waves (SAW) in subwavelength-diameter optical fiber (SDOF) owing to the strong coupling between sound and light at interface between the optical microwire and the surrounding air. These SAW modes or Rayleigh waves cannot be optically excited in standard single-mode fiber (SMF), but they also have recently been evidenced in whispering gallery modes (WGM) microresonators and in silica microspheres [START_REF] Matsko | Optomechanics with Surface-Acoustic-Wave Whispering-Gallery Modes[END_REF][START_REF] Bahl | Stimulated optomechanical excitation of surface acoustic waves in a microdevice[END_REF]. In our experiment, we investigated a 8-cm long silica optical microwire drawn from a commercial single-mode fiber (SMF) by hydrogen flame. It has a waist of 1 µm and linear losses of 0.1 dB/cm at 1.55 µm. We then performed spontaneous Brillouin backscattering measurements at 1.55 µm using a standard heterodyne detection [START_REF] Beugnot | Complete experimental characterization of SBS in photonic crystal fiber[END_REF]. Figure 1(a) shows the spontaneous Brillouin spectrum for an input CW power of 100 mW. Several peak frequencies with different weights and linewidths can be seen in a frequency range from 6 GHz to 11 GHz. The highest frequency peak at 10.86 GHz comes from the standard SBS in the untapered SMF sections of about 2 m long. We can also observe in Fig. 1(a) three other acoustic frequencies centered at 8.33 GHz, 9.3 GHz and 10 GHz, respectively. Like in SMF, these three acoustic modes exhibit a spectral linewidth around 25 MHz in agreement with the acoustic phonon lifetime in fused silica. Two other acoustic resonances at lower frequencies around 6 GHz can also be observed in Fig. 1(a). To identify the physical origin of these elastic resonances, we performed numerical calculations of the phonon wave packets by solving the elastodynamics equation subject to an electrostrictive force (for detailed calculations, See Refs. [START_REF] Beugnot | Electrostriction and guidance of acoustic phonons in optical fibers[END_REF][START_REF] Laude | Generation of phonons from electrostriction in small-core optical waveguides[END_REF]). With this model, we are able to fully characterize the optoacoustic properties of our silica microwire including both the forward and backward Brillouin spectra without the requirement of a full band structure. We also consider in our model the phonon lifetime by taking into account the elastic losses assuming a complex tensor [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF]. This loss model is compatible with the usual assumption that the product of the quality factor Q and the acoustic frequency is a constant for a given material (e.g., for silica, Q×f =5 THz). Figure 1(b) shows the normalized elastic energy distribution versus the acoustic frequency and for varying the waist diameter from 1 µm to 1.35 µm. The elastic energy distribution gives an overview of the theoretical Brillouin spectrum. As it can be seen, we retrieve most of the acoustic resonances experimentally observed in the waist region of the microwire. There is however a slight difference with experiment in frequency position that can be attributed to the nonuniformity of the microwire, but we can clearly identify the two modes around 6 GHz and the 3 other ones around 9 GHz. To go further into detail, we plot in Fig. 1(c) the spatial distribution of the optical mode intensity in the waist region at a wavelength of 1.55 µm. The black circle marks out the interface between the optical field inside the microwire and the strong evanescent field outside the wire. A large portion of this evanescent field enables the diffusion of light by surface acoustic wave and the highly confined optical field into the microwire induces an electrostrictive force 20 times larger than in SMF [START_REF] Laude | Generation of phonons from electrostriction in small-core optical waveguides[END_REF]. Figure 1(b) shows this electrostrictive force distribution, which was numerically calculated from our model. The kinetic energy density or elastic energy of the two surface acoustic waves at 5.382 GHz and 5.772 GHz are plotted in Figs. 1(e,f) for a waist of 1.05 µm. The elastic energy is mainly localized at the surface of the microwire, unlike the longitudinal acoustic waves, shown in Figs. 1(g,h), for which the generated phonon wave packet is guided inside the microwire's core by the optical force. Thus these SAW modes can only be excited in sub-wavelength waveguides exhibiting a strong evanescent field. They propagate along the surface with an acoustic velocity of 3410 m/s, which is lower than both the shear velocity (3740 m/s) and the longitudinal velocity (5960 m/s) [START_REF] Thurston | Elastic waves in rods and clad rods[END_REF].

To summarize, we demonstrated in an optical microfiber with a 1 µm waist the efficient generation of surface acoustic wave Brillouin scattering (SAWBS) at new frequency around 6 GHz. This work should contribute to a further understanding of the intriguing nonlinear opto-acoustic interactions in sub-wavelength optics. In addition, SAWBS could be used in various sensor applications since the SAW frequency should strongly depend on the material surrounding the optical microwire.

This work is supported by the Université de Franche-Comté and the ANR LABEX ACTION.

Figure 1 :

 1 Figure 1: (a) Experimental Brillouin scattering spectrum measured at 1.55 µm in 8-cm long silica optical microfiber with a waist of 1 µm and for a CW pump power of 100 mW. (b) Computed elastic energy of phonon wave packet (i.e., the theoretical Brillouin spectrum) versus the acoustic frequency and for a microwire waist varying from 1 µm to 1.35 µm. The black line shows the SBS phase-matching condition (K=2Kp). (c-g) Numerical results for a 1.05 µm silica rod diameter. (c) Optical power in the fundamental TE-like mode at a wavelength of 1.55 µm (neff = 1.1887). (d) Spatial distribution of the electrostrictive force with 1 W optical power. (e,f) Elastic energy density of resonant surface acoustic waves at 5.382 GHz and 5.772 GHz, respectively. (g,h) Elastic energy density of longitudinal acoustic modes at 8.37 GHz and 9.235 GHz.


Jean-Luc Menet
email: jean-luc.menet@univ-valenciennes.fr

Ion-Cosmin Gruescu

COMPARATIVE LIFE CYCLE ASSESSMENT OF ALCALINE CELLS AND NI-MH RECHARGEABLE BATTERIES

Keywords: Life Cycle Assessment, LCA, batteries, alkaline cells, environmental impacts, climate change indicator

The present study presents some results obtained by applying the LCA methodology to evaluate the environmental footprint of alkaline cells and Ni-MH batteries. The approach is motivated by the increasing number on markets of electronic systems needing local and portable electricity.

The proposed study focuses on a comparison between these two elements on the basis of the same produced energy. The environmental impacts are evaluated by using the Ecoinvent 2.0 data-base, with midpoint indicators associated to the CML method.

The main result of the study let appear that the difference between the environmental impacts of the studied elements is mainly due to the transportation phase. Besides, the emphasis on rechargeable batteries is only justified from an environmental point of view if the consumer uses them according to recommendations of use, i.e. recharges them 250 times.

INTRODUCTION

Rechargeable or not, cells and batteries are known to greatly impact the environment. There are supposed to be sorted following a specific way which depends on the country. For example, in France, a transposition of the EU directive 2006/66/CE has modified the L.541-10 article of the environment code, and lead to the decree 2009-1139 [START_REF]relatif à la mise sur le marché de piles et accumulateurs et à leur élimination modifiant la partie réglementaire le code de l'environnement[END_REF].

Nevertheless, these chemical elements can be often found in classical trash cans, among usual wastes. In France, data for year 2009 are the following [2]: Among the 250 000 tons which were put on the market for all the categories, about 127 000 tons concerned cells and portable batteries; More than 207000 tons of used cells and batteries have been collected (concerning cells and portable batteries, the collection rate is 32%) There are about 40000 places where the collection is organized. The collected cells and batteries 84% are recycled and 90% are recovered. In other words, less than 30% of these elements are really recycled or recovered. Besides, in all cases, the real question is to know their impact on the environment. But such studies are not so numerous. A recent research report [START_REF] Olivetti | Life cycle impacts of alkaline batteries with a focus on end-of-life[END_REF] has studied the environmental impact of alkaline batteries. Another report [START_REF] Nors | Opportunities and challenges in utilization of life cycle data in consumers' portable batteries with focus end of life[END_REF] shows that by identifying portable batteries (with RFID), it is possible to connect a product's end of life datum to their design and manufacturing, which should help sustainability on the whole life cycle. The Life Cycle Assessment methodology, which is a standardized approach [START_REF]Environmental management -Life cycle assessment -Principles and framework[END_REF][START_REF]Environmental management -Life cycle assessment -Requirements and guidelines[END_REF] allows the quantification of the environmental footprint for goods, services and processes. It is then possible to identify some main points allowing the diminution of the environmental impacts since their earlier design stage. This methodology has been recently used in the field of cells and batteries [START_REF] Uniross | Environmental impact of batteries[END_REF] by a manufacturer who concluded that rechargeable batteries have up to 32 times less impact than classical cells.

We have recently used the LCA methodology for performing some case studies concerning building elements [START_REF] Menet | A comparative life cycle assessment of exterior walls constructed using natural insulation materials[END_REF][START_REF] Menet | Environmental footprint of building elements using Life Cycle Analysis methodology, 21è[END_REF]. The idea was to show that the "common sense" generally leads to amazing conclusions, with environment impacts which can be higher for products said to be environmentally friendly. The same conclusion has been recently obtained for a system using renewable energy, especially because of the use of portable batteries [START_REF] Menet | A simplified Life Cycle Assessment applied to a coupled Solar and Eolic street light[END_REF].

In the present work, we aim to do a comparative study of the environmental impacts of Ni-MH rechargeable batteries and alkaline cells. To achieve this aim, a LCA study is made following the standardized methodology reported on figure 1.

GOAL AND SCOPE DEFINITION

The goal of the present study is the comparison of environmental impacts of the two following elements:

-An alkaline cell (AA type made of alkaline manganese, 1.5 V, 2850 mAh) -A portable battery (AA Nickel metal hybrid battery, 1.5 V, 650 mAh) To compare the two elements using the LCA methodology, a Functional Unit must be chosen: in the present study, it corresponds to a delivered power of 1W during 1000 hours, i.e. an energy of 1kWh. This energy is the one which has been used in reference [START_REF] Uniross | Environmental impact of batteries[END_REF]. Let us notice that this energy corresponds to 234 alkaline cells or 5.2 Ni-MH batteries (supposed to be rechargeable 250 times). As the quality of the data is not insured, we admit in the scope definition that the results will only be quantitative in terms of range of order, which means that they can only be significant if the difference between the environmental impacts of the two elements is significant.

The power adapter is not considered in the study.

INVENTORY ANALYSIS (LCI)

Data are not reported in the following. More details about the inventory analysis can be found in ref [START_REF] Wurier | ACV comparative des piles et accumulateurs[END_REF]. Let us notice that the four steps of the LCA methodology are considered: the production phase, the transportation, the use and the end of life. The use phase is only useful for Ni-MH batteries because of the charge.

IMPACT ASSESSMENT (LCIA)

In order to evaluate the environmental impacts from the inventory analysis, the Eco-invent 2.0 data-base is used [12]. The environmental impacts were estimated using midpoint indicators associated with the CML method [13]. Since each indicator has a different scale and unit of measurement, the relative impacts are normalized, according to the specifications of ISO procedures [START_REF]Environmental management -Life cycle assessment -Principles and framework[END_REF][START_REF]Environmental management -Life cycle assessment -Requirements and guidelines[END_REF]. Each indicator is then expressed in points. A point represents the potential impact for a given indicator divided by the value of the same impact for a mean European people during one day. The chosen indicators and the corresponding units are reported in table 1. Because all the data are not known (for example the weights of the materials in the cells or batteries), and because the processes are not all known, some approximations are made according to the scope of the study (step 1). The results of the calculation are reported in the figure 2.

INTERPRETATION

Two main conclusions can be deduced from the environmental impact assessment:

-For the whole life cycle, the alkaline cell seems to have the greatest impact on the environment than the Ni-MH batteries whatever the indicator; -The difference between the two studied elements is essentially due to the transportation phase needed for the alkaline cells. This is mainly due to the different number of elements in the two cases to achieve the same function: 1kWh of energy: 45 times more cells than batteries. The main results of the present study are reported in the table 2: it is shown that, except for photochemical pollution, the present results and those of ref. [START_REF] Uniross | Environmental impact of batteries[END_REF] are of the same range of order. Let us remind that the UNIROSS study takes the adapter into account [START_REF] Uniross | Environmental impact of batteries[END_REF] but does not communicates all the results. Besides, the studied batteries are not exactly the same in the two studies.

CONCLUSION

The UNIROSS study concludes that the main origins of environmental impact are:

-the production and use phases (charge cycles) for rechargeable batteries -the production phase (between 70 and 100%) for disposable batteries On the contrary, the present study shows that it is the transportation phase which creates the greatest impact, and makes consequently the difference between the two elements.

The present study does not take into account the adapter, and the calculations must be refined. However, the Ni-MH battery is supposed to be recharge 250 times. If it is not the case (with for example only 25 recharges), and/or if the adapter is used more than necessary, the previous results could be very different. This shows that the emphasis on rechargeable batteries is only justified from an environmental point of view if the consumer uses them according to recommendations of use.

Figure 1 :

 1 Figure1: Different phases of the LCA methodology[START_REF]Environmental management -Life cycle assessment -Principles and framework[END_REF][START_REF]Environmental management -Life cycle assessment -Requirements and guidelines[END_REF]

Figure 2 :

 2 Figure 2: Environmental impacts for each life phase (1) alkaline cells (2) Ni-MH batteries

Table 1 -

 1 Chosen indicators for the calculation of the environmental impacts

	Letter	Indicator	Unit
	NRE	non-renewable energy consumption	MJ eq.
	RD	resources depletion	kg Sb eq.
	GWP	100 year Global Warning Potential	kg CO 2 eq.
	A E	Acidification Eutrophication	kg SO 2 eq. kg PO4 2-eq.
	PP	photochemical pollution	kg C 2 H 4 eq.
	AT	aquatic toxicity	kg 1.4-DB eq.
	HT	human ecotoxicity	kg 1.4-DB eq.

Table 2 -

 2 Main results and comparison with ref.[START_REF] Uniross | Environmental impact of batteries[END_REF] (n.c.: not communicated)

		(I)	(II)	I/II	I/II
		Alkaline cells	Ni-MH batteries	(present study)	ref [7]
	NRE	25	6.2	4.0	n.c.
	RD	47	2.9	16	23
	GWP	22	1.5	15	28
	A	20	7.0	2.9	9
	E	4	0.23	18	n.c.
	PP	7.8	2.3	3.4	30
	AT	12	4.8	2.4	n.c.
	HT	4.2	0.71	6.0	n.c.

ACKNOWLEDGEMENTS

This study has been partly made using the free Bilan Produit ® software [12] jointly developed by the French Agency for the Environment and the Energy (ADEME) and the Cergy-Pontoise University.

The authors thank the students Thomas Wurier and Raphaël Mondonneix which have worked on the subject during their university projects [11].


J.-P Chehab

P Garnier
email: pierre.garnier@u-picardie.fr

Y Mammeri

Long-time behavior of solutions of a BBM equation with generalized damping

Keywords: BBM equation, dispersion, dissipation. MS Codes. 35B40, 35Q53, 76B03, 76B15

We study the long-time behavior of the solution of a damped BBM equation ut + uxuxxt + uux + Lγ (u) = 0. The proposed dampings Lγ generalize standards ones, as parabolic (Lγ (u) = -∆u) or weak damping (Lγ (u) = γu) and allows us to consider a greater range. After establish the local well-posedness in the energy space, we investigate some numerical properties.

Introduction

The one-way propagation of small amplitude, long wavelength, gravity waves in shallow water can be described by the Korteweg-de Vries equation (KdV) [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal and on a new type of long stationnary waves[END_REF] u t + u x + u xxx + uu x = 0.

or its regularized version, the Benjamin-Bona-Mahony (BBM) equation

u t + u x -u xxt + uu x = 0.
Many works related to the damped KdV equation can be found in the literature ( [START_REF] Amick | Decay of solutions of some nonlinear wave equations[END_REF][START_REF] Cabral | Chaos for a damped and forced KdV equation[END_REF][START_REF] Chehab | Numerical study of a family of dissipative KdV equations[END_REF][START_REF] Chehab | On damping rates of dissipative KdV equations[END_REF][START_REF] Ghidaglia | Weakly damped forced Korteweg-de Vries equations behave as a finite-dimensional dynamical system in the long time[END_REF][START_REF] Ghidaglia | A note on the strong convergence towards attractors of damped forced KdV equations[END_REF][START_REF] Goubet | Asymptotic smoothing effect for weakly damped forced Korteweg-de Vries equations[END_REF][START_REF] Goubet | Asymptotic smoothing and the global attractor of a weakly damped KdV equation on the real line[END_REF][START_REF] Ott | Damping of solitary waves[END_REF][START_REF] Vento | Global well-posedness for dissipative Korteweg-de Vries equations[END_REF][START_REF] Vento | Asymptotic behavior of solutions to dissipative Korteweg-de Vries equations[END_REF][START_REF] Temam | Infinite-dimensional dynamical systems in mechanics and physics[END_REF] and references therein). These account for a wide variety of different results, such as regularizing effect of the damping, asymptotic behavior, existence of attractor or numerical computations. Few literatures are concerned the damped BBM equation [START_REF] Hayashi | Large time asymptotics for the BBM-Burgers equation[END_REF][START_REF] Wang | Strong attractors for the Benjamin-Bona-Mahony equation[END_REF]. This paper consists in establishing theoretical and numerical results for the solution of a generalized damped BBM equation.

We introduce the following equation, called dBBM, for x ∈ T(0, L), L > 0 and t ∈ R u t + u x -u xxt + uu x + L γ (u) = 0, 1 where the operator L γ is defined by its Fourier symbol

L γ (u)(k) := γ k ûk .
Here ûk is the k-th Fourier coefficient of u and (γ k ) k∈Z are positive real numbers chosen such that

T u(x)L γ (u)dµ(x) = k∈Z γ k |û k | 2 ≥ 0.
Standard dampings are included choosing γ k = k 2 as parabolic damping or γ k = γ as weak damping. The proposed sequence (γ k ) k∈Z allows us to consider a greater range of damping. In particular, one may wonder if it is necessary to absorb all frequencies of a long wavelength gravity wave like the solution of KdV or BBM. To obtain this kind of results, we can consider sequences γ k → 0 when |k| → +∞ or even γ k = 0 for large |k|.

In comparison with the standard BBM equation for which the H 1 -norm is preserved, we notice that the H 1 -norm decreases for the solution u of the damped equation dBBM. More precisely, we have for all time t ∈ R

d dt ||u(t)|| 2 H 1 = -|u(t)| 2 γ ≤ 0,
and the natural space to study the well-posedness is the space H γ (T), defined as

H γ (T) :=    u ∈ L 2 (T); k∈Z γ k |û k | 2 < +∞    ,
equipped with the norm

|u| γ := k∈Z γ k |û k | 2 .
In order to simplify the writings, k can denote either an integer as an index, or the value 2πk L . When there is no ambiguity, we denote C the different constants appearing in the following results.

The paper is organised as follows. We establish in Section 2 the local wellposedness of the dBBM equation. We prove some important estimates about the space H γ (T). Some qualitative properties of the solution are also given. In Section 3, numerical schemes are presented to solve the damped equation and to preserve the qualitative properties of the continuous solution. Section 4 deals with the numerical results. In particular, we manage to build a family of dampings weaker than the standard ones. These, for example taking lim |k|→+∞ γ k = 0, yet provide dampen solutions.

Analysis of the problem

After studying the space H γ (T), we establish the well-posedness of the Cauchy problem associated with the dBBM equation.

Proper energy space H γ (T)

Here we take γ k > 0, ∀k ∈ Z. We first state some properties of injection.

Proposition 2.1. Assume that k∈Z 1 γ k < +∞. Then there exists a constant C > 0 such that u ∞ ≤ C|u| γ .
The injection

H γ (T) ֒→ L ∞ (T) is continuous. Proof. Let u ∈ H γ (T). We notice that u(x) = k∈Z û(k)e ikx . Then |u(x)| ≤ k∈Z |û(k)| = k∈Z 1 √ γ k √ γ k |û(k)|.
We assumed that γ k > 0. Hence, the Cauchy-Schwarz inequality implies for all x ∈ T :

|u(x)| ≤   k∈Z 1 γ k   1 2   k∈Z γ k |û(k)| 2   1 2 =   k∈Z 1 γ k   1 2 |u| γ .
This completes the proof.

Remark 2.2. The above result is also a condition to have û l 1 ≤ C|u| γ .

Proposition 2.3. We assume that for all k ∈ Z we have

γ k > β k . Let ρ N = max k≥N β k γ k . Then lim N →+∞ ρ N = 0 if and only if the continuous injection H γ (T) ֒→ H β (T) is compact.
Proof. The condition is necessary, indeed if there exists α > 0 such that ρ N > α, ∀N , then the norms |u| β and |u| γ are equivalent, the injection cannot be compact. Let us prove now that the condition is sufficient. First we have for u ∈ H γ (T) :

|u| β = k∈Z β k |û k | 2 ≤ k∈Z γ k |û k | 2 = |u| γ .
This shows that the injection is continuous. Now we prove that the injection is compact. We use finite rank operators and we take the limit. Let I N be the orthogonal operator on the polynomials of frequencies k such that -N ≤ k ≤ N . We have

I N u = |k|≤N ûk e ikx .
Hence

|(Id -I N )u| 2 β = |k|≥N +1 β k |û k | 2 , ≤ |k|≥N +1 β k γ k γ k |û k | 2 , ≤ ρ N +1 |u| 2 γ -→ N →+∞ 0.
Therefore Id is a compact operator and consequently the injection is compact.

Now we give some conditions on the sequence (γ k ) k∈Z so that H γ (T) is an algebra.

Proposition 2.4. Let u and v be two functions in H

γ (T). Assume that there exists a constant C > 0 such that √ γ k ≤ C √ γ k-j + √ γ j for all k, j ∈ Z.
Then we have the following inequality:

|uv| γ ≤ C (|u| γ v l 1 + |v| γ û l 1 ) , Moreover if k∈Z 1 γ k < +∞ then H γ (T) is an algebra.
Proof. Let u, v ∈ H γ (T). We have

|uv| 2 γ = k∈Z γ k | uv(k)| 2 .
We remind that uv(k) = û ⋆ v(k). We use the inequality

√ γ k ≤ C √ γ k-j + √ γ j .
We obtain for all k ∈ Z

√ γ k | uv(k)| ≤ C   j∈Z √ γ k-j |û(k -j)v(j)| + j∈Z √ γ j |û(k -j)v(j)|   .
Hence

|uv| 2 γ ≤ C 2 k∈Z   j∈Z √ γ k-j |û(k -j)v(j)| + j∈Z √ γ j |û(k -j)v(j)|   2 , ≤ C 2 k∈Z      j∈Z √ γ k-j |û(k -j)v(j)|   2 +   j∈Z √ γ j |û(k -j)v(j)|   2    , ≤ C 2 ( √ γ k |û|) ⋆ |v| 2 l 2 + |û| ⋆ √ γ k |v| 2 l 2 .
We remind that for two functions f and g defined from N to R such that f ∈ l 1 and g ∈ l 2 , we have

|f | ⋆ |g| 2 l 2 ≤ g 2 l 2 f 2 l 1 . Thus |uv| 2 γ ≤ C |u| 2 γ v 2 l 1 + |v| 2 γ û 2 l 1 . We know there exists a constant c > 0 such that û l 1 ≤ c|u| γ if k∈Z 1 γ k < +∞. Then, there exists C > 0 such that |uv| γ ≤ C|u| γ |v| γ .
Remark 2.5. With γ k = k 2s we find the standard Sobolev injection. Indeed, we have the inequality

√ γ k ≤ C √ γ k-j + √ γ j and if s > 1 2 then k∈Z 1 γ k < +∞.

Local well-posedness

We can study the well-posedness of problem dBBM. We consider the Cauchy problem

u t + u x -u xxt + uu x + L γ (u) = 0, x ∈ T, t ∈ [0, T ] (1)
u(x, t = 0) = u 0 (x).

(

) 2
Theorem 2.6. Assume that √ γ k+j ≤ C( √ γ k + √ γ j ) for all k and j in Z, Moreover, for all M > 0 with |u 0 | γ ≤ M and |v 0 | γ ≤ M , there exists a constant C 1 > 0 such that the solutions u and v, associated with the initial data u 0 and v 0 respectively, satisfy for all t ≤ 1

C 0 M |u(•, t) -v(•, t)| γ ≤ C 1 |u 0 -v 0 | γ .
Proof. Let write the initial value problem (1)-(2) as

∂ t u = Au + f (u), (3)
u(t = 0) = u 0 (x), (4)
where

Au = -(1 -∂ xx ) -1 (∂ x u + L γ (u)) , f (u) = -(1 -∂ xx ) -1 ∂ x u 2 2 .
According to the Duhamel's formula, u is solution of the Cauchy problem (3)-( 4) if and only if

u(t) = φ (u(t)) = S t u 0 + t 0 S t-τ f (u(τ )) dτ,
where

S t v = k∈Z e ikx e - (ik+γ k ) 1+k 2 t v(k).
The purpose is to show that u is the unique fixed point of φ. We define the closed ball

B(T ) = {u ∈ C ([0, T ]; H γ (T)) ; |u(t) -u 0 | γ ≤ 3|u 0 | γ } .
Our purpose is to apply the Banach fixed point theorem.

Let u ∈ B(T ). We show that φ (u(t)) ∈ B(T ). From the triangle inequality, we have:

|φ (u(t)) | γ ≤ |S t u 0 | γ + t 0 |S t-τ (f (u)) | γ dτ.
On the one hand, since γ k ≥ 0

|S t u 0 | 2 γ = k∈Z γ k | S t u 0 (k)| 2 , ≤ k∈Z γ k |e - γ k 1+k 2 t u 0 (k)| 2 , ≤ |u 0 | 2 γ ,
and

|S t-τ f (u(τ )) | γ ≤ |f (u(τ )) | γ .
On the other hand, we need to upper bound the term |f (u(τ )) | γ . We have

|f (u(τ )) | 2 γ = | -(1 -∂ xx ) -1 ∂ x u 2 2 (τ ) | 2 γ , = k∈Z γ k F -(1 -∂ xx ) -1 ∂ x u 2 2 (τ ) 2 , = k∈Z γ k ik 1 + k 2 u 2 2 2 . However, since |k| 1+k 2 ≤ 1 and H γ (T) is an algebra, it gets |f (u(τ ))| γ ≤ C |u| 2 γ . The Duhamel's formula implies for 0 ≤ t ≤ T |φ (u(t)) | γ ≤ |u 0 | γ + t 0 |u(τ )| 2 γ dτ, ≤ |u 0 | γ + C sup t∈[0,T ] |u(t)| γ 2 T. But u ∈ B(T ), so |u(t)| γ -|u 0 | γ ≤ |u(t) -u 0 | γ ≤ 3|u 0 | γ , that implies |u(t)| γ ≤ 4|u 0 | γ . We have φ (u(t)) ∈ B(T ) if |φ (u(t)) -u 0 | γ ≤ 3|u 0 | γ . The inequality |φ (u(t) -u 0 ) | γ ≤ 2|u 0 | γ + CT 16|u 0 | 2 γ ≤ 3|u 0 | γ , is true if we choose 0 < T ≤ |u 0 | γ 16C|u 0 | 2 γ = 1 16C 1 |u 0 | γ . Now we show that φ is a contraction mapping on B(T ). Let u, v ∈ B(T ). We have |φ (u(t)) -φ (v(t))| γ = t 0 S t-τ (f (u) -f (v)) dτ γ , ≤ t 0 u 2 2 - v 2 2 γ (τ )dτ. We notice that u 2 -v 2 = (u -v)(u + v). Since H γ (T) is an algebra, we have u 2 -v 2 γ ≤ C |u -v| γ |u + v| γ , ≤ C |u| γ + |v| γ |u -v| γ .
And for u, v ∈ B(T ),

u 2 -v 2 γ ≤ 8C |u 0 | γ |u -v| γ .
We infer for all t ∈ [0, T ]

|φ (u(t)) -φ (v(t))| γ ≤ 8C |u 0 | γ t 0 |u -v| γ (τ )dτ, ≤ 8C |u 0 | γ T sup t∈[0,T ] |(u -v)(t)| γ . Hence sup t∈[0,T ] |φ (u(t)) -φ (v(t))| γ ≤ 8C |u 0 | γ T sup t∈[0,T ] |u(t) -v(t)| γ . Consequently φ is a contraction mapping if 8C |u 0 | γ T < 1, i.e., T < 1 8C|u 0 | γ .
Finally, from the Banach fixed-point theorem, φ has a unique fixed-point solution of u(t) = φ (u(t)). So, there exists a unique solution of the Cauchy problem.

It remains to prove the continuity with the initial data. Let u and v solutions of the Cauchy problem ( 1)-( 2) with initial data u 0 and v 0 respectively, such that

|u 0 | γ ≤ M and |v 0 | γ ≤ M . The Duhamel's formula gives for t ∈ [0, T ], T ≤ 1 C 0 M |u -v| γ ≤ |u 0 -v 0 | γ + t 0 |f (u) -f (v)| γ dτ, ≤ |u 0 -v 0 | γ + C ′ T |u 0 | γ + |v 0 | γ sup t∈[0,T ] |u -v| γ . It implies |u -v| γ ≤ C 1 |u 0 -v 0 | γ .

Behavior of the solution

In this part, we aim to get some estimations of the damping rate. Actually in the case without damping, the H 1 -norm of the solution is invariant during time. Here this norm is decreasing. We adapt the work done on KdV equation [START_REF] Chehab | On damping rates of dissipative KdV equations[END_REF].

Let us begin with the linear equation, that reads

u t -u txx + u x + L γ (u) = 0 x ∈ T, t ∈ [0, T ]. (5)
Let u be valued in L 2 (T), ∀t > 0. We write u as a Fourier series and, due to the orthogonality of the trigonometric polynomials, we obtain

(1 + k 2 ) dû k (t) dt + (γ k + ik)û k (t) = 0. Hence ûk (t) = e - γ k +ik 1+k 2 t ûk (0).
It follows that

|u| 2 H 1 = k∈Z (1 + k 2 )|û k (t)| 2 = k∈Z (1 + k 2 )e - 2γ k 1+k 2 t |û k (0)| 2 . Proposition 2.7. Let γ k > 0, ∀k ∈ Z and u 0 ∈ H δ (T) where δ k = (1+k 2 ) 2 γ k
. Then the unique solution u of (5) verifies

|u| 2 H 1 ≤ min e -1 2t |u 0 | 2 δ , |u 0 | 2 H 1 , ∀t > 0.
More generally, if u 0 ∈ H βδ (T) then

|u| 2 β ′ ≤ e -1 2t |u 0 | 2 βδ where β ′ k = (1 + k 2 )β k .
Proof. On one hand, the scalar product in L 2 (T) of ( 5) with u provides 1 2

d dt |u| 2 H 1 + k∈Z γ k |û k (t)| 2 = 0. Hence d dt |u| 2 H 1 ≤ 0. Consequently |u| 2 H 1 ≤ |u 0 | 2 H 1 .
On the other hand, to obtain the second inequality, we write

|u| 2 H 1 = k∈Z (1 + k 2 )|û k | 2 , = k∈Z (1 + k 2 )e - 2γ k 1+k 2 t |û k (0)| 2 , = k∈Z (1 + k 2 ) 2 γ k 1 + k 2 e - 2γ k 1+k 2 t 1 γ k |û k (0)| 2 .
Since the function β → βe -2βt is uniformly bounded by e -1 2t , we infer that

|u| 2 H 1 ≤ e -1 2t |u 0 | 2 δ .
In the same way we obtain the third inequality from

β k |û k | 2 = β k e - 2γ k 1+k 2 t |û k (0)| 2 = (1 + k 2 ) γ k 1 + k 2 e - 2γ k 1+k 2 t β k γ k |û k (0)| 2 .
We can prove the two following results with similar arguments:

Proposition 2.8. Assume that γ k ∈ [0, 1], ∀k ∈ Z and u 0 ∈ H δ (s)
, where δ

(s) k = (1+k 2 ) s+1 γ s k
. Then for all s > 0, u the solution of (5) verifies ∀t > 0

|u| 2 H 1 ≤ min e -s s 2t s |u 0 | 2 δ (s) , |u 0 | 2 H 1 .
Proposition 2.9. We assume that there exist three positive constants α, β and C such that u the solution of (5) verifies:

i. | u 0k | 2 ≤ Cγ 2δ k with δ = α + β.
ii.

k∈Z

(1 + k 2 ) 2α+1 γ 2β k < +∞.
Then

|u| 2 H 1 ≤ Ce -2a α t 2α k∈Z (1 + k 2 ) 2α+1 γ 2β k = O 1 t 2α .
We can now deal with the non-linear equation. We can find similar kind of decreasing but less explicit than in the linear case. We remind the equation

u t -u txx + u x + uu x + L (u) = 0 (6)
The scalar product in L 2 of ( 6) with u yields 1 2

d dt |u| 2 H 1 + |u| 2 γ = 0.
Proposition 2.10. We assume that i.

γ k > 0 ∀k ∈ Z, ii. u 0 ∈ H γ (T) ∩ H 1 (T) with L 0 u 0 (x)dx = 0, Then lim t→+∞ |u| H 1 = 0.
Proof. Because of γ k > 0, we have

v = 1 L L 0 v(x, t)dx.
Let u be the solution of dBBM valued in H γ (T) ∩ H 2 (T). Then ū(t) = e -γ 0 t ū(0).

Proof. We integrate the equation in space on the interval [0, L] and we obtain

L 0 ∂u ∂t + L γ (u) dx = 0. But L 0 L γ (u)(x, t)dx = L 0 k∈Z γ k ûk e 2iπkx L dx = Lγ 0 û0 = γ 0 L 0 u dx.
Hence we have

d dt   1 L L 0 u(x, t)dx   + 1 L L 0 L γ (u)(x, t)dx = 0, Therefore dū dt + γ 0 ū = 0.
We deduce the desired result.

The solution converges to 0 in H 1 with a rate depending on the H γ -norm. Thus, as in [START_REF] Chehab | On damping rates of dissipative KdV equations[END_REF], we introduce the following ratio function

G : (u, t) → G(u, t) = |u| γ |u| H 1 = k∈Z γ k |û k | 2 k∈Z (1 + k 2 )|û k | 2 . ( 7
)
To simplify the writings, we use G(t) instead of G(u, t) when no confusion is possible.

Proposition 2.12. Let u be the unique solution of dBBM valued in H γ (T)∩ H 1 (T). We assume that G is C 1 in time. Then we have the following equalities.

i.

|u(t)| 2 H 1 = e -2 t 0 G 2 (s)ds |u 0 | 2 H 1 , ii. |u(t)| 2 γ = G 2 (t)e -2 t 0 G 2 (s)ds |u 0 | 2 H 1 . In particular, lim t→+∞ |u| 2 H 1 = 0 if and only if t → G(t) / ∈ L 2 t (0, +∞).
Proof. As previously, we take the scalar product in L 2 (T) of the equation dBBM with u and we obtain 1 2

d dt |u| 2 H 1 + |u| 2 γ = 0. Since |u| 2 γ = G 2 (t)|u| 2 H 1 , we have 1 2 d dt |u| 2 H 1 + G 2 (t)|u| 2 H 1 = 0 thus |u(t)| 2 H 1 = e -2 t 0 G 2 (s)ds |u 0 | 2 H 1 .
We

G 2 (t)e -2 t 0 G 2 (s)ds |u 0 | 2 H 1 = 0.
Remark 2.13. We can establish similar results for the forced equation

u t -u txx + u x + uu x + L γ (u) = f. Proposition 2.14. Assume that f ∈ H 1/γ (T) and u 0 ∈ H 1 (T) ∩ H γ (T).
Then the solution u of the forced equation satisfies:

|u(t)| 2 H 1 ≤ e - t 0 G 2 (s)ds |u 0 | 2 H 1 + L 0 e - t s G 2 (τ )dτ |f | 2 1/γ ds.
Proof. The scalar product of dBBM with u gives 1 2

d dt |u| 2 H 1 + |u| 2 γ =< f, u > . But | f, u | ≤ |f | 1/γ |u| γ and |u| 2 γ = G 2 (t)|u| 2 H 1 .
From the Young inequality, we have 1 2

d dt |u| 2 H 1 + G 2 (t)|u| 2 H 1 ≤ 1 2 |f | 2 1/γ + 1 2 |u| 2 γ . It involves that d dt |u| 2 H 1 + G 2 (t)|u| 2 H 1 ≤ |f | 2 1/γ . Using the Gronwall's lemma, we obtain |u(t)| 2 H 1 (T) ≤ e - t 0 G 2 (s)ds |u 0 | 2 H 1 + L 0 e - t s G 2 (τ )dτ |f | 2 1/γ ds.
Multiplying each term with G 2 (t), it gets

|u(t)| 2 γ (T) ≤ G 2 (t)e - t 0 G 2 (s)ds |u 0 | 2 H 1 + L 0 e - t s G 2 (τ )dτ G 2 (t)|f | 2 1/γ ds.

Numerical schemes

In this part, we study and present some numerical schemes well suited for the long-time behavior of the solution of dBBM equation. We also give details on their implementation

Numerical schemes

Let start with time discretization. We denote by D the operator of derivation in space, by D 2 the operator of second derivation in space and by u n the approximation of u at time t n = n∆t.

We first introduce the fully bacward Euler method. We have the following discretization:

(Id -D 2 ) u n+1 -u n ∆t + Du n+1 + 1 2 D (u n+1 ) 2 + L γ u n+1 = f. ( 8
) Proposition 3.1. If u 0 = u 0 ∈ H 1 (T) and f ∈ H 1/γ (T)
then the sequence (u n ) n∈N generated by the fully backward Euler method is well defined, in H 1 (T) and

|u n+1 | 2 H 1 + |u n+1 -u n | 2 H 1 + ∆t k∈Z γ k |û n+1 k | 2 ≤ |u n | 2 H 1 + ∆t k∈Z 1 γ k | fk | 2 . Moreover, if f = 0 then lim n→+∞ |u n | H 1 = 0.
Proof. We use the following identity

(Id -D 2 )(u n+1 -u n ), u n+1 = - 1 2 |u n | 2 H 1 -|u n+1 -u n | 2 H 1 -|u n+1 | 2 H 1 .
We compute the scalar product of (8) with u n+1 and we obtain

|u n+1 | 2 H 1 + |u n+1 -u n | 2 H 1 -|u n | 2 H 1 + 2∆t L γ u n+1 , u n+1 = 2∆t u n+1 , f , where L γ u n+1 , u n+1 = k∈Z γ k |û n+1 k | 2 ≥ 0.
Hence we have

|u n+1 -u n | 2 H 1 + |u n+1 | 2 H 1 + 2∆t k∈Z γ k |û n+1 k | 2 = |u n | 2 H 1 + u n+1 , f .
The Young inequality on the right hand side provides

u n+1 , f ≤ 1 2 k∈Z γ k |û n+1 k | 2 + 1 2 k∈Z 1 γ k | fk | 2
We infer the expected result

|u n+1 | 2 H 1 + |u n+1 -u n | 2 H 1 + ∆t|u n+1 | 2 γ ≤ |u n | 2 H 1 + ∆t k∈Z 1 γ k | fk | 2
If f = 0, we directly have the equality

|u n+1 | 2 H 1 + |u n+1 -u n | 2 H 1 + ∆t|u n+1 | 2 γ = |u n | 2 H 1
We deduce that the sequence (|u

γ k |û n k | 2 = 0.
As γ k > 0, we have lim n→+∞ ûn k = 0 and consequently C = 0.

Proposition 3.2. Let (u n ) n∈N be the sequence generated by the fully backward Euler method. Assume that f = 0. We set

G (n) = |u n |γ |u n | H 1 . Then we have |u n | 2 H 1 ≤   n j=1 1 1 + 2∆t(G (j) ) 2   |u 0 | 2 H 1 .
Moreover, if ∆t is small enough and if (G

(j) ) j∈Z / ∈ l 2 then lim n→+∞ |u n | H 1 = 0.
Proof. We compute the scalar product of ( 8) with u n+1 and we obtain

1 2∆t |u n+1 | 2 H 1 -|u n | 2 H 1 + |u n+1 -u n | 2 H 1 + (G (n+1) ) 2 |u n+1 | 2 H 1 = 0. Hence 1 + 2∆t(G (n+1) ) 2 |u n+1 | 2 H 1 + |u n+1 -u n | 2 H 1 ≤ |u n | 2 H 1 .
Therefore

|u n+1 | 2 H 1 ≤ |u n | 2 H 1 1 + 2∆t(G (n+1) ) 2 .
We infer by induction that

|u n | 2 H 1 ≤   n j=1 1 1 + 2∆t(G (j) ) 2   |u 0 | 2 H 1 .
Now if ∆t is small enough, we have log

  n j=1 1 1 + 2∆t(G (j) ) 2   = - n j=1 log(1 + 2∆t(G (j) ) 2 ) ≃ -2∆t n j=1 (G (j) ) 2 . If (G (j) ) / ∈ l 2 then lim n→+∞ |u n | H 1 = 0.
Remark 3.3. The forward Euler method can be applied to the BBM equation, contrarily to the KdV equation. Moreover, it is more difficult to establish qualitative properties as Propositions 3.1 and 3.2. The scheme is written as

(Id -D 2 ) u n+1 -u n ∆t + Du n + 1 2 D (u n ) 2 + L γ u n = f.
Let introduce the Sanz-Serna scheme, given by

(Id -D 2 ) u n+1 -u n ∆t + D u n+1 + u n 2 + 1 2 D   u n+1 + u n 2 2   + L γ u n+1 + u n 2 = f. ( 9
)
Proposition 3.4. Assume that u 0 ∈ H 1 (T) and f ∈ L 2 ∩ H 1/γ (T). Then the scheme ( 9) is stable in H 1 (T) for all ∆t > 0.

Proof. We take the scalar product of ( 9) with u n+1 +u n

2

. We obtain

|u n+1 | 2 H 1 -|u n | 2 H 1 2∆t + 1 4 |u n+1 + u n | 2 γ = f, u n+1 + u n 2 .
Using the duality and the Young's inequality, it implies

|u n+1 | 2 H 1 -|u n | 2 H 1 2∆t + 1 4 |u n+1 + u n | 2 γ ≤ 1 2 |f | 2 1/γ + 1 2 u n+1 + u n 2 2 γ . Hence |u n+1 | 2 H 1 + ∆t 4 |u n+1 + u n | 2 γ ≤ |u n | 2 H 1 + ∆t|f | 2 1/γ .
Consequently we have the stability on every time interval [0, T ].

Remark 3.5. The Crank-Nicolson scheme is given by

(Id -D 2 ) u n+1 -u n δt + D u n+1 + u n 2 + 1 4 D (u n+1 ) 2 + (u n ) 2 + L γ u n+1 + u n 2 = f. ( 10
)
We cannot have uniform H 1 -bounds for u n from this scheme because

D (u n+1 ) 2 + (u n ) 2 , u n+1 + u n 2 = 0.

Implementation

The Fourier Transform in space is applied to the dBBM equation. The equation ( 6) becomes

(1 + k 2 ) dû k dt + (ik + γ k )û k + ( uu x ) k = 0.
We remind that the term L γ (u) is written in the Fourier space as

k∈Z γ k ûk e ikx .
Then we use the schemes introduced in the previous subsection for the discretization in time and the FFT for the space discretization. These are implicit methods and we apply a fixed-point method. We solve equations which can be written as u n+1 = Φ(u n , u n+1 ) by Picard iterates.

Algorithm 1 Picard iterates v 0 = u n , m = 0 while Φ(u n , u m ) -u m > ǫ do v m+1 = Φ(u n , v m ) m = m + 1 end while u n+1 = v m
Unfortunately, even if the used scheme is unconditionally stable, this method presents instabilities. The stability of the fixed-point iterate can be improved by applying extrapolation like technic as follows (see also [START_REF] Abounouh | Discrete Schrödinger equations and dissipative dynamical systems[END_REF]). The principle is to change the iterate v m+1 = Φ(u n , v m ) by j) denotes the j-th composition of Φ with itself. The parameter α k m is computed such as minimizing the Euclidean norm of the linearized part of the residual r m+1 = v m+1 -Φ(u n , v m ). We have

v m+1 = v n -(-1) k α k m ∆ k Φ v m , where ∆ k Φ v m = k j=0 C k j (-1) k-j Φ (j) (u n , v m ), Φ (
α k m = (-1) k ∆ 1 Φ v m , ∆ k+1 Φ v m ∆ k+1 Φ v m , ∆ k+1 Φ v m .
In fact, we take k = 1, and we have the following algorithm Algorithm 2 Improved fixed-point algorithm

v 0 = u n , m = 0 while Φ(u n , v m ) -v m > ǫ do ∆ 1 Φ v m = Φ(u n , v m ) -v m ∆ 2 Φ v m = Φ(u n , φ(u n , v m )) -2Φ(u n , v m ) + v m α 1 m = - ∆ 1 Φ v m ,∆ 2 Φ v m ∆ 2 Φ v m ,∆ 2 Φ v m v m+1 = v m + α 1 m ∆ 1 Φ v m m = m + 1 end while u n+1 = v m

Numerical results

Let us begin by stating the parameters chosen for the simulations. We present the results obtained with the Sanz-Serna scheme for the non-forced equation (f = 0) on the interval [-L, L], with L = 50, discretized in 2 10 points. We obtain similar results with fully backaward Euler and Crank-Nicolson schemes. The time step is ∆t = 0.01. Let us remind that k can denote either an integer as an index, or the value 2πk L .

- We consider the soliton as initial datum

u 0 (x) = -3(1 -c) cosh -2 ± 1 2 c -1 c (x -d) ,
with c = 1.5 and d = 0.2L . We begin with a laplacian damping in Figure 1. The figure shows the effect of the damping on the solution. Since γ 0 = 0, the solution tends to the mean value of u 0 . Figure 2 corresponds to the simulation of the following standard dampings

• the constant: γ k = 1, ∀k,
• the laplacian:

γ k = k 2 , ∀k,
• the bilaplacian: γ k = k 4 , ∀k. Here the initial datum is the soliton.

For these dampings, we verify that the H 1 -norm decreases as expected. Moreover the solution converge to the norm value of u 0 for the laplacian and bilpalacian dampings. We can also classify these dampings. Here the constant one is more efficient than the two other. But the laplacian one is more efficient than the bilaplacian one. These results are in agreement with the ones about the KdV equation [START_REF] Chehab | Numerical study of a family of dissipative KdV equations[END_REF][START_REF] Ghidaglia | Weakly damped forced Korteweg-de Vries equations behave as a finite-dimensional dynamical system in the long time[END_REF][START_REF] Ghidaglia | A note on the strong convergence towards attractors of damped forced KdV equations[END_REF][START_REF] Goubet | Asymptotic smoothing effect for weakly damped forced Korteweg-de Vries equations[END_REF][START_REF] Goubet | Asymptotic smoothing and the global attractor of a weakly damped KdV equation on the real line[END_REF][START_REF] Ott | Damping of solitary waves[END_REF].

We extend, in Figure 3, the simulations to the dampings such that γ k > 0 with lim k→+∞ γ k = 0. We consider here three "tend to 0" dampings, with polynomial or exponential decay as follows

• γ k = 1 1+|k| , • γ k = 1 (1+|k|) α , α = 3 • γ k = e -|k| .
These dampings do not fit with the asymptions of the theorem 2.6 but they still damp down the soliton. Indeed we notice that the H1 -norm is decreasing. We can classify the dampings. Precisely, the first one damps more than the two other and the third one is more efficient than the second one. Since these dampings tends to 0, their effects are principally on low frequencies. And for low frequencies, we have

1 1 + |k| > e -|k| > 1 (1 + |k|) α .
Then the classification of dampings in Figure 3 seems licit. We now consider band limited dampings, namely

γ k = 1 if |k| ≤ N 0 if |k| > N . (I(N ))
We define M := max(k), where k is the finite vector of frequencies used for the numerical tests in Figure 4. To keep γ k > 0, ∀k, we choose

γ k = 1 if |k| ≤ N e -a||k|-N | if |k| > N . (EI(N ))
or

γ k = 1 if |k| ≤ N 1 (1+||k|-N |) a if |k| > N . (P I(N ))
We denote the classic band limited damping by I(N ), the one corrected with exponential decay EI(N ) and the third one corrected with polynomial decay P I(N ). Here we set a = 3, we obtain γ k > 0 with lim k→+∞ γ k = 0. The difference between these dampings is when |k| > N . For these frequencies, since γ k = 0, the dampings with exponential or polynomial decay should damp more than the ideal one (I(N )). It is the case when N = M 16 , moreover the damping with polynomial decay is more effective than the one with exponential decay. But we observe similar damping when N = N 4 . We notice the damping depends on the band width. It seems rightful because the equation considered is dispersive. Indeed, for this kind of equation, high frequencies appear and these frequencies cannot be damped especially if the band width is too small.

Let us illustate this by considering now the gaussian as initial datum

u 0 (x) = e -x 2 σ 2 .
Here the standard deviation is fixed equal to σ = max(k) 10 . Then, in order to see the effect of the damping on the low frequencies, we study the influence of the band limited damping size in Figure 5. We observe that the more N increases, the more efficient is the damping. However, it seems that all dampings with N ≥ 3σ 2 behave similarly. Besides, for lower N , the different kind of dampings can be classified. We notice that for all N and all dampings, the decreasing speed is similar at beginning and the difference appears after. The difference appears when the high frequencies does. Figure 6 presents results with L equal to π, 2π and 3π. We notice that if L < 2π, the laplacian damping is less effective than the bilaplacian one. It is the contrary if L > 2π. But, the dampings are similar when L = 2π. We also observe that for a fixed damping, the greater is L, the slower is the damping.

Final comments

The damping operator L γ considered in this article allows to have a large variety of dampings. We first get back the standard dampings (for example -∆u). But we also have less common dampings like band-limited ones or weaker ones, e.g., such that lim |k|→+∞ γ k = 0. Moreover, with this frequential approach, we can adjust the dampings by frequency bands. This is interesting in order to build cheap and efficient dampings. Using these dampings, the numerical observation shows a damping for energy norms like the H 1norm. This is consistent with the results about the asymptotic behavior.

k∈Z 1 γ

 1 k < +∞ and u 0 ∈ H γ (T). Then there exists T = 1 C 0 |u 0 | γ > 0 and there exists a unique u ∈ C ([0, T ], H γ (T)) solution of the Cauchy problem (1)-(2).

Figure 1 :

 1 Figure 1: Solution at times t= 0 (-), t = 5 (• -•), t = 10 (---), t = 20 (• • • ), t = 30 (-) for γ k = k 2 .Here the initial datum is the soliton.

Figure 2 :

 2 Figure 2: Evolution of u H 1 and G(t) with respect to time for standard dampings, the constant (-), the laplacian (---) and the bilaplacian (• • • ).Here the initial datum is the soliton.

Figure 3 :

 3 Figure 3: Evolution of u H 1 and G(t) with respect to time for "tend to 0" dampings, γ k = 1 1+|k| (-), γ k =

16 Figure 4 : 4 (

 1644 Figure 4: Evolution of u H 1 and G(t) with respect to time for band limited type dampings, I M 4 (-), P I M 4 (-), EI M 4 (-), I M 16 (• -•),P I M 16 (---) ,EI M 16 (• • • ).Here the initial datum is the soliton.

2 Figure 5 : 4 (

 254 Figure 5: Evolution of u H 1 and G(t) with respect to time for band limited type dampings, I 3σ 4 (• -•), P I 3σ 4 (---), EI 3σ 4 (• • • ), I 3σ 2 (-),

  = π bilaplacian, L = π laplacian, L = 2π bilaplacian, L = 2π laplacian, L = 3π bilaplacian, laplacian, L = π bilaplacian, L = π laplacian, L = 2π bilaplacian, L = 2π laplacian, L = 3π bilaplacian, L = 3π

Figure 6 :

 6 Figure 6: Evolution of u H 1 and G(t) with respect to time for L = π withγ k = k 2 (• -•) and γ k = k 4 (---), L = 2π with γ k = k 2 (-) and γ k = k 4 (-), L = 3π with γ k = k 2 (• • • ) and γ k = k 4 (• -•).Here the initial datum is a sine.

  We notice that we also have |u|H 1 ≤ |u 0 | H 1 . Since u ∈ H 1 ∩ H γ ,we deduce that lim For all function v smooth enough, we set

	d|u| 2 H 1 dt H 1 is decreasing and consequently lim |u| 2 0. Finally, since γ k > 0, we have lim t→+∞ ûk = 0 and then u = 0, which means ≤ 0 and then the function t → t→+∞ γ = |u(t)| 2
	that C = 0.
	Lemme 2.11.

t→+∞ |u(t)| 2 H 1 = C.

  find lim

	t→+∞ Then, we have |u| 2 H 1 = 0 if and only if t → G(t) / ∈ L 2 t (0, +∞).
		t
	|u(t)| 2 γ = G 2 (t)e -2	0	G 2 (s)ds	|u 0 | 2 H 1 .
	The proposition 2.10 involves that lim t→+∞	|u| γ = 0 and consequently
	lim t→+∞

  n | H 1 ) n∈N is decreasing and lower bounded by 0. Then it is convergent to a positive constant C. It follows that

	lim n→+∞ k∈Z

(1+|k|) 3 (• -•) and γ k = e -|k| (---).Here the initial datum is the soliton.

This way of damping with frequency filters is particulary of interest for its flexibility and its efficiency. It would be a good perspective to perform similar work for the control wave equations.