Lei Mo
email: lei.mo@inria.fr

Member, IEEE Angeliki Kritikakou
email: ange-liki.kritikakou@irisa.fr

Shibo He

Energy-Aware Multiple Mobile Chargers Coordination for Wireless Rechargeable Sensor Networks

Keywords: Wireless rechargeable sensor networks, mobile charger coordination, perpetual operation, mixed-integer linear program, decomposition method

Wireless charging provides dynamic power supply for Wireless Sensor Networks (WSNs). Such systems, are typically considered under the scenario of Wireless Rechargeable Sensor Networks (WRSNs). With the use of mobile chargers (MCs), the flexibility of WRSNs is further enhanced. However, the use of MCs poses several challenges during the system design. The coordination process has to simultaneously optimize the scheduling, the moving time and the charging time of multiple MCs under limited system resources (time and energy). Efficient methods that jointly solve these challenges are generally lacking in the literature. In this paper, we address the multiple MCs coordination problem under multiple system requirements. Firstly, we aim at minimizing the energy consumption of MCs, guaranteeing that every sensor will not run out of energy. We formulate the multiple MCs coordination problem as a mixed-integer linear programming and derive a set of desired network properties. Secondly, we propose a novel decomposition method to optimally solve the problem, as well as to reduce the computation time. Our approach divides the problem into a subproblem for the MC scheduling and a subproblem for the MC moving time and charging time, and solves them iteratively by utilizing the solution of one into the other. The convergence of proposed method is analyzed theoretically. Simulation results demonstrate the effectiveness and scalability of the proposed method in terms of solution quality and computation time.

I. INTRODUCTION

With the rapid development of wireless energy transfer technology [START_REF] Kurs | Wireless power transfer via strongly coupled magnetic resonances[END_REF], Wireless Rechargeable Sensor Networks (WRSNs) have become an important research topic [START_REF] Akan | Internet of hybrid energy harvesting things[END_REF], [START_REF] Huang | Software defined networking for energy harvesting internet of things[END_REF]. Unlike the traditional Wireless Sensor Networks (WSNs), WRSNs are able to circumvent the limitations introduced by the energyconstrained sensor nodes, and, thus, achieving a perpetual network lifetime. To recharge a large-scale sensor network, it is necessary to use multiple mobile chargers (MCs). Compared with single MC dispatch problem, multiple MCs dispatch problem further involves the coordination among the MCs [START_REF] Lu | Wireless charging technologies: fundamentals, standards, and network applications[END_REF], [START_REF] Sangare | Mobile charging in wireless-powered sensor networks: optimal scheduling and experimental implementation[END_REF]. By properly coordinating the MCs, we can finish the complex charging task and also enhance the scalability and the robustness of the system. The multiple MCs coordination problem usually considers WRSNs with three main components: stationary sensor nodes, MCs and one stationary base station (BS). The sensor nodes periodically report their residual energies to the BS. Based on the collected information, the BS makes a charging decision and sends the control commands to the MCs, which will then schedule their movements and adjust their moving time and charging time. During the movement, the MCs stop at the working points and charge the sensor nodes.

The charging process operates in cycles. At each cycle, the MCs start from the BS, move forward the sensor nodes to charge them and return to the BS to charge or replace the battery after finish the assigned energy charging tasks.

The prior work on multiple MCs coordination optimization can be classified according to the following criteria: whether 1) the objective is single, i.e., the problem optimizes only the MC scheduling, or multiple, i.e., the problem jointly optimizes the MC scheduling, the moving and the charging time (time required for the MCs to charge the sensor nodes), 2) the energy of the MC is limited or unlimited, 3) the coordination mechanism is distributed or centralized, and 4) the solution is optimal or non-optimal. Table I provides a summary of several representative works in the literature. The majority of the works in this area deals with the single objective of MC scheduling and utilizes centralized mechanisms to coordinate the MCs. Although the distributed coordination is considered in [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF], [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF], [START_REF] Madhja | Hierarchical, collaborative wireless energy transfer in sensor networks with multiple mobile chargers[END_REF], the coordination process is divided into several phases and a distributed mechanism only exists in some specific phases. However, if multiple system requirements are taken into account, it is difficult to realize a fully distributed mechanism. By jointly considering the scheduling, the moving time and the charging time of the MCs, the replenished energy of the sensor nodes can be optimized. Therefore, we can avoid to fully charge the sensor nodes, leading to a better usage of system resources (time and energy). On the other hand, the number of sensor nodes is usually much larger than the number of MCs, and, thus, one MC requires to charge several sensor nodes in one cycle. If an MC does not have enough energy to finish the assigned tasks (e.g., moving or charging tasks), this MC has to return to the BS for battery charging or replacement. Therefore, the number of available MCs is time-varying, which influences the scheduling, the moving time and the charging time of the MCs. Moreover, the charging optimization problems usually aim to enhance the charging efficiency [START_REF] Rault | Avoiding radiation of on-demand multi-node energy charging with multiple mobile chargers[END_REF], [START_REF] Jiang | Joint charging tour planning and depot positioning for wireless sensor networks using mobile chargers[END_REF], [START_REF] Wang | Recharging schedules for wireless sensor networks with vehicle movement costs and capacity constraints[END_REF], [START_REF] Lin | GTCharge: a game theoretical collaborative charging scheme for wireless rechargeable sensor networks[END_REF] or reduce the energy consumption of the MCs [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF], [START_REF] Wang | NETWRAP: an NDN based realtimewireless recharging framework for wireless sensor networks[END_REF], [START_REF] Xu | Towards perpetual sensor networks via deploying multiple mobile wireless chargers[END_REF]. Therefore, the problem becomes more realistic, when the constraints related to the limited MC energy budget are taken into account.

It is worth noting that the extensions from: 1) the single objective case to the multiple objectives case, or 2) the MC energy unbounded case to the MC energy bounded case, are not straightforward. This is due to the fact that: 1) additional variables and constraints need to be added into the problem so as to satisfy the new requirements, and 2) the coupling between the optimization variables makes the problem hard to solve. Heuristics are popular methods to solve charging optimization problems [START_REF] Wang | NETWRAP: an NDN based realtimewireless recharging framework for wireless sensor networks[END_REF], [START_REF] Xu | Towards perpetual sensor networks via deploying multiple mobile wireless chargers[END_REF], [START_REF] Jiang | Effective on-demand mobile charger scheduling for maximizing coverage in wireless rechargeable sensor networks[END_REF]- [START_REF] Liang | Maintaining sensor networks perpetually via wireless recharging mobile vehicles[END_REF]. They are able to find a feasible (non-optimal) solution in a short time. However, they are sensitive to the modification of the problem structure and the parameters, while it is hard to guarantee the quality of the solution. Finding the optimal solution is very important. By doing so, we can find out how far is the system performance based on a feasible solution from the optimal one, and how to improve heuristic methods based on the quality of the solution. Since the charging behaviors of the MCs are coordinated by the BS (which has more computing resources than the sensor nodes and the MCs), BS is able to run complex MC coordination algorithms, including algorithm for optimal solution.

In this work, we address the multiple MCs coordination problem to perform the sensors charging task. The primary goal is to achieve perpetual operation, i.e., guarantee that the residual energy of each sensor node will never fall below the operational energy level. On this basis, we aim to minimize the total energy consumption of MCs by determining the scheduling of the MCs and adjusting their moving and charging time. According to this work, we answer the following questions: 1) What is the condition to achieve perpetual operation? 2) If this is possible, how to schedule the MCs to perform the charging tasks such that the total energy consumption of MCs is minimized? 3) Is there a way to achieve optimal solution while avoiding high computational complexity? The main distinction of this paper from previous work is that: the scheduling, the charging time and the moving time of multiple MCs are jointly optimized, under the constraint that the energy budget of the MCs is limited and none of the sensors will run out of its energy.

A. Related Work

