Discrete Particle Swarm Optimization Routing Protocol for Wireless Sensor Networks with Multiple Mobile Sinks

Mobile sinks can achieve load-balancing and energy-consumption balancing across the wireless sensor networks (WSNs). However, the frequent change of the paths between source nodes and the sinks caused by sink mobility introduces significant overhead in terms of energy and packet delays. To enhance network performance of WSNs with mobile sinks (MWSNs), we present an efficient routing strategy, which is formulated as an optimization problem and employs the particle swarm optimization algorithm (PSO) to build the optimal routing paths. However, the conventional PSO is insufficient to solve discrete routing optimization problems. Therefore, a novel greedy discrete particle swarm optimization with memory (GMDPSO) is put forward to address this problem. In the GMDPSO, particle’s position and velocity of traditional PSO are redefined under discrete MWSNs scenario. Particle updating rule is also reconsidered based on the subnetwork topology of MWSNs. Besides, by improving the greedy forwarding routing, a greedy search strategy is designed to drive particles to find a better position quickly. Furthermore, searching history is memorized to accelerate convergence. Simulation results demonstrate that our new protocol significantly improves the robustness and adapts to rapid topological changes with multiple mobile sinks, while efficiently reducing the communication overhead and the energy consumption.


Introduction
Recently, wireless sensor networks (WSNs) have gained enormous attention for their wide range of applications [1]. WSNs consist of many tiny sensor nodes and one or multiple sinks, sensor nodes gather data from the sensing environment to sinks by communicating with each other. Energy efficiency is the most important issue in WSNs due to the limited battery capacity of sensor nodes. In WSNs with static sinks, the nodes close to the sinks would become hotspots and die earlier than others, because they are more likely to be the intersection of multihop routes and need to deplete their battery supplies to transmit huge amounts of data for other sensor nodes to the sinks [2]. Node death will result in a series of problems, such as disruptions in the topology, reduction of sensing coverage and packets loss, and so on. Therefore, routing protocols designed for WSNs have to incorporate load-balancing in order to achieve uniformity of energy consumption throughout the network. Mobile sinks are a good solution to avoid the hotspots [3,4]. In mobile based routing protocols, the hotspots around the sink

Related Work
Various routing protocols have been proposed for MWSNs. TTDD (a Two-Tier Data Dissemination) [15] initially builds a rectangular grid structure which divides the network into cells with several dissemination nodes that are used to relay the query and data to/from the proper source. Whenever sinks request data, they query the network by local flooding within the cell and the query packets are relayed to the source nodes through dissemination nodes. A data path from the source to the sink is then established using the reverse of the path taken by the data request. However, TTDD suffers from the high overhead of constructing a separate grid for each source node. Unlike TTDD, Grid-Based Energy-Efficient Routing From Multiple Sources to Multiple Mobile Sinks (GBEER) [16] constructs a single combined rectangular grid structure for all possible sources to eliminate the high overhead of constructing separate grids for each source. However, the frequent grid change due to solving the hotspots introduces extra traffic on numerous nodes residing between the crossing points. In place of a rectangular grid, HPDD (Hexagonal Path Data Dissemination) [17] utilizes a common grid structure composed of hexagons. Which can provide shorter data and sink advertisement routes, but it suffers from the same hotspot problem. SEAD (Scalable Energy-efficient Asynchronous Dissemination) [18] constructs a minimum-cost weighted Steiner tree for the mobile sink by considering the distance and the packet traffic rate among nodes to save communication energy. Like TTDD, separate trees are constructed for each source; thus, the overhead of establishing such trees is very high. Ring Routing [19] designs an easily accessible ring structure to mitigate the anticipated hotspot problem with low overhead and minimize the data reporting delays. However, its questionable scalability is not good. Besides, the overhead of the initial ring construction for large or sparse networks will be high.
In agent-based schemes, one or more agents, which take on the role of representatives for the sources or the sink, are selected to relay the traffic between sources and the sink. These protocols

System Mode
We assume such a MWSN scenario: A wide open area is covered with a large number of homogeneous sensor nodes; there are N source nodes and M mobile sinks. Sinks can move anywhere within the sensor field, at any time with a fixed speed. Each source node sends its sensed data through multi-hops to the nearest sink via one and only one routing path. Several nodes can route packets to the mobile sink. In our network mode, individual node only owns the local information, such as its unique ID, residual energy level, delay and the distances to its neighbors, which can be estimated based on the received signal strength [24,25]. Our model does not need any extra positioning system such as GSP.
We use the same energy model and fault model proposed in [10]. The energy consumption equation of sensing and transferring m bit data is as follows: E tx pm, dq " pβ 1`β2ˆd ispi, jq n qˆm (2) where dis (i, j) is the distance between node i and node j, E sens is the energy consumption of sensing m bits of data, E tx (m, d) and E rx (m) are the energy consumption of sending and receiving m bits of data, n is the channel attenuation index, α 1 , β 1 , β 2 and γ 1 are energy consumption parameters of sensing circuit, sending circuit, sending amplifier and receiving circuit, respectively. Let P fault be the probability of sensor node's failure in the network. We assume that the probability of the mobile sink failure is approximately 0. If any relay node is failed, it is handled by our routing algorithm.

