Jianmin Liu
email: liujianmin18z@ict.ac.cn

Qi Wang
email: wangqi08@ict.ac.cn

Chentao He
email: hechentao@ict.ac.cn

Katia Jaffrès-Runser
email: katia.jaffres-runser@irit.fr

Yida Xu
email: xuyida@ict.ac.cn

Zhenyu Li
email: zyli@ict.ac.cn

Yongjun Xu

Yida Xu

Yongjun Xu

Yida Xu

B Eng

Yongjun Xu

QMR:Q-learning based Multi-objective optimization Routing protocol for Flying Ad Hoc Networks

Keywords: multi-objective routing, FANETs, Q-learning, adaptive parameters

A network with reliable and rapid communication is critical for Unmanned Aerial Vehicles (UAVs). Flying Ad Hoc Networks (FANETs) consisting of UAVs is a new paradigm of wireless communication. However, the highly dynamic topology of FANETs and limited energy of UAVs have brought great challenges to the routing design of FANETs. It is difficult for existing routing protocols for Mobile Ad Hoc Networks (MANETs) and Vehicular Ad Hoc Networks (VANETs) to adapt the high dynamics of FANETs. Moreover, few of existing routing protocols simultaneously meet the requirement of low delay and low energy consumption of FANETs. This paper proposes a novel Q-learning based Multi-objective optimization Routing protocol for FANETs to provide low-delay and low-energy service guarantees. Most of existing Q-learning based protocols use a fixed value for the Q-learning parameters. In contrast, Q-learning parameters can be adaptively adjusted in the proposed protocol to adapt to the high dynamics of FANETs. In addition, a new exploration and exploitation mechanism is also proposed to explore some undiscovered potential optimal routing path while exploiting the acquired knowledge. Instead of using past neighbor relationships, the proposed method re-estimates neighbor relationships in the routing decision process to select the more reliable next hop. Simulation results show that the proposed method can provide higher packet arrival ratio, lower delay and energy consumption than existing good performing Q-learning based routing method.

I. INTRODUCTION

Nowadays, Unmanned Aerial Vehicles (UAVs) have become one of the most important technical areas. To achieve more complex applications which are very difficult for traditional Mobile Ad Hoc Networks (MANETs) or individual UAV, Flying Ad Hoc Networks (FANETs) consisting of UAVs have been intensively studied. FANETs play an important role in the Internet-of-Everything. Significant applications of FANETs are related to the quick deployment of connectivity in situations, such as surveillance [START_REF] Luis | Multi-uav agent-based coordination for persistent surveillance with dynamic priorities[END_REF], emergency communications and rescue [START_REF] Erdelj | Wireless sensor networks and multi-uav systems for natural disaster management[END_REF]. However, FANETs are a kind of wireless networks with dynamic and unsustainable topologies due to the high mobility of nodes. Thus, routing in a network with easily disconnected feature is regarded as one of the main challenges in FANETs.

Many routing protocols have been designed for wireless ad hoc networks: proactive, reactive and hybrid routing. Proactive routing creates routes before packets are forwarded. However, maintaining the routing table information makes larger control overhead. Instead, reactive routing creates routes when packets are forwarded. But it brings larger delay due to discovery delivery paths. Hybrid routing takes a trade-off between proactive and reactive routing. It combines the advantages of low delay for proactive routing and low network control overhead for reactive routing. It is mainly suitable for networks with stable network topology.

Due to the changing network topology, geographic information based routing becomes the primary option for improving routing performance. For example, GPSR [START_REF] Karp | Gpsr:greedy perimeter stateless routing for wireless networks[END_REF] is a typical geographic information based routing protocol, which designs a packet delivery strategy based on the geographic information of the network node. However, when routing hole problem occurs frequently, the hop count increases because the next hop of the decision is random. The majority of algorithms so far have not paid significant attention to the movement pattern of nodes or numerous assumptions about nodes and network topology are considered. However, in FANETs, the high mobility of nodes leads to the highly dynamic and unsustainable topology. An adaptive and highly autonomous protocol is desired, which means that the routing protocol of FANETs should have the ability to find a reliable neighbor to complete the transmission through perceiving the change of the environment adaptively. Q-learning is an adaptive machine learning with environmental feedback as input, which contributes to adaptive routing design. In Q-learning, agents could constantly adjust their action strategies according to the reward of environmental feedback to better adapt to the dynamic and unsustainable topology.

The routing protocol based on Q-learning relies on the local data of the neighboring nodes and it did not make any assumptions about the environment. The existing Q-learning based routing protocols such as QGrid [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF], QLAR [START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF], QGeo [START_REF] Woo | Qgeo: Q-learning based geographic ad-hoc routing protocol for unmanned robotic networks[END_REF], make the best choice among the neighbors at any moment to transmit a packet to the destination. Due to the requirement of real-time data transmission and the limited energy of UAVs, it is crucial for a routing protocol to provide low delay and low energy consumption service guarantees. However, the delay and energy metrics often conflict. The optimal routing solution maximizing for instance delay may not be the one that minimized energy. A multi-objective optimization routing J o u r n a l P r e -p r o o f protocol is desired to concurrently meet the the low delay and low energy consumption requirements. However, most of the existing routing protocols focus on mono-objective optimization, and seldom pay attention to multi-objective optimization. Thus, this paper proposes a multi-objective optimization routing protocol to provide the service guarantees, which cannot be ensured in the existing Q-learning based routing protocols. In addition, most of them use fixed values for Q-learning parameters and lack of a reasonable scheme to balance exploration and exploitation. Due to the high dynamic of UAV networks, if Q-learning parameters including learning rate and discount factor are fixed, the accuracy of action selection declines, and the selected link may have the low probability of connecting to a neighbor node. The faster the network topology changes, the bigger the learning rate should be to pay more attention to new information, and the smaller discount factor should be to reflect unstable future expectations. However, the fixed learning rate and discount factor are not capable of reflecting the dynamic mobile environment, which leads to poor performance. Meanwhile, the reasonable compromise between exploration and exploitation, could make it possible to explore some undiscovered potential optimal routing path while exploiting the acquired knowledge. Due to these limitations, in this paper, we propose a novel Q-learning based Multi-objective optimization routing protocol. The main contributions of this paper are as follows:

• Joint optimization of delay and energy consumption.

Without fixed route vector table, QMR utilizes Q-learning to perform multi-objective optimization routing instead of mono-objective as in [START_REF] Woo | Qgeo: Q-learning based geographic ad-hoc routing protocol for unmanned robotic networks[END_REF]. In the multi-objective optimization, end-to-end delay and energy consumption are simultaneously minimized.

• Adaptively adjust Q-learning parameters. Due to the mobility of nodes, the link quality is extremely unstable. In this method, each link is given a different learning rate, and each node is given a different discount factor. Furthermore, the learning rate and discount factor are adaptively adjusted according to the network condition.

• Re-estimate neighbor relationships. High mobility of nodes leads to unstable neighbor relationships. Therefore, it is incorrect to use the past neighbor relationships to determine the current relationships. In our QMR, neighbor relations are re-estimated in the routing decision process to get the most reliable next hop.

• Improve the exploration and exploitation mechanism.

In FANETs, the balance process of exploration and exploitation should not be simply regulated by learning time, but regulated by the network condition. Different from traditional methods, such as ε greedy strategy, Boltzmann and Upper-Confidence-Bound(UCB), we propose a adaptive mechanism of exploration and exploitation, which balances exploration and exploitation according to the network condition.

This paper is organized as follows. The related work of routing protocol for ad hoc networks is presented in Second II. Section IV proposes a Q-learning based Multi-objective optimization Routing protocol for FANETs. The performance evaluation is then given in Section V. Finally, Section VI makes a conclusion about this paper.

II. RELATED WORK

At present, there is no a proprietary routing protocol for FANETs [START_REF] Yang | An optimization routing protocol for fanets[END_REF]. Most routing protocols of FANETs are modifications of routing protocols of MANETs. MANET routing protocols can be divided into static routing, proactive routing, reactive routing and hybrid routing.

A. Static routing protocols

In static routing protocols, static routing tables of static routing protocols must be computed and loaded before the task starts and cannot be updated during an operation. Due to this limitation, these protocols are not fault tolerant and do not apply to dynamically changing environments.

Load Carry and Delivery Routing (LCAD) [START_REF] Cheng | Maximizing throughput of uav-relaying networks with the load-carryand-deliver paradigm[END_REF] is the first routing protocol in FANET, where a UAV flies and carries data from a ground node through to the destination. The purpose of LCAD is to maximize the throughput and increase security. Although LCAD achieves higher throughput, data delivery delay is longer in LCAD due to use of a single UAV.

Multi Level Hierarchical Routing (MLHR) [START_REF] Koray | Networking models in flying ad-hoc networks (fanets): Concepts and challenges[END_REF] solves the scalability problems of large-scale vehicle networks which performance degrades as the size increases. The size and operation area can be increased by organizing the network as hierarchical structure. Analogously, UAV networks can be grouped into multiple clusters in which only cluster head is connected outside the cluster. However, frequent change of cluster heads leads to large network overhead.

B. Proactive routing protocols

In Proactive routing protocols, all nodes store the routing information from this node to other nodes by a routing table .  When the topology of the network changes, nodes need maintain and update their routing tables by exchanging routing information. However, in FANETs, the high mobility of nodes leads to the frequent change of the network topology. So, such routing protocols are not suitable to be used in FANETs due to bandwidth constraints.