In the literature, the MC scheduling is usually converted into the traveling salesman problem (TSP) [START_REF] Lu | Wireless charging technologies: fundamentals, standards, and network applications[END_REF], where a MC constructs a tour of all nodes only once. For the multiple MCs case, the TSP can be extended to the multiple traveling salesman problem (m-TSP) [START_REF] Wang | NETWRAP: an NDN based realtimewireless recharging framework for wireless sensor networks[END_REF]. Based on the weighted sum of the traveling time and the residual lifetime of nodes, a heuristic algorithm is proposed in [START_REF] Wang | NETWRAP: an NDN based realtimewireless recharging framework for wireless sensor networks[END_REF] to select the nodes to recharge. In [START_REF] Rault | Avoiding radiation of on-demand multi-node energy charging with multiple mobile chargers[END_REF], the point-to-multipoint energy transfer is considered and the aim is to maximize the number of nodes that are charged at each stop to reduce the charging delay. A TSP with multiple time windows (TSPMTW) is formulated and solved using constraint programming model. The aim in [START_REF] Xu | Towards perpetual sensor networks via deploying multiple mobile wireless chargers[END_REF] is to minimize the sum of traveling distances of all chargers. This can be formulated as a q-root TSP problem, i.e., find q closed tours covering all locations such that the total length of the q tours is minimized. Considering WRSNs with negligible charging time, the problem of minimizing the number of MCs for maintaining the operations of the networks is studied in [START_REF] Beigel | On optimal scheduling of multiple mobile chargers in wireless sensor networks[END_REF]. This problem can be optimally solved with linear complexity. The problem of maximizing the network coverage is considered in [START_REF] Jiang | Effective on-demand mobile charger scheduling for maximizing coverage in wireless rechargeable sensor networks[END_REF], where each sensor monitors a circular area. The focus of [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF] is to jointly schedule the MCs and adjust their charging time so that the energy consumption of the MCs is minimized. The authors design a decentralized method to find the optimal solution. However, the energy of the MCs is unlimited in all above studies, compared to our work.

The MC energy is limited in [START_REF] Jiang | Joint charging tour planning and depot positioning for wireless sensor networks using mobile chargers[END_REF]- [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF]. The authors in [START_REF] Jiang | Joint charging tour planning and depot positioning for wireless sensor networks using mobile chargers[END_REF] jointly optimize the BS positioning and the MC scheduling. They maximize the charging efficiency, while reducing the number of required MCs and BSs. The considered problem is solved by a three-step design: 1) scheduling planning, 2) candidate BS identification, and 3) BS deployment and scheduling assignment. In [START_REF] Wang | Recharging schedules for wireless sensor networks with vehicle movement costs and capacity constraints[END_REF], the aim is to maximize energy efficiency, while including no sensor node outage. A profitable TSP considering deadlines is formulated and a greedy algorithm and an adaptive algorithm are used to solve the problem. In [START_REF] Dai | Minimizing the number of mobile chargers for large-scale wireless rechargeable sensor networks[END_REF], the problem of minimizing the number of the chargers is proved to be NP-hard and an approximation algorithm is proposed to solve the problem. Under similar context to [START_REF] Xu | Towards perpetual sensor networks via deploying multiple mobile wireless chargers[END_REF], a q-root TSP is developed in [START_REF] Liang | Maintaining sensor networks perpetually via wireless recharging mobile vehicles[END_REF] to schedule the MCs, while further aiming to minimize the number of MCs. In [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF], based on the energy status of the sensors and the MCs, the network is divided into several regions. Each MC is assigned to a network region and charges the sensors in that region. However, the above studies mainly focus on the heuristic methods, compared to our work.

To enhance the flexibility of the charging process, the MCs are allowed to transfer energy not just to the sensor nodes, but also between themselves [START_REF] Zhang | Collaborative mobile charging[END_REF], [START_REF] Lin | GTCharge: a game theoretical collaborative charging scheme for wireless rechargeable sensor networks[END_REF]. The MCs are divided into two groups in [START_REF] Madhja | Hierarchical, collaborative wireless energy transfer in sensor networks with multiple mobile chargers[END_REF], one that charges the sensor nodes and the other charges the MCs. Such kind of two-layer architecture can be considered as future work. Moreover, the above studies mainly focus on the MC scheduling problem, since the MC charging time optimization is not taken into account except in [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF]. Note that different scheduling schemes lead to different charging time decision. The scheduling problem and the charging time problem should be jointly addressed to find the optimal solution.

With respect to the coordination mechanism, some works propose distributed approach. In [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF], [START_REF] Madhja | Hierarchical, collaborative wireless energy transfer in sensor networks with multiple mobile chargers[END_REF], the distributed mechanism exists among the region partition phase, i.e., the MCs communicate with each other to make the region partition decision. The results in [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF], [START_REF] Madhja | Hierarchical, collaborative wireless energy transfer in sensor networks with multiple mobile chargers[END_REF] show that the performance of the distributed coordination is lower than the centralized one, since distributed coordination utilizes local network information instead of global information to make the decision. In [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF], the MCs communicate with each other to make the charging time decision under the given MC scheduling decision determined by the BS. However, the actions between the sensors and the MCs are still coordinated by the BS in [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF], [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF], [START_REF] Madhja | Hierarchical, collaborative wireless energy transfer in sensor networks with multiple mobile chargers[END_REF]. The BS collects and processes the energy information from the sensors and sends it to the MCs. Different from our preliminary results in [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF] where the MC energy constraints are not taken into account and the solution is based on classical Benders decomposition (BD) approach [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF], this paper considers that the MC energy budget is limited and proposes a novel solution to further reduce the computation time.

B. Contributions

Our main contributions focus on how to formulate the jointdesign problem and how to solve this problem efficiently. More precisely:

1) We propose a novel multiple MCs coordination framework for WRSNs that jointly optimizes the scheduling, the moving time and the charging time of the MCs. Solving these correlated problems separately leads to non-optimal solution. To enhance charging efficiency, based on the energy consumption model of the sensor nodes and the charging model of the MCs, we derive a lifetime-based charging sequence for the sensor nodes and divide the sensor nodes into several groups according to the number of available MCs. We also obtain a set of desired network properties, such as how to determine a sensor node that requires charging in the current charging cycle or not, and how much energy a sensor node should be replenished so as to achieve perpetual operation. The joint-design problem of minimizing the energy consumption of the MCs, while guaranteeing the perpetual operation of the WRSNs, is formulated as an MILP problem. 2) We present an Optimal multiple MCs Coordination algorithm, referred to as OMC, to solve the joint-design problem. Since the moving time and the charging time of the MCs is influenced by the MC scheduling decision, the OMC decomposes the overall problem into an integer linear programming (ILP)-based master problem (MP) for the MC scheduling and a linear programming (LP)-based slave problem (SP) for the MC moving and charging time. We prove that by solving the MP and the SP in each iteration, we obtain a lower bound and an upper bound for the overall problem. Moreover, in each iteration, by constructing a new constraint according to the solution to the SP and adding this new constraint to the MP in next iteration, the gap between the lower bound and the upper bound gradually reduces. Unlike the classical BD approach, we prove that by relaxing the MP to an LP to find a feasible solution and by replacing the optimal solution to the MP with the feasible solution during the iteration between the MP and the SP, the optimality of the solution to the overall problem is still guaranteed. 3) Finally, we conduct extensive simulations results to analyze the quality of the solution, the computation time and the scalability of the proposed OMC algorithm.

The obtained results show that the proposed method is able to achieve the optimal solution with reduced computation time compared to state-of-the-art optimal methods. Moreover, we present how we can further balance the computation time with the quality of the solution by controlling the iteration process between the MP and the SP.

C. Paper Organization

The remainder of this paper is organized as follows. Section II presents the system model and formulates the problem. Then, the multiple MCs coordination mechanism is described in Section III. Finally, Section IV shows the simulation results and Section V concludes this work.

II. SYSTEM MODEL AND PROBLEM FORMULATION

In this section, we first present the sensor energy consumption model and the wireless energy charging model. Then, we formulate the multiple MCs coordination problem as an MILP that takes perpetual operation and energy efficiency into account. The notations followed in this paper are: for a matrix a, a ij is the (i, j) th element of a; for a vector b, b i is the i th element of b; (•) T is the operator for the transpose of a matrix/vector. Let x = [x 1 . . . , x n ] T and y = [y 1 . . . , y n ] T . x y represents x i ≤ y i , (1 ≤ i ≤ n). In the following, node is used in the place of sensor node.

