Ye Liu
email: ye.liu1@alcatel-lucent.com

Chung Shue Chen
email: cs.chen@alcatel-lucent.com

Chi Wan Sung

Distributed Enhanced Inter-Cell Interference Coordination (eICIC) in LTE-Advanced HetNets: A Potential Game Approach

Keywords: LTE-A, enhanced inter-cell interference coordination (eICIC), almost blank subframe (ABS), cell selection bias (CSB), heterogeneous networks

In this paper we propose a distributed algorithm for jointly optimizing almost blank subframe (ABS) and cell selection bias (CSB) patterns in Long Term Evolution-Advanced (LTE-A) heterogeneous networks (HetNets). We formulate the optimization problem as an exact potential game, where a Nash equilibrium point is guaranteed to be achieved within finite number of plays. Through simulations, we are able to demonstrate the fast convergence of the algorithm, an increase in average user rate, and a tremendous improvement on the service fairness of the users.

I. INTRODUCTION

Heterogeneous cellular networks become an attractive solution to the exponentially growing data traffic of today's mobile networks. In frequency reuse-1 settings, however, the strong interference from the macrocells can severely degrade the reception quality of the users attached to the picocells. Besides, users who are located near and can attach to a small cell may be attached to a macrocell because the reference signal received power (RSRP) from the macrocell can be much larger than that from the picocell, which may lead to underutilization of the picocell. To tackle the aforementioned challenges in heterogeneous networks (HetNets) and fully utilize the power of picocells, 3rd Generation Partnership Project (3GPP) proposes enhanced inter-cell interference coordination (eICIC), where almost blank subframe (ABS) and cell selection bias (CSB) patterns can be employed to tune the network and optimize the radio resource utilization [START_REF] Lopez-Perez | Enhanced intercell interference coordination challenges in heterogeneous networks[END_REF].

Several algorithms of optimizing the ABS and/or CSB settings have been proposed. Tall et al. propose an algorithm that optimizes ABS and CSB separately in [START_REF] Tall | Self organizing strategies for enhanced ICIC (eICIC)[END_REF], where the ABS patterns are simplified as fractional numbers. A centralized algorithm that jointly optimizes ABS and CSB patterns is proposed by Deb et al. in [START_REF] Deb | Algorithms for enhanced inter-cell interference coordination (eICIC) in LTE HetNets[END_REF], where the interfering macrocells of a picocell offer ABSs on the same subframes. Pang et al. discuss a distributed algorithm to determine the number of ABSs in [START_REF] Pang | Optimized time-domain resource partitioning for enhanced inter-cell interference coordination in heterogeneous networks[END_REF], yet they do not consider CSB optimization. Simsek et al. propose a learning algorithm that dynamically optimize CSB in frequency domain in [START_REF] Simsek | Enhanced intercell interference coordination in HetNets: Single vs. multiflow approach[END_REF] and extend the idea to optimize CSB in both time and frequency domain in [START_REF]Learning based frequency-and time-domain inter-cell interference coordination in HetNets[END_REF].

In this paper, we provide a distributed algorithm that jointly optimizes the ABS and CSB patterns in HetNets. We formulate the problem as an exact potential game, where the cells in the network are the players and the ABS and/or CSB patterns are the strategies. A decentralized eICIC algorithm is therefore derived, based on the given potential game. We show that such an approach can be used to maximize various objectives, such as the sum of users' throughput and the sum of logarithmic utilities of users' throughput. Through simulation studies, we show that the algorithm can quickly converge, and a 50% improvement on average user throughput can be achieved by optimizing the ABS and CSB patterns. Meanwhile, the service fairness among users can also be greatly improved.

The rest of the paper is organized as follows: Section II presents the system model of Long Term Evolution-Advanced (LTE-A) HetNets and the joint ABS and CSB optimization problem. Section III derives the exact potential game formulation. Section IV discusses the required information to perform the optimization. Section V shows the simulation results and Section VI gives the conclusion.

II. SYSTEM MODEL

Consider a LTE-A HetNet where macrocells and picocells coexist. Let M be the set of all macrocells and P be the set of all picocells. In the sequel a station may either be referred to as a macrocell or a picocell. Let U i be the set of users associated with station i, where i ∈ M ∪ P. The |M||P|by-1 vector γ specifies the CSB values of all macrocells and picocells, where γ i , the CSB value of station i, belongs to a pre-defined set C. Let P i Rx,u be the RSRP of user u from station i, where i ∈ M ∪ P and the received power depends on the transmission power of the station, the distance between the station and the UE, and the shadowing loss. User u will be attached to the following station

g(u, γ) arg max i (P i Rx,u + γ i ). (1)
Assume there are T subframes in the time domain and F resource blocks (RBs) in the frequency domain. Each subframe and RB pair forms a physical resource block (PRB), where a subframe has a duration of 1ms and a RB has a bandwidth of 180 kHz. Let B := T •F be the total number of PRBs available at each station. All macrocells transmit at a fixed power level on PRBs which correspond to non-ABSs, and all picocells transmit at a fixed power level on all PRBs. All subframes used by all stations are assumed to be synchronized.

The achieved rate of user u at PRB b is denoted as r u,b , where b ∈ [1, B]. Due to the small-scale fading and potentially different ABS patterns of the interfering macrocells, a user can get different throughput on different PRBs. Let α m be the Tby-1 binary vector specifying the ABS pattern of macrocell m, where α m,i = 0 means the i-th subframe of macrocell m is an ABS and α m,j = 1 means the j-th subframe of macrocell m is a non-ABS, respectively. The Signal-to-Noiseplus-Interference Ratio (SINR) of a user u associated with an macrocell m on PRB b is given as: is the interference from picocells in P \ {p}. Given the above SINR definitions and user association, the throughput at PRB b of a user u, namely r u,b , can be calculated by Shannon's capacity formula.

SINR m u,b = h m u,b • P m Rx,u • α m,τ (
Assume x u,b is a binary value specifying whether user u occupies the b-th PRB of its serving station, where x u,b = 1 and x u,b = 0 mean user u occupies and does not occupy the b-th PRB, respectively. Let w u be a nonnegative weighting factor of UE u. The optimization problem can be stated as:

maximize i∈M∪P u∈Ui w u • ln B b=1 x u,b • r u,b , (2) subject to u∈Um x u,b = α m,τ (b) , ∀m ∈ M, b ∈ [1, B], α m ∈ A, (3)
u∈Up x u,b = 1, ∀p ∈ P, b ∈ [1, B], (4)
x u,b ∈ {0, 1}, ∀u ∈ U , b ∈ [1, B], (5)
γ(i) ∈ C, ∀i ∈ M ∪ P, (6)
where A is the set of vectors from which macrocells can select their ABS pattern. Note that we choose to maximize the proportional fair (PF) objective, i.e., sum of log of user rates, because it strikes a balance between network throughput and user fairness (e.g., see [START_REF] Hou | Self-organized resource allocation in LTE systems with weighted proportional fairness[END_REF]).

III. DISTRIBUTED SOLUTION USING POTENTIAL GAME

In this section, we aim at developing a distributed solution to the problem described in Section II by formulating an exact potential game.

A. Exact Potential Game Formulation

In an exact potential game, there exists a potential function such that the change in the potential function due to the change of a player's strategy equals to the change of the utility of that player. A Nash equilibrium can be reached within a finite number of plays if a player is randomly chosen to update its strategy so that the player's utility is maximized in each step [START_REF] Neel | Analysis and Design of Cognitive Radio Networks and Distributed Radio Resource Management Algorithms[END_REF]. In [START_REF] Singh | Distributed downlink resource allocation in cellular networks through spatial adaptive play[END_REF], an exact potential game is formulated to optimize parameters in LTE, where the users in the system are the players. In the following, we show that the same idea can be used in the current context where a potential game is formulated by using the stations as the players. After all, it not practical to let the users to control the cellular network.

Let the stations be the players, so that M ∪ P defines the set of players. In the sequel we use station and player interchangeably. Each player i chooses a strategy s i from a set of strategies S i . The underlying assumption is that ABS is allowed in macrocells but not in picocells, and all stations can set their CSB values according to C. Let vector s specify the strategies of all players, where s(i) gives the strategy of station i. Define the strategies selected by all players except player i as s -i (s 1 , ..., s i-1 , s i+1 , ..., s |M|+|P| ), and define

(s ′ i , s -i ) (s 1 , ..., s i-1 , s ′ i , s i+1 , ..., s |M|+|P| )
so that player i adopts strategy s ′ i and other players' strategies are specified by s -i .

Let N i {i} ∪ N Int i ∪ N Att i , where i ∈ M ∪ P, N Int i is the set of stations whose attached users can be interfered by the transmissions of station i, N Att i is the set of stations whose user attachment can be affected by the change of CSB value of station i. Let the utility of station i ∈ M ∪ P be

U i (s) u∈Ui w u • ln B b=1 x u,b • r u,b , (7)
where x u,b is obtained by some scheduling scheme given s so that the constraint in (3), ( 4) and ( 5) are satisfied. The utility of the whole system is

U (s) = i∈M∪P U i (s), (8)
Let the payoff function of player i be

V i (s) j∈Ni U j (s).
We have,

U (s ′ i , s -1 ) -U (s) = j∈M∪P U j (s ′ i , s -1 ) -U j (s) = j∈Ni (U j (s ′ i , s -1 ) -U j (s)) + j∈M∪P\Ni (U j (s ′ i , s -1 ) -U j (s)) = j∈Ni (U j (s ′ i , s -1 ) -U j (s)) (9) = V i ((s ′ i , s -i )) -V i (s), (10)
where ( 9) follows from the fact that the change of player i's strategy will only affect the utilities of players in N i . Eqn. [START_REF] Chen | Round-robin power control for the weighted sum rate maximisation of wireless networks over multiple interfering links[END_REF] shows that U (•) is a potential function of the game (M ∪ P, (S i , i ∈ M ∪ P), (V i , i ∈ M ∪ P)).

Note that the utility function in ( 7) can be modified to other functions without affecting the validity of the potential game formulation, as long as the utility function depends on the users who are attached to a single station. For example, it can be u∈Ui w u • B b=1 x u,b •r u,b , or in other words the weighted sum of user rates [START_REF] Chen | Round-robin power control for the weighted sum rate maximisation of wireless networks over multiple interfering links[END_REF].

The procedure of game play is given as follows:

1) Initialize the ABS and CSB patterns of all macrocells and picocells. 2) Randomly select a player i ∈ M ∪ P. The player then searches through all strategies in S i and chooses a strategy so that V i is maximized. 3) Repeat step 2 until some stopping criterion is met. The utility of a station is dependent not only on its ABS and/or CSB patterns but also on the way it schedules the users. In the following we present the PF scheduler used in later simulation studies to meet the objective in (2).