The Optimized Link State Routing (OLSR) [START_REF] Clausen | Optimized link state routing protocol (olsr)[END_REF] is a wellknown proactive routing protocol where two types of messages including "hello" and "topology" are used to finish routing. The "hello" message is used to find neighbor nodes in the communication rang and maintain neighbor node list. Whereas "topology" message is use for maintaining the topology information in routing tables. However, this protocol has high control overhead due to periodically exchanging messages.

Link-quality and traffic-load aware optimized link state routing protocol (LTA-OLSR) has improved OLSR to apply to FANETs [START_REF] Pu | Link-quality and traffic-load aware routing for uav ad hoc networks[END_REF]. LTA-OLSR integrates the link quality and traffic load schemes with OLSR. The link quality scheme is designed to distinguish the link quality between the node and its neighboring nodes based on the statistical information of received packets. The traffic load scheme can ensure a light-load path by considering MAC layer channel contention J o u r n a l P r e -p r o o f information and the number of packets stored in the buffer. Compared to OLSR, LTA-OLSR can provide reliable and efficient communication in FANETs, because the link quality and the traffic load are considered.

C. Reactive routing protocols

In reactive routing protocols, a route is created only when a packet need to be transmitted from the source to the destination, and nodes do not need to maintain the routing information in real time. So, this type of routing protocols are also called on-demand routing protocols. Although reactive routing protocols reduce overhead problem of proactive routing protocols, they have large delay due to route construction.

Ad-hoc On-demand Distance Vector (AODV) [START_REF] Samir R Das | Ad hoc on-demand distance vector (aodv) routing[END_REF] is a wellknown reactive protocol in mobile ad hoc network. AODV consists of three phases: routing discovery, transmission of packet and route maintenance. In AODV, nodes (source and relay nodes) hold only one entry for each destination and store next hop information corresponding to each data communication.

Link Stability Estimation-based Preemptive Routing (LEP-R) protocol for FANETs is proposed on the basis of AOD-V [START_REF] Li | Lepr: Link stability estimation-based preemptive routing protocol for flying ad hoc networks[END_REF]. LEPR constructs multiple reliable link-disjoint paths with a new link stability metric which takes into account the past, current and future statuses of link stability. In addition, a preemptive route maintenance mechanism is proposed to repair links that may be broken soon.

D. Hybrid routing protocols

Hybrid routing takes a trade-off between proactive and reactive routing. It combines the advantages of low delay for proactive routing and low network control overhead for reactive routing. In Hybrid routing protocols, the network is divided into different zones where proactive protocol and reactive protocol are used for routing in intra zone and inter zone, respectively.

Zone Routing Protocol (ZRP) [START_REF] Haas | The zone routing protocol (zrp) for ad hoc networks[END_REF] is a hybrid routing protocol which works on the concept of zones. Each node has an alternate zone separated by a predefined range called R. Proactive routing is used to route inside the zone. When information needs to be sent outside the zone, reactive strategy is used.

Node Energy Monitoring Algorithm for Zone Head Selection (NEMA) improves ZRP by adding energy constraints to hence the lifetime of MANET [START_REF] Ghode | Nema: Node energy monitoring algorithm for zone head selection in mobile ad-hoc network using residual battery power of node[END_REF]. NEMA consists of two parts. In the first part, a zone head selection algorithm is designed to select a zone head with maximum residual power. In the second part, node energy monitoring algorithm is designed to monitor the change of the residual energy of each node and set different phase for each node according to residual energy level.

E. Q-learning based routing protocol

Complex flight environments and diverse flight tasks have caused FANETs to be in an unpredictable random fluctuation state. Therefore, the above routing protocols are difficult to adapt to the change of the network in real time, which may deteriorate network communication performance for a long time. Hence, an adaptive and highly autonomous protocol which is capable of discovering reliable communication links adaptively and autonomously is desired. Reinforcement learning (RL) is an adaptive learning method that belongs to the category of machine learning. It is a good idea to use reinforcement learning to solve routing problems of FANETs.

For the first time, Boyan and Littman used RL to solve the routing problem in static networks. Based on Q-learning, an adaptive algorithm Q-routing is proposed [START_REF] Boyan | Packet routing in dynamically changing networks: a reinforcement learning approach[END_REF]. Getting the shortest path from the source to the destination may be a good way to achieve fast routing. However, a path with the minimum number of hops cannot be the best route because it may be heavily congested. Thus, Q-routing learns a routing policy that balances minimizing the number of "hops" a packet will take with the possibility of congestion along popular routes. As a result, Q-routing is better than nonadaptive algorithm based on precomputed shortest path. The [START_REF] Oddi | Energy balancing in multi-hop wireless sensor networks: An approach based on reinforcement learning[END_REF] extends the Q-Routing and proposes an energy balancing routing algorithm, designed for Wireless Sensor Networks. The goal of the energy balancing routing algorithm is to optimize the network lifetime by balancing the routing effort among the sensors in consideration of their current remaining batteries. Today, Q-learning based routing protocols for Ad Hoc Networks are gradually emerging. In [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF], a Q-learning Based Routing Protocol (QGrid) for Vehicular Ad Hoc Networks is proposed. QGrid divides the area into different grids to make routing decision from macroscopic and microscopic aspects. The optimal next-hop grid is determined in macroscopic aspect, and then a specific vehicle is selected in the optimal next-hop grid as next-hop vehicle in microscopic aspect. Although the simulation confirmed that QGrid has better performance than the existing position based routing protocols, such as higher packet delivery ratio, end-to-end delay, energy consumption of nodes are not considered in QGrid. The [START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF] proposes a Q-Learning based adaptive routing (QLAR) for MANETs. QLAR develops a new model to detect the mobility level of each node in the network and a new metric to account for the static and dynamic routing. Although end-to-end delay of QLAR is lower than OLSR, QLAR still does not consider energy consumption of nodes.

Recently, Jung et al. proposed a Q-Learning-Based geographic ad hoc routing protocol for Unmanned Robotic Networks (QGeo). QGeo uses a distributed routing decision mechanism based on the geographical location information of nodes. In QGeo, the reward value of the action is related to the packet travel speed. Meanwhile, link condition and location estimation error are considered when calculating the travel time. Experimental results show that QGeo has higher packet delivery ratio and lower end-to-end delay compared with QGrid [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF] in mobile scenarios. However, there are still some shortcomings in QGeo. QGeo does not consider energy consumption, which means it could not simultaneously provide low delay and low energy consumption service guarantees. Limited battery lifetime is considered a major drawback of UAVs [START_REF] Abeywickrama | Empirical power consumption model for uavs[END_REF]. Although carrying an extra battery pack can power the UAV continuously for a period of time, the battery pack is J o u r n a l P r e -p r o o f still limited. Thus, it is very necessary to design a routing protocol with low energy consumption. In addition, QGeo utilizes a fixed learning rate in Q-learning approach. Learning rate is used to control the speed of updating Q-value. A fixed learning rate implies that the speed of updating Q-value is constant. However, links in FANETs are extremely unstable. Therefore, the speed of updating Q-value associated with the link should be adaptively adjusted as network environment changes.

Based on the above analysis, Table 1 compares the characteristics of exiting routing protocols for FANETs and Qlearning based routing protocols. Static, proactive, reactive and hybrid routing protocols do not belong to adaptive routing protocols, which are not able to discover reliable communication links adaptively and autonomously. Although there are some adaptive routing protocols based on Q-learning, they still have some shortcomings. On one hand, existing Q-learning based routing protocols pay little attention on multi-objective optimization to consider end-to-end delay and energy consumption simultaneously. On other hand, Q-learning parameters, learning rate and discount factor, cannot be adaptively adjusted according to the network condition.

For the above motivation, we proposed a multi-objective routing protocol where end-to-end delay and energy consumption are optimized concurrently. Furthermore, we also propose a method of adaptively adjusting Q-learning parameters and a new exploration and exploitation mechanism in Q-learning, in order to adapt to the high dynamics of FANETs.

III. SYSTEM MODEL

In this section, we present a multi-hop FANET model and Q-learning model of the FANET.

A. Multi-hop FANET Model

In this paper, we consider a UAV network comprising multiple UAVs and a ground station, as shown in Fig 1 . The ground station is regarded as the destination node to receive signals from the UAVs, one UAV is regarded as the source node to send signals, and the rest UAVs is regarded as the relay nodes to forward signals. The key issue is how to find an optimal path so that signals transmitted along the path can successfully reach the destination with low delay and low energy consumption.

B. Q-learning Model of FANETs

1) The basic of Q-learning technology: Reinforcement learning is an area of machine learning where agents alter their actions in a specific environment with the goal of maximizing results. In order to achieve the goal, agents constantly adjust their action strategies through the reward of environmental feedback. Reinforcement learning usually considers the reward by estimating value functions, state value functions or action value functions [START_REF] Sutton | Reinforcement learning: An introduction[END_REF]. Q-learning is an optimization of the offpolicy temporal difference(TD) [START_REF] Sutton | Learning to predict by the methods of temporal differences[END_REF]. The foundational idea of the Q-learning is learning from interaction with environment. It uses action value functions to get the feedback from environment. In the technique, the agent chooses an action in a particular state according to the reinforcement (Q-value). Reinforcement consists of direct reward and the future Q-value expectation. Through reinforcement, agents can assess how good an action in the current state and makes a better action at the next step. The goal of the agent is to maximize expectation for cumulative rewards over the long run. The basic iterative formula for the Q-value is as follows:

Q(s t , a t ) ← Q(s t , a t )+α[r t+1 +γmax a Q(s t+1 , a)-Q(s t , a t )]
(1) Where, α is the learning rate, γ is the discount factor, and they are set between 0 and 1. r t+1 is the direct reward value that the agent, at time t, takes action a t in state s t . After the agent taking action a t , the state of the agent changes from the s t state to the s t+1 state. max a Q(s t+1 , a), future Q-value expectation, is maximum Q-value when agent selects possible action a in the next state s t+1 . Fig. 2 shows the interaction between agent and environment in Q-learning. 2) Q-learning Model: In our QMR routing, considering a packet coming from a source and directed to a destination by multi-hop communication, the whole network is considered as an environment, and each packet in the network represents an agent. Each state of the agent indicates a node who holds the packet. For example, when a packet is at node i, the current state associated with this packet is s i . An action a i,j represents the decision of the packet (agent) to be forwarded from node i to neighbor node j. By this action, the state of the agent moves from s i to s j and the agent receives a reward about the action.

IV. Q-LEARNING BASED MULTI-OBJECTIVE

OPTIMIZATION ROUTING In this section, QMR, a Q-learning based Multi-objective optimization Routing protocol for FANETs, is introduced. In QMR, nodes utilize a reinforcement learning algorithm [START_REF] Koray | Networking models in flying ad-hoc networks (fanets): Concepts and challenges[END_REF] No No --Static routing FANETs LTA-OLSR [START_REF] Pu | Link-quality and traffic-load aware routing for uav ad hoc networks[END_REF] Yes No --Proactive routing FANETs LEPR [START_REF] Li | Lepr: Link stability estimation-based preemptive routing protocol for flying ad hoc networks[END_REF] Yes No --Reactive routing FANETs NEMA [START_REF] Ghode | Nema: Node energy monitoring algorithm for zone head selection in mobile ad-hoc network using residual battery power of node[END_REF] No Yes --Hybrid routing FANETs OPT-EQ-Routing [START_REF] Oddi | Energy balancing in multi-hop wireless sensor networks: An approach based on reinforcement learning[END_REF] No Yes Both α and γ are fixed ε greedy Q-learning based routing WSNs QLAR [START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF] Yes No Both α and γ are fixed ε greedy Q-learning based routing MANETs QGrid [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF] No No α is fixed, γ is variable Greedy and Markov Q-learning based routing VANETs QGeo [START_REF] Woo | Qgeo: Q-learning based geographic ad-hoc routing protocol for unmanned robotic networks[END_REF] Yes No α is fixed, γ is selected form two values ε greedy Q-learning based routing FANETs Type of network: (WSNs: Wireless Sensor Networks, MANETs: Mobile Ad Hoc Networks, VANETs: Vehicular Ad Hoc Networks, FANETs:Flying Ad Hoc Networks).

without knowledge of entire networks to make optimal routing decisions considering low-delay and low-energy service.

To solve the routing problem caused by the high mobility within FANETs, adaptive Q-learning parameters and a new exploration and exploitation mechanism for Q-learning are also proposed to enhance the routing performance. The QMR module consists of Routing neighbor discovery, Q-learning, Routing decision and Penalty mechanism. The flowchart of the QMR framework is shown in Fig. 3.

Start

Get location via GPS

Get location via GPS

Is the time to get the location? In QMR, nodes acquire their geographic location information by GPS. Further, the routing neighbor discovery is implemented by sending HELLO packets. When a data packet generated from a source and directed to a destination, Qlearning is the key component of routing decision. In the routing decision, if a routing hole problem [START_REF] Yu | Geographical and energy-aware routing: A recursive data dissemination protocol for wireless sensor networks[END_REF] is encountered (i.e., all the neighbors of a node are distant than the distance from this node to the destination), the penalty mechanism is triggered. Algorithm 1 shows QMR method and the details of the work at each phase of QMR.

A. Routing neighbor discovery

Each node periodically sends HELLO packets that include the node's geographic location, energy, mobility model, queuing delay, and discount factor. When a node receives HELLO packets, it uses the information in HELLO packets to establish and maintain its neighbor table. The Neighbor table contains the important information of each neighbor including the geographic location, energy, mobility model (i.e., moving speed and direction of the neighbor), arrival time, discount factor, queuing delay and Q-value. Moreover, neighbor table also stores some information about links from current node to its neighbors, such as learning rate and MAC delay. Each node utilizes the information of its neighbor table to perceive the network condition. If the information of a neighbor is not refreshed after a ExpireTime [START_REF] He | Speed: A stateless protocol for real-time communication in sensor networks[END_REF], it will be removed from the neighbor table. The time interval of the HELLO packet and the ExpireTime can be adjusted according to the moving speed of the node. The higher the moving speed of the node is, the smaller the time interval of the HELLO packet and the ExpireTime are. On the contrary, the time interval of the HELLO packet and the ExpireTime become larger.

B. Q-learning in QMR

In our QMR, Q-learning is used to make multi-objective optimization routing decisions, where end-to-end delay and energy consumption are considered in the reward function. To adapt to the dynamic mobile environment, Q-learning parameters including learning rate and discount factor are adaptively adjusted according to the delay and neighbor mobility. Furthermore, an adaptive mechanism of the exploration and exploitation in the Q-learning is also proposed to adapt to the dynamic network.

1) Energy metric: In order to balance the energy consumption of nodes, we consider the energy of the node when selecting the next hop. In QMR, we use the ratio of the residual energy of the node to the initial energy (i.e., node residual energy level [START_REF] Li | Enhancing real-time delivery in wireless sensor networks with twohop information[END_REF]) as a energy metric to measure the energy consumption of the node. The energy metric E i of node i can be expressed as follows:

E i = E res i E init i (2)
J o u r n a l P r e -p r o o f if there are neighbors whose actual velocity is greater than 0 then 13:

Select the neighbor associated with the maximum actual velocity; Give minimum reward to related neighbors.

28:

Update Q-value of related neighbors.

29: end if

Where E res i is the residual energy of the node i, and E init i is the initial energy of the node i. The larger the E i is, the lower the energy consumption of the node i is; otherwise, the higher the energy consumption is.

2) Reward function: In QMR, we jointly consider the delay and energy consumption in the reward function. The expression of the reward function f R (s t , a t ) is as follows:

fR(st, at) =    rmax
, when st+1 is destionation rmin , when st is local minimum ω * e -delay i,j + (1ω) * Ej , otherwise

where ω(0 < ω < 1) is the weight for one-hop delay. Suppose that at time t, the state of the data packet (agent) is s i , meaning that the data packet stays at node i. The action a t is to select node j for forwarding data packet. The state of the data packet will be converted from state s i to state s j , i.e., at the next time t+1, the state of the packet is s j . delay i,j is the one-hop delay from node i to node j, E j is the residual energy level of the next hop j. When the next hop j is the destination node, it means that the destination node is the neighbor of node i, so the link from node i to node j gets the maximum reward value r max . When node i is a local minimum, it means that all neighbors of node i are farther away from the destination than node i. If the packet is forwarded to any neighbor of node i, the packet will take longer time to reach the destination node. Therefore the minimum reward value r min is derived.

3) Adaptive Q-learning parameters: In Q-learning, the learning rate determines to what extent the newly acquired information overrides the old information. If the learning rate is higher, the Q-value is updated faster. Most of the existing Q-learning based routing protocols have a fixed learning rate. However, in FANETs, if the link between nodes is more unstable, the one-hop delay will be larger. Then the speed of updating Q-values should be faster. Therefore, we introduce a method of adjusting the learning rate to adapt to the speed of updating Q-value by assigning a corresponding learning rate for each link.

We use the one-hop delay to measure the stability of the link, where a relatively stable link has a smaller delay. The normalized one-hop delay ε i,j can be expressed as follows:

ε i,j = |delay i,j -µ i,j | σ i,j (4)
where delay i,j is the one-hop delay from node i to node j, µ i,j and σ i,j (σ i,j = 0) are the mean value and variance of onehop delay respectively. Since the learning rate ranges from 0 to 1, we introduce an exponential function for adaptive learning rate. The adaptive learning rate α i,j associated with the link from node i to node j can be expressed as follows:

α i,j = 1 -e -εi,j σ i,j = 0 0.3 σ i,j = 0 (5)
In the Q-learning algorithm, discount factor γ represents the stability of future Q-value expectation. A high value of discount factor indicates that future Q-value expectation is stable, while the low value indicates the vulnerable Q-value expectation. Since the routing decision is to find the reliable neighbor to forward the packet, it is important to adjust the discount factor according to the mobility of neighbors in the adjacent time intervals. The faster neighbors of the node move in adjacent time intervals, the more unstable the future Q-value expectation associated with the node is. Then the discount factor of the node is lower. On the contrary, the discount factor is higher. For the node i, the discount factor γ i is defined as follow:

γ i =1- |N i (t -1) N i (t)|-|N i (t -1) N i (t)| |N i (t -1) N i (t)| (6)
where N i (t -1) and N i (t) are the neighbor sets of the node i at time t -1 and time t respectively.

J o u r n a l P r e -p r o o f 4) Exploration and Exploitation: In order to select the most reliable neighbor to forward the packet in FANETs, we also investigate how to balance the relationship between exploration and exploitation in the Q-learning. The exploration is to search for unknown actions (get new knowledge), but too much exploration makes it difficult to retain some better actions. While the exploitation is to take advantage of explored actions which may generate high returns, but too much exploitation makes it difficult to select some undiscovered potential optimal actions.