A. System Model

We consider a WRSN with n nodes {s 1 , . . . , s n } and m MCs {c 1 , . . . , c m }. For each node s i (1 ≤ i ≤ n), its energy is mainly consumed for the data transmitting and receiving. We adopt the following energy consumption model [START_REF] Xie | Multi-node wireless energy charging in sensor networks[END_REF]:

r i = ρ n k=1,k =i f ki + n j=1,j =i C ij f ij + C ib f ib , (1)
where f ki (f ij ) and f ib are the data flow rate from node s k to node s i (from node s i to node s j ) and from node s i to BS, respectively. ρ and C ij (C ib ) are the rate of energy consumption for receiving a unit of data, and transmitting a unit of data from node s i to node s j (BS), respectively. We assume that the system runs with a given routing protocol [START_REF] Pantazis | Energy-efficient routing protocols in wireless sensor networks: a survey[END_REF], and, thus, the energy consumption rate r i of node s i is invariant with time.

On the other hand, we consider the following wireless energy charging model [START_REF] He | Energy provisioning in wireless rechargeable sensor networks[END_REF], [START_REF] Shu | Near-optimal velocity control for mobile charging in wireless rechargeable sensor networks[END_REF]:

p r = G s G r κ L p ω 4π(d + ζ) 2 p 0 = ς (d + ζ) 2 p 0 , (2)
where ς = GsGrκω 2 16Lpπ 2 , p 0 is the source power of the MC, p r is the received power of the node, d is the distance between the MC and the node, G s is the source antenna gain, G r is the receive antenna gain, L p is the polarization loss, ω is the wavelength, κ is the rectifier efficiency and ζ is a parameter to adjust the Friis' free space equation for the short distance transmission [START_REF] Na | Energyefficient mobile charging for wireless power transfer in internet of things networks[END_REF]. We consider point-to-point energy transfer [START_REF] Rault | Avoiding radiation of on-demand multi-node energy charging with multiple mobile chargers[END_REF], where a MC charges only one sensor node at a time by approaching it at a very close distance so the charging process has the maximum efficiency possible. As the experiment shown in [START_REF] Peng | Prolonging sensor network lifetime through wireless charging[END_REF], the charging efficiency η = pr po ≈ 6% when d ≈ 0. Moreover, we assume that all the MCs have the same charging power.

B. Problem Formulation

The parameters and the variables used during the problem formulation are summarized in Table II

i in cycle C k L k i lifetime of node s i v moving speed of MC c j τ s battery charging/replacement time of MC c j moving energy related coefficient of MC c j E j initial energy of MC c j in cycle C k E k j,l residual energy of MC c j at the beginning of round R k l t max j,l maximum charging time of MC c j in round R k l σ
maximum number of rounds in one cycle maximum charging and moving energy consumed by a MC in one round

m k l,1 number of schedulable MCs in the round R k l that estimated in round R k 1 (l ≥ 1) σ k i,1 number of rounds in cycle C k that estimated in round R k 1 with respect to node s i h k l number of schedulable MCs in round R k l p k l number of nodes require charging in round R k l σ k i number of
q k ij,l = 1 if MC c j is scheduled to charge node s i 0 otherwise t k ij,l time of MC c j spends to charge node s i g k ij,l moving time of MC from node s j to node s i TABLE II.
MAIN NOTATIONS and 2) it decides how much energy should be replenished for each node. It is worth noting that different nodes {s 1 , . . . , s n } have different initial energy levels {e k 1 , . . . , e k n } and different energy consumption rates {r 1 , . . . , r n }. To evaluate the energy charging priorities of the nodes, we sort all the nodes according to their lifetimes in an increasing order:

L k = {L k 1 , L k 2 , . . . , L k n }, (3)
where

L k i = e k i -emin ri
is the lifetime of node s i . There is no need to charge all the nodes {s 1 , . . . , s n } in one cycle C k , since some nodes have enough energy to work until a later cycle (e.g., C k+1 ). Therefore, in cycle C k , all the nodes {s 1 , . . . , s n } can be divided into two sets: 1) the serving set S k , and 2) the non-serving set Sk . If s i ∈ S k , this node is charged in the cycle C k , else, it can be charged in a later cycle.

Lemma 2.1: The cycle C k contains at most σ k rounds, where σ satisfies the inequalities:

σ k -1 l=1 m k l,1 < n, (4)
σ k l=1 m k l,1 ≥ n. (5)
Proof: For the MC c j , the worst case happens when 1) all the nodes {s 1 , . . . , s n } require charging in one cycle, 2) the nodes charged by the MC c j in each round require fully charging (i.e., the replenished energy is e max -e min ), and 3) the distances between the nodes that are charged by the MC c j in two adjacent rounds are the maximum (i.e. the moving distance is d max ). Therefore, the maximum charging and moving energy consumed by the MC c j in one round is

= e max -e min η + d max ,
where is an energy related coefficient.

If the initial energy of the MC c j in the cycle C k is E j , it can perform the charging task at least

ρ k j = E j = E j η e max -e min + d max η rounds (i.e., {R k 1 , . . . , R k ρ k j }).
Then, the MC c j has to return to the BS for battery charging/replacement. This process takes at most τ s + 2dmax v time, where v is the moving speed of MC. Since the charging time and the moving time of the MC c j in one round takes at most emax-emin pr + dmax v time, the battery charging/replacement of the MC c j takes at most φ = p r (τ s v + 2d max ) v(e max -e min ) + p r d max [START_REF] Sangare | Mobile charging in wireless-powered sensor networks: optimal scheduling and experimental implementation[END_REF] rounds. Therefore, the MC c j can be scheduled again after the round R k ρj +φ . We assume that the available energy of the MC c j after each battery charging/replacement is the same (i.e., E j ).

Based on the number of rounds that the MCs {c 1 , . . . , c m } can charge the nodes (i.e., {ρ k 1 , . . . , ρ k m }), we can estimate how many MCs are available in the round R k l (e.g., m k l,1

MCs). Note that n nodes require charging in the cycle C k . Based on the number of available MCs in each round (i.e., {m k 1,1 , m k 2,1 , . . .}), we can estimate how many rounds are included in the cycle C k , e.g., σ k rounds, where σ k satisfies (4) and ( 5).

Remark 2.1: In Lemma 2.1, the moving and charging energy consumed by the MC c j in one round is assumed to be . However, the real consumed energy of the MC c j is no more than , and, thus, the real number of rounds in the cycle C k is no more than σ k . Lemma 2.2: If the lifetime of the node s i (1 ≤ i ≤ n) satisfies the inequality:

L k i < (σ k i,1 + σ k+1 -1) + τ d , (7)
s i ∈ S k , else, s i ∈ Sk .
Proof: For the node s i , the worst case happens when 1) all the nodes {s 1 , . . . , s n } require to be charged in the cycles C k and C k+1 , 2) the charging time and the moving time of the MC in each round of the cycles C k and C k+1 is maximum (i.e., ), and 3) this node is charged again in the last round of cycle C k+1 . As illustrated in Fig. 1, if the lifetime of node s 1 satisfies the inequality

L k 1 ≥ (σ k 1,1 + σ k+1 -1) + τ d , s 1 ∈ Sk , else, s 1 ∈ S k
. This inequality implies that if s 1 is a non-serving node, it should have the enough energy to work σ k 1,1 rounds in the current cycle C k and σ k+1 -1 rounds in the next cycle C k+1 , where each round in the cycles C k and C k+1 takes time.

For the i th node in the sequence L k (i.e., s i ), we assume that n i nodes before it do not require charging and all the nodes after it require charging. Therefore, n -n i nodes require charging in the cycle C k . Based on the number of available MCs in each round (i.e., {m k 1,1 , m k 2,1 , . . .}), we can estimate how many rounds exist in the cycle C k (e.g., Remark 2.2: Note that 1) the real number of rounds in the cycle C k+1 is no more than σ, 2) the node s i may not be charged in the last round of cycle C k+1 , and 3) some nodes after the node s i may not require charging. The real waiting time of node s i is no more than (σ k i,1 + σ k+1 -1) + τ d . Let ñk (ñ k ≤ n) denote the number of nodes in the serving set S k . Therefore, the refined charging sequence is