B. PF Scheduler

In this paper, we consider the following PF scheduler [START_REF] Sesia | LTE -The UMTS Long Term Evolution: From Theory to Practice[END_REF]: For station i, the b-th PRB at subframe t, given that the subframe is not an ABS, will be allocated to the following user

u b arg max u∈Ui r u,b r u (t) , where b ∈ [1, B], t = τ (b), r u (t)
gives the long-term average throughput of user u in subframe t and it is calculated as

r u (t) = (1 - 1 t c )r u (t -1) + 1 t c Σ {b|τ (b)=t} r u,b • I{ u b = u}.
In the above equation, t c is the time window which is a design parameter and I{•} is the indicator function.

IV. INFORMATION EXCHANGE REQUIREMENT

It is important to know the information that is required to play the potential game, as it is directly related to implementation complexity. For the selected player i, the required information is: 1) Users who attach to player i when player i adopts one of the CSB values in C. 2) Utilities of stations in N i \ {i}, i.e., {U j |j ∈ N i \ {i}}, when player i adopts different strategies in S i . For stations in N i \ {i}: 1) They need to know which users will attach to them if station i uses one of the CSB values in C. 2) They need to know the interference power from station i to the users attached to stations in N i \ {i}. In case station i is an macrocell, the interference power is dependent on the ABS patterns that station i may adopt.

V. NUMERICAL STUDIES Fig. 1 shows a randomly generated hexagonal HetNet consisting of seven clusters of macrocells, where the clusters' borders are marked with bolded lines. For the center cluster, the macrocells are located at the center of the hexagons and one picocell is randomly placed within each hexagon so that the distance between the picocell and its closest macrocell follows the triangular distribution. 14 users are randomly located within 40 meters of each picocell in the center cluster, and another 6 users are randomly located within each hexagon in the center cluster. The six surrounding clusters of the center cluster are copies of the center clusters. The sets M and P contain macrocells and the picocells, respectively, in the center cluster. The users are assumed to be static and the buffers of the stations are assumed to be fully loaded. All users are equally weighted, i.e., w u = 1 for all u. Each PRB experiences independent Rayleigh fading with variance 1. Shadowing is correlated, and the correlation is achieved by generating a common shadowing value and divide the sum of the common shadowing value and the individual shadowing value by √ 2. Other parameters can be found in Table I. Fig. 2 shows the possible ABS patterns of an macrocell. The characteristic of the set of ABS patterns in Fig. 2 is that subframes which are ABSs in an ABS pattern with lower ABS ratio are also ABSs in an ABS pattern with higher ABS ratio. This kind of design may allow the surrounding macrocells of a picocell to create subframes that are almost interference free even if the macrocells uses different levels of ABS ratios.

On the other hand, the potential game formulation allows a designer to use any kind of ABS and CSB patterns, as the strategy set of a player can be easily modified. Next we specify the set N i for i ∈ M ∪ P. We assume that the transmission of station i will create interference to the users who are located within the same hexagon of station i and the six surrounding hexagons. Also, we assume that a user can only be associated to stations in the same hexagon and its six surrounding hexagons. Therefore, N Int i and N Att i include stations located within the same hexagon of i and the twelve surrounding hexagons. For example, N 1.1 includes all stations in the center cluster in Fig. 1.

Fig. 3 plots the change of global utilities (obtained by ( 8)) using different scheduling schemes as the number of potential game plays increases. For the PF scheduler, three window sizes are evaluated, namely t c = 5, t c = 20, and t c = 40. The plot is obtained by averaging 186 randomly generated network topologies. The global utilities increase as the game is played more times because in each round a player uses a best response strategy, i.e., select a strategy that maximizes its payoff function. Also, the algorithm converges in around 150 iterations. We can also observe that the global utility does not improve much if only the CSB values are optimized, and much of the improvement achieved by ABS and CSB optimization can also be achieved by ABS optimization along. We believe this is reasonable because most of the users are located near the picocell in each hexagon so that they receive fairly good service from the picocell, and each picocell has plenty PRBs to serve its associated users.

Fig. 4 plots the empirical cumulative distribution functions (CDFs) of the rates of the users averaged over their number of allocated PRBs after 300 iterations of game play, where both ABS and CSB patterns are optimized. We can observe that the users' rates per PRB are greatly increased due to the ABS and CSB optimization. Specifically, the median rate after ABS and CSB optimization is increased by roughly 60% compared to the case where no optimization is performed. Fig. 5 and Fig. 6 plot the CDFs of the averaged rates per PRB when ABS optimization is performed and when CSB optimization is performed, respectively. We can see from Fig. 5 that the improvement due to ABS optimization is similar to the improvement due to ABS and CSB optimization. Also, we can observe from Fig. 6 that by performing CSB optimization along, little improvement can be obtained. These observations coincide with the observations from Fig. 3.

Table II summarizes some key performance indicators of various optimization scenarios. We can observe that the ABS and CSB optimization gives about 10% improvement on the mean achieved user rate compared to the case where no optimization is performed, where a user's achieved rate is calculated by adding the capacity of each PRB that is allocated to that user. Also, the Jain's fairness indices on users' throughput when ABS and CSB optimization is performed are greatly improved compared to other cases. Moreover, we can see that the PF scheduler with time window being 5 has the

VI. CONCLUSION

We propose a distributed eICIC optimization algorithm based on potential game formulation. We show that such formulation can be used to optimize different objectives, e.g., maximizing proportional fair and maximizing sum of users throughput. Simulation studies show that the optimization

  b) P M\{m} Int,u,b + P P Int,u,b + N 0 , where h m u,b is the fading gain on PRB b from macrocell m to user u, τ (b) returns the index of subframe of PRB b, P M\{m} Int,b is the interference from macrocells in M\{m} at PRB b, P P Int,b is the interference from picocells in P at PRB b, and N 0 is the Additive-White-Gaussian Noise (AWGN). Note that P M\{m} Int,u,bdepends on the ABS patterns of the macrocells in M \ {m} and the large and small scale fading gains of the channels from the macrocells in M \ {m} to user u. Also, note that P P Int,u,b depends on the large and small scale fading gains of the channels from the picocells in P to user u. The SINR of a user u associated with picocell p ∈ P at PRB b is given as: is the channel gain from picocell p to user u at PRB b, P M Int,u,b is the interference from macrocells in M, and P P\{p} Int,u

  a) ABS and CSB optimization: S i := A × C for i ∈ M. S i := C for i ∈ P. b) ABS optimization: S i := A for i ∈ M. S i := { -→ 1 } for i ∈ P, where -→ 1 means the vector of all 1s. c) CSB optimization: S i := C for i ∈ M ∪ P.

Fig. 1 .

 1 Fig. 1. A randomly generated hexagonal HetNet layout. The triangles represent macrocells and the dots represent picocells. The number 4.1 near one of the triangles indicates that it is an macrocell located in the first (also the center) cluster with reference number 4. The users are not plotted for clean presentation.

Fig. 2 .

 2 Fig.2. A set of ABS patterns, where the filled boxes indicate ABSs. Subframes that are set as ABSs in a pattern with a lower ABS ratio are also set as ABSs in a pattern with a higher ABS ratio.

Fig. 3 .

 3 Fig.3. Averaged global utility against the number of times the potential game is being played using various scheduling schemes.

Fig. 4 .Fig. 5 .

 45 Fig.[START_REF] Pang | Optimized time-domain resource partitioning for enhanced inter-cell interference coordination in heterogeneous networks[END_REF]. CDF plot of throughput of users who are located in the center cluster in all randomly generated topologies when joint optimization (i.e., ABS and CSB joint optimization) is performed.

Fig. 6 .

 6 Fig.6. CDF plot of throughput of users who are located in the center cluster in all randomly generated topologies when CSB optimization is performed.

