Luis Pedro Montejano

Jonathan Chappelon

Jorge Luis

Ramírez Alfonsín

Ramsey for complete graphs

Keywords:

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

Let K n be a complete graph and let r ≥ 2 be an integer. A r-edge coloring of a graph is a surjection from E(G) to {0, . . . , r -1} (and thus each color class is not empty). Let k ≥ t ≥ 1 be positive integers. We denote by K [k,t] the complete graph on k vertices from which a set of edges, induced by a clique of order t, has been dropped, see Figure 1 Let k 1 , . . . , k r and t 1 , . . . , t r be positive integers with

k i ≥ t i for all i ∈ {1, . . . , r}. Let R([k 1 , t 1 ], . . . , [k r , t r ]
) be the smallest integer n such that for any r-edge coloring of K n there always occurs a monochromatic K [k i ,t i ] for some i. In the case when k i = t i for some i, we set

R([k 1 , t 1 ], . . . , [k i-1 , t i-1 ], [t i , t i ], [k i+1 , t i+1 ], . . . , [k r , t r ]) ≤ t i .
We note that equality is reached at min 1≤i≤r {t i |t i = k i }. Since the set of all the edges of K [t i ,t i ] (which is empty) can always be colored with color i. We also notice that the case R([k 1 , 1], . . . , [k r , 1]) is exactly the classical Ramsey number r(k 1 , . . . , k r ) (the smallest integer n such that for any r-edge coloring of K n there always occurs a monochromatic K k i for some i). We refer the reader to the excellent survey [START_REF] Radziszowski | Small Ramsey numbers[END_REF] on Ramsey numbers for small values. In this paper, we investigate R([k 1 , t 1 ], . . . , [k r , t r ]).

General upper bound

In this section we present a recursive formula (Lemma 2.1) that yields to an explicit general upper bound (Theorem 2.2). The latter contains the wellknown explicit general upper bound for R([k 1 , 1], . . . , [k r , 1]) due to Graham and Rödl [START_REF] Graham | Numbers in Ramsey theory[END_REF] (see Equation ( 4)).

The following recursive inequality is classical in Ramsey theory

r(k 1 , k 2 , . . . , k r ) ≤ r(k 1 -1, k 2 , . . . , k r ) + r(k 1 , k 2 -1, . . . , k r ) + • • • + (1) +r(k 1 , k 2 , . . . , k r -1) -(r -2)
In the same spirit, we have the following.

Lemma 2.1 Let r ≥ 2 and let k 1 , . . . , k r and t 1 , . . . , t r be positive integers with

k i ≥ t i + 1 ≥ 2 for all i. Then, R([k 1 , t 1 ], . . . , [k r , t r ]) ≤ R([k 1 -1, t 1 ], [k 2 , t 2 ], . . . , [k r , t r ]) +R([k 1 , t 1 ], [k 2 -1, t 2 ], . . . , [k r , t r ]) . . . +R([k 1 , t 1 ], [k 2 , t 2 ], . . . , [k r -1, t r ]) -(r -2).
A similar recursive inequality has been treated in [START_REF] Shi | A bound for multicolor Ramsey numbers[END_REF] in a much more general setting in which a family of graphs are intrinsically constructed via two operations disjoin unions and joins (see also [START_REF] Huang | New upper bounds for Ramsey numbers[END_REF] for the case r = 2). However, it is not clear how the latter could be used to obtain Lemma 2.1 that allows us to give the following general upper bound for R([k 1 , t 1 ], . . . , [k r , t r ]) (which was not considered in [START_REF] Shi | A bound for multicolor Ramsey numbers[END_REF]). Theorem 2.2 Let r ≥ 2 be a positive integer and let k 1 , . . . , k r and t 1 , . . . , t r be positive integers such that k i ≥ t i for all i ∈ {1, . . . , r}. Then,

R ([k 1 , t 1 ], . . . , [k r , t r ]) ≤ max 1≤i≤r {t i } k 1 + • • • + k r -(t 1 + • • • + t r ) k 1 -t 1 , k 2 -t 2 , . . . . . . , k r -t r
where

n 1 +n 2 +•••+nr n 1 ,n 2 ,......,nr
is the multinomial coefficient defined by

n 1 +n 2 +•••+nr n 1 ,n 2 ,......,nr = (n 1 +•••+nr)! n 1 !n 2 !•••nr!
, for all nonnegative integers n 1 , . . . , n r . Theorem 2.2 is a natural generalization of the well-known explicit upper bound for classical Ramsey numbers. Indeed, an immediate consequence of Theorem 2.2 (by taking t i = 1 for all i) is the following classical upper bound due to Graham and Rödl [3, (2.48

)] R ([k 1 , 1], . . . , [k r , 1]) ≤ k 1 + • • • + k r -r k 1 -1, . . . , k r -1 • (2) Let k ≥ t ≥ 2 and r ≥ 2 be integers and let R r ([k, t]) = R([k, t], . . . , [k, t] r
).

An immediate consequence of Theorem 2.2 (by taking k

= k 1 = • • • = k n and t = t 1 = • • • = t n ) is the following inequality R r ([k, t]) ≤ t r(k -t) k -t, . . . , k -t (3) Moreover, if t = 1 then R r ([k, 1]) ≤ (rk -r)! ((k -1)!) r • (4)

Exact values

By the so-called Chvátal's result [START_REF]Tree-complete Ramsey numbers[END_REF], we know that the exact value of the Ramsey number of K [START_REF] Huang | New upper bounds for Ramsey numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF] (a star) versus cliques is given by R([n, 1], [START_REF] Huang | New upper bounds for Ramsey numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]) = 3n -2 for all n ≥ 1. We then naturally focus our attention to the Ramsey number of K [START_REF] Huang | New upper bounds for Ramsey numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF] versus cliques with either a dropped edge or a dropped triangle, see [START_REF] Chappelon | On Ramsey numbers of complete graphs with dropped stars[END_REF] where R([m, 1], [n, 2]) has been computed for numerous cases. We provide the new following exact values of Ramsey numbers.

Theorem 3.1 Let n ≥ 2 be an integer. Then, • R([n, 2], [4, 3]) = 2 for n = 2, • R([n, 2], [4, 3]) = 5 for n = 3, • R([n, 2], [4, 3]) = 3n -5 for n ≥ 4. Theorem 3.2 Let n ≥ 2 be an integer. Then, • R([n, 3], [4, 3]) = 3 for n = 3, • R([n, 3], [4, 3]) = 6 for n = 4, • R([n, 3], [4, 3]) = 8 for n = 5, • R([n, 3], [4, 3]) = 11 for n = 6, • R([n, 3], [4, 3]) = 3n -8 for n ≥ 7.
3.1 An estimation for R([n, 2], [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]) By considering K [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF] as the book graph B 3 , it was proved in [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Sudakov | Large K r -Free Subgraphs in K s -Free Graphs abd Some Other Ramsey-Type Problems[END_REF] that R([n, 1], [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]) ≤ 3n 2 log(n/e) , for all positive integers n.

The following result is a first estimation for the value R([n, 2], [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]).

Theorem 3.3 Let n ≥ 2 be an integer. Then,

• R([n, 2], [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]) = 2 for n = 2,

• R([n, 2], [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]) = 7 for n = 3,

• R([n, 2], [START_REF] Li | On Book-Complete Graph Ramsey Numbers[END_REF][START_REF] Graham | Numbers in Ramsey theory[END_REF]) ≤ 3 n+1 2 -5n + 4 for n ≥ 4.

  .

Fig. 1 .

 1 Fig. 1. (a) K [5,3] and (b) K [4,2]

Email: lmontejano@cimat.mx

Email: jramirez@um2.fr

Email: jonathan.chappelon@um2.fr


M Mazar
email: mmazar@cesi.fr

M Sahnoun
email: msahnoun@cesi.fr

B Bettayeb
email: bbettayeb@cesi.fr

N Klement

AFROS2018

Keywords: Optimization, Simulation, Multi-Agent System, Agricultural Robot, Mildew, Ultraviolet

Introduction

Robots are used in several economic fields such as industry, hospitalization and agriculture. This latter, knows a real revolution concerning the robotization of its tasks because of their painfulness and/or dangerousness for humans as well as the desertion of human worker from the agricultural sector. Several research works interested to agricultural robots, aim to optimize plants treatment against diseases [START_REF] Oberti | Selective spraying of grapevines for disease control using a modular agricultural robot[END_REF] and harvesting [START_REF] Southall | An autonomous crop treatment robot : Part i. a kalman filter model for localization and crop/weed classification[END_REF] of crops.

Recent researches in horticulture reveal that the use of type-C ultraviolet (UV-C) treatment against powdery mildew is a sustainable alternative to chemical treatment. The robotization of UV-C treatment in the horticulture field become a necessity due to the hardness of this kind of treatment manually. A new robot (called UV-Robot) is developed to ensure this treatment, while guaranteeing the safety of the human operator against the UV-C radiation. It is able to perform several levels of UV-C treatments by varying its speed when passing behind plants. The robot is wireless and uses a battery that limits the number of plants treated before it needs to be recharged.

The scheduling of UV-Robot tasks represent a real challenge regarding the complexity of the system, the constraint concerning its battery capacity, the procedure of treatment and the dynamic level of diseases.

In this abstract, we study an optimization problem aiming to minimize the necessary number of a robot's missions during a limited period of time, while ensuring the treatment of all infected plants in the greenhouse. This period must be less than one day (24 hour) [Danan, 2009] to grantee a stationary level of the disease in the greenhouse. A period of treatment contains several missions of at least 30 minutes, depending on the number of plants to be trated and their levels of disease. In fast, the higher is the level of disease, the more is the time of its treatment, which implies a higher power consumption. Each mission is composed of a set of treatment and movement tasks. The robot starts each mission from the charging station and should come back to this station before the exhaustion of its battery. The robot's battery is recharged during a relatively long time (4 hour in average) between each two successive missions.

In rest of this abstract, we present the mathematical formulation of the optimization problem stated above in section 2. Then, section 3 contains the description of the solving approach that we propose. We end this abstract by summarizing our main contribution and defining the following steps for this work.

Problem formulation

As mentioned in the introduction, the problem studied consists of optimizing the robotized treatment of plants' diseases in a greenhouse. Figure 1 shows a representation of the greenhouse and the undirected connected graph corresponding to the robot routing between the charging station and the row corridors. The greenhouse is composed of N -1 two-sided crops (plants) rows and N rows of dead-end corridors through which the robot passes to treat the plants and return back to the entry point at the beginning of the row. Each color represents a level of disease that affects the plant. An infected plant should be treated on its both sides facing the two surrounding corridors.

We assume that splitting the treatment is not allowed for individual sides of each plant, neither for the whole plants to be treated in the two sides of the same row corridor, i.e. the robot visits each row corridor at most once.

Let i ∈ {0, 1, .., N } be the index representing the charging station (i = 0) and the row corridors (i ≥ 1). We distinguish two types of tasks for the robot : 1) the treatment of the infected plants in a row corridor and ; 2) the motion between the rows and from the charging station to the rows. Let w ii be the energy consumption in Wh (Watt-hour) corresponding to the treatment of the i th row and w ij the energy needed to move from row i to row j (i = j).