σ k i,1 rounds). Therefore, if (7) is satisfied, s i ∈ S k , else, s i ∈ Sk . S 1 ! d R 1 k R 2 k R i k R 1 k+1 R 2 k+1 R " k+1 ! d S i " 1,1 k " i,
Lk = {L k 1 , . . . , L k m k 1,1 R k 1 , L k m k 1,1 +1 , . . . , L k m k 1,1 +m k 2,1 R k 2 , . . . , L k ñk }, (8)
and the refined charging cycle is defined as Definition 2.3: One charging cycle means that the nodes in the serving set S k are charged once.

Based on the refined charging sequence given by ( 8), we can schedule the MCs to charge the nodes.

2) Optimize Charging Behavior: Without loss of generality, we assume that the current round is R k l . If the residual energy of the MC c j satisfies the inequality:

E k j,l ≥ e max -e min η + 2d max , (9)
the MC c j is available in the round R k l , else, it has to return to the BS for battery charging/replacement. [START_REF] Xu | Towards perpetual sensor networks via deploying multiple mobile wireless chargers[END_REF] shows that the available MC should have enough energy for moving, charging a node and returning to the BS. We assume that h k l (h k l ≤ m) MCs are available in the round R k l .

Lemma 2.3: If the replenished energy of node s i satisfies the constraint:

E k i ≤ h k l j=1 p r t k ij,l ≤ E k i , (10)
where

E k i = [(σ k i + σ k+1 -1) + τ d ]r i -(e k i -τ k l-1 r i ), (11)
E k i = e max -(e k i -τ k l-1 r i ), (12)
no matter in which round of the next cycle this node is placed, it will never deplete its residual energy before being charged again.

Proof: If the residual energy of a node is larger than the energy consumed during the waiting time to be charged, this node can work perpetually. For the i th node in the sequence Lk (i.e., s i ), the worst case happens when 1) all the nodes after it (i.e., {s i+1 , . . . , s ñ}) require charging in the current cycle C k , 2) all the nodes (i.e., {s 1 , . . . , s n }) require charging in the next cycle C k+1 , 3) each round in the cycles C k and C k+1 takes time, and 4) this node is charged again in the last round of cycle C k+1 .

With the residual energy E k j,l , the MC c j is able to perform the charging task from round R k l to round R k l+ρ j,l . Then, it is available again after the round R k l+ρ k j,l +φ , where

ρ k j,l = E k j,l = E k j,l η e max -e min + d max η ,
and φ is given by ( 6). Therefore, we can estimate how many MCs can be scheduled in the rounds {R k l , R k l+1 , . . .} (e.g., {m k l,l , m k l+1,l , . . .}). Note that the number of uncharged nodes is known when charging the node s i (i.e., ñ -i + 1 sensors). Based on the number of available MCs in the later rounds (i.e., {m k l,l , m k l+1,l , . . .}), we can estimate how many rounds are left in the cycle C k (e.g., σ k i rounds). Since the cycle C k+1 contains at most σ k+1 rounds, the maximum waiting time of node s i to be charged again is (σ k i + σ k+1 -1) + τ d . On the other hand, since each node s i has a maximum energy level e max , the replenished energy of node s i should not exceed e max -(e i -τ k l-1 r i ), where e i -τ k l-1 r i is the residual energy of node s i in the current round R k l . Therefore, we have [START_REF] Beigel | On optimal scheduling of multiple mobile chargers in wireless sensor networks[END_REF].

When determining the sequence (8), we consider the worst case, which means s i ∈ Sk even if [START_REF] Wang | NETWRAP: an NDN based realtimewireless recharging framework for wireless sensor networks[END_REF] is satisfied. Substituting (7) into (11), we have

E k i > [τ k l-1 + (σ k i -σ k i,1 ) ]r i -e min .
Based on different values of the parameters, the lower bound of replenished energy E k i can be either positive or negative. If E k i < 0, s i ∈ S, since this node has enough energy to work until at the end of next cycle, and, thus, it can be removed from the charging sequence. Note that in the last round of a cycle, the number of sensors is no more than the number of MCs. We assume that p k l (p k l ≤ h k l ) nodes require charging in the round R k l . Since one MC is responsible for at most one node and each node is charged by one MC, the MC scheduling variable q k ij,l is bounded by

p k l i=1 q k ij,l ≤ 1, 1 ≤ j ≤ h k l , (13)
h k l j=1 q k ij,l = 1, 1 ≤ i ≤ p k l . (14)
On the other hand, since the charging time of the MC c j in the round R k l is limited by its maximum charging time t max j,l , we have

0 ≤ t k ij,l ≤ t max j,l q k ij,l , 1 ≤ i ≤ p k l , 1 ≤ j ≤ h k l . (15)
Note that in Lemma 2.3, the maximum waiting time of a node s i is assumed to be (σ k i + σ k+1 -1) + τ d . It means the interval between two adjacent rounds R k l and R k l+1 should not exceed . Therefore, we introduce a (continuous) variable g k ij,l to adjust the moving time of the MCs, and let the variables q k ij,l , t k ij,l and g k ij,l satisfy the constraints:

g k ij,l ≥ q k ij,l d ij v , 1 ≤ i ≤ p k l , 1 ≤ j ≤ h k l , (16)
0 ≤ t k ij,l + g k ij,l ≤ θ k l , 1 ≤ i ≤ p k l , 1 ≤ j ≤ h k l , (17)
where the time threshold θ k l is given by [START_REF] Lin | GTCharge: a game theoretical collaborative charging scheme for wireless rechargeable sensor networks[END_REF]. For the nodes that are charged in the round R k l , the number of the nodes in the later rounds is the same, and, thus, the number of remaining residual rounds is the same, which is denoted as σ k l . We assume that n k l-1 nodes require charging and nk l-1 nodes do not require charging in the previous rounds {R k 1 , . . . , R k l-1 }. Lemma 2.4: To ensure that the interval between two adjacent rounds is smaller than , we can set

θ k l = min{ , k l+1 , . . . , k l+σ k l }, (18)
where

k j = L ι k j -τ k l-1 -(j -l) -dmax v and ι k j = n k l-1 + nk l-1 + j-1 i=l m k i,l + 1 (l + 1 ≤ j ≤ l + σ k l )
. Proof: To determine the value of time threshold θ k l , we need to know how long the nodes in the later rounds {R k l+1 , . . . , R k l+σ l } can survive. Note that m k l,l MCs can be scheduled in the round R k l and the previous rounds

{R k 1 , . . . , R k l-1 } contain n k l-1 + nk l-1 nodes. The node s ι k l+1
in the next round R k l+1 has the shortest lifetime, where

ι k l+1 = n k l-1 + nk l-1 + m k l,l + 1. Since q k ij,l , t k ij,l , and g k ij,l
are variables, the exact charging and moving time in the round R k l is unknown. Therefore, we set τ k l = τ k l-1 + , since is an upper bound of the charging and moving time in one round. To guarantee that the nodes in the round R k l+1 can be charged in time, the maximum waiting time of the nodes in the round

R k l+1 should not exceed k l+1 = L ι k l+1 -τ k l -dmax v = L ι k l+1 -τ k l-1 --dmax v . For the nodes in the later round R k j (l + 1 ≤ j ≤ l + σ k l ), the maximum waiting time should not exceed k j = L ι k j - τ k j-1 -dmax v = L ι k j -τ k l-1 -(j -l) -dmax v , where ι k j = n l-1 + nl-1 + j-1 i=l m k i,l + 1.
Note that the time threshold θ k l should not exceed the maximum charging and moving time . Therefore, we have [START_REF] Lin | GTCharge: a game theoretical collaborative charging scheme for wireless rechargeable sensor networks[END_REF]. Theorem 2.1: The sufficient condition to achieve perpetual operation is that the constraints (10), ( 16) -( 18) are satisfied.