TABLE I PARAMETERS

 I FOR GENERATING HETNET TOPOLOGIES

	Parameter	Value
	Inter-macrocell distance	1000 m
	Minimum distance from macro to user	35 m
	Minimum distance from pico to user	10 m
	Minimum distance from macro to pico	75 m
	Antenna per site	Omnidirectional × 1
	Macrocell power	40 W
	Picocell Power	1 W
	Noise density	-174 dBm/Hz
	Noise figure	9 dB
	Bandwidth per RB	180 kHz
	CSB values	0

dB, 5 dB, 10 dB, 15 dB Log-normal shadowing SD 10 dB Path loss from macrocell to user 128.1 + 37.6 log 10 d, d in km Path loss from picocell to user 140.7 + 36.7 log 10 d, d in km

TABLE II FAIRNESS

 II INDICES AND AVERAGE ACHIEVED USER RATE can lead to 10% increase in average achieved user rate while at the same time greatly increase the service fairness of the network.

	Optimization	Scheduler	Jain's Index	Mean UE Rate
				(bit/s/Hz)
		PF tc = 5	0.8263	186.89
	ABS+CSB	PF tc = 20 PF tc = 40	0.8276 0.8281	192.93 193.66
		RR	0.7904	154.09
		PF tc = 5	0.7799	183.25
	ABS	PF tc = 20 PF tc = 40	0.7798 0.7796	189.49 190.30
		RR	0.7152	149.34
		PF tc = 5	0.7213	167.33
	CSB	PF tc = 20 PF tc = 40	0.7246 0.7250	171.27 171.77
		RR	0.6517	124.93
		PF tc = 5	0.7002	167.86
	No optimization	PF tc = 20 PF tc = 40	0.7057 0.7063	171.79 172.27
		RR	0.6343	125.41

algorithm

ACKNOWLEDGEMENT

This work was supported by ANR project IDEFIX under grant number ANR-13-INFR-0006 and in part by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China, under Project CityU 121713. A part of the work was carried out at LINCS (www.lincs.fr).


Arnaud Delaplace

Ste ´phane Roux

Gilles Pijaudier-Cabot

AVALANCHE STATISTICS OF INTERFACE CRACK PROPAGATION IN FIBER BUNDLE MODEL: CHARACTERIZATION OF COHESIVE CRACK

This paper considers a model of crack propagation taking place at the interface between a rigid support and an elastic plate. The interface is modeled using a fiber bundle model (i.e., describing a damage behavior using a discrete set of elastic brittle elements having a random strength). This paper studies the fluctuations of the force required to propagate the crack along the interface. The statistics of avalanches, defined as a series of elements that are broken simultaneously under a load that decreases with the crack advance, are studied numerically and analytically. Local fiber breakage kinetics is related to a correlation length, which sets the size of a fracture process zone. FIG. 1. Response of Four-Point Bending Test on Fiber-Reinforced Concrete (50-cm-long Specimens with 8-cm Square Section; Points Represent Force Peak for Three Curves; Curves Show Variability of Response for Nominally Identical Specimens and Local Fluctuations of Force Displacement Curves

INTRODUCTION

Damage in quasi-brittle materials such as concrete is due to the accumulation of microscopic crack nucleation, propagation, and most importantly, arrest. Arrest of microcracks is an essential feature to produce progressive damage; otherwise, a brittle behavior would result. To allow for microcrack arrest, it is therefore important to deal with material heterogeneities. In spite of the importance of this factor, most damage approaches in continuum mechanics treat the solid as homogeneous and little attention is paid to the detailed material microstructure producing the progressive damage behavior. Only the mean effects of the microcrack nucleation, propagation, and arrest are considered. The effects that are neglected (i.e., the intrinsic heterogeneity of the solid) may, however, play a role, and this study illustrates some of these effects in terms of the fluctuations of the macroscopic response due to the existence of an underlying heterogeneity.

An experimental example illustrates this discussion. Fig. 1 shows the measured force-displacement relation for a fourpoint bending test performed with a steel-fiber reinforcedconcrete specimen. Different curves corresponding to different samples of nominally identical material and geometry are shown. One observes a rather large variability of the response. This noise is not an artifact of the measurement device but really reflects the material behavior. It renders difficult the identification of a constitutive behavior law from such a test. Moreover, considering one single curve, one observes a noisy response consisting in abrupt drops of the force at fixed strain. Elucidating the statistical feature of such a ''noise'' is one of the objectives of this work.

Another effect that may result from this heterogeneity is the size effect (i.e., size dependence of the macroscopic strength), a crucial property when scaling of structural response has to be performed. However, in the particular model considered hereafter, such a size effect does not exist.

This paper will mainly focus on the fluctuations of the macroscopic response, in a geometry such that a usual macroscopic modeling will produce a steady response. This particular choice will allow one to study the fluctuations directly and to emphasize the effective noise, which would be obtained experimentally, in the form of ''avalanches'' (i.e., a series of microcrack events that are produced quasi-simultaneously under constant conditions of loading to be specified more precisely below). The notion of avalanches is here borrowed from a large body of works in statistical physics where this notion was proven to be useful in the analysis of various critical systems [START_REF] Hemmer | The distribution of simultaneous fiber failures in fiber bundles[END_REF][START_REF] Paczuski | Avalanche dynamics in evolution, growth and depinning models[END_REF]Delaplace et al. 1999a). Those avalanches depend not only on the total statistical distribution of critical loading but, more importantly, on time (or displacement) correlations.

There are few models that can address such questions without having to rely heavily on numerical simulations. Among these, one of the most intensively studied is the fiber bundle model introduced by [START_REF] Daniels | The statistical theory of the strength of bundles of threads[END_REF]. It consists of a series of elastic brittle ''fibers'' loaded in parallel between two rigid bars. The advantage of this simple model is that almost all of its properties can be solved for analytically. Moreover, it has been extended to account for a wide variety of more complex phenomena, such as time-dependent failure. The fact that the fibers are loaded through two rigid parallel bars implies that all fibers are subjected to the same displacement; thus, the load is always equally shared by all the surviving fibers. This is called the global load sharing rule. The opposite limit of a very compliant interface has been introduced by [START_REF] Harlow | Approximations for the strength distribution and size effect in an idealized lattice model of material breakdown[END_REF] and gives rise to a very different behavior where the system behaves as macroscopically brittle. An intermediate case, taking into account a coupling of the fibers to a semi-infinite elastic medium, has more recently been studied by Delaplace et al. (1999a). This last case gives rise to a homogeneous damage up to a localization point where the most unstable case has a wavelength equal to the system size, provided the elastic block is not too compliant. A major conclusion from this study was that, at the inception of interface cracking, there is no such thing as a finite size fracture process zone. If a process zone of finite size has to be expected, it must develop after the inception of cracking, in the propagation regime that is of concern in this paper.

Studying the original Daniels' model, [START_REF] Hemmer | The distribution of simultaneous fiber failures in fiber bundles[END_REF] have worked out analytically the statistics of avalanches for this system. This is one of the few works that really extracted information about the fluctuation of the force-displacement characteristics for this system. The same results can immediately be transferred to the case of a more realistic load transfer mechanism, such as the case of an elastic block. The statistics of avalanches has also been extended to the local load sharing rule (limit of extremely compliant coupling) in [START_REF] Kloster | Burst avalanches in solvable models of fibrous materials[END_REF]. The aim of the present study is to extend the study of these avalanches to a system derived from the fiber bundle model, where fibers are connected to a deformable media and where the geometry is suited to an inhomogeneous, but steady, state of loading. The geometry mimics the propagation of a crack front at the interface between a rigid interface and a deformable plate (or rather, a beam in two dimensions). It can be looked at as a schematic model for mode I crack propagation in a quasi-brittle material (the rigid side of the interface representing an axis of symmetry). The purpose of introducing a nonlinear interface is to concentrate cracking onto a specific region, a line in the present approach, the same as in a cohesive crack model.

In the context of a steady-state crack propagation, this problem has been studied with different approaches, motivated by an experimental study of the roughness of crack fronts by [START_REF] Schmittbuhl | Direct observation of a selfaffine crack propagation[END_REF] and [START_REF] Delaplace | High resolution description of a crack front in an heterogeneous plexiglass block[END_REF]. Discrete fuse networks have been studied by [START_REF] Zapperi | Planar cracks in the fuse model[END_REF] to mimic this situation, considering the steady-state regime. The problem of nucleation of such interfacial cracks using a variety of load sharing rules has also been developed. However these studies were focused on either the global average behavior or the morphology of the crack front. The focus here will be different in the sense that one concentrates on the characterization of the fluctuations of the force-displacement curve.

This paper is organized as follows. First, the model is defined and a mean-field version of it is introduced in which the profile of imposed displacement along the interface is fixed. The statistical distribution of the loading is obtained simply in this limit and is shown to converge to a Gaussian as the size of the ''process zone,'' as compared to the microstructural size, diverges. Then, ''avalanches'' are introduced and their statistics studied numerically and analytically in the framework of the simplified model. Finally, the numerical study of the complete model is reported, in which the displacement at the interface results from beam deformation. The results are shown to be similar to those obtained for the simplified model.

ZIP MODEL

General Case of Elastic Coupling

First recall the definition of the Daniels' model (1945). A collection of elastic brittle fibers are loaded in parallel between two rigid bars. Each fiber has a perfectly brittle behavior with the same stiffness . However, the strength of the fibers are random uncorrelated variables. The loading is an imposed displacement, and the response is the global force-displacement curve.