Matrix (1) represents all possible robot's tasks in the greenhouse. The treatment tasks are placed in the diagonal (w ii energy consumption for the treatment of the row i). In a given mission the robot goes from a row i to the row j by executing the motion task consuming w ij energy.

W =       w 00 w 01 • • • w 0N w 10 w 11 • • • w 1N . . . . . . . . . . . . w N 0 w N 1 • • • w N N       (1)
The objective is to define the optimal treatment and charging schedule, which minimizes the total duration of a treatment. This optimization problem can be formulated as MILP as follows :

minimize Z(Y ) = k∈K y k (2)
Subject to :

j∈N i∈N w ij x k ij ≤ Cy k ∀k ∈ K (3) i∈N x k ij = 2 ∀j ∈ N (4) j∈N x k ij = 2 ∀i ∈ N (5) x k 00 = 1 ∀k ∈ K (6) j∈N x k 0j = 2 ∀k ∈ K (7) y k ≥ y k+1 ∀k ∈ K (8)
Where :

-C : Capacity of the battery at the beginning (in Watt) of the mission w ij : Power consumption (in Watt) of the task ij x k ij : a binary decision variable to select tasks

x k ii =
1 if the robot is treating the row i during the mission k 0 otherwise.

x k ij =
1 if the robot goes from i to j ∀ i = j 0 otherwise. y k : a binary decision variable to select the the k th mission 2) is the objective function, which correspond to the minimization of the number of missions. The total energy consumption to perform the tasks of each mission k must not exceed the battery's capacity of the robot (constraint (3)). The two inequalities (4) and ( 5) express the constraint that the robot should execute a moving task between each two plant treatment tasks, i.e. when two tasks are selected in the diagonal of the matrix (1), the task corresponding to their line-column intersection should be also selected. The robot must start all missions from the charging station, this is expressed by the constraint (6). Constraint (6) means that the robot must return after the end of each mission to the charging station.

y k = 1 if the mission k is planned 0 otherwise. Equation (
This problem as formulated above is equivalent to the one-dimensional bin packing problem. This problem is known to be NP-hard [START_REF] Bernhard | Combinatorial optimization : Theory and algorithms[END_REF]. This formulation assumes that the dynamic behavior of the diseases in the greenhouse is negligible and one robot is able to treat all the infected plants in a stationary period of time. However, in the case of big instances (bigger or several greenhouses), these assumptions can not be verified. It becomes mandatory to consider the dynamic evolution of the diseases and the possibility to use several robots, which makes the problem much more complex to solve. However, the developed model can be used with in a simulation-optimization approach where the optimization algorithm in run several times during the simulation process. For that, a simulation model is developed, which is presented in the following section.

Simulation model

A farm composed of several greenhouses adopting robotized treatment present a high level of system dynamics and interaction between its components. In order to reduce its modeling task, a Multi-Agent System (MAS) model is then proposed. As shown in Figure 2, this system contains seven agents interacting with each other. Each agent's behavior depends on its own characteristics and the evolution of its environment.

The greenhouse agent is the place that contains the plants, the robot and the charging station. each plant agent grows in the greenhouse and can be treated by the robot. The robot agent treats plants in the greenhouse with UV-C. The charging station agent is installed in the greenhouse and charges of the robot's battery when it is placed on. The UV-C lamps are a class of agents installed on and controlled by the robot. The monitoring agent schedules missions of the robot and exchanges

Figure 1 -

 1 Figure 1 -Diagram of the greenhouse with the tasks assigned to the robot

Acknowledgment

This research was made possible thanks to €1.35 million financial support from the European Regional Development Fund provided by the Interreg North-West Europe Programme in context of UV-ROBOT

Plants

Robots

Greenhous

Conclusion

In this abstract, we explained a plant treatment system in horticulture, which represents a complex system because it is composed of several interacting entities. Two models are presented : the first is MILP optimization model equivalent to the bin packing problem, and the second is MAS based simulation model. We will implement our two models in a simulator able to simulate and optimize the process of UV-C robotized treatment. the first trials of the simulation will be done on relatively small instances to tune and validate the approach. After, the simulation prototype will be implemented in several real cases, using the robot developed by Octinium (a project partner), to prove the effectiveness of this innovative solution.


Ismail Bensekrane
email: ismail.bensekrane@polytech-lille.fr

Pushpendra Kumar

Achille Melingui
email: achillemelingui@gmail.com

Vincent Coelen
email: vincent.coelen@polytech-lille.netjahan.pushp@gmail.com

Yacine Amara
email: amara.yacine@gmail.com

Rochdi Merzouki
email: rochdi.merzouki@polytech-lille.fr

Energy Planning for Unmanned Over-Actuated Road Vehicle

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.

Energy Planning for Unmanned Over-actuated Road Vehicle

Ismail Bensekrane 1 Pushpendra Kumar 1 Achille Melingui 2 Vincent Coelen 1 Yacine Amara 3 Rochdi Merzouki 1 Abstract-Energy planning for unmanned road vehicles (URV)s is an important step for the management of the autonomous driving. This energy planning depends on the model for power consumption estimation related to URVs. Generally URVs are over-actuated, and this property leads to multiple kinematic configurations for driving. Consequently, it adds more constraints and more complexity for energy planning. In this paper, a Neuro-Fuzzy model is developed for power consumption estimation for different driving modes configurations of URV. Furthermore, a dynamic programming algorithm is applied to find the optimal velocity profile, and the optimal configuration mode in each segment of the road for an over-actuated URV called RobuCAR, used for experimental validation.

Index Terms-Power consumption estimation, Energy planning management, Neuro-Fuzzy modeling, Dynamic programming.

I. INTRODUCTION A. Power Consumption by URV

The power consumption estimation is mainly related to the URV's dynamic parameters. Battery Management System (BMS) [START_REF] Szumanowski | Battery management system based on battery nonlinear dynamics modeling[END_REF] is used to monitor instantly the electric power of batteries for vehicles. Energy Management Strategy (EMS) is used to identify the energy optimal control of an overactuated vehicle [START_REF] Brembeck | Energy optimal control of an over actuated robotic electric vehicle using enhanced control allocation approaches[END_REF], to maximize the traveling distance [START_REF] Chen | Energy management and driving strategy for in-wheel motor electric ground vehicles with terrain profile preview[END_REF]. Both BMS and EMS cannot estimate the instantaneous power consumption of an URV during the trajectory tracking under kinematic and dynamic conditions of driving. The kinematic parameters are essentially represented by the degree of actuation redundancy, degree of steerability, and the road geometry in terms of path curvature [START_REF] Bensekrane | Towards adaptive power consumption estimation for over-actuated unmanned vehicles[END_REF].

The development of a model of power consumption is often elaborated using a linear regression model or exponential logarithmic model [5], [START_REF] Chuy | Power modeling of a skid steered wheeled robotic ground vehicle[END_REF]. Generally, those models can not represent with accuracy due to presence of non-linearities and not considering the driving modes of navigation and road curvature. In [START_REF] Bensekrane | Towards adaptive power consumption estimation for over-actuated unmanned vehicles[END_REF], a model is presented taking into account the path curvature C, degree of steerability δ S , and degree of actuation redundancy δ R , to estimate analytically power consumption. But, analytical models cannot model all the system non-linearities and complexities. Therefore, it is required to look for new methods with better accuracy as Neural network and Adaptive Neuro Fuzzy Inference System (ANFIS) [START_REF] Bilgehan | Comparison of anfis and nn modelswith a study in critical buckling load estimation[END_REF].

Neural network modeling is among many types of modeling methods with better accuracy, and often used to characterize a nonlinear system. Many models have been developed for the power consumption estimation for electric or hybrid vehicles using neural network. In [START_REF] Felipe | Energy consumption estimation in electric vehicles considering driving style[END_REF], a neural network based prediction of energy consumption is developed for an electric vehicle considering the driving style of driver. Sakayori et al. [START_REF] Sakayori | Energy efficient slope traversability planning for mobile robot in loose soil[END_REF] used multilayer perception (MLP) neural network to predict power consumption for an autonomous mobile robot, including the slope, heading angle, and velocity profile as inputs. Masikos et al. [START_REF] Masikos | Reliable vehicular consumption prediction based on machine learning[END_REF] presented a MLP neural network for an electric vehicle to predict the energy consumption in urban area with a 10% of average accuracy. In [START_REF] Togun | Prediction of torque and specific fuel consumption of a gasoline engine by using artificial neural networks[END_REF], neural network is used to estimate torque and brake fuel specific consumption (BFSC) for the gasoline engine. The limitation of these works is the high percentage of average error, which can not be appropriate for autonomous navigation of an over-actuated URV. Furthermore, these works have not considered many possible driving modes configurations of an over-actuated URV. [START_REF] Bilgehan | Comparison of anfis and nn modelswith a study in critical buckling load estimation[END_REF] find that ANFIS present better performances than Neural network modeling. In this work ANFIS is considered to estimate the instantaneous power consumption for different driving mode configurations of URV.

B. Energy Planning Management

In addition to power consumption estimation modeling, it is important to have a good strategy for energy planning management. The main objective of energy planning management is to optimize energy consumption by URV for a given trip. Many researchers developed energy optimization strategies using dynamic programming ( [START_REF] Ozatay | Cloud-based velocity profile optimization for everyday driving: A dynamic-programming-based solution[END_REF], [START_REF] Zhang | Role of terrain preview in energy management of hybrid electric vehicles[END_REF]). Moreover, authors developed energy optimization methods using model predictive control (MPC) and compared with dynamic programming based methods( [START_REF] Chen | Energy management and driving strategy for in-wheel motor electric ground vehicles with terrain profile preview[END_REF], [START_REF] Quaglia | Communication-aware bandwidthoptimized predictive control of motor drives in electric vehicles[END_REF]). Furthermore, the power management has been studied by using control allocation for over-actuated URV based on parametric torque distribution [START_REF] Dizqah | A fast and parametric torque distribution strategy for fourwheel-drive energy-efficient electric vehicles[END_REF]. But, this approach considered only dynamics of the vehicle and not the road profile. Chen et al. [START_REF] Chen | Energy management and driving strategy for in-wheel motor electric ground vehicles with terrain profile preview[END_REF] estimated the optimal velocity profile and torque using dynamic programming method for an URV considering its front steering. In [START_REF] Zhang | Power management comparison for a dual-motor-propulsion system used in a battery electric bus[END_REF], neural network and MPC based method is developed for energy optimization and compared with dynamic programming algorithm. In [START_REF] Liu | Integrated torque vectoring control for a three-axle electric bus based on holistic cornering control method[END_REF], quadratic programming is used to develop energy optimization for URV. In this paper, we are applying a dynamic programming algorithm to find the optimal velocity profile and the optimal driving mode in each segment of the trajectory.

C. Problem Statements and Contribution

Based on the above literature review, we found that not much models are available for power consumption estimation by an URV, in addition, analytical models can not provide precise results due to non-linearities of the system and parametric uncertainties. Moreover, existing methods of energy planning consider optimization of velocity profile for a given trajectory, but optimization of different driving modes have not been considered. In this paper, a model for power consumption estimation is developed using Neuro-fuzzy approach for minimizing the modeling errors. Furthermore, a dynamic programming is developed for energy planning management of an URV, which optimizes the velocity profile and the driving mode for each section of a given trajectory (Fig. 1). In Fig. 1, power is estimated for an URV using a qualitative model (ANFIS), which is compared with the experimental measurements of power consumption obtained from a professional real-time simulation software for vehicle dynamics used by automotive industries namely, SCANeR Studio [START_REF]Oktal website for scaner studio[END_REF]. The inputs to model are driving mode, road profile (curvature of the trajectory), velocity profile, and acceleration for a given trajectory. In the next step, a directed graph is generated to apply the optimization method. For this, road sampling is performed to obtain different segments of it based on curvature of the road considering three types of road segments i.e., line, clothoid, and arc. Then, a 3D grid is constructed with three axes of velocity, arclength of the segment, and driving mode; the purpose of using 3D grid is to consider different driving modes on a third axis, which is not possible with a 2D grid [START_REF] Ozatay | Cloud-based velocity profile optimization for everyday driving: A dynamic-programming-based solution[END_REF]. In this grid, all the feasible points for each driving mode are given. Then, total energy consumption in each segment of the road are calculated, which are weights of digraph edges. Finally, digraph is constructed in the 3D grid, where nodes represent starting and ending points for the different segments of the road profile. Dynamic programming is applied on the digraph to find the optimal global solution in terms of optimal velocity Fig. 2: RobuCAR and optimal driving mode for each segment of the road profile.

D. Paper Organization

The remaining paper is organized in the following sections. In Section II, power consumption model of an URV is developed. In section III, the energy planning management algorithm is developed. Section IV discusses about simulation results. Section VI concludes this paper with a summary of our observations.

II. POWER CONSUMPTION ESTIMATION

In this section, Neuro-fuzzy model is developed for the power consumption estimation of a four wheeled overactuated URV. The model considers the degree of actuation redundancy δ R , degree of steerability δ S , and the path curvature C; for definition of these terms one can refer to [START_REF] Bensekrane | Towards adaptive power consumption estimation for over-actuated unmanned vehicles[END_REF].

In this paper, we analyze a four wheeled URV called RobuCAR (Fig. 2). RobuCAR is a four wheeled autonomous vehicle with on-board computer and electronics. The four wheels are independently driven by DC motors. There are two steering control systems for the steering of front and rear wheels. There are total six input actuators (traction and steering) to control its motion in a plane with three degrees of freedom (DoF) i.e., longitudinal, lateral, and yaw motions; therefore, RobuCAR is an over-actuated URV. Moreover, it is a non-holonomic URV because it cannot move in the three DoF without steering. According to δ S and δ R , RobuCAR can have different possible driving modes configurations as shown in Fig. 3.

4,

2, 3

T S R N N     2 S   4, 1, 2 T S R N N     1 S   4, 0, 1 T S R N N     0 S  

Fig. 3: Different types of configuration of RobuCAR

Refer to Fig. 3, N T and N s denote total number of active traction and steering actuators in a configuration, respectively. The Dual configuration (δ S = 2) has two steering control systems for the steering of front and rear wheels. In the Single configuration (δ S = 1), there is only one steering control system for the steering of front wheels. The Skid configuration (δ S = 0) has no steering system, and all the four wheels are the fixed tractive wheels. We represented three configuration with a number 4, which denotes that all the four wheels are tractive. In addition, δ R for each configuration is shown in Fig. 3, it can be observed that δ R depends on δ S .

A. Adaptive neuro fuzzy inference system (ANFIS) ANFIS have been successfully used for approximation tasks [START_REF] Jang | Anfis: adaptive-network-based fuzzy inference system[END_REF], [START_REF] Buragohain | A novel approach for anfis modelling based on full factorial design[END_REF]. In this paper, ANFIS is used to estimate power consumption of an URV. The ideal may be to build an ANFIS model which emulates the different URV actuation modes. However, as we will see in the simulation results, the URV behavior from one actuation mode to another is so different that three separate ANFIS models, one for each actuation mode, have been built to improve regression performance. In the following subsections the ANFIS structure and the learning algorithm are presented.

1) ANFIS structure: The ANFIS structure considered in this work is shown in Fig. 4. The inputs to ANFIS model are curvature of the road C, velocity V and acceleration A of the URV, while power consumed by URV is the output of the ANFIS model. Takagi and Sugeno's type rules are considered [START_REF] Jang | Anfis: adaptive-network-based fuzzy inference system[END_REF] as follows:

Rule1: If C is C 1 and V is V 1 and A is A 1 , then P ω is p 1 C + q 1 V + r 1 A + s 1 .
Where, p 1 , q 1 , r 1 , and s 1 are constants, while, V 1 , A 1 , and C 1 , are the fuzzy sets in the antecedents, and P ω = f (C,V, A) is a crisp function in the consequent.

N N C 1 V 1 C n V n A 1 A n . . . . . . . . .

Fig. 4: ANFIS structure for one driving mode

The output of each ANFIS layer is evaluated as follows: Layer 1: Each i th node in this layer is given by:

O 1 i = µ A i (x) (1)
where, x is the input to i th node, A i is the linguistic variable associated with this node function, and µ A i is the membership function of A i . The generalized bell-shaped membership function is chosen for µ A i (x) :

µ A i (x) = 1 1 + (x -c i /a i ) 2 b i (2)
where, {a i , b i , c i } is the premise parameter set. Layer 2: Each node in this layer calculates the firing strength w i of a rule. The output of each node is the product of all the incoming signals to it:

O 2 i = w i = µ C i (C) × µ V i (V ) × µ A i (A), i = 1, 2, ..., n.
(3) Layer 3: The output of each node is the normalized firing strength given by:

O 3 i = w i = w i w 1 + w 2 + • • • + w n (4)
Layer 4: Every i th node in this layer is an adaptive node with a node function given by:

O 4 i = w i f i = w i (p i C + q i V + r i A + s i ) (5)
where, w i is the output of Layer 3 and {p i , q i , r i , s i } is the consequent parameter set.

Layer 5: This layer comprises one node, and computes the overall output as the summation of all incoming signals, i.e,

O 5 i = overallout put = ∑ i w i f i = ∑ i w i f i ∑ i w i (6)
2) Learning algorithm: By observing the ANFIS structure, it can be observed that the values of the premise parameters and the output of the system can be expressed as a linear combination of the consequent parameters

(p 1 , q 1 , r 1 , s 1 , . . . , p n , q n , r n , s n ). f = w 1 w 1 + w 2 + • • • + w n f 1 + w 2 w 1 + w 2 + • • • + w n f 2 + • • • + w n w 1 + w 2 + • • • + w n f n = w 1 f 1 + w 2 f 2 + • • • + w n f n = (w 1 C) p 1 + (w 1 V ) q 1 + (w 1 A) r 1 + (w 1 ) s 1 + (w 2 C) p 2 + (w 2 V ) q 2 + (w 2 A) r 2 + (w 2 ) s 1 + • • • + (w n C) p n + (w n V ) q n + (w n A) r n + (w n ) s n (7)
In the learning process, the consequent parameters are identified by using the least squares estimate in the forward pass. In the backward pass, the premise parameters are updated by mean of the gradient descent algorithm using the error signals [START_REF] Jang | Anfis: adaptive-network-based fuzzy inference system[END_REF].

3) Learning and validation of ANFIS models: The database for ANFIS model learning is built by using the engineering software SCANeR studio.

Different training paths containing data of the curvature points in the vehicle's trajectory are considered. In SCANeR studio environment, the velocity profile in range [START_REF] Szumanowski | Battery management system based on battery nonlinear dynamics modeling[END_REF][START_REF] Zhang | Power management comparison for a dual-motor-propulsion system used in a battery electric bus[END_REF] km/h with a step size of 1 km/h is used for each driving mode; the velocity profile, acceleration, path curvature, and the consumed power are recorded. We obtained the databases of 262079, 264634, and 400965 samples respectively for, Dual mode, Single mode, and Skid mode. Each database is divided into three subsets: training, validation and test sets. The training set is used during the learning phase, and the test set is only used to evaluate the performance of ANFIS models. For a good generalization and to avoid over-fitting, the validation set is used during the training phase and the early-stopping method is applied for training. The earlystopping method requires that after a period of training (an epoch) using the training set, the weight matrices of ANFIS are fixed, and the ANFIS operates in the forward mode using the validation set. The process is reiterated until the mean square error (MSE) on the validation set reaches its minimum value.

For ANFIS parameters initialization, three membership functions are considered for each linguistic variable, and the grid partitioning method is used to initialize the fuzzy rules, the premise and consequent parameters. The different ANFIS models converge approximately after 2000 epochs. The following MSE, 2.56e -3 , 2.27e -3 , 3.03e -2 respectively for, Dual mode, Single mode, and Skid mode are obtained on the corresponding test sets.

In Fig. 5(a,c,e), a comparison of power consumption of the URV using ANFIS and SCANeR is given for Skid, Single, and Dual configurations, respectively. The errors between these models are shown in Fig. 5(b,d,f ). From the results, it can be observed that ANFIS model follows almost the same profile as SCANeR. However, there are peaks of errors at some points due to turns in the trajectory.

III. ENERGY PLANNING MANAGEMENT

The energy planning strategy considers the dynamics of URV and the road geometry for each segment of the planned trajectory. To reach this goal, a cost function is given to minimize. min

V,M J = n ∑ i=1 E k (t i ,V i ,C i , M i,k ) (8)
s.t. V min ≤ V i ≤ V max T min ≤ T i ≤ T max a min ≤ a i ≤ a max (9) E k (s i ,V i ,C i , M i,k ) =
Where, E k and P k represent energy consumed and instantaneous power consumption for the kth configuration of URV, respectively. The parameters t i , V i , and C i represent traveling time, actual velocity, and curvature in ith segment, respectively. T i represents the total torque for the four wheel drive configuration with lower and upper limits T min and T max , respectively. a i represents the total acceleration for the four wheel drive configuration with lower and upper limits a min and a max , respectively.