Proof: The constraint [START_REF] Beigel | On optimal scheduling of multiple mobile chargers in wireless sensor networks[END_REF] makes sure that all nodes will not run out of their residual energies before being charged again, under the condition that the interval between two adjacent rounds is smaller than a predefined time threshold, which is guaranteed by the constraints ( 16) - [START_REF] Lin | GTCharge: a game theoretical collaborative charging scheme for wireless rechargeable sensor networks[END_REF].

For the objective function, we consider minimizing the total energy consumption of the MCs. Based on this objective and all the aforementioned constraints, the Primal Problem (PP) is formulated as

PP : min Q k l ,T k l ,G k l Φ = p k l i=1 h k l j=1 (q k ij,l d ij + p 0 t k ij,l ) (19)
s.t. ( 10) - [START_REF] Lin | GTCharge: a game theoretical collaborative charging scheme for wireless rechargeable sensor networks[END_REF].

where

q k l = [q k ij,l ], t k l = [t k ij,l ], and g k l = [g k ij,l ].
Since the binary and the continuous variables are coupled with each other linearly, the PP is an MILP problem.

III. MULTIPLE MCS COORDINATION SCHEME

In this section, we present an optimal multiple MCs coordination (OMC) algorithm to solve PP. The structure of OMC is shown in Fig. 2

A. Formulations of MP and SP

For convenience, the matrices and the vectors are used to denote the constraints and the variables. Therefore, the PP is reformulated as

PP1 : min x,y Φ(x, y) = g T x + f T y (20) s.t. Ax b 1 , Cx + Dy b 2 ,
where x and y are the vectors of binary and continuous variables, respectively. g and f are the vectors of the objective function coefficients. A, C and D are the matrices of the coefficients in the constraints. b 1 is an u-dimensional vector and b 2 is an v-dimensional vector. Note that the objective function of the PP1 contains the binary variables x as well as the continuous variables y, while the MP only considers the binary variables x. To facilitate the iteration between the MP and the SP, we introduce an auxiliary (continuous) variable Φ into the MP as the objective function, where Φ and Φ have the same physical meaning. Based on the structure of the PP1, the corresponding MP is

MP : Φ L = min x, Φ Φ (21) s.t.    Ax b 1 , C 1 : Φ ≥ g T x + λ(i) T (Cx -b 2 ), ∀i ∈ A, C 2 : 0 ≥ ϕ(j) T (Cx -b 2 ), ∀j ∈ B,
where C 1 and C 2 are the sets of feasibility constraints (FCs) and infeasibility constraints (ICs), respectively. They are generated from the solution to the dual slave problem (DSP) [START_REF] He | Energy provisioning in wireless rechargeable sensor networks[END_REF].

A and B are the sets of iterations where the DSP has bounded and unbounded solutions, respectively. λ(i) is the solution to the DSP at the i th iteration, while ϕ(j) is the solution to the dual feasibility check problem (DFCP) (26) at the j th iteration. Let (x(l), Φ(l)) denote the solution to the MP at the l th iteration. Therefore, the corresponding MP is

SP : Φ U = min y 0 Φ(x(l), y) = g T x(l) + f T y (22) s.t. Cx(l) + Dy b 2 ,
Comparing the SP with the PP1, we observe that their formulations are the same, except that the binary variables x in the SP are fixed. Let (x * , y * ) denote the optimal solution to the PP1. We have Φ * = Φ(x * , y * ). Compared with the PP1, the MP only considers the MC scheduling variables x, whereas the constraints with respect to the charging and moving time variables y are relaxed. Solving the MP yields a lower bound Φ L of Φ * . On the other hand, since the MC scheduling decision x(l) may be just a feasible solution (not optimal yet), solving the SP with x(l) yields an upper bound Φ U of Φ * . Therefore, we have Φ L ≤ Φ * ≤ Φ U . To reduce the gap between Φ L and Φ U , a new constraint FC (or IC) is added into C 1 (or C 2 ) at each iteration. When the gap is smaller than a predefined threshold ε, the optimal solution (x * , y * ) is found.

B. Iterations between MP and SP 1)

Step 1 -Initialization: Initialize the iteration counter l = 0, the MP solution x(0), the lower bound Φ L = -∞, and the upper bound Φ U = ∞. The sets C 1 and C 2 are set to null. The initial solution x(0) can be given arbitrarily, as long as it satisfies the constraint Ax(0) b 1 .

2) Step 2 -Solving SP: In this paper, rather than solving the SP directly, we solve its dual problem. This is because the SP and the DSP are equivalent due to the strong duality [START_REF] Boyd | Convex optimization[END_REF], and the new constraints can be constructed according to the solution of the DSP.

To construct the dual of the SP, we introduce Lagrange multipliers λ

[λ i ] (1 ≤ i ≤ v) to the SP. Therefore, the DSP is DSP : max λ 0 g T x(l) + λ T (Cx(l) -b 2 ) (23) s.t. f + D T λ 0.
Since the DSP is an LP, it can be solved very fast using standard algorithms, such as simplex method or interior point method [START_REF] Mo | Distributed node coordination for real-time energy-constrained control in wireless sensor and actuator networks[END_REF].

3) Step 3 -Solving MP: Based on the solution to the DSP, we have:

1) If the DSP is infeasible, the SP has an unbounded solution. Therefore, the PP1 is infeasible. 2) If the DSP has a bounded solution λ(l), A ← {m} ∪ A.

We have g T x(l) + f T y(l) = g T x(l) + λ(l) T (Cx(l)b 2 ) due to the strong duality, where y(l) is the solution to the SP at the l th iteration. At the same time, the upper bound Φ U (l) is updated by

Φ U (l) = min{Φ U (l-1), g T x(l)+λ(l) T (Cx(l)-b 2 )}.
Since we have Φ(l) < g T x(l)+f T y(l) = g T x(l)+λ(l) T (Cx(l)-b 2 ), the stopping criterion is not met, and, thus, x(l) is a nonoptimal solution to the PP1. To avoid selecting again the non-optimal solution x(l), a new FC

Φ ≥ g T x + λ(l) T (Cx -b 2 ) ( 24
)
is generated and added into C 1 at the (l + 1) th iteration. 3) If the DSP has an unbounded solution (i.e., g T x(l) + λ(l) T (Cx(l) -b 2 ) = +∞), B ← {l} ∪ B. Due to the strong duality, the SP has no feasible solution under the given x(l). Note that the feasibility of the SP is related to the constraints rather than the objective function. This problem may be feasible if the positive (continuous) variables ξ Since the FCP is an LP, the strong duality exists between the FCP and its dual problem. Instead of solving FCP, we solve its dual problem. To construct the dual of the FCP, we introducing Lagrange multipliers ϕ = [ϕ i ] (1 ≤ i ≤ v) to the FCP. Therefore, the dual of the FCP (DFCP) is

[ξ i ] (1 ≤ i ≤ v)
DFCP : max ϕ 0 ϕ T (Cx(l) -b 2 ) (26) s.t. 1 -ϕ 0, D T ϕ 0.
The DFCP can be solved by a similar method used to solve the DSP. Let ξ(l) and ϕ(l) denote the solutions to the FCP and the DFCP at the l th iteration, respectively. If SP exists infeasible constraints, the corresponding variables are non-zero, while the others are zero. Due to the strong duality, we have

1 T ξ(l) = ϕ(l) T (Cx(l) -b 2 ) > 0.
To avoid selecting again the infeasible solution

x(l), a new IC 0 ≥ ϕ(l) T (Cx -b 2 ) ( 27
)
is generated and added into C 2 at the (l + 1) th iteration. Since solving DFCP and solving FCP are equivalent, and ϕ T (Cx -b 2 ) is a function with respect to x but not 1 T ξ (i.e., 0 ≥ 1 T ξ is an invalid constraint for the MP), we construct the IC through the solution to the DFCP rather than the solution to the FCP.

For the new generated FC [START_REF] Shu | Near-optimal velocity control for mobile charging in wireless rechargeable sensor networks[END_REF] or IC [START_REF] Boyd | Convex optimization[END_REF], all the parameters are constant except Φ and x. They are the variables to the MP. When the MP is solved, the iteration counter l increases, and Step 2 to Step 3 is repeated. The iteration stops when Φ U (l) -Φ L (l) ≤ ε is satisfied.