In the model, the fibers are connected to a semi-infinite 1D beam on one side and a rigid substratum on the opposite side. One could also consider another elastic beam instead of the rigid part (even having a different stiffness) without any change in the formulation of the problem. A normal displacement is imposed at one point of the elastic beam, which may move along the interface, as if a wedge was pushed in a double cantilever geometry. Fig. 2 illustrates the model. In the notation, the position of fiber i is x i and its extensional displacement is denoted by y i . The spacing between fibers is set to 1 (i.e., x iϩ1 Ϫ x i = 1), thus defining a fixed microstructural size. The maximum extension a fiber can stand before breaking y c is chosen randomly, as mentioned above. This distribution is uniform in the interval [0:1] in the analytical part of this study, whereas other distributions could also be considered in the simulations. A fiber that is strained up to its critical extension is broken instantaneously. The normal displacement imposed on the beam is such that all fibers under this extension are broken with a probability of 1.

The interface consists basically in three parts. Under the wedge, all fibers are broken. This constitutes the crack by itself. Ahead there is a ''process'' zone where a finite fraction of fibers are broken. Finally, moving away from the crack tip, broken fibers become more scarce and their influence can be neglected. There is no clear-cut separation in the last two regions; however, the probability of fiber rupture ahead of the crack tip decays exponentially and thus can be used to define the extension of the process zone. The terminology for the fracture process zone is borrowed from [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF], who studied fracture in concrete with the help of a cohesive crack model. They assumed that the fracture process zone was a region of intense damage with the material, although aggregates bridge the crack and allow for stress transfer across the crack.

Imposed Displacement Profile

A simplified model where the displacement profile along the upper beam is imposed will be implemented in the analytical study. The interface opening has an exponential shape, which captures some of the features of a beam deflection (as will be seen later). For any abscissa x, the profile y is given by where the length scale is considered as a fixed parameter; and U = time-dependent horizontal displacement of the edge.

The major advantage of this simplifying assumption is that it allows derivations of closed-form expressions for the loading while it is not too far away from the exact beam deflection.

SIMPLIFIED ANALYTICAL MODEL Mean Behavior and Fluctuations

The total force exerted on the wedge is simply the sum of all fiber contributions

F = f (U ) (2) i i
where each individual force f i is

f (U ) = y(x , U ), if y(x , U ) < y (x) ( 3 a) i i i c f (U ) = 0, if y(x , U ) > y (x) ( 3 b) i 1 c
In the case of a uniform distribution in the interval [0:1] of the breaking extension of individual fibers, one may compute the average force in fiber i

U Ϫ i UϪ i ͗ f (U )͘ = 1 Ϫ exp exp (4) i ͫ ͩ ͪͬ ͩ ͪ
for the ith fiber at a location such that i Ն U (the fiber spacing is set to 1) and ͗ f i ͘ = 0 otherwise. Hence, the mean total force is

͗ F(U )͘ = exp(d(U )/) exp(Ϫi/) Ϫ exp(2d(U )/) i ␣ и exp(Ϫ2i/) = [exp(d(U )/)(1 ϩ ␣) 2 1 Ϫ ␣ i Ϫ exp(2d(U )/)␣] (5)
where d(U ) = U Ϫ int(U ) (i.e., the noninteger part of U); and ␣ ϵ exp(Ϫ1/). Fig. 3 shows the resulting force as a function of U. One observes a small amplitude component due to the discrete nature of the fiber distribution. However, the amplitude of the latter modulation is very small and typically it is negligible when compared to the fluctuation of the force, when the loading length scale is large enough. For integer valued displacement d(U ) = 0, and for large , the mean force takes its minimum value and can be approximated by the following expression [from ( 5

)] 1 ͗ F͘ Ϸ 1 Ϫ (6) ͩ ͪ 2 2 6
The maximum force is obtained for displacement such that

d͗ F(U )͘ = 0 dd(U) and leads to exp(1/) ϩ 1 1 1 d(U ) = log Ϸ ϩ (7) ͩ ͪ 2 2 8
Substituting the last expression in (5) results in the peak mean force

1 max͗ F(U )͘ Ϸ 1 ϩ (8) ͩ ͪ 2 2 1 2
and thus relative peak-to-valley difference vanishes as 1/(4 2 ) for large . Therefore in the following, one may disregard this systematic variation and only consider the case of integer displacement U.

The variation of the total force can also be obtained in a similar way, thanks to the fact that the threshold forces for different fibers are uncorrelated. The variance of the total force is simply the sum of the variances of the f i , and thus using

2 2 22 2 23 ͗ f ͘ Ϫ ͗ f ͘ = y [(1 Ϫ y) Ϫ (1 Ϫ y) ] = y [1 Ϫ y] (9) i i
one obtains for the total force variance 2 (F )

ϱ 2 2 2 2 (F ) ϵ ͗F ͘ Ϫ ͗F͘ = exp(3d(U )/) exp(Ϫ3i/) ͫ i=1 ϱ Ϫ exp(4d(U )/) exp(Ϫ4i/) ͬ i=1 (10) hence 3 4 ␣ ␣ 2 2 (F ) = exp(3d(U )/) Ϫ exp(4d(U )/) ͫ ͬ 3 4 1 Ϫ ␣ 1 Ϫ ␣ (11)
Considering only integer value

U 3 ␣ 2 2 (F ) = ͫ ͬ 2 3 4 5 (1 Ϫ ␣)(1 ϩ 2␣ ϩ 3␣ ϩ 3␣ ϩ 2␣ ϩ ␣ ) (12)
In the limit of a large , the above expression reduces to

2 1 2 (F ) Ϸ 1 Ϫ (13) ͩ ͪ 2 12
Therefore the relative standard deviation of the total force will approach (F )/͗F͘ → for >> 1. 1/ 3

͙

The entire distribution of the total force can be obtained analytically in the context of this simplified model. In the limit of a large length , one can simply observe that F is given by a sum of statistically independent random variables and the law of large numbers applies. Consequently, F will have a Gaussian distribution. Thus the above computed average and variance are sufficient to specify entirely the distribution. Keeping only the leading terms, one can write

2 6 3(2F Ϫ ) ͙ p(F ) = exp Ϫ (14) ͩ ͪ 2 2 ͙

Correlations

The above section gave some information concerning the total force distribution. To characterize the evolution of the signal F(U), one needs to work out the correlations in this function. More precisely, knowing the force F(U ), the question is, What can be inferred for F(U ϩ ⌬U )?

One can get some insight into this question by computing the expectation value of dF(⌬U ) = F(U ϩ ⌬U ) Ϫ F(U ). The force difference can again be written as a sum of independent statistical variables

dF = df (15) j j
where j extends over the unbroken fibers after the crack tip U.

The random variables df j assumes the following values: where the last condition has been written so that the sum can extend over all fibers for j > U ϩ ⌬U. One has to consider also the fibers in the range U < j < U ϩ ⌬U, which are broken with probability 1, but were surviving under the displacement U with probability 1 Ϫ y j (U ). At the steady-state value of the force, ͗dF͘ = 0 by definition. The expectation value of ͗dF 2 ͘ can thus be obtained by summing up the variances of the df j

2 2 ͗dF ͘ = (1 Ϫ exp(Ϫ⌬U/)) (F ) (17)
Thus the squared force difference increases first linearly with ⌬U and saturates to a constant value equal to twice the variance of the force. The interpretation of this property is straightforward: the total force is the sum of the order of uncorrelated random variables. Thus over this length scale, F(u) behaves as a random walk. However, for larger distances, the fluctuations of F become uncorrelated. As a consequence of this observation, note that the typical variation of the force over a short displacement ⌬U scales as Consequently, the fluctuating part of the signal be-⌬U.

͙

comes nondifferentiable when the microstructural size goes to 0, keeping fixed (note that this construction implies a redefinition of the physical scale because one chose here to measure distances in terms of the microscopic distance, the fiber separation distance). It is important to stress that such a conclusion has already been shown for the fiber bundle model in the case where a global load sharing rule was used and thus for a homogeneous macroscopic loading.

If the limit of a bundle with an infinite number of fibers is expected to represent a continuum response, the fact that the constitutive response is not differentiable is a striking departure from traditional assumptions. In fact, it means that when this continuum limit is considered, a smoothing of the constitutive response is performed at the same time so that the response becomes differentiable. It follows that upon taking this limit, the information contained in the fluctuations of the response for large-size systems is lost. As will be seen in the next section, this information enlightens the crack propagation regime and yields a parallel between the fiber bundle model and cohesive crack models in the context of this study.

Avalanches

There exists various definitions of avalanches. The first one is the most natural, consisting of fixing a level of force and computing the distance ⌬U over which the crack can propagate. The avalanches are characterized by their statistical distribution, p 1 (⌬U, F ).

To give a global characterization of the signal without considering a specific value of the force, one may consider the above avalanches for any crack length U such that a fiber would break and then average over all breaking events. These avalanches are the ''forward'' avalanches. A similar construction can be performed for avalanches ending at each fiber breaking event. They correspond to a similar construction as the previous ones after reversing the arrow of time. They are called ''backward'' avalanches. In some cases such as invasive percolation, these two statistics have very different properties. Fig. 4 shows the computed forward avalanche distribution for the present model where a fixed displacement profile is set. Three values of have been used. The first one is = 100 (⅜) with 10 8 broken fibers, the second one is = 1,000 (ϫ) with 2.10 7 broken fibers, and the third one is = 10,000 (᭟) with 5.10 6 broke fibers.