A. Road Profile

In a Geographical Information System (GIS), a road is represented by polylines [START_REF] Jakkula | Extraction of road network topology for transportation and gis applications[END_REF], describing the central axis of its surface. The road path is decomposed into horizontal and vertical curves representing respectively the roads mapped onto the plane (Fig. 6(a)). In this paper, we use the cubic spline interpolation approach for polylines describing the considered trajectory [START_REF] Ahlberg | The theory of splines and their application[END_REF].

Fig. 1 :

 1 Fig.1: The proposed methodology for power consumption estimation and energy planning management.

  Fig. 5: Power estimation

  0, δ R = 1 Skid(Over-actutaed) δ s = 1, δ R = 2 Single(Over-actutaed) δ s = 2, δ R = 3 Dual(Over-actutaed)

Fig. 6 :

 6 Fig. 6: (a) Campus trajectory of the University of Lille and (b) the filtered and real curvature

Fig. 11 :

 11 Fig. 11: Percentage of modes use in campus trajectory

978-1-5386-6203-8/18/$31.00 ©2018 IEEE

B. Trajectory Filtering

The GIS data for real road is often tainted by errors of acquisition systems, which explains the fluctuations presented in the estimated curvature of the 2D trajectory. The curvature of the trajectory can take three forms:

where, s represents the arc length; while a and b are constant parameters of clothoide in Cs plane. The estimated and filtered curvature for the campus of university of Lille 1 is shown in Fig. 6 (b) (Cs plane), with total of 138 segments.

In each segment of the trajectory, a unique drive pattern is applied (acceleration, constant velocity, or deceleration) as shown in Fig. 7. In addition, it is considered that each segment (line, clothoid, and arc) has constraint of driving mode, which are detailed in Table I.

C. Dynamic Programing

To solve the non-linear problem of energy planning, we apply dynamic programming algorithm on weighted digraph.

1) 3D grid construction: The grid consists of 3 axes, arclength which is represented with a set of samples of distance d = {d 0 , d 1 , d 2 , ..., d n }, velocity v = {v min : ∆v : v max }, and driving mode m = {Skid, Dual, Single}. Each point in the grid is obtained with respect to all constraints namely, geometrical, kinematic, and dynamic such as velocity, acceleration, and torque; as shown in Fig. 8. In this approach, arclength is segmented according to the geometrical criteria, which is usually non uniform segmentation.

2) Digraph construction: A digraph is a weighted directed graph, where the edges have a direction associated with them. A directed graph is an ordered pair G = (Ve, Ed), where Ve is a set of vertices representing end points of each segment and Ed is a set of edges representing connection between vertices as a function of velocity, mode, and arclength. The weight on each edge represents energy consumed to travel the corresponding segment for each mode.

Each vertex is denoted with a specific number say, jkl, where j, k, and l represent jth row of velocity vector, kth row of mode vector, and lth row of segment vector, respectively. In addition, a staring and an ending vertex is added on the top and the bottom of the digraph. The starting vertex is connected to the first set of vertices with zero weights on edges. Similarly, the ending vertex is connected to the last set of vertices with zero weights on edges. The developed digraph is used to find shortest path between the starting and the ending points by applying dynamic programming.

IV. SIMULATION BASED ON EXPERIMENTAL DATA

A set of experimental data from embedded sensors of the RoboCAR during its autonomous driving along the trajectory of Fig. 6(a) are collected for three kinematic configurations. These data are then used for power consumption estimation and energy planning for the driving management.

In Fig. 9 and 10, the optimal velocity profile (4-16 km/h) and the optimal configuration modes are presented for each segment of the filtered trajectory in s -V and X -Y planes, respectively. The algorithm shows that Single mode is the most optimal for the given optimal velocity profile as shown in Fig. 11. The Dual mode should be used in sharp turns and Skid mode in some line segments along the trajectory. Furthermore, the energy consumption for each driving mode is given in Table II, while considering only one mode for the whole trajectory tracking. From the Table II, it can be observed that the proposed optimal driving strategy with changing driving modes consumes less energy in comparison to the one mode driving.

TABLE II: Comparison of energy consumption

Energy consumption (J) Skid for optimal velocity 3.1573e+06 Single for optimal velocity 3.0279e+05 Dual for optimal velocity 4.5242e+05 Optimal modes and optimal velocity 1.6912e+04

V. CONCLUSION

This paper presents a method of energy planning of an over-actuated URV based on power consumption estimation model. The power consumption is estimated using a qualitative ANFIS model for different configuration modes of an URV called RobuCAR. Dynamic programming algorithm is applied to obtain the global optimal velocity profile and the optimal driving mode associated with each segment of the trajectory in off-line mode. The combinations of optimal driving modes and optimal velocity gives optimal energy consumption for the whole trajectory tracking in comparison to one mode driving. For future work, an extension to the case of faulty actuator will be considered for the energy planning strategy.


Wei Wang

Bernard Bonello
email: bernard.bonello@insp.jussieu.fr

Bahram Djafari-Rouhani

Yan Pennec

Jinfeng Zhao

Double-negative Pillared Elastic Metamaterial

We present a single-phase double-sided pillared metamaterial with effective mass density and elastic modulus simultaneously negative. The negative effective mass density is achieved owing to the combination of the symmetric bending resonance and the antisymmetric compressional resonance of both pillars, whereas the negative effective elastic modulus results from their symmetric compressional resonance. Comparisons with the single-sided pillared structure are made. The dependence of the width of the double-negative band against the geometric parameters of the unit cell is also investigated. In addition, the negative refraction and the zero-index refraction of the symmetric Lamb wave in the deepsubwavelength scale are observed, while at the same frequency, the propagation of the antisymmetric Lamb wave is forbidden.

I. INTRODUCTION

Elastic metamaterials are artificial composites featuring frequency dependent constitutive properties that are not present in nature. Owing to their locally resonant substructures, their effective parameters, including the mass density, the bulk, and the shear modulus, can be dynamically set to a positive, zero or negative value [START_REF] Wu | Evidence of complete band gap and resonances in a plate with periodic stubbed surface[END_REF][START_REF] Yao | Experimental study on negative effective mass in a 1D mass-spring system[END_REF][START_REF] Huang | Wave attenuation mechanism in an acoustic metamaterial with negative effective mass density[END_REF][START_REF] Zhu | Microstructural design and experimental validation of elastic metamaterial plates with anisotropic mass density[END_REF], enabling thus to manipulate the propagation of elastic waves in the subwavelength scale. In analogy to their electromagnetic counterparts, elastic metamaterials with simultaneously negative effective mass density (NMD) and elastic modulus (NEM) have attracted considerable attention notably because of their great potential for the negative refraction of elastic waves or the over-the-diffraction-limit imaging [START_REF] Fang | Ultrasonic metamaterials with negative modulus[END_REF][START_REF] Lee | Composite acoustic medium with simultaneously negative density and modulus[END_REF][START_REF] Yang | Coupled membranes with doubly negative mass density and bulk modulus[END_REF][START_REF] Zhu | Negative refraction of elastic waves at the deep-subwavelength scale in a single-phase metamaterial[END_REF][START_REF] Dong | Topology optimization of anisotropic broadband double-negative elastic metamaterials[END_REF][START_REF] Zhang | Investigation of surface acoustic wave propagation in composite pillar based phononic crystals within both local resonance and Bragg scattering mechanism regimes[END_REF]. The double negativity can be realized either by combining two different substructures, each supporting a different resonant mode or by constructing a single structure where two resonances occur at a single frequency [START_REF] Ding | Metamaterial with simultaneously negative bulk modulus and mass density[END_REF][START_REF] Jin | Pillar-type acoustic metasurface[END_REF]. In the wake of the seminal work by Liu et al. [START_REF] Liu | Locally resonant sonic materials[END_REF], a three-constituents chiral metamaterial made of soft-coated heavy cylinders embedded into a polyethylene matrix that exhibits the double negativity thanks to the combination of the translational and the rotational resonances of the heavy core has been proposed [START_REF] Liu | An elastic metamaterial with simultaneously negative mass density and bulk modulus[END_REF]. Using the same two resonances but this time in a single-phase metamaterial, Zhu et al. [START_REF] Zhu | Negative refraction of elastic waves at the deep-subwavelength scale in a single-phase metamaterial[END_REF] observed experimentally the subwavelength negative refraction of longitudinal elastic waves and demonstrated in this way the occurrence of simultaneous NMD and NEM in this system.

A single-phase solid metamaterial which unit cell comprises two independent resonators, each of them supporting a distinct resonant mode, was proposed by Oh et al. [START_REF] Oh | Elastic metamaterials for independent realization of negativity in density and stiffness[END_REF] to independently tune the effective mass density and stiffness. However, in most of these aforementioned systems, the double negativity only occurs in a narrow frequency band which may constitute a drawback for many applications. To overcome this difficulty, Dong et al. [START_REF] Dong | Topology optimization of anisotropic broadband double-negative elastic metamaterials[END_REF] developed a topology optimization scheme to design a two-dimensional (2D) single-phase anisotropic elastic metamaterial with broadband doubly negative effective property. They further demonstrated numerically the superlensing effect in the deep-subwavelength regime as well as the cloaking effect for both the longitudinal and transverse waves.

Beside these designs that mainly aim at manipulating bulk waves, doubly negative elastic metamaterials involving other types of elastic waves, such as Lamb waves, have attracted since a few years more and more interest. Actually, a decade ago a new type of phononic crystal consisting of an array of pillars on a plate was simultaneously proposed by Pennec et al. [START_REF] Pennec | Low-frequency gaps in a phononic crystal constituted of cylindrical dots deposited on a thin homogeneous plate[END_REF] and Wu et al. [START_REF] Wu | Evidence of complete band gap and resonances in a plate with periodic stubbed surface[END_REF]. These structures have the ability of exhibiting both Bragg gaps and hybridization gaps associated to local resonances of the pillars, where the latter can be pushed into low frequency (subwavelength) range. These structures can therefore serve both as phononic crystals or metamaterials and provide a useful and practical platform for several physical properties such as waveguiding, filtering, acoustic lensing or liquid sensing. A few years ago, Gusev et al. [START_REF] Gusev | Double-negative flexural acoustic metamaterial[END_REF] introduced a doubly negative metamaterial for flexural waves based on resonant elements attached to a plate. In this system, negative density is achieved by normal-force interactions, whereas negative modulus is obtained thanks to the introduction of rotationally resonant mechanical elements. More recently, a class of 2D single-phase metamaterials have shown noteworthy dispersion properties including a Dirac-like cone at the center of the Brillouin zone (BZ), zero-group velocity in the fundamental Lamb modes [START_REF] Zhu | Phononic thin plates with embedded acoustic black holes[END_REF], and anomalous refraction on transmitted antisymmetric Lamb mode involving mode conversion of a symmetric mode [START_REF] Zhu | Anomalous Refraction of Acoustic Guided Waves in Solids with Geometrically Tapered Metasurfaces[END_REF]. Among the most promising dynamic systems for applications requiring the full control over the phase are the zero-index elastic/acoustic metamaterials [START_REF] Zhu | Double-Zero-Index Structural Phononic Waveguides[END_REF]. Actually, the wave transmitted through a medium exhibiting this dynamic property remains unaffected since it is not subjected to any phase change. As a consequence of this remarkable property, these structures could serve as a basis to shape an "inaudibility cloak". Let us also mention a few other relevant works using the platform of phononic plates. For instance, Assouar et al. [START_REF] Assouar | Enlargement of a locally resonant sonic band gap by using doublesides stubbed phononic plates[END_REF] reported in a theoretical work the enlargement of the locally resonant band gap in 2D sonic crystals based on a double-sided stubbed plate.