C. Convergence Analysis

From the MP ( 21), we observe that the real constraints are Ax b 1 and C 2 . C 1 can be treated as the objective function. Therefore, the MP can be solved by only considering the binary variables x. Let λ(k) denote the bounded solution to the DSP at the k th iteration (k ∈ A). Comparing the following ILP problem

Φr (k) = min x g T x + λ T (k)(Cx -b 2 ) (28) s.t. Ax b 1 , C 2 : 0 ≥ ϕ(j) T (Cx -b 2 ), ∀j ∈ B,
with the MP, we have Φ(l) = max ∀k∈A { Φr (k)}. Although the MP includes a continuous variable Φ, it can be solved by only considering the binary variables x. Lemma 3.1: The lower bound Φ L (l) and the upper bound Φ U (l) on the optimal objective function value Φ * are derived from the solution to the MP and the SP, respectively.

Proof: Without loss of generality, we assume that Φ(l) = Φr (k) (k ∈ A). Therefore, we have

Φ(l) = Φr (k) = min x g T x + λ T (k)(Cx -b 2 ) ≤g T x * + λ T (k)(Cx * -b 2 ) (29a) ≤ max λ 0 g T x * + λ T (Cx * -b 2 ) (29b) =Φ * , (29c)
where (29a) holds since x * is not the optimal solution to the problem [START_REF] Mo | Distributed node coordination for real-time energy-constrained control in wireless sensor and actuator networks[END_REF] with λ(k). (29b) holds since λ(k) is not the optimal solution to the problem ( 23) with x * . (29c) holds since solving the problem ( 23) with x * we obtain Φ * . From ( 29), we observe that Φ(l) is a lower bound of Φ * . Depending on the solution of the DSP, the value of the objective function of the DSP can be either finite or infinite. If the DSP has an unbounded solution, i.e., g T x(l) + λ(l) T (Cx(l)b 2 ) = +∞, it is obvious that +∞ is an upper bound of Φ * . Therefore, we focus on the case when the DSP has a bounded solution. Note that

Φ U (l) = min{Φ U (l -1), g T x(l) + λ(l) T (Cx(l) -b 2 )} = min 1≤i≤l {g T x(i) + λ(i) T (Cx(i) -b 2 )}, (30) and g
T x(i) + λ(i) T (Cx(i) -b 2 ) = min y 0 Φ(x(i), y) ≥ min y 0 Φ(x * , y) = Φ * , (31)
where [START_REF] Albert | Solving mixed integer linear programs using branch and cut algorithm[END_REF] holds since x * is the optimal solution to the problem [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF]. From ( 30) and ( 31), we observe that Φ U (l) is an upper bound of Φ * . Remark 3.1: Note that the SP has the same formulation as the PP1 except that the binary variables x in the SP are fixed. In addition, the SP and the DSP are equivalent due to the strong duality. From (29), we observe that the auxiliary variable Φ has the same physical meaning as the objective function of DSP, and, thus, Φ has the same physical meaning as the objective function of PP1.

Lemma 3.2: The lower bound sequence {Φ L (0), . . . Φ L (l)} is increasing, while the upper bound sequence {Φ U (0), . . . Φ U (l)} is decreasing.

Proof: With iteration l increasing, more constraints are added into the MP. Therefore, the feasible region of the MP shrinks. Since the MP is a minimization problem, the nonoptimal values of Φ * (i.e., { Φ(0), . . . , Φ(l)}) are excluded by the constraints C 1 and C 2 , and, thus, Φ L (l + 1) = Φ(l + 1) is larger than the previous lower bounds {Φ L (0), . . . , Φ L (l)}.

We assume that Φ U (l+1) > Φ U (l). This contradicts the fact that Φ U (l + 1) = min{Φ U (l), λ(l + 1) T (Cx(l + 1) -b 2 )}. Therefore, Φ U (l +1) is smaller than the previous upper bounds {Φ U (0), . . . , Φ U (l)}.

Theorem 3.1: With the FC [START_REF] Shu | Near-optimal velocity control for mobile charging in wireless rechargeable sensor networks[END_REF] and the IC ( 27) added into the MP, the algorithm converges.

Proof: Based on Lemma 3.1 and Lemma 3.2, as well as the fact that the non-optimal and the infeasible solutions of the binary variables x are excluded by the FCs and the ICs, and the dimension of the binary variables x is finite, the algorithm converges in a finite number of iterations.

Theorem 3.2: The FC and the IC generated by solving the DSP with x(l) do not exclude the optimal solution (x * , y * ), where x(l) is an arbitrary feasible solution to the MP at the l th iteration.

Proof: If the DSP has a bounded solution λ(l) with x(l), the corresponding FC is

Φ ≥ g T x + λ(l) T (Cx -b 2 ). (32)
On the other hand, if the DSP has an unbounded solution with x(l), the corresponding IC is

0 ≥ φ(l) T (Cx -b 2 ), (33)
where φ(l) is solution to the DFCP with x(l). In the following, we prove that the optimal solution (x * , y * ) to the PP1 does not violate the constraints [START_REF] Rothberg | An evolutionary algorithm for polishing mixed integer programming solutions[END_REF] and [START_REF] Randazzo | A comparison of optimal methods for local access uncapacitated network design[END_REF].

If the DSP has a bounded solution λ(l) with x(l), suppose that x * and Φ * violate the FC [START_REF] Rothberg | An evolutionary algorithm for polishing mixed integer programming solutions[END_REF]. Therefore, we have

Φ * < g T x * + λ(l) T (Cx * -b 2 ). (34)
However, (34) contradicts the fact that Φ * is the optimal value to the DSP with x * :

Φ * = max λ 0 g T x * + λ(Cx * -b 2 ) ≥ g T x * + λ(l) T (Cx * -b 2 ).
Therefore, the FC (32) does not exclude the optimal solution x * . On the other hand, if the DSP has an unbounded solution with x(l), x(l) is excluded by the IC [START_REF] Randazzo | A comparison of optimal methods for local access uncapacitated network design[END_REF]. Since x * = x(l), x * does not violate the IC (33). Theorem 3.3: With the FC (32) and the IC (33) added into the MP, the algorithm converges.

Proof: Since the dimension of binary variables x is finite, according to Theorem 3.2, the solution converges to the global optimal value within a finite number of iterations.

Remark 3.2: Since MP is an ILP, this problem is still hard to solve directly compared with the LP-based SP. Moreover, the size of MP will increase with the number of iterations, since at each iteration a new FC or IC is added into the MP. Therefore, the computational complexity of OMC is dominated by the cost of solving the MP at each iteration. Based on Theorem 3.3, we can circumvent the above difficulties by replacing the optimal solution to the MP x(k) with the feasible solution x(k) during the iteration between the MP and the SP. Such a feasible solution can be efficiently found by using the heuristics, such as feasibility pump (FP) method [START_REF] Fischetti | The feasibility pump[END_REF]. Algorithm 1 summarizes the implementation details of the OMC algorithm. Output: Optimal scheduling decision x * , optimal charging and moving time decision y * ; Set initial values:

l = 0, Φ L (0) = -∞, Φ L (0) = +∞, {x(0)|Ax(0) b 1 }, ε; C 1 and C 2 are set to null; while Φ U (l) -Φ L (l) > ε do
Solve MP [START_REF] Xie | Multi-node wireless energy charging in sensor networks[END_REF] to obtain temporary scheduling decision x(k); Update lower bound Φ L (l) = max ∀k∈A { Φr (k)}; Solve DSP [START_REF] He | Energy provisioning in wireless rechargeable sensor networks[END_REF] [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] is infeasible; end l ← l + 1; end Optimal scheduling decision x * = x(l); Submit x * into PP1 and solve it to obtain optimal charging and moving time decision y * ;