The main observation is that the distribution exhibits two distinct behaviors according to the value of the avalanche size ⌬ with respect to a crossover value ⌬*. The first regime, ⌬ < ⌬*, is a power law p 1 (⌬) ϰ with an exponent

Ϫ 1 ⌬ = 1.50 Ϯ 0.05 (18) 1
The second regime, ⌬ > ⌬*, is also a power law but with a different exponent 2 = 2.05 Ϯ 0.10 ( 19) 2 Fits to both of these power laws are plotted in Fig. 4. Finally, the crossover scale ⌬* scales as : the graph used the scaled variables ⌬/ and scaled distribution to show that the 1 p (⌬) 1 three curves collapse onto a single master curve. This data collapse shows that, indeed, the crossover scale ⌬* is proportional to .

One can easily understand the value of the two exponents. The first one corresponds to a regime where the force versus crack length U displays correlations similar to a random walk. The forward avalanche, in this case, can be interpreted as the time required for a random walk to return to the origin. This well-known statistical problem is indeed a power law of exponent 3/2 in agreement with the first regime. Note that this behavior is exactly the one that has been established for the global load sharing fiber bundle with rigid boundaries [START_REF] Hemmer | The distribution of simultaneous fiber failures in fiber bundles[END_REF][START_REF] Hansen | Criticality in fracture: The burst distribution[END_REF], using also a mapping onto a random walk problem.

For large avalanches, ⌬ > ⌬*, the forces are uncorrelated. One thus may resort to this simple case to work out the avalanche statistics: let (t) be a random uncorrelated noise, with a distribution p() and cumulative distribution P() = ͐ p(x) dx. Starting at a given value of = 0 , the probability that an avalanche is >⌬ is Q(⌬, 0 ) = P( 0 ) ⌬ because the different values are uncorrelated. The cumulative distribution P 1 of forward avalanches is obtained from the integration of the above Q distribution over all starting points of distribution p(), hence

1 1 ⌬ ⌬ P (⌬) = P( ) p( ) d = u du = (20) 1 0 0 0 ͵ ͵ ⌬ 0
where one has used ⌬ >> 1. The avalanche distribution p 1 is obtained from the derivative of the cumulative distribution and leads to the power law p 1 (⌬) = ⌬ Ϫ2 for all distributions p().

In this problem, the large avalanches correspond to this case and, indeed, one observes an exponent 2 = 2. One sees in this particular example that a simple statistical analysis performed on the force signal allows one to extract the correlation length ϰ ⌬* without knowing it beforehand. This correlation length defines the size of the fracture process zone.

There is another definition of avalanches that might be easier to handle in experimental studies. In the above construction, one needs to know the entire signal F(U ). This requires an infinitely stiff mechanical testing device. In contrast, if a finite stiffness K is used, the test would not be able to follow the characteristic of the system unconditionally. After a local maximum (U 0 , F 0 ), U would jump directly to a new position

U 1 such that F(U ) < F 0 Ϫ K(U Ϫ U 0 ), for U 0 < U < U 1 , and F(U 1 ) > F 0 Ϫ K(U 1 Ϫ U 0 )
. Such a construction defines biased avalanches characterized by the value of the stiffness K. K , K 2 . The avalanches are power-law distributed up to ⌬ c with an exponent 3 Ϸ 1.5. One has seen above that the fluctuating part of the force signal had a self-affine character with a typical variation over a scale ⌬U proportional to dF = A⌬U 1/2 . Thus, at small scales, the fluctuations dominate whereas at a large avalanche size, the bias takes over the fluctuations. Thus there is a characteristic scale for the largest avalanche size, such that K⌬ c = A(⌬ c ) 1/2 ; hence

c Ϫ2 ⌬ ϰ K (21)
For smaller sizes, the bias can be neglected; hence, the biased statistics are again similar to a first return time distribution for a random walk. Hence a power law of exponent 3/2 is ex-pected consistently with the numerical results. This regime is, however, valid only for the case where K Ϫ2 << . Otherwise, above , avalanches cannot be distributed in a similar way because the force signal becomes uncorrelated. Solving for the case of an uncorrelated signal with a Gaussian distribution, one finds that the avalanches are uniformly distributed up to a maximum scale depending on K as ⌬ c ϰ K Ϫ1/2 . Fig. 7 illustrates this scaling behavior for a small value of .

Finally, the distribution of backward avalanches has been checked. One finds the same distribution as for forward avalanches, with the same two power laws. Superposition of forward and backward avalanche distributions is shown in Fig. 8. This result may appear as natural; however, in other models [START_REF] Paczuski | Avalanche dynamics in evolution, growth and depinning models[END_REF]) (e.g., invasive percolation), a very different behavior is obtained for the two definitions. In this case though, the random walk analogy as well as the case of uncorrelated noise signal both support such a symmetrical behavior.

MODEL WITH DEFORMABLE BEAM Mean Deflection of Beam

One now comes back to the complete model where the beam is not assumed to have a constant displacement profile. Rather one wishes to compute here the average shape of the deformed beam. Introducing E and I, the young modulus and transverse geometrical inertia of the beam, respectively, one can write, in the spirit of a continuum modeling (i.e., for length scales much larger than the fiber separation), an equation for the mean deflection of the beam y(x) This equation holds for a uniform distribution of critical fiber extension between 0 and 1, and for y < 1, whereas d 4 y/dx 4 = 0 for larger y. The boundary conditions are y(ϱ) = dy(ϱ)/dx = 0, y(0) = 1 (imposed displacement), and d 2 y(0)/dx 2 = 0 (no torque being applied at the loading point). The reason one has to distinguish between y and Y is that the deflection is not a monotonous function of x. Because damage is irreversible, one has to compute the maximum damage having been met by the corresponding section of fibers. One does not know the analytical solution to this problem; however, one can see that the quadratic nonlinear term becomes important at a large distance from the origin. Thus the asymptotic shape will have the following expression: ͙ ͩ ͪ

The oscillatory component is the one that makes the deflection nonmonotonous and thus requires the distinction between y and Y. Omitting this oscillatory component and the nonlinearity, one retrieves the exponential decay that was imposed in the simplified version of the model. A numerical integration based on a relaxation method has been done, and the results are shown in Fig. 9.

Avalanches

The problem requires a much longer computation time compared to the simplified model because of the larger number of degrees of freedom. For numerical convenience, only the process zone is dealt with, thus neglecting the presence of a few broken fibers ahead of this region. However, this zone is taken into account in the computation by introducing a boundary condition at the end of the process zone that represents an infinite beam connected to the substratum through intact fibers. This involves two relations between the derivatives of order 0-3 of the deflection function y(x). The length of the domain considered numerically is set equal to .

Fig. 10 shows forward avalanche distributions for two stiffnesses ratios and hence two values of . In both cases the previous results are recovered (i.e., two power laws with exponents 1 = 1.50 Ϯ 0.05 and 2 = 2.1 Ϯ 0.15). One checks the reliability of the measure by plotting and evaluating the mean shape of the beam during crack propagation. The results are very close to the simplified model. Hence, the fluctuation analysis from the force-displacement response, an accessible experimental information, provides information on the existence of a well-defined length scale related to the fracture process zone in the interface.

CONCLUSIONS

The microstructure of heterogeneous materials plays a crucial role in the rupture process. This paper introduces a model of a 1D crack propagation that includes randomness. Fluctuations that are encountered on the force-displacement response are studied through an avalanche distribution analysis. Two cases are considered. First, a simplified model with a prescribed displacement opening profile is considered along the interface and it could be solved analytically. Second, a more complete model, where the coupling is done between a flexible beam and the fiber bundle, is considered. The latter was studied only through numerical simulations. Analytical predictions are supported by the numerical result of the simplified model. The main conclusions are • Although the response shows large fluctuations, an analysis of the avalanche distributions provides information on a redistribution length that scales the fracture process zone. • In the limit of an interface of an infinite size, the analytical model yields a response that is nondifferentiable. It cannot be compared with the smooth response expected in the context of a continuous interface. Therefore, the fluctuations contain important information on the interface fracture process, which is lost in a continuum approach. • Two regimes are found on the avalanche distributions.

The first one is a power law with an exponent 1.5, which can be mapped onto a random walk. It concerns avalanches whose sizes are smaller than the fracture process zone. The second one is also a power law with an exponent 2, which corresponds to an uncorrelated noise. It governs large avalanche sizes. • A ''biased avalanche'' has been introduced, which takes into account a finite stiffness of the loading device or of the system itself. Similar to the previous definition of avalanches, it exhibits two different regimes depending on the stiffness and on the size of avalanches. • The crossover avalanche size between the two regimes gives natural access to the extent of the fracture process zone. This information may prove to be a precious experimental tool to provide access to this a priori unknown scale.

1

  Asst. Prof., Laboratoire de Me ´canique et Technologie, Ecole Normale Supe ´rieure de Cachan, 61 Ave. du Pre ´sident Wilson, F-94235 Cachan, France. 2 Prof., SVI, unite ´mixte Centre National de la Recherche Scientifique-Saint Gobain, 39 quai Lucien, Lefranc, 93303 Aubervilliers cedex, France. 3 Prof., Laboratoire de Ge ´nie Civil de Nantes-Saint Nazaire, Ecole Centrale de Nantes, 1 rue de la Noe, F-44321 Nantes cedex 3, France.