The relative bandwidth was significantly increased by a factor of two compared to the classical one-sided stubbed plate. The elastic energy carried by SAW can also be localized and confined within micropillars, as it has been demonstrated lately by Benchabane et al. [START_REF] Benchabane | Surface-Wave Coupling to Single Phononic Subwavelength Resonators[END_REF]. Bilal et al. [START_REF] Bilal | Trampoline metamaterial: Local resonance enhancement by springboards[END_REF] proposed a metamaterial in which the pillars stand in between the holes regularly patterned in a plate. This structure, sometimes referred as a trampoline, allows to broaden the subwavelength band gap by a factor of four. Furthermore, it has been shown that such pillared metamaterials exhibit NMD at low frequencies where homogenization theories apply [START_REF] Oudich | Negative effective mass density of acoustic metamaterial plate decorated with low frequency resonant pillars[END_REF][START_REF] Oudich | Rayleigh Waves in Phononic Crystal Made of Multilayered Pillars: Confined Modes, Fano Resonances, and Acoustically Induced Transparency[END_REF]. It is noteworthy that, owing to the coupling between the Lamb waves and the resonant modes, a strong anisotropy of the effective mass density may be observed around the resonant frequencies. This makes the pillared structures suitable for the manipulation and the subwavelength focusing of elastic wave as well as excellent candidates for numerous applications ranging from biomedical imaging to energy harvesting [START_REF] Ramahi | Metamaterial particles for electromagnetic energy harvesting[END_REF][START_REF] Carrara | Metamaterialinspired structures and concepts for elastoacoustic wave energy harvesting[END_REF].

In this work, we have investigated a single-phase double-sided pillared metamaterial. We show that this structure exhibits both NMD and NEM allowing for the propagation of the zero-order symmetric Lamb wave in a frequency interval where the propagation of the antisymmetric Lamb wave keeps forbidden. For comparison, we firstly have investigated a single-sided pillared metamaterial featuring NMD. The mechanism responsible for this single negativity, and in turn for the formation of the low frequency band gap, is described as being the combination of a bending mode and a compressional mode of the pillar. We then have investigated a double-sided pillared metamaterial that displays the double negativity. On the basis of the preceding analysis, it is shown that NMD is achieved within a frequency range where occur both the symmetric bending resonance and the antisymmetric compressional resonance of the pillars, whereas NEM only involves the symmetric compressional resonance. We develop a simplified massspring model to theoretically describe this behavior. Additionally, the double negativity of the metamaterial is further established through both numerical simulations of the negative refraction of symmetric Lamb mode and the zero-index refraction in the deep-subwavelength regime. Finally, the width of frequency range where the metamaterial is doubly negative against the geometric parameters of the unit cell is presented.

II. SINGLE-SIDED PILLARED METAMATERIAL

A. Lamb waves in a periodic structure

Before analyzing the band structure and the dynamical properties of a pillared metamaterial, we briefly recall in this section the main features of Lamb waves in periodically structured plates [START_REF] Royer | Elastic Waves in Solids[END_REF].

Lamb waves are the solutions of the equations of motion for a harmonic wave propagating in a medium with a finite thickness along one direction of space and infinite along the other two, with the boundary conditions that the components of the stress in the direction normal to the free surfaces are zero. Owing to the symmetry with respect to the mid-plane of the plate, these modes can be distinguished as symmetric (S) and antisymmetric (A) modes. In a homogeneous isotropic medium, the three lowest branches are respectively the antisymmetric flexural branch and the symmetric shear horizontal and longitudinal branches. The higher branches have a cut-off frequency as the wave vector k approaches 0. All these modes are highly dispersive [START_REF] Royer | Elastic Waves in Solids[END_REF]. The phononic plate is created by inserting periodic inclusions (such as holes) in the plate or by depositing pillars on it. The dispersion curves of this periodic inhomogeneous medium are now described in the corresponding 2D BZ with their dynamical properties having specific characteristics including band gaps, the folding of the dispersion curves and the mode conversion. In particular, the strong resonances associated with pillars can be very efficient in opening band gaps.

Although the branches have more complex character and coupling between the displacement components than in a homogeneous plate, some branches still keep a significant flexural, longitudinal or shear horizontal character. In the case of one-sided pillars, the symmetry with respect to the mid-plane of the plane is strictly speaking broken, while this symmetry property can be restored in the case of double-sided pillars. This allows to study separately the symmetric and antisymmetric modes and in particular the design of their corresponding band gaps. We shall see that this property will be very helpful in the design of the double negative material which is the objective of this paper.

B. Band structure and localized modes

We analyze in this section the dynamical properties of a metamaterial made of a square lattice with pillars erected on one face of a thin homogeneous plate. The elementary unit cell is displayed in Fig. 1(a) where a denotes the lattice constant, h and d stand for the height and the diameter of the pillar respectively, and e is the thickness of the matrix plate. The periodicity is in the plane (x, y) and z-axis is chosen perpendicular to the plate. In order to obtain dynamic properties in the MHz range, we assumed these values (in µm) a = 200, h = 325, d = 120, and e = 100. The complete structure was made of isotropic steel whose physical parameters are the elastic modulus E = 200GPa, the Poisson's ratio v = 0.3, and the mass density ρ = 7850kg/m 3 . Both the plate and the pillars being made from the same material, a strong coupling between the Lamb waves propagating in the plate and the vibrations of the pillars may be expected close to the resonant frequencies [START_REF] Jin | Pillar-type acoustic metasurface[END_REF]. This configuration much differs from composite heterostructures where a soft rubber provides the bonding between the resonators and the stiff plate [START_REF] Li | Expansion of lower-frequency locally resonant band gaps using a double-sided stubbed composite phononic crystals plate with composite stubs[END_REF][START_REF] Assouar | Broadband plate-type acoustic metamaterial for low-frequency sound attenuation[END_REF][START_REF] Hu | Enlargement of the locally resonant Lamb wave band gap of the phononic crystal plate at the deep sub-wavelength scale[END_REF]

∫∫∫ ∫∫∫

. A complete band gap opens up in between 3.17MHz to 3.48MHz. In this frequency range, the wavelengths of the lowest order antisymmetric (A 0 ) and symmetric (S 0 ) Lamb modes in the matrix plate are about 500 μm and 1800 μm respectively i.e. more than twice and nine times the lattice constant a. Therefore, this gap cannot result from Bragg interferences. This rather suggests that this complete band gap should be ascribed to local resonances of the pillars. To support this assertion, we have computed the displacement field distribution and the deformation of the unit cell at point M of the BZ for the three branches converging to f~3.1MHz (see inset in Fig. 1(b)). The result is displayed in Fig. 1(c).

Clearly, the deformation concentrates in the pillar whereas it is almost zero in the plate. This is further supported by the flat branches in the low frequency edge of the gap that correspond to zero group velocities and are evidences of local resonances as widely documented in literature [START_REF] Assouar | Enlargement of a locally resonant sonic band gap by using doublesides stubbed phononic plates[END_REF][START_REF] Oudich | A sonic band gap based on the locally resonant phononic plates with stubs[END_REF][START_REF] Pourabolghasem | Physics of band-gap formation and its evolution in the pillar-based phononic crystal structures[END_REF][START_REF] Lee | Origin of negative density and modulus in acoustic metamaterials[END_REF].

Furthermore, the deformations shown in Fig. 1(c) allow to identify two second-order bending modes, labelled as A and B in Fig. 1(c), and one compressional mode, labelled as C. Given the square symmetry of the unit cell, modes A and B are polarized along orthogonal directions and occur at the same frequency.

C. Formation of the low frequency complete band gap

We investigate here the formation of the low frequency band gap using FEM with frequency response analysis. The method, which is only valid within the long wavelengths limit, consists of applying an external displacement field U on the four lateral boundaries of the unit cell while leaving free both the face supporting the pillar and the rear face of the plate [START_REF] Zhu | Microstructural design and experimental validation of elastic metamaterial plates with anisotropic mass density[END_REF][START_REF] Oudich | Negative effective mass density of acoustic metamaterial plate decorated with low frequency resonant pillars[END_REF][START_REF] Zhu | 10126 -Revised version 27 in a metamaterial plate with anisotropic mass density[END_REF]. The induced force F exerted on these boundaries is then derived from the stress average over the four boundaries of the unit cell considered as a homogeneous media. In the harmonic regime at frequency 2 ω π , F and U are related through [ ]

2 eff V ω ρ = - F U ,
0 exp i U i t ω = U
x (i = 1,2,3) along one of the axes in the principal coordinate system while the other two components are set to zero. In addition, symmetry arguments allow to reduce the number of independent matrix elements. For the structure we consider here, z-axis being a four-fold symmetry axis the matrix elements are invariant through the exchange of indices 1 and 2 and therefore 11 This is in good agreement with the band gap shown in Fig. 1(b) that goes from 3.17MHz to 3.48MHz. The small discrepancy of about 2.8% on the frequency at the lower edge of the gap may be readily ascribed to the phase change across the unit cell, not accounted for in the calculations since the preceding theory is only valid within the long wavelengths limit.

Beside this comprehensive approach where the structure is considered as an effective medium, one can also explain the formation of the low frequency band gap owing to a local analysis of the motions of the pillars. In this framework, the resonator is regarded as a mass-spring system coupled to the plate where a wave at frequency 2 ω π and wave vector k propagates. Let us first consider the first-order compressional resonance shown in Fig. 1(c) for which the pillar vibrates along the z-direction and exerts a vertical traction force on the matrix plate. This can be modeled as a mass-spring subsystem with mass m and spring stiffness k C depicted by the right panel in Fig. 3(a). Gusev et al. [START_REF] Gusev | Double-negative flexural acoustic metamaterial[END_REF] already set the equations of motion of this corresponding mechanical system and obtained an analytical expression for the effective mass density component 33 ρ . The governing equations of this simplified model in the frequency domain can be written as:

( ) ( ) 4 2 4 2 C C ˆˆˆẑ z z z D u e u nk u u x m u k u u ρ ω ω ∂ - = - ∂ = - ⎧ ⎪ ⎨ ⎪ ⎩ . ( 1
)
In Eq. ( 1), Fourier transformations of both the vertical displacement of the matrix plate z u and mass u allow to eliminate the time dependence and lead to the variables ˆz u and û written in the frequency domain;

( ) 1), namely:

1 2 1 2 2 2 2 C 1 e k m e D n ω ρ ω ω ⎡ ⎤ ⎢ ⎥ ⎛ ⎞ ⎢ ⎥ = ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ - ⎢ ⎥ ⎣ + ⎦ ⎛ ⎞ ⎜ ⎟ ⎝ ⎠ . ( 2
)
Equation [START_REF] Yao | Experimental study on negative effective mass in a 1D mass-spring system[END_REF] shows that k 2 turns imaginary if ω takes a value in between C ω and ( )

1 2 C 1 nm e ω ρ +
, hence preventing the propagation of the elastic waves in the plate. Moreover, Eq. ( 2) allows to define the effective mass density as:

( )

2 2 C 1 eff nm e ρ ω ρ ω ω - = + ⎛ ⎞ ⎜ ⎟ ⎝ ⎠ , (3)
which is negative when the excitation frequency is in the same interval. This explains the occurrence of the negative value of component 33 ρ . This simplified model can be further extended to the case of a 2D phononic crystal by setting periodic conditions and implementing a plane wave expansion method to derive the dispersion curves [START_REF] Xiao | Flexural wave band gaps in locally resonant thin plates with periodically attached spring-mass resonators[END_REF]. The solution is represented as black solid lines in Fig. 1(b).

It is less intuitive to establish the governing equations that accurately describe the complex motions of both the pillar and the matrix plate for the second-order bending resonance. However, the problem can be simplified by equating the single-sided pillared unit cell to the classical cantilevered beam in which the rigid body is replaced by a deformable medium and where the oscillation of the pillar is excited by the elastic wave that propagates in the plate. In this simplified model depicted in Fig. 3(b), the overall motion is broken down into two motions that can be studied independently, namely the motion of the pillar modeled by an internal moment and the motion in the plate driven by a traction force parallel to the surface. This decoupling allows for a comprehensive analyze of the effective mass density and a better understanding of the occurrence of NMD in this type of structures.

Let us first consider the motion of the pillar in the sagittal plane depicted by the mid panel in Fig. 3(b).

Since no traction force is considered, the motion reduces to a pendulum motion and can be well modeled by a mass-torsional spring with mass m, rotational inertia J and torsional stiffness k B rigidly attached to the plate. Then, the governing equations for waves propagation along x-direction in the plate-resonator coupled system can be written as [START_REF] Qian | Bandgap properties in simplified model of composite locally resonant phononic crystal plate[END_REF]:

4 2 2 B 4 2 2 B ˆẑ z z z u u D e u n k x x x u J k x θ ρ ω ω θ θ ⎧ ⎛ ⎞ ∂ ∂ ∂ - = - ⎪ ⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ ⎪ ⎝ ⎠ ⎨ ∂ ⎪ ⎛ ⎞ = - ⎜ ⎟ ⎪ ∂ ⎝ ⎠ ⎩ , (4)
where ˆz u and θ are normal displacement of the matrix plate and rotation angle of torsional spring in frequency domain; k B and J determine the resonant frequency of the mass-torsional spring subsystem through

1 2 B B k J ω ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠
. By solving Eq. ( 4), we can derive the wave number k, namely:

1 2 2 2 2 2 2 B B 2 2 B B 2 2 1 4 2 1 1 nJ nJ k D e D ω ω ω ω ω ω ω ρ ⎧ ⎫ ⎡ ⎤ ⎛ ⎞ ⎪ ⎪ ⎢ ⎥ ⎜ ⎟ ⎪ ⎪ = - ⎢ ⎥ ⎨ ⎬ ⎜ ⎟ ⎢ ⎥ ⎪ ⎪ ⎜ ⎟ - - ⎢ ⎝ ⎠⎥ ⎪ ⎪ ⎣ ⎦ ⎭ + ⎩ . ( 5
)
Equation ( 5) allows to define the effective mass density at frequency 2 ω π as being:

( )

2 2 2 2 2 2 B 1 4 1 eff n J De ω ρ ω ρ ω ω ⎛ ⎞ ⎜ ⎟ = + ⎜ ⎟ - ⎜ ⎟ ⎝ ⎠ . ( 6
)
Equation [START_REF] Lee | Composite acoustic medium with simultaneously negative density and modulus[END_REF] shows that eff ρ is always positive whatever the frequency. Therefore, the resonance of the torsional spring cannot explain on its own that 11 ρ becomes negative in a specific frequency range, and in turn the occurrence of NMD. Therefore, it is essential to consider the traction force acting on the foot of the pillar to understand the occurrence of NMD evidenced by the simulations above.

This force excites the plate into an oscillation and the pillar can be regarded as an additional mass M attached to the unit cell by compression springs with stiffness k H . For waves propagating along x-direction, the problem reduces to the classic 1D mass-spring metamaterial depicted in the right panel in Fig. 3(b).

One can show that the effective mass density can be written as [START_REF] Yao | Experimental study on negative effective mass in a 1D mass-spring system[END_REF]:

0 2 2 H 1 ( ) eff M m ρ ω ω = + - , (7)
where m 0 stands for the mass of the matrix in the unit cell and ( ) ( )

1 2 H H 0 k M m ω = +
is the angular frequency of the assembly "unit cell+additional mass". Clearly, the effective mass density turns negative

III. DOUBLE-SIDED PILLARED METAMATERIAL

As discussed in the previous section, the single-sided pillared metamaterial can achieve NMD but not NEM. In order to simultaneously realize both negative properties, we have considered a double-sided pillared metamaterial that, under certain conditions [START_REF] Assouar | Enlargement of a locally resonant sonic band gap by using doublesides stubbed phononic plates[END_REF][START_REF] Li | Expansion of lower-frequency locally resonant band gaps using a double-sided stubbed composite phononic crystals plate with composite stubs[END_REF][START_REF] Hu | Enlargement of the locally resonant Lamb wave band gap of the phononic crystal plate at the deep sub-wavelength scale[END_REF][START_REF] Ma | Extending and lowing band gaps in one-dimensional phononic crystal strip with pillars and holes[END_REF], enables enlargement of the width of the low frequency band gap. We show in this section that the double negativity can be obtained through an appropriate choice of the dimensions of the pillar and matrix plate, both of which being made of steel. The unit cell is shown in Fig. 4(a). Two identical pillars are symmetrically arranged on both sides of the matrix plate. The geometric parameters of the unit cell are the same as those of the above single-sided pillared metamaterial. The corresponding band structure is shown in Fig. 4(b). In contrast to the band structure of the single-sided pillared metamaterial, an isolated propagative negative-slope branch arises inside the complete band gap that opens between 3.25MHz and 3.76MHz. This is almost twice the width of the band gap obtained with the single-sided pillared metamaterial. This enlargement is due to the strong coupling between resonances of the double-sided pillars and Lamb modes in the matrix plate [START_REF] Assouar | Enlargement of a locally resonant sonic band gap by using doublesides stubbed phononic plates[END_REF].