In order to balance the relationship between exploration and exploitation, the ε greedy strategy, boltzmann mechanism [START_REF] Bloembergen | Evolutionary dynamics of multi-agent learning: A survey[END_REF] and Upper-Confidence-Bound (UCB) [START_REF] Yves | Regret Bounds and Minimax Policies under Partial Monitoring[END_REF] are usually adopted in the field of reinforcement learning. The ε greedy strategy is that the quantitative allocation method is used to explore with a small probability ε, which cannot be adaptively adjusted. Although UCB and boltzmann mechanism can change the degree of exploration, the degree of exploration is regulated according to time. This method of regulating exploration through time may be more adaptive to static environments. However, for dynamic environments such as FANETs, the balance between exploration and exploitation should be regulated by the network condition rather than simply by time.

In QMR, we propose a new mechanism of exploration and exploitation in Q-learning routing. The mechanism balances exploration and exploitation according to network condition, including actual velocity of a data packet traveling over the link, link quality and intimacy of neighbor relationship. The balance between exploration and exploitation in this method is mainly reflected in the following two aspects:

• When selecting the next hop, instead of directly selecting the link with the largest Q-value, the actual velocities of a data packet traveling over links from the current node to its neighbors are evaluated. These neighbors associated with links that meet velocity requirements are selected (more details are discussed in IV-C2). By filtering the neighbors with actual velocities, the data packet may arrive at the destination node with smaller delay. Meanwhile, neighbors that have not been selected in the past but meet the velocity requirement may also be used as the next hop, which means the agent wants to explore new actions.

• Among the neighbors that meet the velocity requirement, the κ-weighted Q-value of each link is calculated.

The neighbor with the maximum κ-weighted Q-value is selected. The weight κ that depends on link quality and intimacy of neighbor relationship (see IV-C3 for details), is used to evaluate the current condition of the link. Furthermore, Q-value is used to measure the past condition of the link. The Q-value can be regarded as the empirical value of learning link, which reflects that the agent is exploiting the learned knowledge. On the basis of Q-value, the weight is introduced, which actually reflects that the agent is balancing the relationship between exploration and exploitation.

C. Routing decision

Considering the real-time transmission, each data packet has a corresponding deadline to constrain the velocity of the data packet during the transmission. Neighbors whose actual velocities are not less than the constraint velocity will be selected as the candidate neighbors for routing decision. Furthermore, Q-learning is utilized to perform the optimal routing decision upon these candidate neighbors.

1) Delay constraint: For the real-time application, data packets generated by source have to be transmitted to the destination before a certain deadline. In QMR, when a relay node decides to forward a data packet, the latest deadline for this data packet will be updated by subtracting its passed time.

Assuming that node i selects node j as the next hop to forward data packets, node i calculates the one-hop delay delay i,j from node i and node j using the medium access delay (MAC delay) and the queuing delay recorded in the neighbor table. Since data packets travel at the speed of light in wireless media, propagation delay is approximately nanoseconds in the communication range of the order of hundred meters. Therefore, propagation delay is negligible compared to the MAC delay and queuing delay. The one-hop delay delay i,j is expressed as follows:

delay i,j = D mac i,j + D que i,j (7)
where D mac i,j is the MAC delay, which is the time needed by the medium access protocol to either successfully deliver the packet or drop it in case of repeated failures. D que i,j is the queuing delay, which is the time for the packet to reach the head of the transmission queue.

The MAC delay t m is estimated by ACK packets and is calculated by

t m = t ACK -t send ( 8
)
where t ACK is the moment when a node receives the ACK packet from its neighbor, and t send is the moment when the node sends the data packet to its neighbor. The method of window mean with exponentially weighted moving average (WMEWMA) is used to update MAC delay. For a node i with m neighbors, it always maintains m sliding windows with length n. Each window records the MAC delay of the last n data packets sent by node i to node j. The formula for the lth updated MAC delay is as follows:

D mac i,j (l) = (1 -β) l-1 k=l-n D mac i,j (k) n + βt m ( 9
)
where β(0 < β < 1) is the tunable weighting coefficient. The MAC delay t m is measured from the time of node i sends the packet to the time node i gets an acknowledgement (ACK) from the node j, which is expressed as follows: Similarly, the queuing delay is updated by

D que i,j (l) = (1 -β) l-1 k=l-n D que i,j (k) n + βt q ( 10
)
where t q is the waiting time for the data packet to reach the head of the transmission queue.

J o u r n a l P r e -p r o o f

Assume that node i sends a data packet to node j, and the deadline of the data packet at nodes i and j are respectively deadline i and deadline j . The deadline is updated by:

deadline j = deadline i -delay i,j (11)
2) Velocity Constraint: To meet delay constraint, we define the requested velocity and actual velocity. The requested velocity to transmit the data packet from node i to destination D has to meet the constraint of the deadline for end-to-end packet delivery. When node i will forward a data packet, it calculates the requested velocity according to the current deadline of the data packet. The requested velocity to transmit the data packet at node i can be expressed as follows:

V i = d i,D deadline i ( 12
)
where d i,D is the distance from node i to destination D, and deadline i is the deadline of the data packet at node i.

In QMR, node i considers the mobility of neighbor j when calculating the actual velocity of a data packet from node i to its neighbor j. Here, we assume that the node's moving speed is fixed for a certain period of time. The position of the neighbor node is predicted by the moving speed and moving direction of the neighbor node recorded in the neighbor table. Let node i add node j to the neighbor table at time t 1 . The location of node j at time t 1 is (x j (t 1 ), y j (t 1 )). The moving speed is V m j and the moving direction is angle xy j . t 2 is the current moment, that is, node i makes a routing decision at time t 2 . Assuming that node i selects node j as the next hop and the data packet reaches node j at time t 3 , the position of node j at time t 3 can be estimated as:

x j (t 3 ) = x j (t 1 ) + V m j * cos(angle j ) * (t 3 -t 1 ) (13)
ŷ j (t 3 ) = y j (t 1 ) + V m j * sin(angle j ) * (t 3t 1 ) (14)

t 3 = t 2 + delay i,j (15)
According to the current position of node i and the predicted position of neighbor node j, the actual velocity v i,j from node i to node j can be obtained.

v i,j = d i,D -d j,D delay i,j (16)
where d i,D is the distance between the real position of the node i at the time t 2 and the destination D. d j,D is the distance between the predicted position of the node j at the time t 3 and the destination D.

To meet the deadline, the actual velocity to forward the packet should be not less than the requested velocity. It means that the candidate node must satisfy v i,j no less than V i .

3) The selection of the optimal forwarding node: Due to the high mobility of nodes, the neighbor relationship and link quality between two nodes are extremely unstable. To solve this problem, we introduce a weighted Q-value as an indicator for selecting the next hop. The weight κ of the Qvalue depends on the neighbor relationship coefficient M i,j and the link quality LQ i,j . The expression of κ is as follows:

κ=M i,j *LQ i,j (17)
The link quality LQ i,j is calculated using the forward delivery ratios df i,j and reverse delivery ratios dr i,j of the link. df i,j represents the probability that a data packet successfully arrives at the recipient (i.e., node j); dr i,j represents the probability of sender (i.e., node i) successfully receiving ACK packets. In [START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF], df i,j and dr i,j are measured using hello message. In this paper, we also use hello message to measure df i,j and dr i,j . The expression of LQ i,j is as follows:

LQ i,j = df i,j * dr i,j (18)
The neighbor relationship coefficient M i,j indicates the intimacy between node i and node j. Due to the high dynamics of nodes, the neighbor relationships may vary constantly. Therefore, it is not a good way to determine the current neighbor relationships by past neighbors. We need to reestimate the neighbor relationships between nodes. Let node i add node j to its neighbor table at time t 1 , and node i makes routing decision at time t 2 . Here assuming that node i selects node j as the next hop, the data packet arrives at node j at time t 3 , then we can use ( 13) and ( 14) to estimate the position of node j at time t 3 . According to the estimated position of node j at time t 3 and the actual position of node i at time t 2 , the intimacy of node i and node j within (t 3t 2 ) time can be estimated. The coefficient of neighbor relationship M i,j as:

M i,j = 1 - di,j R , d i,j R 0, d i,j > R (19)
d i,j = ( x j (t 3 ) -x i (t 2 )) 2 + ( ŷ j (t 3 ) -y i (t 2 )) 2 (20)
where d i,j is the distance between the position of node i at time t 2 and the estimated position of node j at time t 3 . R is the propagation range of the node. If d i,j is greater than R, the data packet cannot reach node j from node i, indicating that node i has a weak neighbor relationship with node j. So, M i,j is equal to 0. On the contrary, if d i,j is smaller than R, M i,j becomes lager as the decrease of d i,j , which indicates that node i has strong intimacy with neighbor node j.

When the node i makes a routing decision, a candidate neighbor associated with the largest κ-weighted Q-value is selected as the next hop. It can be expressed as following:

max κ * Q(s i , a i,j ) s.t. v i,j ≥ V j (21)
where s i is the state of the data packet (agent) at current time.

One possible action a i,j is to select a neighbor node j as the next hop.

If the set of candidate neighbors is empty, while there are neighbors whose actual velocities are greater than 0, the neighbor associated with the maximum actual velocity will be selected as the next hop. Moreover, if all neighbors whose actual velocities are no greater than 0 (i.e., the routing hole problem is encountered), the penalty mechanism is triggered.

D. Penalty mechanism

The emergence of routing holes will increase the delay of the data packet. In order to reduce the routing hole problem, penalty mechanism is proposed in QMR. Penalty mechanism J o u r n a P r -p o o is triggered in the following two cases. One case is that when node j makes routing decision, if the routing hole problem is encountered, it gives previous hop node i a feedback. After receiving the feedback, the previous node i will give a minimum reward r min to the link from node i to node j and use Q-learning to update the Q-value of the link again. Another case is that forwarding node i does not receive the ACK packet from the next hop node j, node i will give the minimum reward r min of the link from node i to node j, and update the corresponding Q-value.