FIG. 2 .

 2 FIG. 2. Schematic Representation of ''Zip'' Model (Point Where Displacement Is Imposed May Move Along x-Axis)

  FIG. 3. Evolution of Mean Total Force (Bold Curve) Scaled by and Shifted by Ϫ1/2 as Function of Displacement U ( Has Been Set to 10, to 1; Long Dashed Lines Are Boundary of ͗F ͘-Values That Are Obtained from Eqs. (6) and (8); Thin Curve Represents Evolution of Mean Force for = 20 and Shows How Fast Variation with Decreases)

FIG. 4 .

 4 FIG. 4. Avalanche Distribution for Zip Model with Constant Beam Shape (Active Area Size Are 100 (⅜), 1,000 (ϫ), and 10,000 (᭟); Dashed Line Is Guideline with Slope of Ϫ1.5; Continuous Lines Are Guideline with Slope of Ϫ2.05)

  df = ( y (U ϩ ⌬U) Ϫ y (U)), with probability (1 Ϫ y (U ϩ ⌬U)) df = Ϫy (U ), with probability (y (U ϩ ⌬U ) Ϫ y (U )) (16b)

FIG. 7 .

 7 FIG. 7. Statistical Distribution

  FIG. 5. Force Signal (Thin Line), When Followed with Loading Device Having Finite Stiffness K, Gives Access to Bold Curve Envelope, Where Decreasing Parts Have Slope ϪK

  Fig. 5 illustrates these avalanches. Fig. 6 shows the statistical distribution p b (⌬, K) obtained for different values of the stiffness K. The distribution has been scaled by whereas the avalanche size has been scaled by 2 1

FIG. 8 .

 8 FIG. 8. Comparison of Forward (⅜) and Backward (ϫ) Avalanches (Computations Have Been Done with = 100)


João Peterson

Vinicius Lopes

Americo Cunha

João Victor Ligier

Lopes Peterson

Vinicius Gonçalves Lopes

Maximization of the electrical power generated by a piezo-magneto-elastic energy harvesting device

Keywords: nonlinear dynamics, energy harvesting, chaos, optimization, test 0-1 for chaos

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Energy harvesting devices are mechanisms able to collect energy from external sources, and store it as electric potential, for a possible use as electric energy [START_REF]Energy Harvesting Technologies[END_REF]. There is a wide range of possible applications for such devices, that goes from electrical sources for equipment displays to indicators lights mounted on embarked dispositives, such as cars'. They can also be useful to support small equipment placed far from electrical distribution system or in dangerous (restricted access) areas, where the best option may be an independent supply of energy [START_REF] Spies | Handbook of Energy Harvesting Power Supplies and Applications[END_REF]. Furthermore, in general these equipment have an extremely rich non-linear behavior [START_REF]Energy Harvesting Technologies[END_REF].

In order to fully understand the potential of energy harvesting, and to improve its efficiency, it is necessary to understand how their parameters' values on its nonlinear dynamic behavior. In this sense, this work aims to study the nonlinear dynamics of the energy harvesting device proposed by [START_REF] Erturk | A piezomagnetoelastic structure for broadband vibration energy harvesting[END_REF], which is composed by the classical Moon and Holmes magneto-elastic beam [START_REF] Moon | A magnetoelastic strange attractor[END_REF], coupled with laminae of piezoelectric material. The paper also aims to solve an optimization problem, in order to maximize the mean electrical power in the device. Therein, in section 2 presents the modeling of the harvesting device. Numerical experiments with this model are discussed in section 3. Finally, the section 4 summarizes the main conclusions of this work.

2 Modeling of the nonlinear dynamic system

Physical model

For the nonlinear dynamics study proposed here, the mathematical and theoretical development are focused on the harvesting device proposed by [START_REF] Erturk | A piezomagnetoelastic structure for broadband vibration energy harvesting[END_REF], which is sketched in Figure 1. The structure supports a slim beam made of a ferromagnetic material, and subjected to two magnetic fields, provided by magnets placed in the lower part of the structure. The excitation of an external force creates vibration and, once this movement is perceived by the piezoelectric laminae, placed on highest part of blade, the mechanical energy is converted into voltage.

Mathematical model

As shown in [START_REF] Erturk | A piezomagnetoelastic structure for broadband vibration energy harvesting[END_REF], the dynamics of the harvesting device of Figure 1 is described by the following set of ordinary differential equations

ẍ + 2ξ ẋ - 1 2 x(1 -x 2 ) -χυ = f cos Ωt, (1) υ
+ λυ + κ ẋ = 0, ( 2
)
where x is the beam's extreme displacement, ξ is the mechanical damping ratio, χ represents a piezoelectric coupling term in the mechanical equation, f is the amplitude of excitation, Ω is the forced excitation frequency, υ means the voltage by the dispositive, κ is a piezoelectric coupling term in the electrical circuit equation and, finally, λ is a reciprocal time constant. These differential equations are supplemented by the initial conditions x 0 , ẋ0 , and v 0 . All of these parameters are dimensionless, assuming the following nominal values Ω = 0.8, ξ = 0.01, χ = 0.05, κ = 0.5, λ = 0.05, x 0 = 1, ẋ0 = 0, and v 0 = 0. The value of f will be specified soon.

3 Numerical experiments

Exploration of nonlinear dynamics

The voltage signal may have regular (non chaotic) or chaotic behavior [START_REF] Strogatz | Nonlinear Dynamics and Chaos: with applications to physics, biology, chemistry and engineering[END_REF], depending on the values of certain parameters of the harvesting device underlying dynamic system, such as f . This fact is illustrated in Figure 2, through the voltage time series, and in Figure 3 through Poincaré sections in x ẋ plane and the corresponding phase space trajectory projection. Similar results are observed if x 0 is changed. Therefore, depending on the chosen value for the pair x 0 and f , it is possible to conclude that the dynamics presents different characteristics. This fact was observed for the first time by [START_REF] Erturk | A piezomagnetoelastic structure for broadband vibration energy harvesting[END_REF], which did a similar analysis using the pair ( ẋ0 , f ).

Optimization of mean power in the harvesting device

Aiming to find a pair (x 0 , f ), in a given operating window x 0 min ≤ x 0 ≤ x 0max and f min ≤ f ≤ f max , that maximizes the mean power in the harvesting device, an optimization problem is solved in this work. The objective function is the mean power

P = 1 τ t+τ t P (t ) dt , (3)
being P(t) = υ 2 (t) the instantaneous power in the device, and τ a representative interval of (dimensionless) time in the steady state regime. For purposes of energy utilization, only the voltage presenting regular (not chaotic) behavior is of interest. This condition is introduced in the optimization problem as a constraint. In what follows, the strategy to solve the optimization problem is described. First, the operating window is discretized using a rectangular grid of points, where each point corresponds to a pair (x 0 , f ). Then, for each grid point, the system dynamics is integrated. Now the problem turns into separate interesting (non chaotic) time series of voltages, to compute the maximum mean power. At this point, the test 0-1 for chaos [START_REF] Gottwald | The 0-1 Test for Chaos: A review[END_REF] is applied. The pairs (x 0 , f ) associated with the non chaotic dynamics define the admissible region for the optimization problem. The mean power is computed in each point of the admissible region, and the maximum value is determined. The contour map of P , defined over -3 ≤ x 0 ≤ 3 and 0.01 ≤ f ≤ 1, is shown in Figure 4. For f 0.7 the initial displacement interference is negligible. As the force amplitude increases, the contour map reveals variations in P value across the x 0 direction. In fact, it is possible to note an asymmetry with respect to x 0 = 0. The highest (dimensionless) mean power is approximately 0.45, obtained at the pair (x 0 , f ) = (-3, 0.89), which is indicated in Figure 4 with a blue cross, in the border of the domain.

Final remarks

This work analyzed the non-linear dynamics of a piezoelectric harvesting device, and addressed an optimization problem seeking to maximize the electrical power in this device. Results presented by time series of voltage demonstrated that system dynamics may be chaotic or regular, depending on its initial displacement and force amplitude. Such behavior was also observed through an analysis using Poincaré maps. In the optimization problem solution, test 0-1 for chaos was used to identify the admissible region of parameters, and an optimal configuration (x 0 , f ) was obtained numerically. The main contribution of this work is the solution of the optimization problem reported above, once, to the best of authors' knowledge, there is no work in the open literature doing this type of analysis. In a future work, the authors' want to analyze how other parameters' influence over the non-linear dynamic behavior and see which pair of parameters are the most effective to maximize the mean power in the device.

Figure 1 :

 1 Figure 1: Schematic representation of the energy harvesting device proposed by [1].

Figure 2 :

 2 Figure 2: Time series of the voltage in the energy harvesting device, showing regular (left) and chaotic (right) behaviors.

Figure 3 :

 3 Figure 3: Poincaré sections in x ẋ plane (blue dots) and the corresponding projection of the phase space trajectory (green line), showing regular (left) and chaotic (right) behaviors.

Figure 4 :

 4 Figure 4: Contour map of the mean (dimensionless) power in the harvesting device. The maximum value is approximately 0.45, and occurs at (x 0 , f ) = (-3, 0.89) (blue cross).

joao.peterson@hotmail.com

vinigolop@hotmail.com

americo@ime.uerj.br

Acknowledgments

The authors are indebted to the Brazilian agencies CNPq, CAPES, and FAPERJ for the financial support given to this research. They also thank Dr. Jorge Gripp (ITA) for the enlightening comments.