Φ ≥ g T x + λ(l) T (Cx -b 2 ) to C 1 ; Update upper bound Φ U (l) = min{Φ u (l -1), λ(l) T (C xT (l) -b 2 )}; else B ← {l} ∪ B; Add IC: 0 ≥ φ(l) T (Cx -b 2 ) to C 2 ; end else PP1 (
IV. SIMULATION We consider a WRSN as result case study. The mobile robots and the sensor nodes (e.g., Mica2) that equipped with the Powercast chargers and receivers are employed as MCs and rechargeable sensor nodes, respectively. The nodes and the MCs are randomly deployed in a 100 m × 100 m area to monitor the environment and charge the nodes. The BS is at (50, 50) m. The system parameters are summarized in Table III and are adopted from [START_REF] Xie | Multi-node wireless energy charging in sensor networks[END_REF], [START_REF] Peng | Prolonging sensor network lifetime through wireless charging[END_REF]. For the energy model given by (1), λ 1 and λ 2 are the distance-independent and the distancedependent constant terms, respectively. d ij is the distance between the nodes s i and s j . θ is the path loss index. For a regular AA battery, its nominal cell voltage and the quantity of electricity is 1.2 V/2.5 Ah. Since two AA batteries provide an average voltage 2.4 V for the Mica2 node and the operating limit is 2.1 V, we set e min = 2.1 × 2.5 × 3600 = 18900 J and e max = 2.4 × 2.5 × 3600 = 21600 J. Note that different MCs and sensor nodes only change the values of the problem parameters {A, C, D, g, f , b 1 , b 2 }. They do not affect the structure of our problem formulation. The simulations are performed on a PC with quad-core 2.5 GHz Intel i7 processor and 16 GB RAM, and the algorithms are implemented in Matlab 2016a.

Sensor node s i characteristics We present the following evaluation results: 1) the system performance (the energy status of the nodes) with the proposed OMC method; 2) the convergence iteration and the computation time of OMC and the trade-off between the solution quality and the computation time; 3) the algorithm performance (the solution quality and the computation time) comparison of OMC with: i) optimal approaches: decentralized Benders decomposition (DBD) [START_REF] Mo | Decentralized multicharger coordination for wireless rechargeable sensor networks[END_REF], branch and bound method (B&B) [START_REF] Boyd | Branch and bound methods[END_REF], and branch and cut method (B&C) [START_REF] Albert | Solving mixed integer linear programs using branch and cut algorithm[END_REF], and ii) evolutionary approach: genetic algorithm (GA) [START_REF] Rothberg | An evolutionary algorithm for polishing mixed integer programming solutions[END_REF]; and 4) the system performance (the energy consumption of the MCs) using OMC, m-TSP [START_REF] Wang | NETWRAP: an NDN based realtimewireless recharging framework for wireless sensor networks[END_REF] and region partition [START_REF] Madhja | Distributed wireless power transfer in sensor networks with multiple mobile chargers[END_REF] methods to schedule the MCs and charge the nodes.

r i = ρ n k=1,k =i f ki + n j=1,j =i C ij f ij + C ib f ib ρ = 50 nJ/b f ij ∈ [1, 10] kb/s C ij = λ 1 + λ 2 (d ij ) θ λ 1 =
Let g l (k) = si∈R 1 l (e * i -e i (k)) 2 /p 1 l denote the charging error of the l th round at step k, where the charging error is defined as the mean squared error (MSE) between the desired energy levels {e * i } and the residual energy levels {e i (k)} of the nodes. Fig. 3 shows the charging error in the first cycle C 1 , with n = 25 and m = 5. Based on the refined charging sequence L1 , the serving set S 1 contains 22 nodes (i.e., ñ1 = 22). The numbers of charged nodes p 1 l and available MCs h 1 l in each round R 1 l are listed in Table IV. From it, we observe that the cycle C 1 contains five rounds

{R 1 1 , R 1 2 , R 1 3 , R 1 4 , R 1
5 } and the total number of charged nodes in this cycle is less than 22. For some nodes in the charging sequence L1 , the lower bounds of the replenished energy {E 1 i } are negative, and, thus, these nodes can be removed from the charging sequence. Fig. 3 shows that there is no gap between two adjacent rounds, since the residual energy of the nodes decreases gradually without charging. The earlier the MCs start the next charging round, the less energy each node requires to be replenished.

Note that the nodes in the sequence L1 are sorted according to their lifetimes in an increasing order. The desired energy levels of the nodes in the previous rounds are higher than the desired energy levels of the nodes in the later rounds. Moreover, at the beginning of charging process, the charging errors of previous rounds are larger than the charging errors of later rounds. Fig. 4 shows the MSE between the minimum optional energy e min and the residual energy levels {e i (k)} of the nodes (i.e., g e (k) = n i=1 (e i (k) -e min ) 2 /n) in five cycles {C 1 , C 2 , C 3 , C 4 , C 5 }, with n = 25 and m = 5. Due to the introduction of the constraints [START_REF] Beigel | On optimal scheduling of multiple mobile chargers in wireless sensor networks[END_REF], ( 16) -( 18), the nodes will never run out their residual energies before being charged again. Step k Step k We define the OMC convergence iteration as the number of iterations required by the OMC to converge. The convergence iteration of OMC under different n (the number of nodes) and m (the number of MCs) parameters in one round is compared in Fig. 5. Since the number of charged nodes in one round is no more than m, we set n = m and change the value of n from 5 to 25 with a step of 5. With the value of n increasing, more variables and constraints are involved into the PP, and, thus, more iterations are required to find the optimal solution. Fig. 6 compares the computation time of OMC under different n and m parameters in one cycle (C 1 ), where the OMC computation time is defined as the time required to solve PP. Fig. 6 shows that under the given number of MCs (m), the computation time of OMC is almost linearly increased with the number of nodes (n). This is because in each round the number of charged nodes, as well as the number of available MCs, is no more than m, while the number of rounds will increase with n in one cycle. For the rest of the rounds, we obtain the similar results. Under the given ε, the iterations between the MP and the SP stop when the gap between the bounds Φ U (l) and Φ L (l) is smaller than ε and the DSP is feasible under the solution to the MP (i.e., x(l)). Accordingly, by adjusting the value of ε, we can control the quality of the solution. From Fig. 8 and Fig. 7, we observe that the solution quality and the computation time decrease, when ε increases. This characteristic inspires us a feasible PP solution can be found by using the following way: we stop the iterations between the MP and the SP when the SP is feasible for the first time (assume that at the l th iteration). Therefore, the corresponding MC moving time and charging time decision y(l) is obtained under the given MC scheduling decision x(l). The solutions to the PP found by the OMC, DBD, B&B, B&C and GA are compared in Fig. 9. definition of the data sets is similar to "DBD vs OMC". On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25 th and 75 th percentiles, respectively. The whiskers extend to the most extreme data points that are not considered outliers, while the outliers are plotted individually using the '+' symbol. Fig. 9 shows that the solutions given by the OMC, DBD, B&B and B&C are the same. Therefore, the OMC also finds the optimal solution, since B&B and B&C are able to find the optimal solution for the MILP problem. Note that PP is a minimization problem and there is no guarantee of convergence to a global optimum for GA. The solution provided by the OMC has a lower objective function value than GA. The convergence of GA is sensitive to the choice of the genetic operators, the mutation probability and the selection criteria, while fine-tuning of these parameters is often required. for all tuned m and n parameters. Similarly, we construct the data sets for the rest of algorithm comparisons. With the values of n and m increasing, the computation time of all compared algorithms grows. However, OMC has a shorter computation time than DBD, B&B, B&C and GA. Although the SP of DBD is solved in a distributed manner, the reduction of computation time is limited, since the SP is a LP problem. Moreover, the MP of DBD is still an ILP problem, which is more difficult to solve than the LP-based MP of OMC. B&C, which combines the benefits of B&B and Gomory cutting scheme, can better balance optimality, efficiency and stability. Usually, B&C has a faster convergence speed than B&B. Note that the compu-tational complexity of an optimization problem highly relates to the number of variables and constraints. Solving the smaller problems with less variables and constraints (i.e., MP and SP) iteratively is more efficient than solving a single large problem. This result is in line with the comparison in [START_REF] Randazzo | A comparison of optimal methods for local access uncapacitated network design[END_REF]. Compared with OMC, the structure of GA is more complex, since in each iteration GA needs to generate new populations through several procedures, such as selection, reproduction, mutation and crossover. Therefore, GA has a longer computation time. The basic idea of the m-TSP method is the use of a graph to model the charging process, where the vertexes and the edges represent the sensor nodes and the charging costs (i.e., the moving and charging time of the MCs), respectively. The aim of the m-TSP method is to minimize the charging costs under the constraint of visiting all the vertexes. On the other hand, the region partition method divides the sensors into several clusters and each cluster is assigned one MC, based on the energy status of the sensors and the MCs. The MC with higher energy level is assigned to the cluster that contains the sensors requiring to be replenished with higher amount of energy. For all methods, a MC will go back to the BS for battery charging/replacement when it has no enough energy to finish the assigned tasks, event if it is in the middle of charging phase. Fig. 11 compares the total energy consumption of the MCs using the OMC, m-TSP and region partition methods for one cycle. We set m = 5 and change the value of n from 25 to 50 with a step of 5. In contrast to the m-TSP and the region partition methods, the moving and the charging energy of the MCs are jointly optimized by the OMC method. This implies that there is no need to fully charge all the sensors. The replenished energy of a sensor is as much as it is required so as to make sure that this node has the enough energy to work until being charged again. As shown in Fig. 11, the OMC achieves lower energy consumption than m-TSP and region partition. Achieving a lower MC energy consumption, the times of the MCs back to the BS during the charging process can be reduced. In this way, the charging efficiency is further enhanced and the constraints with respect to the network perpetual operation can be more easily satisfied.