E. Illustration for QMR routing decision

Fig. 4 shows a simple network topology diagram. In the network, there are source node, destination node and 5 relay nodes. Suppose that at the current moment t, there are three data packets in the network, i.e., three agents. Packet 1, packet 2 and packet 3 (i.e., agent 1, agent 2, agent 3) stay on node 1, node 4 and node S respectively. Therefore, at time t, the state of packet 1, packet 2, and packet 3 are s 1 , s 4 , and s s , respectively, as shown in Fig. 4(a). Suppose that at time t, the distances between node S, node 1, node 4 and destination node D are 500m, 360m and 245m respectively. The deadline of packet 1 at node 1 is 45 ms; the deadline of packet 2 at node 4 is 35 ms; the deadline of packet 3 at node S is 50 ms. Using [START_REF] Samir R Das | Ad hoc on-demand distance vector (aodv) routing[END_REF], according to the distance between the node and the destination and the deadline of the packet on the corresponding node, the requested velocities of the node S, node 1 and node 4 are 10km/s, 8km/s and 7km/s, respectively. As can be seen from the Fig. 4(a), at the current moment t, the sets of neighbors of node S, node 1 and node 4 are {node 1, node 2}, {node S, node 3, node 4}, {node 1, node 2, node 3, node 5} respectively. The node S, node 1 and node 4 need to select one of their neighbors as the next hop to forward the packet. It is assumed that at time t, the actual velocities for one-hop forwarding are shown in Fig. 4(b). According to the requested velocities of source node S, node 1 and node 4 and the actual velocities of their neighbors, we can obtain the set of candidate neighbors of each node. The sets of candidate neighbors of node S, node 1 and node 4 are {node 1, node 2}, {node 3, node 4}, {node 3, node 5} respectively, according to the velocity requirement that the actual velocity is not less than the requested velocity. Suppose that at time t, the Q-value, weight κ and weighted Q-value of the link corresponding to the candidate neighbors of node S, node 1 and node 4 are shown in Table II. In table II, (i, j) represents the link from node i to node j. At the beginning of the learning phase, the initial Q-value of each link is 0.5. As the weighted Q-value of the link (1,3) is greater than the weighted Q-value of the link [START_REF] Luis | Multi-uav agent-based coordination for persistent surveillance with dynamic priorities[END_REF][START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF], packet 1 (agent 1) will select the link (1,3) to forward. In fact, the Q-value of the link (1,3) is also greater than the Q-value of the link [START_REF] Luis | Multi-uav agent-based coordination for persistent surveillance with dynamic priorities[END_REF][START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF], which implies link (1,3) is better than link [START_REF] Luis | Multi-uav agent-based coordination for persistent surveillance with dynamic priorities[END_REF][START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF] in the past. Packet 1 (agent 1) selects the best link in the past, which means it mainly exploits the knowledge (experience) that has been learned in the past.

For packet 2 (agent 2), although the Q-value of the link (4,5) is smaller than the Q-value of the link (4,3), the weighted Qvalue of the link [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF][START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF] is greater than the weighted Q-value of the link [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF][START_REF] Karp | Gpsr:greedy perimeter stateless routing for wireless networks[END_REF]. Therefore, packet 2 (agent 2) will select the link [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF][START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF]. Since the Q-value of the link [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF][START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF] is not equal to the initial value 0.5, the link (4,5) has been explored before time t. However, packet 2 (agent 2) still selects the link (4,5) to forward, which indicates that links that have been explored in the past but in poor condition will have the opportunity to be re-explored.

For packet 3 (agent 3), although the Q-value of the link (S,2) is smaller than the Q-value of the link (S,1), the weighted Qvalue of the link (S,2) is greater than the weighted Q-value of the link (S,1). Therefore, packet 3 (agent 3) will select the link (S,2). The link (S,2) has not been explored before time t, because the Q-value of the link is equal to the initial value of 0.5. However, packet 3 (agent 3) still selects link (S,2) to forward, which implies that packet 3 (agent 3) decides to explore a new link (i.e., explore new action) due to the link (S,2) with bigger weight (i.e., the link in better condition at present). This means that previously undiscovered links might be explored.

In a word, packet 1 selects the link that was regarded as the best in the past (exploitation), while packet 2 selects the link that was regarded as the poor in the past but may be better at present (exploration). Packet 3 selects the link not previously explored (exploration). Then at the next time t + 1, packet 1, packet 2 and packet 3 respectively stay on node 3, node 5 and node 2. Therefore, at the time t + 1, the states of packet 1, packet 2 and packet 3 (i.e., agent 1, agent 2 and agent 3) are s 3 , s 5 and s 2 respectively, as shown in Fig. 4(b). Meanwhile, these actions (i.e., link [START_REF] Luis | Multi-uav agent-based coordination for persistent surveillance with dynamic priorities[END_REF][START_REF] Karp | Gpsr:greedy perimeter stateless routing for wireless networks[END_REF], link [START_REF] Li | Qgrid: Q-learning based routing protocol for vehicular ad hoc networks[END_REF][START_REF] Serhani | Qlar: A qlearning based adaptive routing for manets[END_REF] and link (S,2)) will receive corresponding rewards.

From this example, we can see that both explored links and unexplored links are likely to be explored when forwarding J o u r n a l P r e -p r o f packets. Due to the frequent change of topology in FANETs, the quality of link is extremely unstable, so the exploration and exploitation mechanism mentioned in this paper can better adapt to the situation of link instability to find a better forwarding path.

V. PERFORMANCE EVALUATION

In this section, our QMR algorithm is implemented and compared with the existing good performing QGeo [START_REF] Woo | Qgeo: Q-learning based geographic ad-hoc routing protocol for unmanned robotic networks[END_REF], using an event-driven wireless networks simulator WSNet 1 .

In QGeo, the packet travel speed is considered as the reward to select optimal next hop. The discount factor selected from two constant values and a fixed learning rate are used to update Q-value in QGeo. The packet travel speed is only considered in reward function, However, the energy consumption is ignored, which makes QGeo difficult to provide low energy consumption service. In contrast, discount factor can be adjusted adaptively according to the mobility of neighbors in the adjacent time intervals, and learning rate can be adjusted adaptively according to one-hop delay in our QMR. In addition, QMR simultaneously considers delay and energy consumption to provide low delay and low energy consumption service.

For the considered scenario, 25 nodes are evenly distributed in an area of 500m x 500m, and the coordinates of destination node are (500, 500). We randomly select one node as the source node to transmit data to the destination node, and the rest of nodes except the destination node are relay nodes. The source emits a periodic flow of data packets whose data interval is set differently for comparison. Initially, Q-value of each link is 0.5. The parameters for the scenario are shown in Table III. The performance metrics including average end-to-end delay, maximum end-to-end delay, packet arrival ratio and energy consumption are considered:

• Average end-to-end delay: The average delay of the data packet from the source node to the destination node.

• Maximum end-to-end delay: The maximum delay of the data packet from the source node to the destination node.