Terminologies
Throughout this paper, we use the following terminologies: ‚ Source j the jth sensor node. Sink i is the ith mobile sink. ‚ PTree i is the routing tree of Sink i . Path j P PTree i and RN k j P PTree i . ‚ l f pRN k j q means the lifetime of the RN k j . l f pPTree i q means the lifetime of the l f pPTree i q. ‚ DelaypRN k j ) is the delay of the RN k j . Delay pPTree i q means the overall delay of the PTree i . ‚ Len(Path j q means the length of the Path j . Len(PTree i q means the overall length of the PTree i . ‚ E cons pRN k j q is the energy consumed by RN k j due to data forwarding. ‚ E Path cons`P ath j˘i s the total energy consumption of routing path Path j . ‚ G(V,E) or G means the connected graph containing all candidate relay nodes.
‚ Set(V) is the set of all sensor nodes of G(V,E), i.e., the set of all candidate relay nodes.
‚ N V is the number of sensor nodes in Set(V). N 1 V is the number of sensor nodes in Set(V) except the agent node, N 1 V " N´1, i.e., the number of Set(V)-{Agent i }. ‚ Neig(node i ) is the set of all those sensor nodes within the communication range of node i , therefore, Neig(node i ) = {node j |0 < dis(node i , node j ) ď R}.
‚ NextRelay(node d ) is the sensor node which is selected as the next relay node for node d . Therefore, NextRelay (node d ) = {node j |@node j dis(node j , Agent) < dis(node d , Agent)}.

New Routing Protocol for MWSNs
In this section, we design an efficient routing strategy, based on the routing strategy, a new routing protocol is presented.

Our Efficient Routing Strategy in MWSNs
In our routing strategy, mobile sinks employ agent-based data gathering scheme [20,26,27] to collect sensed data. Our routing strategy can be described as follows.

Design of the Neighbor Table
Above of all, we design the neighbor table that is the key component of our protocol. It is illustrated in Table 1. self_ID means a node itself, neighbor_ID means its neighbor node, Dis is the distance between seft_ID and neighbor_ID. isRelay is a Boolean flag indicating whether the node pair <self_IDÑneighbor_ID > within a routing path, isRelay = 1 means true and data packet is routed from self_ID to neighbor_ID, its initial value is set to 0. Table 1 is the neighbor table for node 7 and illustrates that node 7 has two neighbor nodes: node 3 and node 9, the distance between node 7 and node 3 is 20, the distance between node 7 and node 9 is 15, in addition, <7Ñ3> within a routing path.
Notably, our neighbor table is a wonderful issue. Some advantages can be summarized as follows: Firstly, neighbor table is the basic component of our protocol. Its basic function is used to record neighbor information which is the key information used to build the optimal routing path. Secondly, neighbor table enhances the MWSN scalability and the protocol robustness. When some nodes are failed, they will be deleted from its neighbors' neighbor table on demand, and never be used to build the new optimal routing path. If the failed nodes are repaired or replaced, or new nodes are added in the MWSNs, corresponding neighbor information will be inserted into their neighbors' neighbor table, thus, they can be used to build the new optimal routing path. Beyond that, when the routing path is changed, it is only to reset the isRelay flag of the relevant nodes, and this feature is well suited to resolve the routing problem in MWSNs. Thirdly, the isRelay flag indicates the routing information. It serves two purposes: routing tree establishment and locating the failed nodes. The routing tree is established by resetting the isRelay flag. In addition, if a relay node has not received any packet from its previous relay node after a reasonable time interval, it implies that the previous relay node has failed. With the help of the isRelay flag, the relay node can quickly find failed node ID from its neighbor table.

Building Routing Tree
Here, we detailed the procedure of building the routing tree from source nodes to sink. This procedure consists of the following steps: Step 1 Network initiation and neighbor discovery. After all sensor nodes are deployed in the sensed field, all nodes can update their neighbor table using RSSI-based distant estimated scheme.
Step 2 Collect network statuses. When a sink move in a target area, the first thing has to be completed is to find the appropriate agent by local-broadcasting 1-hop Agent Query packet (AQ). The sink selected the closest node as the agent, and then broadcasts a Data Query (DQ) packet using unicast local flooding. Each node that receives the DQ packet firstly forwards the DQ to its neighbors and then updates its neighbor table if necessary. When a source matches DQ packet, it replies a Success Response (SR) packet to the sink. Each node that receives the SR packet firstly forwards the SR to its neighbors, then broadcast Data Query Response (DQR) packet that includes its newest state information such as residual energy, delay, and neighbor table data and so on along with its ID to the sink. Once the SR packet is received by sink, the data query procedure is finished immediately or finished after a reasonable time from then. This procedure as illustrates as Figure 1a. There are two special operations in this process are noteworthy. One is that all query packets (i.e., AQ and DQ) and their reply packets (i.e., DQR) are unicast by local flooding to reduce the overall communication and energy overhead. The other one is that, when a source node receivs multi DQ packets from multi sinks, it only replies to the first received DQ packet, which can assure it only sends its sense data to the nearest mobile sink to save energy.  Step 3 Build optimal routing tree. In most situations, multi-source nodes send data packets to a same sink, therefore, we need to build a routing tree. Root node is the agent, and the leaf nodes mean its source nodes. An example of routing tree is shown in Figure 2a. There are two routing paths, i.e., {Path1, Path2}, in the routing tree, the root node is Agenti and the leaf nodes are Source1 and Source2. The routing tree contains 5 relay nodes. Once finished Step 2. Sink has collected a number of neighbor tables along with their remaining energy information and delay information which are wrapped in DQR packets. As illustrated in Figure 1b, the network topology of candidate relay nodes that extracted from the neighbor tables, i.e., G(V, E), can be constructed by sink. Where ei∈E indicates that two adjacent sensor nodes can communicate with each other. G(V, E) is stored as adjacency matrix A, A[i, j] can be calculated as follows: where i, j are two adjacent nodes, and dis(i, j) denote the length of eij∈E. c. Each optimal routing path from each source node towards the sink can be built by using the improved greedy forwarding mechanism from candidates based on the network topology. Notably, even though each routing path (or branch) of a sink is optimal, the routing tree of the sink may not be the optimal due to the intersection relay nodes shared by multi paths. Therefore, the optimal algorithm needs to be employed to achieve the global optimum (the optimal routing tree). We designed GMDPSO , which is detailed in Section 5, to solve the problem.
Step 4 Establish routing tree and begin data transmission. After a mobile sink finished building the optimal routing tree in Step 3, it sends Routing Result (RR) packet, in which the routing tree can be stored as array, singly linked list or other forms, to its source nodes to establish the routing tree. The routing path establishment process can be detailed as follows: Assuming that Pathj = {Sourcej→11→9 →7→3→ → } is a routing path of the routing tree. sends RR packet to its Step 3 Build optimal routing tree. In most situations, multi-source nodes send data packets to a same sink, therefore, we need to build a routing tree. Root node is the agent, and the leaf nodes mean its source nodes. An example of routing tree is shown in Figure 2a. There are two routing paths, i.e., {Path 1 , Path 2 }, in the routing tree, the root node is Agent i and the leaf nodes are Source 1 and Source 2 . The routing tree contains 5 relay nodes. Once finished Step 2. Sink has collected a number of neighbor tables along with their remaining energy information and delay information which are wrapped in DQR packets. As illustrated in Figure 1b, the network topology of candidate relay nodes that extracted from the neighbor tables, i.e., G(V, E), can be constructed by sink. Where e i PE indicates that two adjacent sensor nodes can communicate with each other. G(V, E) is stored as adjacency matrix A, A[i, j] can be calculated as follows: where i, j are two adjacent nodes, and dis(i, j) denote the length of e ij PE.  Step 3 Build optimal routing tree. In most situations, multi-source nodes send data packets to a same sink, therefore, we need to build a routing tree. Root node is the agent, and the leaf nodes mean its source nodes. An example of routing tree is shown in Figure 2a. There are two routing paths, i.e., {Path1, Path2}, in the routing tree, the root node is Agenti and the leaf nodes are Source1 and Source2. The routing tree contains 5 relay nodes. Once finished Step 2. Sink has collected a number of neighbor tables along with their remaining energy information and delay information which are wrapped in DQR packets. As illustrated in Figure 1b, the network topology of candidate relay nodes that extracted from the neighbor tables, i.e., G(V, E), can be constructed by sink. Where ei∈E indicates that two adjacent sensor nodes can communicate with each other. G(V, E) is stored as adjacency matrix A, A[i, j] can be calculated as follows: where i, j are two adjacent nodes, and dis(i, j) denote the length of eij∈E. c. Each optimal routing path from each source node towards the sink can be built by using the improved greedy forwarding mechanism from candidates based on the network topology. Notably, even though each routing path (or branch) of a sink is optimal, the routing tree of the sink may not be the optimal due to the intersection relay nodes shared by multi paths. Therefore, the optimal algorithm needs to be employed to achieve the global optimum (the optimal routing tree). We designed GMDPSO , which is detailed in Section 5, to solve the problem.
Step 4 Establish routing tree and begin data transmission. After a mobile sink finished building the optimal routing tree in Step 3, it sends Routing Result (RR) packet, in which the routing tree can be stored as array, singly linked list or other forms, to its source nodes to establish the routing tree. The routing path establishment process can be detailed as follows: Assuming that Pathj = {Sourcej→11→9 →7→3→ → } is a routing path of the routing tree. sends RR packet to its Each optimal routing path from each source node towards the sink can be built by using the improved greedy forwarding mechanism from candidates based on the network topology. Notably, even though each routing path (or branch) of a sink is optimal, the routing tree of the sink may not be the optimal due to the intersection relay nodes shared by multi paths. Therefore, the optimal algorithm needs to be employed to achieve the global optimum (the optimal routing tree). We designed GMDPSO , which is detailed in Section 5, to solve the problem.
Step 4 Establish routing tree and begin data transmission. After a mobile sink finished building the optimal routing tree in Step 3, it sends Routing Result (RR) packet, in which the routing tree can be stored as array, singly linked list or other forms, to its source nodes to establish the routing tree. The routing path establishment process can be detailed as follows: Assuming that Path j = {Source j Ñ11Ñ9Ñ7Ñ3ÑAgent i ÑSink i } is a routing path of the routing tree. Sink i sends RR packet to its agent. When Agent i received the RR packet, it firstly search its previous relay node from the received packet according its ID, i.e., node 3. Then it sets isRelay = 1 for node pair {Agent i Ñnode3} in its neighbor table. Next, Agent i replaces the destination address of RR packet with node3 and forwards it to node3. After received the RR packet, node 3 forwards it to its previous relay node7 as the same way, and so on for each relay node in the routing path. Finally, the source Source j received the RR packet, which means that the routing path is established completely and source can begin the data packet transmitting via the reverse path of RS packet towards source node. The detailed data forward process is described as follows: Source j generates a data packet and queries the neighbor_ID from its neighbor table as its next relay node using the condition isRelay = 1, i.e., node11, then the data packet is send to node11. After received the data packet, node 11 immediately queries its next relay node (i.e., node 9) and forwards data packet to node 9 in the same way, and so on for each relay node in the routing path. Finally, all data packets from Source j are received by Agent i , and then Sink i receives data packet from Agent i directly.

Routing Path Recovery
When sink moves out of the radio rang of agent or relay node fails, the original routing path need to be recovered to continue the data gathering. The routing recover procedure is detailed in two cases: (1) When sink moves out of the radio rang of agent.
Sink i cannot receive any packet from Agent i , then it needs to select a new agent and build a new routing path to continue the data gathering. To ensure no previous data packet is lost, Temporary Agent node (TA) is introduced in our routing strategy. This routing recovery result is shown in Figure 2b. In order to help the sink to judge whether or not it is within the radio rang of an agent, the agent and relay nodes are required to send at least one packet to sink for a time interval T. Therefore, even though they have nothing to send, they also need to send a NULL packet to sink. Therefore, if not any packet has been received from Agent i for time interval T, Sink i can knows that it moves away from Agent i . To continue the data gathering, the following steps need to be completed: Step 1 Select the TA.
Sink selects the closest node as the TA. If the sink moves faster, then TA is farther away from Agent i ,which means that more relay nodes are required between TA and Agent i .
Step 2 Build the temporary routing path After TA is selected, Sink i sends Temporary Routing Path Setup (TRPS) packet to Agent i via TA. Once received the TRPS packet, Agent i sends its cached data packets which come from source Source j to TA along the reverse path of TRPS packet, the reverse path is named Temporary Path (TP). Because some data packets will be lost after the sink moved away from Agent i and before TP is established, the Agent i needs to cache the data packets in this time interval to avoid data loss. Once the Agent i received the TRPS packet, these cached data packets are routed to TA via TP.
The procedure is the same to Step 2 in Section 4.1.2. The only difference is that the sink sends Routing Path Reset (RPRS) packet to source Source j to collect the newly network status information.
Step 4 Build optimal routing path.
As same as Step 3 in Section 4.1.2. The only difference is that the optimal routing tree is built for the nearest source nodes using those nodes whose response packet (RPRSR) are received by sink.
Step 5 Establish routing tree and begin data transmission.
The procedure is the same as Step 4 in Section 4.1.2. Step 6 Clear original routing path.
Notably, when the new optimal routing path established, the original routing path is still work on. Before source node begin routing data packet via new optimal routing path, it sends Routing Path Clear (RPC) packet to TA along the original routing path and TP, each relay node on the path that received RPC packet will reset its isRelay = 0 to remove routing state information. Once the original routing path and TP is cleared, TA becomes the new agent Agent i .
(2) When the relay node fails Now we consider the relay node fault tolerant. This routing recovery procedure is shown in Figure 2c.
Step 1 Locate the failed node When a relay node (e.g., . . , RN n j , Agent i , Sink i } has failed, e.g., exhausted its battery power, then Path j is broken, its next relay node RN k`1 j can realize this situation if it cannot receive any packet from RN k j for a predetermined time interval T. Once detecting routing path is broken, RN k`1 j immediately queries its neighbor table to obtain the failure node ID (i.e., RN k j ) and sends a relay node failure (RNF) packet, which includes the failure node ID to Agent i via the original path. The relay node in the original path that received the RNF packet, e.g., RN k`2 j , forward it to the next relay node (RN k`3 j ) and reset its isRelay = 0 to remove routing information.
Step 2 Build new routing path without failed node After received the RNF packet, Agent i removes the failure node RN k j from the set of all candidates Set(V) and executes our new protocol to re-build new optimal routing path with Set(V)-{RN k j }. The subsequent steps are the same as routing recover for mobile moves away the current agent, which are described in Section 4.1.2.

Overview of New Protocol
Based on our efficient routing strategy, a new centralized discrete PSO routing protocol is presented to solve the problem of routing in MWSNs.
The new protocol consists of two phases, the various steps are depicted in the flowchart as shown in Figure 3. When a mobile sink moves in the sense area, the routing tree from nearest source nodes towards it will be established in phase 1. When a mobile sink moves away its agent, or its relay node fails, its original routing path is recovered in phase 2. The detailed description can be seen in Section 4.1. Step 6 Clear original routing path.
Notably, when the new optimal routing path established, the original routing path is still work on. Before source node begin routing data packet via new optimal routing path, it sends Routing Path Clear (RPC) packet to TA along the original routing path and TP, each relay node on the path that received RPC packet will reset its isRelay = 0 to remove routing state information. Once the original routing path and TP is cleared, TA becomes the new agent .
(2) When the relay node fails Now we consider the relay node fault tolerant. This routing recovery procedure is shown in Figure 2c.
Step 1 Locate the failed node When a relay node (e.g., } has failed, e.g., exhausted its battery power, then Pathj is broken, its next relay node +1 can realize this situation if it cannot receive any packet from for a predetermined time interval T. Once detecting routing path is broken, +1 immediately queries its neighbor table to obtain the failure node ID (i.e., ) and sends a relay node failure (RNF) packet, which includes the failure node ID to i Agent via the original path. The relay node in the original path that received the RNF packet, e.g., +2 , forward it to the next relay node ( +3 ) and reset its isRelay = 0 to remove routing information.
Step 2 Build new routing path without failed node After received the RNF packet, removes the failure node from the set of all candidates Set(V) and executes our new protocol to re-build new optimal routing path with The subsequent steps are the same as routing recover for mobile moves away the current agent, which are described in Section 4.1.2.

Overview of New Protocol
Based on our efficient routing strategy, a new centralized discrete PSO routing protocol is presented to solve the problem of routing in MWSNs.
The new protocol consists of two phases, the various steps are depicted in the flowchart as shown in Figure 3. When a mobile sink moves in the sense area, the routing tree from nearest source nodes towards it will be established in phase 1. When a mobile sink moves away its agent, or its relay node fails, its original routing path is recovered in phase 2. The detailed description can be seen in Section 4.1.

Step 1. Network initiation and neighbor discovery
Step 2.Collect network status information Step 4.Build routing path using GMDPSO

YES
Step 1. Select closest node Step 3 Select agent Step 2. Build temporary routing path Step 3. Collect network statuses Step 4. Build new routing path

GMDPSO
Step 5. Establish routing tree Step 6. Clear original routing path 2.2 Rout recover for nodes failed

NO
Step 1. Locate the failed node Step 2. Build new routing path without failed node

GMDPSO
Step 3. Establish routing tree Step 4. Clear original routing path   It can be seen from Figure 3 that the GMDSPO is the core algorithm of our new protocol, which is used to build the optimal tree from the nearest source nodes to the mobile sink in both phase 1 and phase 2.

PSO
PSO is a population-based stochastic searching algorithm, which is inspired by social behavior of bird flocking and fish schooling. The easy implementation, concise framework and fast computational convergence make PSO a popular optimization technique for solving continuous optimization problems.
PSO works with a swarm of a predefined size (say N p ) of particles. Each particle has a position and a velocity vector. The position vector gives a complete candidate solution to the optimization problem, and the velocity vector denotes the position-changing tendency. Each particle is evaluated by a fitness function to judge the quality of the solution to the problem .To search for the optimal solution, a particle iteratively updates its flying velocity and current position according to its own flying experience, i.e., personal best called Pbest i and according to the other particles' flying experiences, i.e., global best called Gbest.
In canonical PSO, a particle updates its position and velocity using the following simple rules: . . , x D i } are the ith particle's (say P i ) velocity and position vector. w is the inertial weight. t means the current iteration and t + 1 is the next itreration. c 1 and c 2 are acceleration factors term as cognitive and socail componeents. r 1 and r 2 are two different uniformly distributed random numbers in the range [0,1].
. . , pbest D i } is the personal best position of P i and Gbest = {gbest 1 , gbest 1 , . . . , gbest D } is the global best position of the whole swarm.
Conventional PSO was designed to optimize continuous problems. However, its fast convergence and easy implementation yet effectiveness have driven us to extend continuous PSO to solve the discrete routing optimization problem in MWSNs.

GMDPSO: Greedy Discrete PSO with Memory
The GMDPSO makes use of the network topology of candidate relay nodes to guide particle's position and velocity updates. The greedy strategy originated from greedy forwarding mechanism is proposed to direct particles to move towards better position. Some small operators, such as position initialization and memory mechanism, are introduced to speed up convergence. This section will describe the proposed algorithm in detail. The whole framework of GMDPSO is given in Algorithm 1.
It is worthy noted that the Prev variable is delicately designed in GMDPSO, it memorizes the search history to avoid the same position be updated more than once during the iterative process, and thus can avoid lots of reduplicative search and promote particle to find a better position faster. In other words, it can speed up the convergence of the algorithm. The function of Prev is called memory mechanism.
The positions are iteratively updated until the termination condition is satisfied. In our approach, there are two termination conditions: predefined iteration number and accuracy requirements. Once one of the two conditions is satisfied, the algorithm stops and the particle Gbest is found, Gbest.X represents the clustering result or routing result.

Algorithm 1 Framework of GMDPSO Algorithm
Parameters: particle swarm size N p , number of iterations Gen, inertia weight w, learning factors c 1 and c 2 ; Input: (1) adjacency matrix A for Sink i ; (2) the residual energy vector E res = {E res (1), E res (1) Step 5: output Gbest.X

Fitness Function Derivation
Fitness function is an important issue because it directly affects the final results. Now, we construct our fitness function (fitness) to evaluate the individual particle of population, in GMDPSO, the optimal tree is built such that it minimizes the fitness value. We have three objectives when build the optimal routing tree for Sink i (i.e., PTree i ): The first one is to maximize the lifetime of the routing tree to achieve energy-consumption balancing; The second one is to minimize the length of the routing tree to achieve energy conservation and enhance reliability; and the third one is to minimize the communication delay to enhance network throughout. In order to tune the three objectives to build optimal routing tree, we need to normalize the above three objectives and use the weight sum approach (WSA) [28] to construct our multi-objective fitness function fitness. Therefore, our fitness function depends on the following parameters: (1) Lifetime of relay node Our first objective is to maximize the routing tree lifetime, which is defined as the time interval from the establishment of the routing tree till its first relay node dies (FND). This is possible if we can maximize the lifetime of relay node that has least lifetime. Let RN k j be the kth relay node on the Pathj. The lifetime of RN k j is defined as follows: l f´RN k j¯" E res´R N k j¯{ E cons´R N k j¯( 7) where Eres(RN k j ) is the residual energy of relay node RN k j , and E cons pRN k j q is its consumed energy. It can be seen from Equation (7) that, even though RN k j keeps more E res (RN k j ), it maybe die faster if its E cons pRN k j q is bigger than others at the same time. Now, we begin the energy consumption analysis for relay node. RN k j must consume energy to forward the incoming data packets, which come from source node whose routing path to sink goes through it. Before calculating the routing energy consumption, it is needed to calculate the total number of incoming packets, which come from other relay node toward sink. Since multi routing paths may share some relay nodes, the number of incoming packets can be recursively calculated as follows: The relay node RN k j will consume its energy to receive N in (RN k j ) incoming packets and forward them. Therefore, the total data forwarding energy consumption of RN k j can be calculated as follows: where m is the size of data packet. Let minL f PTree i be the minimum lifetime of the routing tree PTree i . Therefore, our first objective is: Using Equation (9) to substitute the E cons pRN k j q in Equation (7), then, l f´RN k j¯" E res´R N k j¯{ E cons´R N k j¯" E res´R N k j¯{ N in´R N k j¯ˆmˆ´γ 1`β1`β2ˆd is´RN k j , RN k`1 j¯n¯ ( 11) It can be seen that, if E res (RN k j ) and N in (RN k j ) are fixed, then lf (RN k j ) is negatively correlated with dis(RN k j , RN k`1 j ). Therefore, it can be deduced that relay node with minimum lifetime is the relay node with the longest transition distance. That is, where argmaxf (x) returns the value of x that maximizes f (x). Bigger minL f PTree i means smaller fitness value. Therefore, fitness is inversely proportional to the minL f PTree i , i.e., f itness 9 1{minL f PTree i where 0 < 1/minL f PTree i ď 1. minL f PTree i is used to avoid node with lower residual energy being selected as sharing relay node whose load burden is heavy.
(2) Routing Path Length In our network mode, the energy is mainly consumed to collect network statuses and forward data to sink. The unicast local flooding mechanism is adopted to save energy for collecting network statuses. It is well known that the energy consumption of sensor node in WSNs is subject to the transmission distance-shorter data dissemination paths lead to decreased energy consumption together with increased throughput and reliability, which also can be deduced from Equation (9). Therefore, our second objective is to minimize the length of the routing tree to minimize the forwarding data energy consumption. Let Path j be routing tree of Agent i . Now, we calculate the total energy consumption of PTree i (i.e., E PTree i cons ). According to Equation (9), the energy consumption of each routing path Path j (Path j Ď PTree i ), i.e., E Path j cons , can be calculated as follows: Therefore, E PTree i cons can be calculated as follow: For each Path j assume that no packet is lost during routing data to agent, then N in (RN k j ) is the same for each relay node in Path j which can be considered as a constant value. Let N in (RN k j ) = N in . Then Equation (14) is improved as follows: where ř N Path j i"1 dis´RN k j , RN k`1 j¯i s length of the routing path Path j . It is can be seen that E Path j cons is only positively correlated with routing path.
Therefore, to minimize E PTree i cons , the total routing tree length Len(PTree i ) needs to be minimized. That is, our second objective can be described as follows: ( 17) Smaller Len (PTree i ) is better. Therefore, fitness is proportional to the Len (PTree i ) i.e., f itness 9Len pPTree i q Now, we normalize Len (PTree i ) by allLen, i.e., allLen can be calculated as follows: allLen " # 1 i f node j within the communicaiton rang o f node i , @i, j : where N V is the number of candidate relay nodes. Therefore, It can be deduced from Equation (21) that 0 < Len (PTree i ) ď 1. (

3) Communication Delay
End-to-end delay is another important performance metric of MWSNs. Given a fixed channel bandwidth, less the delay, higher the throughout. Let Delay pPTree i q be the total communication delay of PTree i . Then Delay pPTree i q can be calculated as following: Therefore, our third objective is: ( 23) Clearly, smaller Delay pPTree i q is better. Then our fitness is proportional to the Delay pPTree i q i.e., f itness 9Delay pPTree i q Similar with Equation (19), Delay pPTree i q also can be normalized as follows: After normalizing the above three objectives, the final fitness function for PTree i is: where ω 1 , ω 2 and ω 3 are three control parameters, 0 < ω i ď 1 and ω 1 + ω 2 + ω 3 = 1. In the paper, we give the same weight to them, that is, ω 1 = ω 2 = ω 3 = 0.33.

Particles Representation
How to encode the routing tree is very critical. To make GMDPSO feasible for discrete MWSNs scenario, the position and velocity of particle in our protocol are redefined.

Definition 1 (Position).
The position vector provides a routing tree. The position vector of ith particle is defined as is dimension of X i which means the number of candidate relay nodes excluding the agent, the d component x d i P {0, N} is the sensor ID, which maps node k (k = x d i ) as the next relay node of node d . That is to say, x d i indicates that node d forwards data to node k .
Notably, X i means a routing path tree that includes multi routing paths from multi source nodes to the same root node. Suppose mobile sink Sink i has N S source nodes, and then X i includes N S routing paths. A graphical illustration of particle representation can be seen in  Notably, Xi means a routing path tree that includes multi routing paths from multi source nodes to the same root node. Suppose mobile sink Sinki has NS source nodes, and then Xi includes NS routing paths. A graphical illustration of particle representation can be seen in Figure 4.  9,10,5,9,5,5,13,13,13,4, 3}, in Xi, for example, x 11 I = 4 indicates that source node node11 sends data packet to node4. A routing tree is constructed by decoding the encoded particle, in which each route path from a source node to sink can be built by appending relay nodes one by one till the agent node is selected as the end.

Illustration 1.
Consider mobile sink Sink i masters a subnet of MWSN G (V, E) with 13 sensor nodes and one mobile node as shown in Figure 4a. V = {node 1 , node 2 , . . . , node 13 }. node 13 is selected as agent, node 11 and node 12 are two source nodes, which means two routing paths will be built. The dimension of the particle position vector is N 1 = N´1 = 13´1 = 12. As shown in Figure 4b. Particle i is encoded as X i = {x 1 i , x 2 i , x 3 I , . . . , x 12 i } = {6, 9, 10, 5, 9, 5, 5, 13, 13, 13, 4, 3}, in X i , for example, x 11 I = 4 indicates that source node node 11 sends data packet to node 4 . A routing tree is constructed by decoding the encoded particle, in which each route path from a source node to sink can be built by appending relay nodes one by one till the agent node is selected as the end. Figure 4c, the path tree can be decoded as follows: firstly, the routing path of Source 1 (i.e., node 11 ) is built as Path 1 : Source 1 Ñnode 4 Ñnode 5 Ñnode9ÑAgent 1 .Then, routing path of Source 2 (i.e., node 12 ) is built as Path 2 :Source 2 Ñnode 3 Ñnode 10 ÑAgent 1 , and at last, the whole routing tree can be constructed by combining Path 1 and Path 2 .

As shown in
It can be seen from Figure 4 that our discrete position definition is straightforward and easy to decode and will lower the computational complexity, especially in the case of large-scale MWSNs, because the dimensions of the fitness function is equal to the size of candidate relay nodes collected by sink which is smaller than the entire MWSN size.

Definition 2 (Velocity).
Velocity is a very crucial component in PSO, by working on the position vector, it guides a particle and determines whether it can reach its destination and by how fast it could. Our discrete velocity of particle i is defined as In canonical PSO, velocity is used to learn knowledge from itself and swarm and finally leads the particle to a better position. In addition, a threshold V max is used to inhibit particles from flying out of the boundaries because there is a situation whereby when the speed of a particle is substantial. Unlike the continuous optimization, we have known that, in our discrete MWSNs scenarios, to compare two different routing paths from the same source node to the same agent, we only take care about whether their relay nodes are the same. Furthermore, we only need to compare the two relay nodes ID value in the same position of two different routing paths. There are only two results: is equal or no, therefore, the velocity can be encoded binary, and we defined 0 means two relay nodes are the same, 1 means they are different.
The first motivation of the velocity definition is to actually reflect the differences between two position vectors. The second motivation of the velocity definition is to prevent particles from flying away, because our velocity is binary-coded, we no longer need V max parameter.

Particle Swarm Initialization
A good initialization mechanism can reduce the searching space to reach global optima faster and promote diversity. Conventional random initialization method for PSO based algorithm is not applicable for our algorithm. The main reason is that random sequence of edges usually results in invalid routing tree that does not terminate on the agent node or that have loops. Therefore, we need to design a more efficient initialization method for our protocol. Based on our particle representation for discrete MWSNs, the position vectors initialization focus on how to map the next relay node, in other words, how to select the next relay node for the current one, for example, how to map another node as the next relay node for current node node d . Our main idea to solve this problem is that the next relay node for the current one is randomly selected from its neighbors. The mapping is done as follows: Let Neig(node d ) = {node 1 , node 2 , . . . , be the neighbors of node d , and A is the adjacent matrix for G (V, E), Then, node k " Index pNeig pnode d q , nq , 1 ď n ď |Neig pnode d q| (27) where Index (Neig(node d ), n) is an indexing function that indexes nth node of Neig(node d ) as the next relay node, and n is a randomly generated uniformly distributed integer number. Table 2 shows the nodes and their neighbors. Besides, Table 1 also illustrates how the next relay node is chosen. {node 6 , node 8 , node 10 , node 12 } -node 10 node 4 {node 2 , node 5 , node 7 , node 11 } 2 node 5 node 5 {node 2 , node 4 , node 6 , node 7 , node 9 , node 10 } -node 9 node 6 {node 1 , node 3 , node 5 , node 7 , node 10 } -node 10 node 7 {node 4 , node 5 , node 6 , node 11 } 2 node 5 node 8 {node 3 , node 10 , node 13 } -node 13 node 9 {node 2 , node 5 , node 10 , node 13 } -node 13 node 10 {node 3 , node 5 , node 6 , node 8 , node 9 , node 13 } -node 13 node 11 {node 4 , node 7 } 1 node 4 node 12 {node 1 , node 3 } 2 node 3 node 13 {node 8 , node 9 , node 10 , sink} sink Blue color means the corresponding node are mapped in step1. Green color means the corresponding node are mapped in step 2.Black color means the corresponding node are mapped in step 3.
In order to reduce the randomicity and blindness of swarm initialization, and at the same time speed up the convergence of our algorithm, the position vectors are initialized in such way: Step 1 Agent node is forced to be mapped as the relay node of its neighbor nodes (seeing the blue number in Figure 4b. In detail, firstly, position vector of i particle is empty, that is, Xi = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0). Then assign agent node (i.e., node 13 ) to its neighbor nodes, as shown in Table 2, then Xi = ({0, 0, 0, 0, 0, 0, 0, 13, 13, 13, 0, 0}. Marker '-' in Table 2 means that the relay node of the corresponding node is specified directly instead of choosing from its neighbors using Equation (27).
Step 2 Maps agent's neighbors to be relay node of the agent's neighbor's neighbor. Notably, if some of the agent's neighbors have the same neighbor, which cause more than one node will be mapped as the relay node of the same neighbor at the same time, and then we chose the nearest neighbor as the relay node. For example, as shown in Table 2, the neighbor of the agent (i.e., node 13 ) is Neig (node 13 ) = {node 8 , node 9 , node 10 }, then node 13 is the common next relay node of Neig (node 13 ). Next, we assign node 8 , node 9 , and node 10 to their neighbors respectively. Neig(node 8 )YNeig(node 9 )YNeig(node 10 ) = {node 2 , node 3 , node 5 , node 6 }, which means that node 8 , node 9 , and node 10 should be assigned to node 2 , node 3 , node 5 and node 6 ,that is to say, the third column of Table 2 (i.e., column n) of the corresponding nodes, i.e., node 8 , node 9 , node 10 , node 2 , node 3 , node 5 and node 6 should be marked with '-'. It is easy to be seen form Figure 4a that Neig (node 2 ) = node 9 , Neig (node 6 ) = node 10 . For node 5 , because Neig (node 9 )XNeig(node 10 ) = {node 5 }, which means that node 9 can be chosen as the next relay node for node 5 , and so do node 10 . We chose node 9 as the relay node of node 5 with assumption that node 9 is the nearest neighbor of node 5 . For node 6 , we choose node 10 as its next relay node in the same way. After finished Step 2, then the position vector Xi will look as following: Xi = {0, 9, 10, 0, 9, 10, 0, 13, 13, 13, 0, 0}.
Step 3 The relay nodes of remaining nodes are mapped by randomly choosing a node from their neighbor. For example, for Source 1 (i.e., node 11 ), let random n = 1, then its first neighbor is chosen as the next relay node.
Step 4 Finally, after finished the above three initiation steps, position vector Xi is initiated as Xi = {6, 9, 10, 5, 9, 10, 5, 13, 13, 13, 4, 3}. It is worth noting that the invalid routing tree, in which contains one or more invalid routing path (that does not terminate on the agent node or that have loops), will be punished with a very high fitness.
The velocity vectors are initialized as all-zero vectors. The P best vectors are initialized in the same manner as the position vectors, and the vector is set as the best position vector in the original population.
Comparing to the conventional random initialization, the advantages of our initiation method are follows: (1) Since node can and only can select its next relay node from its neighbors based on the network topology of G(V, E), our method takes advantage of this feature to reduce the vast search space significantly. (2) In our method, agent node and its neighbors are mapped firstly in step1 and step3, and this reduces the likelihood of invalid routing path. In addition, it further reduces search space and drives particle to find its personal best position faster. That is, it speeds the algorithm convergence.

Velocity and Position Update
In canonical PSO, velocity gives a particle the moving direction and tendency. After updating its velocity, the particle builds its new position using the new velocity. However, in the proposed algorithm, the particle position and velocity vectors have been redefined in a discrete integer form, and thus, the mathematical updating rules in the canonical continuous PSO no longer fit the discrete case. In order to meet the requirements of building routing tree in discrete MWSNs, the particle's velocity and position updating rules have been redefined as follows: where w is the dimensional inertial weight vector, c 1 and c 2 are two N 1 V dimensional cognitive and social components, r 1 and r 2 are two N 1 V dimensional random vectors in rang [0,1]. Equation (28) is used to update the old velocity, and Equation (29) is the position updating rule. It can be observed that the new updating rules take the same form in canonical PSO but different the key operator. The following will detail our new operators in discrete PSO. Operator). a. Position a Position builds a velocity vector. Assume that we are given two position vector

Definition 3 (Position
Definition 3 is inspired from the following two aspects: First, it is well known that, in canonical PSO, a particle adjusts its velocity by learning from its old velocity (V i (t)), personal best position (Pbest i -X i (t)) and the swarm global best position (Gbest-X i (t)). It can be seen that the learning process is actually a comparison between the positions, that is to say, for the current node from all its neighbors. However, unlike greedy mechanism that only chooses the nearest neighbor as the best next relay node, our approach chooses the neighbor node, which achieves the best balance of lifetime, distance and delay as the best next relay node by using our fitness function. Therefore, with the beginning from Source j . Path j can be built by selecting the best relay node one after one until the agent node is chosen, and Path j that is built in this way is the optimal routing path. For example, let first bestRN 1 be RN 1 j , then Path 1 j " p Source j Ñ RN 1 j q is the current optimal routing path. Let Path 2 j be the bestRN 2 of RN 1 j , Path 2 j " p RN 1 j Ñ RN 2 j q is the current optimal routing path. Then, where Path j is also the optimal routing. Using the same way, the other routing branches of X t can be built. Next, the fitness of the X t can be calculated using Equation (26).
This search mechanism of bestRN i (i.e., Equation (34)) actually act as our particle searching strategy in GMDPSO, in which a particle update its position by selecting the next relay node that can generate the largest decrement of the fitness value, so it can be regard as a greedy local search strategy.
Our motivation of searching bestRN i in this way is based on the divide-conquer strategy, in order to build an optimal routing tree of X t , we firstly build each optimal routing path branch of X i respectively. When all routing path branches are built completely, then, the whole routing path tree is finished. Similarly, in order to build the optimal routing path branch Path j , we choose the best next relay node for the previous relay node one after one, until the whole Pathj is completed.
However, Equation (34) is only used to build the single optimal routing path branch, rather than build the routing path tree (i.e., X i ). Because when we choose the best relay node for the current node of Path j , we do not consider the effect of the other routing path branch. Therefore, this method may neglect such a case, each routing branch is optimal, but the whole routing tree are weak due to the intersection relay nodes of multi routing paths. In other words, optimal routing path branches do not mean the optimal routing tree. As shown in Figure 5, many source nodes send their data packets to same next relay node, and these relay nodes will consume more energy to forward more packets. For example, node 6 need forward data packets of Source 1 , Source 2 and Source 3 to node 10 , while, node 10 needs forward data packets of Source 1 , Source 2 , Source 3 , Source 4 and Source 5 to Agent 1 . The worst result is that these sharing relay nodes will die quickly because they deplete their energy of forwarding too many packets. This will lead to the expensive routing recover. is finished. Similarly, in order to build the optimal routing path branch Pathj, we choose the best next relay node for the previous relay node one after one, until the whole Pathj is completed. However, Equation (34) is only used to build the single optimal routing path branch, rather than build the routing path tree (i.e., Xi). Because when we choose the best relay node for the current node of Pathj, we do not consider the effect of the other routing path branch. Therefore, this method may neglect such a case, each routing branch is optimal, but the whole routing tree are weak due to the intersection relay nodes of multi routing paths. In other words, optimal routing path branches do not mean the optimal routing tree. As shown in Figure 5, many source nodes send their data packets to same next relay node, and these relay nodes will consume more energy to forward more packets. For example, node6 need forward data packets of Source1, Source2 and Source3 to node10, while, node10 needs forward data packets of Source1, Source2, Source3, Source4 and Source5 to Agent1. The worst result is that these sharing relay nodes will die quickly because they deplete their energy of forwarding too many packets. This will lead to the expensive routing recover. Therefore, after all of individual optimal routing path branches are built, the corresponding routing tree Xi should be evaluated using fitness function.
From the above description of our new GMDPSO algorithm, the proposed algorithm has the following features: (1) with a concise framework; (2) the newly defined particle position and velocity are direct and easy to decode; and (3) redefined updating rules based on the new operators are easy to realize and will significantly reduce the computational complexity, especially for large-scale WSNs.

The Redesigned GMDPSO Seems to be Very Suitable for Solving Routing Problem in MWSNsSimulations and Results
In this section, we test our protocol against several well-known protocols: ECPSOAR, IAR, and TTDD, which all can deal with routing problem of multi mobile sinks. The performance is compared Therefore, after all of individual optimal routing path branches are built, the corresponding routing tree X i should be evaluated using fitness function.
From the above description of our new GMDPSO algorithm, the proposed algorithm has the following features: (1) with a concise framework; (2) the newly defined particle position and velocity are direct and easy to decode; and (3) redefined updating rules based on the new operators are easy to realize and will significantly reduce the computational complexity, especially for large-scale WSNs.

The Redesigned GMDPSO Seems to be Very Suitable for Solving Routing Problem in MWSNsSimulations and Results
In this section, we test our protocol against several well-known protocols: ECPSOAR, IAR, and TTDD, which all can deal with routing problem of multi mobile sinks. The performance is compared in terms of the following metrics: average packet delivery ratio (PDR, measured as the average number of successfully delivered packets versus required packets per round), average end-to-end delay (EED) [29], and average energy consumption ratio (ECR, measured as the average energy consumption from source to sink versus the initial energy per round). All simulations are performed using MATLAB R2012b on Windows 7 with Intel core i5-2520M Dual-Core CPU (2.50 GHz) and 8 G RAM. For ease of description for the comparison results of the above metrics, we define another notation, which is called metric comparative advantage for our new protocol as defined below: Metric comparative advantage " For example, PDR comparative advantage to ECSPOAR can be calculated as follow:

Simulation Setting
Simulations are performed on the MWSN, which consists of diverse number of homogenous sensor nodes ranging from 50 to 450. Each sensor node is assumed to have initial energy of 120 J and the mobile sink is assumed to have sufficient energy and cannot be fault.
To build a level playing field, the characteristics of the networks and communication models are configured as illustrated in [10] and as shown in Table 3. The extra PSO parameters used for ECPSOA are fixed to: particle updating energy consumption E PU = 80 pJ the endocrine selection energy consumption E ES = 50 pJ at per iteration, function dimension D = 30, division factory k = 6, and maximum iteration P Gen = 800 .The extra PSO parameters used for GMDPSO are fixed to: ω = 0.7968, c1 = c2 = 1.4926 In addition, the population size of ECPSOA and DPSORR is set to be 60; both of the two algorithms are run for maximum of 800. For the weight sum approach, in our proposed algorithms, we give equal weight to each sub-objective. That to say, we set w 1 = w 2 = w 3 = 0.33.

Performance of GMDPSO
First, we compare the performance of the proposed GMDPSO with ECPSOA, standard PSO [30] and the CPSOA [31]. In order to ensure a fair comparison, we configure these tree algorithms based on the same fitness function in Equation (26) with 450 sensor nodes, all simulation parameters are set to the same value, and the iterated generation for three protocols is fixed to 300.
The test result is shown in Figure 6. It can be seen that: GMDPSO outperforms the other PSOAs in term of convergence rate and the minimum fitness value. This is mainly due to the greedy search strategy based on the MWSNs topology, which avoid the blind search of the particle; another reason is that the memory mechanism reduces the repeated and invalid searching. It also can be seen that GMDPSO has the best initial fitness value, which is achieved by our special particle initial mechanism. For the weight sum approach, in our proposed algorithms, we give equal weight to each sub-objective. That to say, we set w1 = w2 = w3 = 0.33.

Results and Analysis
6.2.1. Performance of GMDPSO First, we compare the performance of the proposed GMDPSO with ECPSOA, standard PSO [30] and the CPSOA [31]. In order to ensure a fair comparison, we configure these tree algorithms based on the same fitness function in Equation (26) with 450 sensor nodes, all simulation parameters are set to the same value, and the iterated generation for three protocols is fixed to 300.

Packet Delivery Ratio
A routing protocol's reliability depends on PDR to the sink. Here, we run these protocols for comparing average PDR. Figure 7 and Table 4 show the comparison results when the speed of sinks is 5 m/s. Figure 7 shows that our protocol is more reliable and robust. In detail, firstly, average PDR of all protocols are decreasing as the number of nodes increases; however, our protocol still outperforms the others and the advantage (i.e., Equation (17)) is becoming more and more obvious as the number of sensor nodes increase in Figure 7d. Secondly, PDR reduces with the increase of the node failure probability (Figure 7a-c). Still, our protocol can deliver more packets than the other protocols with the same network size. Thirdly, with each same P f ault , t , our protocol always keeps the maximum average PDR (Figure 7a-c), which means that our protocol can find a better solution than others respect to different sensor nodes. Meanwhile, Table 4 illustrated that, with each same sensor node, our protocol keeps the maximum average PDR, which means that our protocol can find a better solution than others respect to different P f ault . That is to say, no matter with different sensor nodes or different node failure probability, our protocol can find the best path tree among the 4 protocols. Therefore, based on the above two results, we conclude that our protocol is more robust in PDR.
Our protocol improved the average PDR noticeably due to the following three reasons: The first reason is, the temporary path (TP) is designed to continue the old data packets transmission when sink moves away the old agent, which can avoid the packets loss before the new optimal routing path coming into services. The second reason is: our protocol can quickly locate the failure node and quickly build an alternative optimal path to recover the broken path link by its fast convergence feature and without any flooding when relay nodes fail, which also can reduce the data packets loss. The third reason is: the proposed GMDPSO can build better alternative routing path than others due to its ability of achieving global optimum, such as shorter total transmission distance and small communication delays, which can enhance the communication reliability and PDR. Moreover, our fitness function minimizes the energy consumption of the relay nodes to reduce the premature death probability of the relay node, which also can reduce the data packets loss due to the broken routing path.
outperforms the others and the advantage (i.e., Equation (17)) is becoming more and more obvious as the number of sensor nodes increase in Figure 7d. Secondly, PDR reduces with the increase of the node failure probability (Figure 7a-c). Still, our protocol can deliver more packets than the other protocols with the same network size. Thirdly, with each same ,t, our protocol always keeps the maximum average PDR (Figure 7a-c), which means that our protocol can find a better solution than others respect to different sensor nodes. Meanwhile, Table 4 illustrated that, with each same sensor node, our protocol keeps the maximum average PDR, which means that our protocol can find a better solution than others respect to different . That is to say, no matter with different sensor nodes or different node failure probability, our protocol can find the best path tree among the 4 protocols. Therefore, based on the above two results, we conclude that our protocol is more robust in PDR. Our protocol improved the average PDR noticeably due to the following three reasons: The first reason is, the temporary path (TP) is designed to continue the old data packets transmission when sink moves away the old agent, which can avoid the packets loss before the new optimal routing path coming into services. The second reason is: our protocol can quickly locate the failure node and   Next, we compare the average EED of the proposed protocol on the same experiment environment as Section 6.2.1. Figure 8 and Table 5 show the comparison result. It can be observed from Figure 8 that the proposed protocol has smaller end-to-end delay than the existing protocols.
More specifically, firstly, with the increase of the number of nodes, EED of all these protocols is increased. Again, our protocol performs best. This is because our protocol adopts the GMDPSO algorithm, whose better performance in reaching the global optimum allows it to build the optimal routing tree with shorter transmission distance. Furthermore, its faster convergence feature make it can build the optimal routing tree more quickly. Moreover, by using our neighbor table, the routing information is stored in each node to improve the speed of route establishment. All these advantages of GMDPSO can help to decrease the end-to-end delay. It is worthy to note that, as illustrated in Figure 8d, the average EED comparative advantage (i.e., Equation (17)) decreases as the number of sensor nodes increases, this is because our greedy search rule spends more time to select the best relay nodes from larger scale nodes to build the optimal routing tree. Nevertheless, the fast response of routing recovery and less communication control overhead by our unicast flooding mechanism also make our protocol end-to-end delay lower than the other protocols. Secondly, the end-to-end delay increases with the increase of the node failure probability by comparing Figure 8a-c, as the less stable topology causes more route recovery operation, which consumes more time for the protocols to maintain the network and prolong the delay. However, our protocol can still achieve the optimal delay; this is because we have designed the quick routing recovery mechanism for failure relay nodes. Thirdly, similar to average PDR, it can be observed from Figure 8a-c and Table 5 that our protocol achieves best average EED with respect to different sensor nodes and different P f ault , which means that our protocol can find a better solution than others. That to say, our protocol keeps better robustness of EED, this is because that the GMDPSO adopted in our protocol can also build the global optimum routing tree in different network size, however, others maybe build different suboptimal routing solutions for different network sizes, which increase the volatility of delay. Moreover, once the existing routing path is broken due to failure nodes, an alternative optimal path can be quickly established for the source node in our protocol. In addition, the other failure nodes (i.e., not the relay nodes) never be selected as the relay nodes due to our fitness function. quickly build an alternative optimal path to recover the broken path link by its fast convergence feature and without any flooding when relay nodes fail, which also can reduce the data packets loss. The third reason is: the proposed GMDPSO can build better alternative routing path than others due to its ability of achieving global optimum, such as shorter total transmission distance and small communication delays, which can enhance the communication reliability and PDR. Moreover, our fitness function minimizes the energy consumption of the relay nodes to reduce the premature death probability of the relay node, which also can reduce the data packets loss due to the broken routing path. Next, we compare the average EED of the proposed protocol on the same experiment environment as Section 6.2.1. Figure 8 and Table 5 show the comparison result. It can be observed from Figure 8 that the proposed protocol has smaller end-to-end delay than the existing protocols. More specifically, firstly, with the increase of the number of nodes, EED of all these protocols is increased. Again, our protocol performs best. This is because our protocol adopts the GMDPSO algorithm, whose better performance in reaching the global optimum allows it to build the optimal routing tree with shorter transmission distance. Furthermore, its faster convergence feature make it can build the optimal routing tree more quickly. Moreover, by using our neighbor   Although mobile sink protocol can alleviate hotspots implicitly by changing the possible high energy consumption zones around the sinks as the sinks move. However, these operations may cause the overall energy consumption in the network to increase. Now, we compare the average ECR of these protocols, which is used to measure the influence of node failure probability and mobile sink speed to the network. Simulations are performed on different network sizes with different sink speeds (speed of 5 m/s, 10 m/s and 20 m/s). Here, the node failure probability is set to 0.1. Figure 9 and Table 6 illustrate that the energy consumption of our protocol is the smallest among these protocols, and the comparative advantage (i.e., Equation (17)) becomes larger as the number of sensor nodes increases in Figure 9d. It worth to note that the average ECR of all these protocols increases as the mobile sinks move faster, because the change of the frequent topology results in frequent routing recovery which introduces heavier communication and energy overhead due to flooding operation. In this situation, our protocol can still consume less energy than the other protocols with the same network size. Besides, similar to average PDR, we can also conclude that our protocol achieves better robustness of EED than others from Figure 9a-c and Table 6.

Conclusions
In this paper, the routing of MWSNs is formulated as an optimization problem and we employ PSO to design an efficient routing protocol to achieve higher energy efficiency and lower communication delay. However, conventional PSO was originally designed for continuous optimization problems, which limits its application in discrete optimization domains. In addition, conventional PSO suffers from the curse of dimensionality, i.e., its performance deteriorates quickly as the dimensionality of the search space increases exponentially. To address these problems, we design a novel GMDPSO to   Our protocol outperforms others in term of ECR. The main reasons are as follows: the improved greedy forwarding mechanism is used ensure each routing branch has the mining energy consumption by selecting the relay node with optimized QoS parameters (energy, delay, energy consumption and so on). In addition, the unicast local flooding mechanism reduces the communication overhead in the network, which can minimize and indirectly reduce energy consumption.
In summary, the overall performance of our protocol outperforms LEACH, SEP, ERP and TPSO-CR in terms of the PDR, EED and EDR, while maintaining the best robust.

Conclusions
In this paper, the routing of MWSNs is formulated as an optimization problem and we employ PSO to design an efficient routing protocol to achieve higher energy efficiency and lower communication delay. However, conventional PSO was originally designed for continuous optimization problems, which limits its application in discrete optimization domains. In addition, conventional PSO suffers from the curse of dimensionality, i.e., its performance deteriorates quickly as the dimensionality of the search space increases exponentially. To address these problems, we design a novel GMDPSO to build the optimal route tree. In GMDPSO, we first deduced a new more suitable fitness function, then redefined the particle position and velocity in a discrete form and subsequently redesigned the particle update rules based on the network topology; consequently a discrete PSO framework was established. When applying the proposed discrete PSO framework to solve the mobile sink route problem, to alleviate prematurity, a greedy local search based method was specially designed for the particle position update rule by improving the greedy forwarding mechanism. Simulations demonstrated that the proposed protocol is effective and promising.