Eber Dantas

Michel Tosin

João Peterson

Vinicius Lopes

Americo Cunha

A mathematical analysis about Zika virus outbreak in Rio de Janeiro

Keywords: nonlinear dynamics, mathematical biology, SEIR epidemic model, model calibration, Zika outbreak in Rio de Janeiro

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Zika virus was first isolated in primates from Zika forest at Uganda about 1947 [START_REF] Dick | Isolations and Serological Specificity[END_REF]. Evidences of the virus in humans were reported in Nigeria at 1968 [START_REF] Moore | Arthropod-borne viral infections of man in Nigeria, 19641970[END_REF]. It's first observation in south America occurred in Easter island, about 2014 [START_REF] Daniele | Zika Virus: Implications for Public Health[END_REF]. First cases in Rio de Janeiro, Brazil, were reported in 2015 [START_REF] Platonow | Estado do Rio Registra Primeiro Caso de Zika Virus[END_REF], quickly evolving to an outbreak in 2016. This scenario caused concern in medical community, health authorities and local population, specially due to a probable relation among Zika virus and other diseases, such as microcephaly [START_REF] Schuler-Faccini | Possible Association Between Zika Virus Infection and Microcephaly Brazil, 2015[END_REF][START_REF]Rapid risk assessment: microcephaly in Brazil potentially linked to the Zika virus epidemic[END_REF] and Guillain-Barré syndrome [START_REF] Cao-Lormeau | Guillain-Barré Syndrome outbreak associated with Zika virus infection in French Polynesia: a Case-Control Study[END_REF], and the realization of the 2016 Summer Olympics at Rio de Janeiro, which, in reason of the large flux of people, could spread Zika virus to several parts of the world.

In this (near) epidemic scenario, developing strategies of control and prevention of the disease is a critical issue. In this sense, a mathematical model able to predict the number of infected people, during the virus outbreak, is an useful tool, that can be used to identify vulnerable aspects of the disease control strategies.

This paper deals with the study of Zika virus outbreak at Rio de Janeiro in 2016. For this purpose, a mathematical model that was used to describe the outbreak in French Polynesia [START_REF] Kucharski | Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 201314 French Polynesia Outbreak[END_REF] is adapted to Rio de Janeiro scenario.

The nominal value of model parameters are estimated based on real data about the outbreak. Numerical simulation results are compared to the outbreak official numbers, in order to calibrate and validate the adapted model.

The rest of this paper is organized as follows. In section 2, the mathematical model is described, as well as its parameters. In section 3, the process of estimation of model parameters is discussed and numerical simulation results are reported, being followed by the comparison of numerical predictions and experimental data in section 4. Finally, in section 5, the main contributions of this work are emphasized and paths for future works are suggested.

Mathematical model

The transmission of Zika virus occurs, mainly, by Aedes genus of mosquitoes, but recently researches indicates that Culex genus is able to act as vector too [START_REF]Fiocruz Identifica Culex no Recife com Potencial para Transmitir o Virus Zika[END_REF]. Medical literature also relates sexual contact and blood transfusions as potential ways of transmission [START_REF] Musso | Potential sexual transmission of Zika virus[END_REF]. But, for purposes of modelling, in this work only the transmission by Aedes mosquitoes is considered.

To describe the Zika virus outbreak in Rio de Janeiro, this work employs an epidemic model of SEIR (susceptible, exposed, infected and recovered) type [START_REF] Keeling | Modeling Infectious Diseases in Humans and Animals[END_REF], adapted from the one used in French Polynesia [START_REF] Kucharski | Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 201314 French Polynesia Outbreak[END_REF]. In this model, the total population is distributed in four groups (susceptible, exposed, infected and recovered) according to the health condition of the individual at time t. The number of susceptible individuals at time t is denoted by S(t), and comprehends healthy individuals which does not have any immunological resistance against Zika virus. The label of exposed individuals, denoted by E(t), applies to the cases which are in an incubation period (latent period), i.e. a time gap while, despite of pathogen exposure, individuals are not infectious (asymptomatic period). The number of infected individuals, dubbed I(t), means affected but not latent anymore. This group is, however, capable of transmitting the disease. The recovered or immune group, represented by R(t), contains the individuals that were exposed and are recovered.

S v S h E h I h R h I v E v γ δ δ δ (βvI h ) 1 N α v β h Iv
The evolution of the diseased vector population obeys almost the same rules, in exception that there is no recovered group and it's lifetime is considered, since it is small compared to human lifetime. The individuals stream trough the model groups is depicted in Figure 2, both for affected and vector populations.

The stream (evolution) of individuals trough the groups follows the rules defined by the following set of differential equations

dS h dt = -β h S h I v , (1)
dE h dt = β h S h I v -α h E h , (2)
dI h dt = α h E h -γ I h , (3)
dR h dt = γ I h , (4)
dC dt = α h E h , (5)
dS v dt = δ -β v S v I h N -δ S v , (6)
dE v dt = β v S v I h N -(δ + α v ) E v , (7)
dI v dt = α v E v -δ I v , (8)
where β represents the transmission rates and γ, the recovery rate, defined as the inverse of infection period; α means disease's incubation ratio, calculated as the inverse of incubation period and δ, the vector lifespan, and N represents human total population. The h subscripts refers to the human parameters, while v, the mosquitoes (vector).

3 Numerical experiments

Model parameters

The following values are used, at first, for the model parameters: N = 16.5 × 10 6 humans, a recent estimation for Rio de Janeiro population [START_REF] Ibge | Estimativas da População Residente no Brasil e Unidades da Federação com Data de Referência em 1 o de Julho de 2015[END_REF]; extrinsic incubation ratio α v = 1/14 days -1 [START_REF] Edward | [END_REF]; intrinsic incubation ratio α h = 1/7.5 days -1 [START_REF] Lessler | Times to Key Events in The Course of Zika Infection and Their Implications: a Systematic Review and Pooled Analysis[END_REF]; human infectious ratio γ h = 1/7 days -1 [START_REF] Gregory | Zika Virus-Associated Microcephaly and Eye Lesions in the Newborn[END_REF]; mosquito lifespan ratio δ v = 1/25 days -1 [START_REF] Zettel | Linnaeus)(Insecta: Diptera: Culicidae)[END_REF]; while the vector-to-human transmission rate β h = 400 days -1 , and human-to-vector transmission rate β v = 40 days -1 were arbitrarily chosen through model experimentation, starting from the Tahiti values proposed for the French Polynesia outbreak [START_REF] Kucharski | Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 201314 French Polynesia Outbreak[END_REF].

Regarding the initial conditions, a fully susceptible population is considered, meaning S h (0) = N for humans and S v (0) = 1 for the proportion of mosquitoes. The hypothesis made for the French-Polynesian model, that initial number of infected (humans or vector) was equal to their respective incubating group is adopted: [START_REF] Kucharski | Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 201314 French Polynesia Outbreak[END_REF]. As such, the value of I v (0) = 0.014 from the Tahiti outbreak is also used [START_REF] Kucharski | Transmission Dynamics of Zika Virus in Island Populations: A Modelling Analysis of the 201314 French Polynesia Outbreak[END_REF]. The initial number of infected humans considered is I h (0) = 412, such number is a sum of sparse reported cases of Zika infection through some Rio de Janeiro cities in 2015 [START_REF] Faria | Zika Virus in the Americas: Early Epidemiological and Genetic Gindings[END_REF]. Of course, their sum is not a accurately representation of the quantity I h at such time (beginning of 2016), but a arbitrary choice to condensate the 2015 information due to the lack of more precise data.

I h (0) = E h (0) and I v (0) = E v (0)
Lastly, R(0) is arbitrarily chosen to be zero, signifying no recovered individual at the beginning of the outbreak, and C(0) = I h (0) because of the system dynamics implications.

Model predictions

In order to verify the representativeness of the proposed model, it is necessary to compare its predictions with field data. In this sense, real data about Zika outbreak, extracted from the epidemiological reports of Brazil's Ministry of Health [START_REF] Ministry | Monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika[END_REF] is used as reference. The report catalog data about the zika infection in the form of suspected accumulated cases since the thirteenth epidemiological week (03/27/2016 until 04/02/2016). The precise number of confirmed cases is actually unknown by the health authorities. Thus, this work took this value as 35%, as a first estimate.

Unfortunately, only six points of significant data are available through these reports so far (depicted on Table 1), and Zika cases along with the associated neurological and congenital syndrome related were not given notifiable conditions by the Brazil Ministry [START_REF] Ministry | Monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika[END_REF].

The comparison between the number of new infected humans per week, predicted by the model, and the real data is shown in Figure 2. Clearly, the reader can see that model response (blue curve) is not in good agreement with the real data (red circles).

Being 2015 the year of the official first confirmed case of Zika in Brazil and the probable epidemic peek that happened in mid-July of the same year [START_REF] Faria | Zika Virus in the Americas: Early Epidemiological and Genetic Gindings[END_REF], the lack of information about the spread of the disease during that year proposed a major difficulty on judging the consistency of the model. The choosing of parameters by ways of physical and biological context was not enough to establish the system on a basis of reasonable comparison with empirical data, and the lack of such data also hindered proper selection of Figure 2: Comparison between the model prediction for the number of new infected humans per week (blue curve), and real data (red circles) obtained from [START_REF] Ministry | Monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika[END_REF]. parameters as β h , β v and initial conditions as IV (0). Therefore, a calibration process is required to analyze the necessary parameters to fit the system.