• Packet arrival ratio: The ratio of the number of data packets received by the destination node (excluding redundant 1 http://wsnet.gforge.inria.fr/ data packets) to the number of data packets transmitted by the source node.

• Energy consumption: We consider as a first approximation that the main energy consumption factor is due to the emission and reception of a data packet. Thus, the energy consumption is defined as the average number of emissions and reception operations performed by all nodes (sources and relays). This value is normalized by the number of data packets sent by the sources. The complete definition is given in [START_REF] Qi | Tdma versus csma/ca for wireless multihop communications: a stochastic worst-case delay analysis[END_REF] In the following simulations, all mobile nodes (including source and relay nodes) utilize the Random Waypoint Mobility Model. In this mobility model, a mobile node moves from its current location to a new location by choosing a direction and speed [START_REF] Camp | A survey of mobility models for ad hoc network research[END_REF]. The new location is randomly chosen in the range of simulation area, while the new speed is evenly chosen from [minspeed, maxspeed]. When the mobile node moves to the newly selected destination at the selected speed, the mobile node pauses for a specified period and then starts to move to another new location. In all simulations, the pause time of mobile nodes set to 0.

A. Evalutation metrics for different data interval

In this simulation, the minspeed and the maxspeed of mobile nodes are 0m/s and 15m/s respectively. Source node sends one thousand data packets at different interval. For each interval, we repeated the simulation 100 times. From Fig. 5 to Fig. 8, the performance of our QMR for different data interval of source is compared with QGeo, considering average end-toend delay, maximum end-to-end delay, packet arrival ratio and energy consumption. From Fig. 5 and Fig. 6, we can see that the average end-toend delay and max end-to-end delay of our algorithm QMR are lower than QGeo. Compared to the QGeo, the average end-to-end delay and max end-to-end delay are averagely reduced by 42% and 47%, respectively. The main reason is that our algorithm constrains the velocity of data packets in the transmission process. QMR requires that the actual velocity of the data packets in the transmission process is not less than the requested velocity within the delay deadline. Thus, the routing path with low delay from source to destination is selected for data transmission. In addition, compared with QGeo, our algorithm not only considers the MAC delay, but Fig. 7: Energy consumption for different data interval also considers the average queuing delay in one-hop delay, which is more realistic end-to-end delay. Fig. 7 shows that energy consumption of our QMR is lower. Compared with QGeo, energy consumption is reduced by 18%. One of the major reasons is that our algorithm considers the energy consumption of the node in the Q-learning reward function. By considering the initial energy and the residual energy of the nodes in a comprehensive manner, the nodes with smaller energy consumption are selected as the next hop, thereby balancing the energy consumption of all nodes in the network. Besides, QMR has good capability to find a low delay path, which results in fewer retransmissions and higher energy utilization efficiency.

Figure 8 represents packet arrival ratio of our QMR. Compared with QGeo, the packet arrival rate of our QMR has increased by 6%. The main reason is that considering the high mobility of nodes, our algorithm re-estimates neighbor In addition, learning rate and discount factor in Q-learning are adaptively adjusted according to the one-hop delay and the mobility of neighbors to predict the status of the link. According to the prediction results, the most stable link is selected for data transmission, which in turn increases the packet arrival ratio. To be concluded, our QMR performs better than the existing good performing QGeo for different data interval.

B. Evalutation metrics for different moving speed

In this simulation, source node sends one thousand data packets at the interval of 30ms. The minspeed of mobile nodes is 0m/s, and the maxspeed of mobile nodes varies from 10m/s to 100m/s. For each maximum speed of mobile nodes, we repeated the simulation 100 times. From Fig. 9 to Fig. 12, the performance of our QMR for different moving speed is compared with QGeo. Figure 9 shows that packet arrival ratio is higher than Qgeo, which is increased by 10% on average. We can see in figure 10 that energy consumption of the node is lower. Compared with Q-Geo, energy consumption is reduced by 19% on average. From figure 9 and figure 12, we can see that in our algorithm, the average end-to-end delay and max end-to-end delay are lower, which is averagely reduced by 48% and 44%, respectively. We can include that our QMR performs better than the existing good performing QGeo for different moving speed.

C. Evalutation metrics for a scenario with faulty relay nodes

In this simulation, we randomly select 10 nodes among 25 nodes to stop working by power off at 1 second after the start of the simulation, to achieve the scenario where numbers of relay nodes are faulty. The performance of the proposed method QMR the existing QGeo under different data interval of source are compared. The minimum speed of mobile nodes is 0m/s, and the maximum speed of mobile nodes is 15m/s. We repeated the simulation 100 times. From Fig. 13 to Fig. 16, we can see that our method still have better performance than the existing method QGeo, even if numbers of intermediate nodes are faulty. From Fig. 13 and Fig. 14, it can be seen that our method QMR still has lower average and max end-to-end delay. Compared to QGeo, the average end-to-end delay and max end-to-end delay are averagely reduced by 40% and 46%, respectively. Meanwhile, energy consumption is reduced by 17% on average, and packet arrival ratio is increased by 6% in comparison to QGeo.

D. Comparisons of exploration and exploitation mechanism

In this simulation, we use QMR to compare the exploration and exploitation mechanism proposed in this paper with traditional methods (ε greedy strategy, UCB, boltzmann mechanism) by average end-to-end delay and packet arrival ratio. The minspeed and the maxspeed of mobile nodes are 0m/s and 15m/s respectively. The data interval of the source to transmit one packet is 10ms and 40ms. The experimental results are shown in the Fig. 17.

From Fig. 17, we can see that our mechanism can enhance the routing performance, where the average end-to-end delay is smallest and the packet arrival ratio is the highest. This is because the regulation method is adaptive to the variation of network condition. Therefore, it can better balance the relationship between exploration and exploitation, which is more suitable for FANETs and other dynamic environments. The drastically changing topology of FANETs has brought great challenges to routing protocols for FANETs. Existing routing protocols for MANETs and VANETs cannot be directly applied in FANETs. In this paper, we propose a novel Q-learning based routing protocol for FANETs. Considering the requirement of low delay and low energy consumption of FANETs, a Q-learning based Multi-objective optimization Routing protocol (QMR) is proposed aiming directly at simultaneously providing low-delay and low-energy service guarantees. In addition, a method of adaptively adjusting Qlearning parameters is proposed to adapt to the high dynamics of FANETs. In this method, learning rate is adaptively adjusted according to one-hop delay, and discount factor is adjusted according to the mobility of neighbors in the adjacent time intervals. The results have demonstrated outstanding performance of our QMR in comparison with the QGeo. The transmission multi-flows with different requirements of QoS will be investigated in the future. Moreover, the implementation of QMR in a physical UAVs would be a challenge and could provide many useful insights.

Fig. 1 :

 1 Fig. 1: A Multi-hop FANET with multi-UAVs

Fig. 2 :

 2 Fig. 2: The interaction between agent and environment in Q-learning

  J o u r n a l P r e -p r o o f

Fig. 3 :

 3 Fig. 3: Flowchart of QMR framework

1 Fig. 4 :

 14 Fig. 4: The network status at time t and t + 1

Fig. 5 :

 5 Fig. 5: Average end-to-end delay for different data interval

Fig. 6 :

 6 Fig. 6: Max end-to-end delay for different data interval

Fig. 8 :

 8 Fig. 8: Packet arrival ratio for different data interval

Fig. 9 :Fig. 10 :Fig. 11 :Fig. 12 :

 9101112 Fig. 9: Packet arrival ratio for different moving speed

Fig. 13 :

 13 Fig. 13: Average end-to-end delay for a scenario with faulty relay nodes

Fig. 14 :Fig. 15 :

 1415 Fig. 14: Max end-to-end delay for a scenario with faulty relay nodes

Fig. 16 :Fig. 17 :

 1617 Fig. 16: Energy consumption for a scenario with faulty relay nodes

TABLE I :

 I Comparison between routing protocols

	Routing protocol	Consider end-to-end delay	Consider energy consumption	Q-learning parameters	Exploration and utilization policy	Type of routing	Type of network
	LCAD [8]	No	No	-	-	Static routing	FANETs
	MLHR

Update neighbor table Update neighbor table Forward data packet? Forward data packet? Select the neighbor with the largest weighted Q-value Select the neighbor with the largest weighted Q-value No candidate neighbor? No candidate neighbor? Update Q-value of the neighbor Update Q-value of the neighbor No neighbor with actual speed than 0 ? No neighbor with actual speed than 0 ? Give minimum reward to related neighbors Give minimum reward to related neighbors Select the neighbor with the maximum actual velocity Select the neighbor with the maximum actual velocity End

	Is the time to get
	the location?
	Send HELLO Send HELLO
	packet ? packet ?
	Send HELLO packet Send HELLO packet

TABLE II

 II

	: Link information
	link weight κ Q-value weight Q-value 0.372 0.39 0.48 0.468 0.504 0.578 (s,1) (s,2) (1,3) (1,4) (4,3) (4,5) 0.6 0.78 0.6 0.65 0.63 0.85 0.62 0.5 0.8 0.72 0.8 0.68

TABLE III :

 III Parameter configuration

	Parameters	Settings
	Area Size	500m x 500m
	Number of Nodes	25
	MAC	802.11 DCF
	Radio propagation	propagation range, rang = 180m
	Interferences	interferences orthogonal
	Modulation	modulation bpsk
	Antenna	antenna omnidirectionnal
	Battery	energy linear
	HELLO Interval	100ms
	ExpireTime	300ms
	update interval of γ	100ms
	size of data packet	127Bytes
	β	0.5
	ω	0.6

VII. ACKNOWLEDGMENTS

This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grant No. 61602447.

She received a three-year Marie-Curie OIF fellowship from the European Union (2007-2010). Her research interest includes the performance evaluation of wireless networks in general, with a special focus on real-time guaranties provision.

J o u r n a l P r e -p r o o f

Conflict of interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

J o u r n a l P r e -p r o o f


Olivier Bernard
email: olivier.bernard@inria.fr

Liudi Lu
email: liudi.lu@inria.fr

Jacques Sainte
email: jacques.sainte-marie@inria.fr

Julien Salomon
email: julien.salomon@inria.fr

Controlling the bottom topography of a microalgal pond to optimize productivity

We present a coupled model describing growth of microalgae in a raceway cultivation process, accounting for hydrodynamics. Our approach combines a biological model (based on the Han model) and shallow water dynamics equations that model the fluid into the raceway. We then describe an optimization procedure dealing with the topography to maximize the biomass production over one cycle (one lap of the raceway). The results show that non-flat topographies enhance microalgal productivity.

I. INTRODUCTION

Research on production of biotechnological microalgae has been booming in recent years. The potential of these techniques finds interests for the cosmetics, pharmaceutical, food and -in the longer term -green chemistry and energy applications [START_REF] René | An outlook on microalgal biofuels[END_REF]. Production is carried out in biophotoreactors that often take the form of a raceway, i.e. a circular basin exposed to solar radiation where the water is set in motion by a paddle wheel [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF]. On top of homogenizing the medium for ensuring an equidistribution of the nutrients necessary for algal growth, the main interest of mixing is to guarantee that each cell will have regularly access to light [START_REF] Demory | How do microalgae perceive light in a high-rate pond? towards more realistic lagrangian experiments[END_REF]. The algae are regularly harvested, and their concentration is maintained around an optimal value [START_REF] Muñoz-Tamayo | Optimizing microalgal production in raceway systems[END_REF], [12]. The algal concentration stays generally below 1% [START_REF] Bernard | Modelling of microalgae culture systems with applications to control and optimization[END_REF]. Above this value, the light extinction is so high that a large fraction of the population is in the dark and does not grow anymore.

Many phenomena have to be taken into account to represent the entire photoproduction process. In this paper, we develop a coupled model to describe the growth of algae in a raceway, accounting for the light that they receive. More precisely, prolonging the study of [START_REF] Bernard | A 2d model for hydrodynamics and biology coupling applied to algae growth simulations[END_REF] with a simpler model, we combine a photosynthesis model, the Han model, and a hydrodynamic model based on the Saint-Venant equations for shallow water flows. This approach enables us to formulate an optimization problem where the raceway topography is designed to maximize the productivity. For this problem, we present an adjoint-based optimization scheme including constraints to respect the shallow water regime. Our numerical tests show that non-trivial topographies can be obtained.

The outline of the paper is as follows: in Section 2, we present the biological and hydrodynamic models underlying our coupled model. In Section 3, we describe the optimization problem and a corresponding numerical optimization procedure. Section 4 is devoted to the numerical results obtained with our approach. We conclude in Section 5 with some perspectives opened by this work.

II. COUPLING HYDRODYNAMICAL AND BIOLOGICAL MODELS

Our approach is based on a coupling between the hydrodynamic behavior of the particles and the evolution of the photosystems driven by the light intensity they received when traveling across the raceway pond.

A. Modelling the photosystems dynamics

The photosystems we consider are cell units that harvest photons and transfer their energy to metabolism. The photosystems dynamics can be described by the Han model [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF]. In this compartmental model, the photosystems can be described by three different states: open and ready to harvest a photon (A), closed while processing the absorbed photon energy (B), or inhibited if several photons have been absorbed simultaneously (C). Their evolution satisfy the following dynamical system

C A B

Photo Photo

   Ȧ = -σIA + B τ , Ḃ = σIA -B τ + k r C -k d σIB, Ċ = -k r C + k d σIB.
Here A, B and C are the relative frequencies of the three possible states with

A + B + C = 1, (1)
and I is the photon flux density, a continuous time-varying signal. Besides, σ stands for the specific photon absorption, τ is the turnover rate, k r represents the photosystems repair rate and k d is the damage rate. Following [START_REF] Lamare | Gradient-based optimization of a rotating algal biofilm process[END_REF] and using (1), we reduce this system to one single evolution equation:

Ċ = -α(I)C + β(I), (2)
where

α(I) = k d τ (σI) 2 τ σI + 1 + k r , β(I) = α(I) -k r .
The net specific growth rate is obtained by balancing photosynthesis and respiration, which gives

µ(C, I) = -γ(I)C + ζ(I), (3)
where

γ(I) = kσI τ σI + 1 , ζ(I) = γ(I) -R.
Here k is a factor that links received energy and growth rate.

The term R represents the respiration rate.

Remark 1: The dynamics of the biomass X is derived from µ:

Ẋ = µ(C, I)X. (4)

B. Steady 1D Shallow Water Equation

The shallow water equations are one of the most popular model for describing geophysical flows, which is derived from the free surface incompressible Navier-Stokes equations (see for instance [START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water; numerical validation[END_REF]).

We assume that the hydrodynamics in the raceway pond corresponds to a laminar flow where the viscosity is neglected. As a consequence, it can be described by the steady Shallow Water equation which is given by

∂ x (hu) = 0, (5)
∂ x (hu 2 + g h 2 2 ) = -gh∂ x z b , (6)
where u is the horizontal averaged velocity of the water, h is the water elevation, the constant g stands for the gravitational acceleration, and z b defines the topography of the raceway pond. The free surface η is given by η = h + z b . This system is presented in direction and the x axis represents the horizontal direction. Besides, I s represents the light intensity at the free surface (assumed to be constant).

Considering Equations ( 5)-( 6), and integrating the former, we get

hu = Q 0 , (7)
for a fixed positive constant Q 0 . On the other hand, since hu is a constant, Equation ( 6) can be rewritten as

hu∂ x u + h∂ x gh + h∂ x gz b = 0.
Assuming that h > 0 and using [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF], we can eliminate u to get

∂ x ( Q 2 0 2h 2 + g(h + z b )) = 0.
Then for two given constants

h(0), z b (0) ∈ R, we have for all x ∈ [0, L] Q 2 0 2h 2 + g(h + z b ) = Q 2 0 2h 2 (0) + g(h(0) + z b (0)) =: M 0 . ( 8
)
Finally we find that the topography z b satisfies

z b = M 0 g - Q 2 0 2gh 2 -h. (9)
Remark 2: This model is relevant if the flow is subcritical and the water elevation is positive, meaning that the following relation shall hold

0 < u(x) gh(x) < 1 x ∈ (0, L),
which can also be written as

h(x) > ( Q 2 0 g ) 1 3
x ∈ (0, L). [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF] This leads us to introduce the threshold value h c := (

Q 2 0 g ) 1 3
, which guarantees that the system is in a shallow water regime. Note that h c is the threshold value of h for which the Froude number equals to one.

From [START_REF] Martínez | Theory of turbid microalgae cultures[END_REF], h is the solution of a third degree equation. Given a smooth topography z b , there exists a unique positive smooth solution of ( 9) which satisfies the subcritical flow condition (see [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]Lemma 1]).

C. Lagrangian trajectories of the algae and captured light intensity

Let z be the vertical position of a particle in the raceway. We first determine the Lagrangian trajectory of an algal cell which starts at a given position z(0) at time 0.

From the incompressibility of the flow, we have

∇ • u = 0,
with u = (u(x), w(x, z)). Here u(x) is the horizontal velocity and w(x, z) is the vertical velocity. This implies that

∂ x u + ∂ z w = 0. (11)
Integrating [START_REF] Muñoz-Tamayo | Optimizing microalgal production in raceway systems[END_REF] from z b to z gives:

0 = z z b ∂ x u(x) + ∂ z w(x, z) dz, = ∂ x z z b u(x)dz + z z b ∂ z w(x, z)dz, = ∂ x (z -z b )u(x) + w(x, z) -w(x, z b ), = (z -z b )∂ x u(x) -u(x)∂ x z b + w(x, z).
where we have used the non-penetration condition w(x, z b ) = 0. It then follows from (9) that

w(x, z) = ( M 0 g - 3u 2 (x) 2g -z)u (x).
The Lagrangian trajectory is consequently characterized by the system

ẋ(t) ż(t) = u(x(t)) w(x(t), z(t)). ( 12
)
To obtain the light intensity observed on this trajectory, we assume the turbidity to be constant over the considered time scale and use Beer-Lambert law:

I(x, z) = I s exp (-ε(η(x) -z)), ( 13
)
where ε is the light extinction coefficient. In doing so, we suppose that the system is perfectly mixed so that the concentration of the biomass X defined in ( 4) is homogeneous and ε is constant. Remark 3: We see that in this approach, the light intensity I couples the hydrodynamic model and Han model: the trajectories of the algae define the received light intensity, which is used in the photosystem dynamics.

III. OPTIMIZATION PROBLEM

In this section, we define the optimization problem associated with the biological-hydrodynamical model. We assume a constant volume for the raceway system. We first introduce our procedure in the case of one single layer, and then extend this to a multiple layers system. For simplicity, we omit t in the notations.

A. One layer problem 1) Optimization functional: The average net specific growth rate for one trajectory is defined by