The isolated branch that ranges from 3.53MHz to 3.57MHz divides the band gap into two narrower ones ranging from (in MHz X is that of a symmetric bending mode (point E in Fig. 4(c)). We show in the next subsection that the former is responsible for NEM whereas the latter, similar to the single-pillared metamaterial, leads to NMD.

A. Formation of the isolated negative-slope branch

Before analyzing the formation of the isolated branch, it is necessary to identify the resonances in the frequency range of interest. This can be done by computing the displacement field in the unit cell upon the global harmonic excitation ( ) As it is the case in the single-sided pillared metamaterial, the formation of the low frequency band gap should be the signature of NMD that results in turn from the combination of the bending and compressional resonances of the two pillars. For verifying this hypothesis, the same formalism as the one applied in Sec. II-B was implemented to calculate the effective mass density matrix against the frequency.

The normalized components 11 ρ and 33 ρ displayed in Fig. 5 the resultant force is zero. Although the resultant moment is not equal to zero, it does not couple with the motion in the plate and this bending resonance cannot induce NMD: 11 ρ and 22 ρ do not turn negative around the antisymmetric bending resonance, except at the exact resonance frequency. The reverse phenomenon occurs at the frequency of B S . Actually at this frequency the resultant moment is equal to zero, but not the resultant traction force that consequently contributes to the negative value of 11 ρ in between 3.26MHz and 3.77MHz. By symmetry, 22 ρ gets also negative in the same interval.

The same argument applies for both compressional modes. Indeed, the resultant force at the frequency of C S mode (i.e. 3.57MHz) is null and this mode does not trigger negative values of the component 33 ρ .

Conversely, the situation at the frequency of C A is similar to the one observed at the compressional resonance in the single-sided pillared metamaterial: the axial forces in the pillars combine and 33 ρ takes negative values in the range 3.04MHz to 5.29MHz. Although the isolated branch is relatively flat, especially along XM in the BZ, the group velocity is nonnull and the propagation of elastic energy at the corresponding frequencies is allowed. Therefore, the dynamic elastic modulus is necessarily negative for any vibrational mode on this branch. Figure 5(b) clearly shows that the narrow interval around 3.57MHz where propagation is allowed corresponds to the first-order symmetric compressional resonance C S which is thus causing NEM and in turn the doubly negative property. It should be noticed however that this resonance involves a purely out-of-plane motion only at points Γ, X and M of the BZ. In between these main points in the BZ, in-plane motions contribute also to the deformation of the pillars, as it can be observed in Fig. 4(c), point E.

B. Occurrence of the negative elastic modulus

Similar to the analytical study developed above to explain NMD achieved by the single-sided pillared metamaterial, one can model the symmetric compressional vibration of the double-sided pillared system as two mass-spring subsystems, each featuring a mass m and a stiffness k C , vibrating in phase as shown in To simplify the derivation of the effective elastic modulus E eff of the unit cell, we assume that the reactive forces F z caused by the vibration of the mass-spring subsystems are uniformly distributed on the matrix plate. Under this hypothesis, the dynamic stresses tensor σ and strains tensor ε in the matrix plate submitted to a harmonic excitation are given by Hooke's law and can be expressed as:

x z y m m C k C k e z A x A 1 x 1 z 2 z x F x F x z z F z F ( ) ( ) x x x y z z z 1 x x z 1 z z x 0 2 1 2 1 F A F A x a E z e E σ σ σ ε σ νσ ε σ νσ ⎧ ⎪ = ⎪ ⎪ ⎪ = ⎪ ⎪ = ⎨ ⎪ ⎪ = = - ⎪ ⎪ ⎪ = = - ⎪ ⎩ (7)
As previously, E is the elastic modulus and v the Poisson's ratio. On the other hand, the equation of motion of the mass-spring subsystems is:

( ) 2 2 z 2 1 2 C z m F k z z t ∂ = -= - - ∂ . ( 8
)
By combining Eqs. ( 7) and ( 8) one easily obtains:

1 2 x 1 2 2 0 x M 2 C 2 1 1 F x E A k a k ν ω ω ω - ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = - ⎟ - + ⎜ ⎟ ⎝ ⎠ , (9)
where

1 2 C 0 k m ω ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠
is the natural frequency of the mass-spring subsystem,

z M 2EA k e =
is the effective stiffness of the plate. Equation ( 9) is nothing else but the Hook's law applied to an effective medium which effective elastic

modulus 1 2 2 2 0 2 1 1 1 eff E E ν ω ω η ω - ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = - ⎜ ⎟ - + ⎜ ⎟ ⎝ ⎠
depends on the frequency 2 ω π and takes into account the internal motions of the unit cell through the parameters

C M k k η =
and ω 0 .

The effective elastic modulus eff E can also be written evidently as:

( )

2 2 2 2 2 0 1 1 1 1 eff E E ν ω ν ω ω η η ⎛ ⎞ ⎜ ⎟ = + ⎜ ⎟ - + - ⎜ ⎟ ⎝ ⎠ . ( 10
)
The effective elastic modulus diverges if ( )

1 2 0 2 1 1 1 ω ω ν η ⎛ ⎞ ⎜ ⎟ = ⎜ ⎟ - + ⎝ ⎠
, which thus must be regarded as the resonant frequency of the effective material. It is interesting to notice here that the resonant frequency of the effective material is not that of the mass-spring subsystem. Furthermore,

0 eff E = when ( ) 1 2 0 1 1 ω ω η = +
and eff E is negative if:

( )

1 1 2 2 0 0 2 1 1 1 1 1 ω ω ω η ν η ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ < < ⎜ ⎟ ⎜ ⎟ + - + ⎝ ⎠ ⎝ ⎠ . ( 11
)
This interval is all the more broad that the Poisson's ratio ν is large. Although this simplified model hardly allows predicting the position and width of the isolated double-negative branch, it provides an intuitive comprehension of the occurrence of NEM and in turn of the double negativity of the double-sided pillared metamaterials.

C. Evolution of the two band gaps and isolated branch vs the geometric parameters of the unit cell

Being the consequence of a resonant phenomenon, the double-negative branch of the proposed structure is relatively narrow which may be a drawback for some applications. However, the width of the frequency band where double negativity occurs can be increased through a proper choice of the geometrical parameters of the unit cell [START_REF] Jin | Pillar-type acoustic metasurface[END_REF]. We investigate in this section the influence of the height and diameter of the pillar as well as the thickness of the matrix plate, on both the width of the forbidden band gap and the negative-slope branch.

Both Figs. 7(a) and 7(b) show the effect of the dimensions of the pillars on both the low frequency band gap and the negative-slope branch. Increasing the height of the pillar leads to the decrease of the central frequency of the band gap, as well as to the decrease of the range where the double negativity occurs. This should be related to the decrease of the compressional resonance frequency for increasing height of the pillars [START_REF] Ramahi | Metamaterial particles for electromagnetic energy harvesting[END_REF]. Moreover, the lower part of the band gap broadens as the height of the pillar increases. In contrast, the propagative branch moves closer to the upper limit of the band gap as the height of the pillar increases leading to the closure of the upper part when the height of each pillar is more than about 350µm.

Remembering that the effective mass density tends towards zero while keeping negative values when the frequency approaches the upper limit of the band gap (see Fig. 5(b)), it is expected that this structure may behave as a zero-index elastic metamaterial. Actually, in that case the phase velocity in the metamaterial tends to infinity and therefore the refractive index (i.e. ratio of the velocity in the background to the velocity in the metamaterial) goes to zero. This point is further developed below.

On the other hand, the bending resonance is very sensitive to the diameter d of the pillars [START_REF] Ramahi | Metamaterial particles for electromagnetic energy harvesting[END_REF] and consequently this parameter has a large impact on the width of the band gap that broadens as the diameter increases (Fig. 7(b)). However, the width of the double-negative band is very little affected by this parameter and remains equal to about 8.5% of the width of the forbidden band whatever the diameter of the pillars is. However, the parameter that most efficiently affects the width of both the band gap and the propagative branch is the thickness of the plate (Fig. 8). Actually, both parts of the band gap slightly decrease as the thickness increases while at the same time the width Δω of the double negativity band linearly increases with increasing thickness of the plate. The relative width of the band exhibiting NEM reaches the maximum value of ~20% of the total band gap when e = 180µm. A simple model allows to summarize these results. Actually, at C S resonance there is almost no displacement at the foot of the pillars (see Fig. 5(a)) that can be modelled as a mass-spring system with a steady attachment point and for which the displacement on top of the pillar is proportional to the axial force. With these assumptions, Hook's law leads to . With the typical dimensions of the unit cell, η is always much less than unity and Eq. ( 11) can be expanded as:

( )

1 2 2 2 2 2
16

d e E h a π ω ν ρ Δ = , ( 12
)
where we have used

1 1 2 2 C 0 2 k E m h ω ρ ⎛ ⎞ ⎛ ⎞ = = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠
. Equation [START_REF] Jin | Pillar-type acoustic metasurface[END_REF] gives full account of the behaviors displayed in Fig. 7 and Fig. 8.

Although relatively narrow even when it is maximum, the width of the doubly negative band remains sufficient for an experimental validation based onto conventional techniques as for instance the excitation of Lamb waves with interdigitated transducers (IDT) on a piezoelectric film. Actually, in the frequency range that we have considered, elastic waves featuring a spectral purity δf better than the width of the NEM band [START_REF]The bandwidth at 3dB of a transducer comprising N pairs of electrodes is given by 0 0.885 f f N δ =[END_REF] can be easily obtained with an IDT comprising only a few tens pairs of electrodes [START_REF] Monkhouse | Flexible interdigital PVDF transducers for the generation of Lamb waves in structures[END_REF][START_REF] Castaings | The performance of flexible interdigital PVDF lamb wave transducers[END_REF],

allowing in turn to tailor the transition from single negativity to double negativity. It is worth noting, however, that our results are not specific to the MHz range and may be transposed either to the GHz or the kHz range by rescaling both the thickness of the plate and the lattice constant while keeping the ratio λ/a = const (λ is the wavelength at the operating frequency) and adapting the dimensions of the pillars accordingly. The high frequency domain fits with nanostructures whose typical dimensions (thickness of the plate, high and diameter of the pillars, and lattice constant) are of a few micrometers or less, with possible applications such as the vibration isolation of sensitive components, the super-focusing, the highresolution imaging, or the cloaking of elastic waves (see below). At the other end of the scale, doublesided pillared metamaterials with typical dimensions in the millimeter range are excellent candidates for applications for noise reduction and filtering in the audio range. However, a narrow band is not necessarily a drawback and, even it becomes an advantage for all applications and devices where a high selectivity in frequency is a requirement.

D. Negative refraction of Lamb waves and cloaking

The negative refraction (NR) of elastic waves through the interface between a phononic crystal and the adjacent medium, basically relates to the bands folding resulting from the periodicity, and to the negative slope of some dispersion branches [START_REF] Pierre | Negative refraction of zero order flexural Lamb waves through a two-dimensional phononic crystal[END_REF][START_REF] Sukhovich | Negative refraction and focusing of ultrasound in twodimensional phononic crystals[END_REF][START_REF] Bucay | Positive, negative, zero refraction, and beam splitting in a solid/air phononic crystal: Theoretical and experimental study[END_REF]. As a result NR in these heterostructures only occurs when the wavelength is of the same order of magnitude as the period and double-negative elastic properties are not relevant to explain this behavior. In contrast, in an elastic metamaterial the negative refraction is the direct consequence of simultaneous NMD and NEM and as such, occurs in the long wavelength limit [START_REF] Oudich | Negative effective mass density of acoustic metamaterial plate decorated with low frequency resonant pillars[END_REF]. The occurrence of NR would provide therefore further evidence of the double negativity in the double-sided pillared metamaterial. For this purpose, we have considered the FEM model shown in Fig. 9 should also notice that, as a resonant process, even the small out-of-plane component of the Lamb mode is able to excite the symmetric compressional mode C S . When NR results from the band folding at frequencies larger than the Bragg gap, it is generally accompanied with a high level of reflection at each interface between the phononic crystal and the background, which may constitutes a severe drawback in many applications. It is worth mentioning that this is not the case here. To show this, we have computed the transmission coefficient of a symmetric Lamb wave impinging at normal incidence a structure made of eleven unit cells along x-axis and infinite along the y-direction. The result displayed in Fig. 9(c) shows that the transmission is equal to unity in the whole band where double negativity occurs.

The result is totally different with the antisymmetric mode, as can be seen from Fig. 9(d). Actually, for this polarization the wave cannot propagate into the metamaterial and is totally reflected at the inlet interface. This comes from the symmetry of this mode that cannot couple neither to the symmetric compressional mode C S nor to the symmetric bending mode B S of the double-sided pillars and therefore prevents the occurrence of the necessary double negativity for the propagation to be allowed. Another current topic is the implementation of acoustic cloaking based on a zero-index metamaterial and much work has been devoted to this research since a few years [START_REF] Dong | Topology optimization of anisotropic broadband double-negative elastic metamaterials[END_REF][START_REF] Graciá-Salgado | Negative mass density and ρ -near-zero quasi-two-dimensional metamaterials: Design and applications[END_REF][START_REF] Edwards | Experimental Verification of Epsilon-Near-Zero Metamaterial Coupling and Energy Squeezing Using a Microwave Waveguide[END_REF][START_REF] Silveirinha | Tunneling of Electromagnetic Energy through Subwavelength Channels and Bends using ε -Near-Zero Materials[END_REF][START_REF] Wu | Total reflection and cloaking by zero index metamaterials loaded with rectangular dielectric defects[END_REF]. This exciting phenomenon can be ascribed to the effective phase velocity that tends towards infinity if the effective mass density is close to zero [START_REF] Graciá-Salgado | Negative mass density and ρ -near-zero quasi-two-dimensional metamaterials: Design and applications[END_REF]. Therefore, an alternative method to obtain a refractive index with a null value consists in combining an infinite effective elastic modulus to a finite effective mass density. As discussed above, the symmetric compressional resonance C S in the double-sided pillared metamaterial causes NEM with an infinite value when the excitation frequency is equal or quite close to 3.57MHz, while the effective density keeps a negative but still finite value at this frequency; zero-index behavior should therefore be observed at this frequency. To check this assumption, we have computed the displacement field in the numerical sample displayed in the top panel in Fig. 10. The metasurface is made of 21×11 unit cells and features a 7a×3a void in its center. A zero-order symmetric Lamb wave is excited at a distance of 1 mm from the left edge of the metamaterial and PML's are implemented on each side of the sample to eliminate any reflection from the boundaries. Periodic boundary conditions are applied on the other two edges. The displacement field when the frequency is tuned to 3.57MHz is shown in the mid panel of Fig. 10. At this frequency, the wave front keeps plane upon transmission through the sample, except around the void where scattering effects are observable. This is a clear confirmation that, as a result of the nearly infinite phase velocity in the metamaterial, the phase of the symmetric Lamb wave keeps a constant value during propagation, allowing for a cloaking effect in this system. In contrast, when the working frequency is tuned to 4MHz, i.e. a frequency where the elastic modulus is positive, the Lamb wave undergoes strong scattering on the void, giving rise to the distorted wave front observable in the bottom panel of Fig. 10.

This simple analysis of the transmission through the metasurface unambiguously shows that the shielding of substructures at specific frequencies may be achieved with this geometry.

IV. CONCLUSION

In this work, we have numerically investigated the negative properties of a single-phase double-sided pillared metamaterial. In contrast to the single-sided pillared structure where NEM never occurs, an isolated propagative branch with negative slope appears in the band structure of the proposed system.

Symmetric Lamb wave at a frequency in this branch can propagate within the metamaterial whereas the propagation of the antisymmetric Lamb mode is forbidden. Comparisons with the case of the single-sided pillared system well explains the formation of the low frequency band gap and the occurrence of the double negativity. Actually, in the single-sided pillared structure the negative effective mass density results from the combination of both the bending and the compressional resonances of the pillar which together allow for the formation of the band gap. At resonance, the bending motion may be broken down into a rotational and a translational motion. The latter causes the effective mass density to be negative, whereas the moment applied to the plate by the rotational motion is ineffective. Afterwards, the doublesided pillared metamaterial is investigated. The mechanism responsible for the occurrence of the effective mass density and elastic modulus simultaneously negative is discussed in details. It is shown that the negative effective mass density is achieved owing to the combination of the symmetric bending resonance and the antisymmetric compressional resonance of the pillars. It is further shown that the negative effective elastic modulus relates to their symmetric compressional resonance. In other words, the doublesided pillared structure has the advantage to allow to separate the symmetric and antisymmetric modes which is not achievable with the single-sided pillared structure. The symmetric band is created by the close vicinity of a bending resonance and a compressional resonance that together allow for the opening of a narrow pass band.

Furthermore, the influence of the geometrical parameters of the unit cell on the width of both the band gap and the double-negative band is studied. The band gap shifts downwards when increasing the height of the pillar or decreasing the diameter of the pillar, while the thickness of the plate weakly affects the position of the band gap. More importantly, both the height of the pillars and the thickness of the matrix plate have a strong impact on the width of the double-negative band. By carefully designing the unit cell, the width can be significantly enlarged. In particular, equation [START_REF] Ding | Metamaterial with simultaneously negative bulk modulus and mass density[END_REF] demonstrates that materials with large Poisson's ratio allow for larger width. Finally, two properties which both directly ensue from the double negativity of the structure, namely the negative refraction and zero-index refraction of Lamb waves in the long wavelength regime, are shown to take place when the frequency lies in the double-negative band. While the central focus of this study was the demonstration of the simultaneous negativity of both the effective mass density and elastic modulus in the pillared metamaterials, the occurrence of a refractive index with a zero value, even in the absence of a Dirac-like cone at the center of the BZ [START_REF] Zhu | Double-Zero-Index Structural Phononic Waveguides[END_REF], constitutes a significant fallout of this work.

  generally giving rise to sharp resonances. The band structure shown in Fig. 1(b) was obtained by applying periodic conditions in both xand ydirections and solving the eigenfrequencies equation with the help of a finite element method (FEM). The color scale in this figure represents the relative amplitude of the out-of-plane component

Figure 1 :

 1 Figure 1: (a) Unit cell of the single-sided pillared metamaterial and first irreducible BZ of the square lattice. (b) Band structure in the low frequency range along the high symmetry axes of the BZ; the color scale corresponds to the relative amplitude of the out-of-plane component, integrated over its volume (see text). For comparison, the band structure of the simplified coupled system described by Eq. (1) is shown as black solid lines. Inset: enlargement of the band gap around 3.1MHz. (c) Normalized total displacement field amplitude and deformation in the unit cell for modes denoted as A, B, and C in the band structure.

  where V is the volume of the unit cell and [ ] eff ρ the 3×3 dynamic effective mass density matrix to be determined. This can be done by successively applying the harmonic displacement

  it has been shown in Ref.[START_REF] Assouar | Enlargement of a locally resonant sonic band gap by using doublesides stubbed phononic plates[END_REF] that in a pillared metamaterial the offdiagonal terms in [ ] eff ρ are null except for values very close to the resonant frequencies where a singularity occurs; therefore the only elements to consider are 11 ρ edge of the band gap are displayed in Fig. 2. The component 11 ρ that mainly relates to in-plane motions (modes A and B in Fig. 1(c)), is negative from 3.21MHz to 3.50MHz, whereas 33 ρ that should be ascribed to out-of-plane displacements (mode C in Fig. 1(c)), becomes negative from 3.26MHz to 4.63MHz. A complete band gap opens when all the components of the effective mass density matrix turn negative, i.e. in between 3.26MHz and 3.50MHz in the present case.

Figure 2 :

 2 Figure 2: Components 11 ρ and 33 ρ of the effective mass density normalized to the mass density of steel, as a function of the frequency.

  the flexural rigidity of the matrix plate that depends on both the elastic modulus E and the Poisson's ratio ν ; n is the number of mass-springs per unit area; frequency of the mass-spring subsystem. Remembering that ˆz u reads derive the wave number k by solving Eq. (

.Figure 3 :

 3 Figure 3: (a) Vertical traction force (left panel) and simplified mass-spring subsystem (right panel) for the compressional resonance C. (b) Scheme of the traction force and moment (left panel), simplified mass-torsional spring subsystem (mid panel) and mass-spring subsystem modelling the bending resonances denoted as A and B in Fig. 1.

  ) 3.25 to 3.53, and from 3.57 to 3.76 respectively. A zoom in on this propagative branch is displayed in the inset in Fig. 4(b).