Model calibration

To calibrate the model, a non-linear curve fitting procedure, based on a trust-region-reflective algorithm [START_REF] Coleman | An Interior, Trust Region Approach for Nonlinear Minimization Subject to Bounds[END_REF], is employed. The algorithm is performed first varying a single parameter, while the others remained fixed. This procedure allows to identify the parameters that most affect the model response. The parameters α h and α v show up to be unreliable, since they required very high initial guesses and usually would deform the quality of the fitting curve. Besides that, the model is very sensible to variations on β h and β v , being more affected by the values of β h . The other parameters do not bring remarkable conclusions.

Next, a two varying parameters attempt is conducted through the method. Most combinations do not bring satisfactory results, to the extent that some pairs of parameters could not even be computed in the adopted tolerance levels, e.g (α h , α v ), probably because of inefficient initial guesses or the low number of data points available. Furthermore, it sounds that some parameters control the quality of the fitting curve. Such control parameters are β h and β v and, unsurprisingly, the pair (β h , β v ) proved to be the best one for fitting purposes. Hence, the pair (β h , β v ) is chosen to be focus of the process of model calibration.

The parameter N is held fixed in all attempts. A change in the initial conditions is also performed, considering in particular I h and I v . The initial number of infected humans proved unresponsive, threatening significant changes on the curve only when varying his order of magnitude, what is beyond of our assumptions. But I v is noteworthy. The value of I v = 0.008 is chosen when testing the (β h , β v ) pair due to compatible general results. Some considerations on this process of investigation to find a good fitting curve regard the model curve peak position (local maxima). Brazil's 2016 Zika peak infection happened on February [19], but precise data for the state of Rio de Janeiro can only be estimated thus far, as discussed. This implicates that the six data points of Table 1 give very few information on the general shape the curve should have around the peek, because of their distance from such event. Additionally, the lack of reliable data about the infection in the beginning of year 2016 and end of 2015 allows a multitude of different shape-wise curves to fit the midyear data provided by the Ministry of Health. These problems manifested themselves while trying the parameter testing method, as a great number of low peak or wrong peak-positioned curves would fit the precise empirical numbers for the region (Table 1) but, of course, contradict the other information about the general distribution of the infection along the year.

The adopted solution for this lack of information is to perform (educated) guesses for the the peak region values, around the month of February. An imprecise number of 3500 new cases per week on February for Rio de Janeiro that circulated around brazilian media [START_REF] Prada | Redução de Casos de Zika Deve Diminuir Medo de Quem Vem para Olimpíadas[END_REF] is used on the iterative algorithm in search of a β h and β v that would better represent a peak like curve around such month. Figure 3.a summarizes the result of only using the data from Table 1, where the results of the estimation algorithm are β h = 0.0017 days -1 and β v = 4.0599 days -1 . Figure 3.b represents the curve when a single unofficial data point (marked by a asterisk) of 3500 new cases per week is assumed on the epidemiological week 7. For this case the fitting algorithm returned in β h = 0.0055 days -1 and β v = 0.6479 days -1 . Finally, Figure 3.c depicts the best fitted curved with respect to the empirical data at disposal, considering three extra unofficial data points: in addition to the 3500 new cases per week on the seventh week, a 2000 new cases/week on the second and 1300 weekly new cases on the third are included. The resulting estimated parameters are β h = 0.0045 days -1 and β v = 0.8296 days -1 .

Figure 3.a indicates a general shape of the curve by only considering the data of Table 1: a early peak on the year 2016 and a decay like appearance around the data points. The introduction of an arbitrary point in Figure 3.b moved the peak closer to the desired week and defined better the shape. However, without comparison data around the maximum region, the algorithm could not direct the parameters to values that would make 3500 be about the highest new infected per week around the seventh week. Using three arbitrary points as disposed on Figure 3.c, intentionally chosen to do the shape of the peak, the curve better distributed itself through the data and the magnitude stayed closer to the expected maximum value. Figure 3: Camparision between the model prediction (blue curve) and data from epidemiological bulletin (red circles). In the top only bulletin data is presented, while in the bottom 'ghost data' (black stars) are also used.

Although graphical results above are satisfactory, the estimated parameters are not good representations of reality. A β h = 0.0039 translates to a vector only contacting a human being in a period of 256 days, which surpass extensively the lifespan of the mosquito. This may indicate this SEIR model is unrealistic to represent Rio de Janeiro outbreak, once it fits the observation data only when unreal parameters are used.

Final remarks

This work adapted a SEIR epidemic model, used to predict French Polynesia Zika spread, to describe the 2016 Zika virus outbreak in Rio de Janeiro, Brazil. Nominal quantities for the parameters were select through available experimental data in conjunction with arbitrary adopted values. A model calibration procedure, that uses a nonlinear curve fitting method, was employed to pick the best parameter values that would fit the model response into the empirical data regarding the evolution of the infection. Due to the low number of real data about the Zika outbreak, the curve fitting procedure was not effective at first. The effects of the lack of information have been solved by incorporating the reference data ansatz for the outbreak peak values. In this second step, the calibration procedure could get a good fit to describe the outbreak. Nevertheless, the adjust parameters obtained this way are unrealistic, which may suggest that this SEIR model is not a realistic representation of the Rio de Janeiro outbreak scenario. In future work the authors intend to construct a probabilistic model of uncertainties, to better describe the model parameters variability, and employ a Bayesian procedure for model calibration, in order to to verify if the model calibration can only be done via unrealistic data.

Figure 1 :

 1 Figure 1: Schematic representation of SEIR epidemic model used to describe Zika virus outbreak.

Table 1 :

 1 of Health until February 17th of 2016 [19], thus leaving 2015 data unreliable. Number of suspected cases of infected humans in 2016

	week number of cases
	13	25930
	16	32312
	18	38196
	20	43516
	21	46027
	32	60176

eber.paiva@uerj.br

michel.tosin@uerj.br

joao.peterson@uerj.br

vinicius.lopes@uerj.br

americo@ime.uerj.br

Acknowledgments

The authors are indebted to the Brazilian agencies CNPq, CAPES, and FAPERJ for the financial support given to this research.


Gilbert Le Goff

Cyrille Lebon

Frédéric Jean

Sébastien Boyer

Louis-Clément Gouagna

Technique de l'insecte stérile à la Réunion : comprendre le fonctionnement de la cible : Aedes albopictus [poster]

teaching and research institutions in France or abroad, or from public or private research centers.

TECHNIQUE DE L'INSECTE STERILE A LA REUNION

Comprendre le comportement de la cible: Aedes albopictus Gilbert Le Goff, Cyrille Lebon, Frédéric Jean, Sébastien Boyer, Louis Clément Gouagna.

Contexte

• Le principe de la Technique de l'insecte stérile (TIS) est de lâcher des moustiques mâles élevés en masse et rendus stériles par irradiation.

• L'accouplement des moustiques mâles stériles avec les moustiques femelles ne donne pas de descendant.

• La TIS présente de nombreux avantages dont notamment le fait qu'elle est écologique et ne cible que l'espèce qu'on veut combattre (pas d'effets collatéraux sur les autres espèces ou organismes non cibles); elle permet de s'affranchir à terme de l'utilisation d'insecticides et peut conduire à la suppression de certaines populations de moustiques.

Etude du comportement sexuel des mâles stériles et non stériles

Objectif: Evaluer l'effet de l'irradiation sur 1) la capacité d'accouplement (Fig. 1), la longévité (Fig. 2) le cycle d'activité (Fig. 3.) la compétition pour l'accouplement (Fig. 4)

Etude de la dynamique des populations d'Aedes albopictus à la Réunion

Objectifs scientifiques

• Chez les moustiques, seules les femelles piquent. Dans une intervention par la TIS, un grand nombre de mâles (qui ne piquent pas) doit être relâché.

• La réalisation de la TIS nécessite au préalable de connaître la biologie et surtout le comportement sexuel des moustiques dans la natures, mais également l'effet de la stérilisation sur le comportement des mâles stérilisés. -moustiques mâles stériles sexuellement aussi compétitifs que les mâles non stériles;

Partenaires

-la TIS est plus efficace quand le nombre de mâles stériles est 5 à 10 fois supérieur au nombre mâles non stériles

Conclusion générale

Les preuves scientifiques de l'efficacité de la TIS sont réunies, mais des efforts sont nécessaires pour démontrer son efficacité sur le terrain à l'île de la Réunion.

•Les informations recueillies sur la biologie du moustique tigre permettent de cibler le moment de l'intervention et d'évaluer l'efficacité de la TIS

Fig. 1 :

 1 Fig. 1: Capacité sexuelle des mâles stériles et non stériles

Fig. 4 :Fig. 3 :

 43 Fig. 4: Compétitivité sexuelle entre mâles stériles et non stériles

  Suivi vidéo au laboratoire: deux pics d'activité même rythme d'activité chez les mâles irradiés et non irradiés La stérilisation des moustiques se fait par irradiation aux rayons gamma (Caesium 137 ) à l'aide de piège à odeur La dynamique des moustiques tigre, vecteurs potentiels d'arbovirus, montre une très nette saisonnalité, avec un pic d'abondance en décembre-janvier Objectif: surveillance du niveau des populations de moustique tigre à La Réunion ♀ ♂