μ1 := 1 T T 0 µ(C, I(x, z))dt. ( 14
)
As we have mentioned in the previous section, in the subcritical case, a given topography z b corresponds to a unique water height h. On the other hand, the volume of our system is defined by

V = L 0 h(x)dx. (15)
Therefore, we choose to parameterize h by a vector a ∈ R N , which will be the variable to be optimized, in order to handle the volume of our system. Given a vector a and the associated h, the optimal topography can be obtained by using [START_REF] Martínez | Theory of turbid microalgae cultures[END_REF].

Our goal here is to optimize the topography to maximize μ1 . In this way, we define the functional with the help of ( 3

The optimal control problem then reads: Find a solving the maximization problem: where p 1 , p 2 and p 3 are the Lagrange multipliers associated with the constraints (16).

max a∈R N μ1 (a). (17
The optimality system is obtained by cancelling all the partial derivatives of L. Differentiating L with respect to p 1 , p 2 , p 3 and equating the resulting terms to zero gives the corrected model equations (16). Integrating the terms p 1 Ċdt, p 2 żdt and p 3 ẋdt on the interval [0, T ] by parts and differentiating L with respect to C, z, x, C(T ), z(T ), x(T ) gives rise to

                                           ∂ C L = -1 T γ(I(x, z; a)) + ṗ1 -α(I(x, z; a))p 1 ∂ z L = 1 T -γ (I(x, z; a))C + ζ (I(x, z; a)) + p 1 -α (I(x, z; a))C + β (I(x, z; a)) ∂ z I(x, z; a) + ṗ2 + p 2 ∂ z w(x, z; a) ∂ x L = 1 T -γ (I(x, z; a))C + ζ (I(x, z; a)) + p 1 -α (I(x, z; a))C + β (I(x, z; a)) ∂ x I(x, z; a) + p 2 ∂ x w(x, z; a) + ṗ3 + p 3 ∂ x u(x; a) ∂ C(T ) L = p 1 (T ) ∂ z(T ) L = p 2 (T ) ∂ x(T ) L = p 3 (T ).
(18) Given a vector a, let us still denote by C, x, z, p 1 , p 2 , p 3 the corresponding solutions of ( 16) and (18). The gradient ∇μ 1 (a) is obtained by

∇μ 1 (a) = ∂ a L,

B. Multiple Layers Problem

We now extend the previous procedure to deal with multiple layers. Let us denote N z the number of layers and C i (resp. z i ) the photo-inhibition state (resp. the trajectory position) associated with the i-th layer. As the optimization functional, we consider the semi-discrete average net specific growth rate over the domain, namely:

μNz (a) = 1 N z Nz i=1 1 T T 0 µ(C i , I(x, z i ; a))dt, (20)
where C i , x, z i verify the constraints (16

) for i = 1, • • • , N z .
Remark 4: Note that the average net specific growth rate over the domain is defined by:

μ∞ := 1 L L 0 1 h(x) η(x) z b (x) µ(C, I(x, z))dzdx.
Our approach consequently consists in considering a vertical discretization of μ∞ , which gives (20). This discretization should not give rise to any problem and is left to a future contribution.

Similar computations as in the previous section give rise to N z systems similar to (18), where C and z are replaced by C i and z i respectively. Denoting by p 1,i , p 2,i the associated Lagrange multipliers, the partial derivatives ∂ x L is a little different from the previous section, more precisely, this can be computed by:

∂ x L = 1 N z Nz i=1 1 T -γ (I(x, z i ; a))C i + ζ (I(x, z i ; a)) ∂ x I(x, z i ; a) + Nz i=1 p 1,i -α (I(x, z i ; a))C i + β (I(x, z i ; a)) ∂ x I(x, z i ; a) + Nz i=1 p 2,i ∂ x w(x, z i ; a) + ṗ3 + p 3 ∂ x u(x; a).
Finally, the gradient ∇μ Nz (a) is given by

∂ a L = 1 N z Nz i=1 1 T T 0 -γ (I(x, z i ; a))C i + ζ (I(x, z i ; a)) ∂ a I(x, z i ; a)dt + Nz i=1 T 0 p 1,i -α (I(x, z i ; a))C i + β (I(x, z i ; a)) ∂ a I(x, z i ; a)dt + Nz i=1 T 0 p 2,i ∂ a w(x, z i ; a)dt + T 0 p 3 ∂ a u(x; a)dt.

IV. NUMERICAL EXPERIMENTS

A. Numerical method 1) Gradient Algorithm : In order to tackle the optimization problem (17), we consider a gradient-based optimization method. The complete procedure is detailed in Algorithm 1. Note that in addition to a numerical tolerance criterion on

Algorithm 1 Gradient-based optimization algorithm

Input: Tol> 0, ρ > 0. Initial guess: a. Output: a Set err := Tol+1 and define h by (21) using the input data.

while err >Tol and h ∞ > h c do Compute u by [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF]. Set x, z as the solutions of the last two equations of (16). Compute I by [START_REF] René | An outlook on microalgal biofuels[END_REF].

Set C as the solution of the first equation of (16). Set p 1 , p 2 , p 3 as the solutions of (18). Compute the gradient ∇μ 1 by (19). a = a + ρ∇μ 1 , Set err := ∇μ 1 . end while the magnitude of the gradient, we have added a constraint on the water height h. The latter guarantees that we remain in the framework of subcritical flows [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF] (and in the range of industrial constraints, see [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF]).

A similar algorithm can be considered to tackle the multiple layers Problem (20). Remark that since there is no interaction between layers, the gradient computation can be partially parallelized when computing z i , C i , p 1,i , p 2,i .

2) Numerical Solvers: We introduce a supplementary space discretization with respect to x to solve our optimization problem numerically, in this way, we set a time step ∆t, and use Heun scheme to compute the trajectories described by (12). The integration is stopped when x reaches L. Denoting by N T the final number of the time steps, we then have N T ∆t ≈ T .