Figure 4 :

 4 Figure 4: (a) Unit cell of the double-sided pillared metamaterial and first irreducible BZ of the square lattice. (b) Band structure in the low frequency range. Inset: enlargement of the isolated branch. (c) Normalized total displacement field amplitude and deformation in the unit cell at the points denoted as D and E in the inset.

.

  On the whole, there are four resonances involving a motion in the plane (x, z), namely: B A (3.17MHz) is the second-order antisymmetric bending resonance; B S (3.26MHz) is the second-order symmetric bending resonance; C A (3.04MHz) is the first-order antisymmetric compressional resonance, and C S (3.57MHz) is the first-order symmetric compressional resonance. The corresponding displacement fields are displayed in Fig. 5(a). Because of the squared symmetry of the unit cell one should also add to this list two bending resonances equivalent to modes B S and B A respectively, but involving motions along y-axis.

  (b), highlight the role played by each resonance. Actually, at the frequencies of modes B A and C S the changes in 11 ρ and 33 ρ are very sharp and cannot cause NMD. On the contrary, both 11 ρ and 33 ρ turn negative at the frequencies of B S and C A which thus directly relate to NMD. The simplified traction forces and moments modeling the pillars at resonances provide a comprehensive understanding of these divergent behaviors. Both are schematically described in Fig. 5(a) by black and red arrows respectively. At B A resonance, the moments induced by the two identical pillars are the same whereas the traction forces are pointing in opposite directions and hence

Figure 5 :

 5 Figure 5: (a) Normalized total displacement field amplitude and deformation of the unit cell at four representative resonance frequencies upon excitation ( ) , 0, 0 i t Ae ω . (b) Components of the effective mass density 11 ρ (red line) and 33 ρ (black line) against the excitation frequency.

Fig. 6 (Figure 6 :

 66 Fig.6(a). In this figure, e is the thickness of the matrix plate; A x = e×a and A z = a 2 are the areas of boundaries parallel to the planes (y, z) and (x, y) respectively. To calculate the traction forces exerted on the plate when the pillars are on a vertical motion, let us consider two harmonic tension forces F x with the same amplitude, both applied along x-axis on the lateral boundaries parallel to the plane (y, z) and pointing in opposite directions (Fig.6(a)). Displacements ±x 1 of the lateral boundaries, ±z 1 of the free surfaces and ±z 2 of mass m result from the forces F x (Fig.6(b)). It is expected that the harmonic deformation induced by F x allows to excite the resonances of the mass-spring subsystems if the frequency is properly tuned.

Figure 7 :

 7 Figure 7: Edges of the band gap (lower panel) and relative width of the double-negative band (upper panel) against (a) the height and (b) the diameter of the pillars. The black (resp. red) solid lines stand for the upper (resp. lower) part of the band gap. In both panels the blue curve is the relative width of the band where the double negativity occurs.

Figure 8 :

 8 Figure 8: Edges of the band gap (lower panel) and relative width of the double-negative band (upper panel) against the thickness of the plate. The black (resp. red) solid lines stand for the upper (resp. lower) part of the band gap. The blue curve is the relative width of the band where double negativity occurs.

  figure displays a 2D prism shaped metamaterial constituted by 120 unit cells arranged in a 45° isosceles triangle within a circular steel plate. The plate was surrounded by a perfectly matched layer (PML) layer to eliminate reflections from the boundaries. The two perpendicular sides of the prism were set parallel to the lattice directions ΓX and XM. A line source 1600µm long was placed 100µm away from the inlet interface perpendicularly to ΓX direction. Lamb waves either symmetric or antisymmetric, were normally incident from the left and were propagating along ΓX direction. The frequency of the exciting waves was set to the mid value of the peak in the isolated branch, namely 3.55MHz. At this frequency, the

Figure 9 :

 9 Figure 9: (a) Finite element model used to demonstrate the negative refraction of Lamb waves; (b) Inplane velocity plot under symmetric excitation at 3.55MHz; (c) Transmission coefficient of a symmetric Lamb wave impinging the phononic crystal at normal incidence; (d) In-plane velocity for an antisymmetric excitation at 3.55MHz.

Figure 10 :

 10 Figure 10: Finite element model (top panel) used to bear out the cloaking effect and total displacement field upon symmetric excitation at 3.57MHz (middle panel) and 4MHz (bottom panel)

Acknowledgements

Wei Wang acknowledges the research scholarship No. 201608310130 provided by the China Scholarship Council.