Let E o (m, n), E d (m, n), E b (m, n), E c (m, n

V. CONCLUSION

In this paper, we studied the problem of coordinating the multiple MCs to charge the sensor nodes, with the aim of enhancing the energy efficiency of MCs, while keeping the network operating perpetually. This problem has been formulated as an MILP problem, which jointly optimized the scheduling, the moving time and the charging time of the MCs. To find an optimal solution with reduced computation time, we proposed an OMC algorithm. It decomposes the multiple MCs coordination problem into two correlated subproblems: an MP for the MC scheduling and a SP for the MC moving time and charging time, and solves the subproblems iteratively. We proved that the proposed OMC algorithm converges to the optimal solution through a finite number of iterations. We also proved that the optimality of the solution is still guaranteed when the optimal solution to the MP is replaced with an arbitrary solution during the iteration between the MP and the SP. The results show that the desired system requirements are satisfied using the proposed OMC algorithm. We can also achieve a trade-off between the solution quality and the computation time by adjusting the stopping criteria of the OMC algorithm.

  rounds left in cycle C k when charging node s i τ k l total moving and charging time in rounds {R k 1 , . . . , R k l } τ d time interval between two adjacent rounds θ k l predefined time threshold k l maximum waiting time of nodes in round R k l ι k l index of node in round R k l that has the shortest lifetime Variables

#

  moving and changing time # moving and changing time Cycle C k Cycle C k+1

Fig. 1 .

 1 Fig.1. An example to determine if node s i requires charging in the current cycle C k or not.

Fig. 2 .

 2 Fig. 2. The structure of OMC algorithm.

  are introduced to relax the constraints. Based on this idea, we construct a Feasibility Check Problem (FCP) . Cx(l) + Dy b 2 + ξ.

Algorithm 1 :

 1 Optimal Multiple MCs Coordination (OMC) Algorithm Input: System parameters A, C, D, f , b 1 , b 2 ;

  50 nJ/b λ 2 = 0.0013 pJ/(b.m 4 ) θ = 4 e min = 18900 J e max = 21580 J ek i ∈ [19000, 20000] J ∆ s = 1 s Mobile charger c j characteristics η = 6 % p 0 = 5 W v = 1 m/s = 1 J E k j ∈ [1000, 5000] J τ s =10000 s, τ d = 50000 s Tuned parameters n m Min/Max/Step 25/50/5 5/15/5

Fig. 3 .

 3 Fig. 3. The charging error of the nodes in the first cycle C 1 .

Fig. 4 .

 4 Fig. 4. The energy status of the nodes in five cycles.

Fig. 7 and

 7 Fig. 8 evaluate the algorithm performance (the solution quality and the computation time) of OMC with the stopping criteria ε varying. We set n = 50 and m = 15 and change the value of ε between the range of [0.1, 1, 10, 100].

Fig. 5 .

 5 Fig. 5. Convergence iteration of OMC with n varying (n = m).

Fig. 6 .

 6 Fig. 6. Computation time of OMC and with m and n varying.

4 Fig. 7 .

 47 Fig. 7. Solution quality of OMC with ε varying.

Fig. 8 .

 8 Fig. 8. Computation time of OMC with ε varying.

Fig. 9 .

 9 Fig. 9. Energy gain of OMC, DBD, B&B, B&C and GA.

Fig. 10 compares

 10 Fig. 10 compares the computation time of the OMC, DBD, B&B, B&C and GA under different n and m parameters. Let T o (m, n), T d (m, n), T b (m, n), T c (m, n) and T g (m, n) denote the computation time of OMC, DBD, B&B, B&C and GA required to solve the PP under the given m and n parameters, respectively. The box plot of "DBD vs OMC" shows the statistical property of the data set {(T d (m, n)-T o (m, n))/T d (m, n)}for all tuned m and n parameters. Similarly, we construct the data sets for the rest of algorithm comparisons. With the values of n and m increasing, the computation time of all compared algorithms grows. However, OMC has a shorter computation time than DBD, B&B, B&C and GA. Although the SP of DBD is solved in a distributed manner, the reduction of computation time is limited, since the SP is a LP problem. Moreover, the MP of DBD is still an ILP problem, which is more difficult to solve than the LP-based MP of OMC. B&C, which combines the benefits of B&B and Gomory cutting scheme, can better balance optimality, efficiency and stability. Usually, B&C has a faster convergence speed than B&B. Note that the compu-

Fig. 10 .

 10 Fig. 10. Time gain of OMC, DBD, B&B, B&C and GA.

Fig. 11 .

 11 Fig. 11. MC energy consumption with OMC, m-TSP and region partition.

TABLE I

 I

	.	CLASSIFICATION OF SOME MULTIPLE MCS CHARGING
		OPTIMIZATION APPROACHES

  max maximum distance between two nodes e max maximum capacity of node battery e min minimum energy for a node to be operational

	Parameters
	n	number of nodes
	m	number of MCs
	d p 0	source power of MC
	p r	received power of node
	η	charging efficiency of MC
	r i	energy consumption rate of node s i
	e k i	initial energy of node s

. 1) Determine Charging Sequence: Definition 2.1: If n nodes are charged once, this process is called one charging cycle. Definition 2.2: If m MCs are scheduled to charge nodes, at most m nodes can be charged simultaneously. This process is called one charging round.

Since the number of nodes is usually much larger than the number of MCs (n m), one charging cycle contains several charging rounds. To achieve perpetual operation, the charging process contains two phases: 1) it determines the sets of the nodes that are charged in each cycle and in each round,

  . To solve PP, the most important step is to find a proper MC scheduling decision q k l . If the value of q k

				l
	is determined, PP is reduced to a LP problem, which has a
	simpler structure and it is easier to solve. Based on this idea,
	OMC tackles the problem by iteratively solving an ILP-based
	master problem (MP) to determine the MC scheduling and a
	LP-based slave problem (SP) to determine the MC charging
	and moving time.
	Feasibility Constraints
	Exclude non-optimal solutions
	for MC scheduling problem	Master Problem (MP)
				MC scheduling
	Infeasibility Constraints		problem
	Exclude infeasible solutions
	for MC scheduling problem
	SP Solution	Slave Problem (SP)		MP Solution
	Corresponding MC moving	MC moving and		Temporary MC
	and charging time decision	charging time problem	scheduling decision

TABLE III .

 III SYSTEM PARAMETERS

TABLE IV .

 IV THE REAL NUMBER OF CHARGED NODES IN CYCLE C 1 .

ACKNOWLEDGMENT

This research was partly supported by the National Natural Science foundation of China (Grant No. 61403340 and 61672458), the ANR ARTEFACT (AppRoximaTivE Flexible Circuits and Computing for IoT) project (Grant No. ANR-15-CE25-0015), and the Zhejiang Natural Science Foundation (Grant No. LR16F020001).