We use the Heun scheme again for computing C via (2). We use a first-discretize-then-optimize strategy, meaning that the Lagrange multipliers p 1 (resp. p 2 , p 3 ) are also computed by a (backward) Heun's type scheme. Note that this scheme is still explicit since it solves a backward dynamics starting from p 1 (T ) = 0 (resp. p 2 (T ) = 0, p 3 (T ) = 0).

B. Parameter settings

1) Parameterization: In order to describe the bottom of an optimized raceway pond, we choose to parameterize h by a truncated Fourier series for our numerical tests so that h(0) = h(L), which means we have accomplished one lap. More precisely, h reads:

h(x) = a 0 + N n=1 a n sin(2nπ x L ). (21)
The parameter to be optimized is the Fourier coefficients

a := [a 1 , • • • , a N ]
. Note that we choose to fix a 0 , since it is related to the volume V of the raceway. Indeed, under this parameterization we have

V = L 0 h(x)dx = a 0 L.
From ( 7) and ( 9), the velocity u and the topography z b read also as functions of a. Once we find the vector a maximizing the functional μNz , we then find the optimal topography of our system.

2) Parameters of the system: The time step is set to ∆t = 0.1 s which corresponds to the numerical convergence of the Heun's scheme, and we take Q 0 = 0.04 m 2 s -1 , a 0 = 0.4 m, z b (0) = -0.4 m to stay in standard ranges for a raceway. We take here L = 10 m and the free-fall acceleration g = 9.81 m s -2 . The initial state of C is set to be its steady state. All the numerical parameters values for Han's model are taken from [START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF] and recalled in Tab. I. In order to determine the light extinction coefficient ε, let us assume that only 10% light can be captured by the cells at the bottom of the raceway, meaning that I b = 0.1I s , we choose I s = 2050 µmol m -2 s -1 which corresponds to the average light intensity at summer (see [START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF]). Then ε can be computed by ε = (1/a 0 ) ln(I s /I b ).

C. Numerical results

1) Convergence test: The first test consists in studying the influence of the vertical discretization number N z . We fix N = 5 and take 100 random initial guesses of a. Note that the choice of a should respect the subcritical condition [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]. For N z varying from 1 to 80, we compute the average value of μNz for each N z . The results are shown in Fig. 3. We observe numerical convergence when N z grows, showing the convergence towards continuous model. In view of these results, we take N z = 40 for the successive studies.

2) Optimization example: We focus now on the shape of the optimal topography. Let us set the numerical tolerance Tol= 10 -10 , we choose N = 5 terms in the truncated Fourier series as an example to show the shape of the optimal topography. As an initial guess, we consider the flat topography, meaning that a is set to 0. The optimal shape of the topography is shown in Fig. 4, and the a for the final iteration reads a = [0.1043, 0.0503, 0.0333, 0.0250, 0.0201]. The resulting optimal topography is not flat, which is not a common knowledge in the community.

3) Study of the influence of N : The last test is given to study the influence of the order of the truncation of the Fourier series. Set N = [0, 5, 10, 15, 20] and keep all the other parameter settings. Tab. II shows the optimal value of our functional μNz for different values of N . Note that, even the value of μNz is small, it is not the numerical error. There is a slight increase of the optimal value of the functional μNz when N becomes larger. However, corresponding values of μNz remain close to the one associated with a flat topography. As for the optimal shape, for instance, we give in Fig. 4 another optimal topography when N = 20. The optimal topography seems to converge to a slope. Here, the topography varies significantly so that viscosity of the fluid can give rise to turbulence (and irreversibility) which is not taken into account in [START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water; numerical validation[END_REF][START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF].

V. CONCLUSIONS AND FUTURE WORKS

A non flat topography slightly enhances the average growth rate. However the gain stays very limited and it is not clear if the practical difficulty to generate such pattern would be compensated by the increase in the process productivity.

Here the class of functions representing the topography has been chosen so that the topographic functions are zero on average. As a consequence, their choice does not affect the liquid volume. Further studies could consist in considering more general class of functions, also considering a non constant volume. Moreover, the optimal algal biomass concentration is itself related to the growth rate ( [START_REF] Martínez | Theory of turbid microalgae cultures[END_REF]) so that the light extinction coefficient is no more a constant, and depends on biomass supported by the parameterization of the system. The computation of the gradient should then be revisited to manage this more complicated case.

Fig. 1 .

 1 Fig. 1. Scheme of the Han model, representing the probability to go from one state to another, as a function of the photon flux density.

Fig 2 .Fig. 2 .

 22 Fig. 2. Representation of the hydrodynamic model.

  I(x, z; a))C + ζ(I(x, z; a))dt, where C, x, z satisfy    Ċ = -α(I(x, z; a))C + β(I(x, z; a)) ẋ = u(x; a) ż = w(x, z; a).

) 2 ) 0 - 0 p 1 Ċ 0 p 2 ż 0 p 3 ẋ

 20010203 Optimality System: Define the Lagrangian of Problem (17) by L(C, z,x, p 1 , p 2 , p 3 , a) γ(I(x, z; a))C + ζ(I(x, z; a)) dt -T + α(I(x, z; a))C -β(I(x, z; a)) dt -T -w(x, z; a) dt -T -u(x; a) dt.

0 - 0 p 1 - 0 p 2 ∂ 0 p 3 ∂

 0010203 γ (I(x, z; a))C + ζ (I(x, z; a)) ∂ a I(x, z; a)dt + T α (I(x, z; a))C +β (I(x, z; a)) ∂ a I(x, z; a)dt + T a w(x, z; a)dt + T a u(x; a)dt.

5 Fig. 3 .

 53 Fig. 3. The value of μNz for Nz = [1, 80].

Fig. 4 .

 4 Fig.4. The optimal topography at the final iteration for N = 5 (up) and for N = 20 (down). The red thick curve represents the topography (z b ), the blue thick curve represents the free surface (η), and all the other curves between represent the trajectories for different layers.

TABLE II THE

 II VALUE OF μNz FOR DIFFERENT

	N 0	Iter 0	μNz 1.232270 10 -5	log 10 ( ∇μ Nz ) -
	5	39	1.250805 10 -5	-10.077503
	10	39	1.251945 10 -5	-10.091920
	15	39	1.252354 10 -5	-10.097072
	20	39	1.252565 10 -5	-10.099389

VI. ACKNOWLEDGEMENTS

This research benefited from the support of the FMJH Program PGMO funded by EDF-THALES-ORANGE.