# Deployment of Wireless Sensor Networks for Oilfield Monitoring by Multiobjective Discrete Binary Particle Swarm Optimization

Authors: Zhen-Lun Yang; Angus Wu; Hua-Qing Min
Journal: Journal of Sensors
 (2016)

---

## Abstract

The deployment problem of wireless sensor networks for real time oilfield monitoring is studied. As a characteristic of oilfield monitoring system, all sensor nodes have to be installed on designated spots. For the energy efficiency, some relay nodes and sink nodes are deployed as a delivery subsystem. The major concern of the construction of the monitoring system is the optimum placement of data delivery subsystem to ensure the full connectivity of the sensor nodes while keeping the construction cost as low as possible, with least construction and maintenance complexity. Due to the complicated landform of oilfields, in general, it is rather difficult to satisfy these requirements simultaneously. The deployment problem is formulated as a constrained multiobjective optimization problem and solved through a novel scheme based on multiobjective discrete binary particle swarm optimization to produce optimal solutions from the minimum financial cost to the minimum complexity of construction and maintenance. Simulation results validated that comparing to the three existing state-of-the-art algorithms, that is, NSGA-II, JGGA, and SPEA2, the proposed scheme is superior in locating the Pareto-optimal front and maintaining the diversity of the solutions, thus providing superior candidate solutions for the design of real time monitoring systems in oilfields.


---

## Body

## 1. Introduction

In recent years, with the advent of industrial wireless sensor networks (IWSNs), the fusion of wireless communication and distributed sensing technologies and real time remote monitoring with IWSNs have been widely adopted in oil, gas, and other resources industries. As reported in [1], IWSN is the only feasible sensing solution for the oil and gas industry benefited from lowering costs and expanding the possibilities for new and future improvements. In fact, oil and gas companies are migrating to standards based wireless mesh sensor networking according to the survey in [2]. In China, at present, the domestic digital oilfield has been in rapid development. Daqing Oilfield Company is playing the leading role while Xinjiang Oilfield and Tarim Oilfield are establishing their digital oilfields. The effective use of wireless technology is essential to the success of digital oilfield management [3]. Wireless real time monitoring system provides a simpler access to real time data and insight into operations in oilfield at acceptable cost. Besides, the maintenance costs, downtime, and energy consumption will be reduced and equipment performance will be improved.To the best of our knowledge, all existing applied studies on IWSNs for oil well remote monitoring are on the system design such as the main structure, hardware component, and software design [4, 5], aiming at extending the application of IWSNs in oil and gas industry. However, the inclusion of further optimization of IWSNs in oilfield monitoring system is essential towards the efficiency of design and system operation. There is still no research focusing on such topic yet. In this paper, a systematic framework providing problem formulations and a novel optimization strategy is proposed and it can be the reference for further attempts to optimize the IWSNs in oilfield management.Typical structure of IWSN for oil well remote monitoring is shown in Figure1 regardless of the fact that it is for production wells or water injection wells. Sensors are installed on the operation components of the wells. Pressure, voltage, and stream temperature of each well are acquired and transmitted to the control center via the wireless networks. Sensors are usually powered by batteries and do not support long distance transmission since energy consumption is a superlinear function of the transmission distance. Sensor nodes are in general configured as an end-device with less communication capability only allowed communicating with relay node or sink node. To guarantee the reliability of data communication and the ease of system setup, a relay node called well relay is deployed next to the operation unit in every well to collect sensor data and forward to the sink node. In this model, sensors at each well are capsulated as a point represented by the position of the corresponding well relay. The deployment of well relays is determined by the condition and landform of each well and should be decided in advance of the placement of sink and relay nodes. In wireless networks with fixed infrastructure, it is sufficient for each node having a wireless link connected to at least one sink node acting as the gateway to collect all the data in its dominating area and deliver the data towards the remote server. The remote communication protocol for sink node can be either McWiLL, WiMAX, or 3G/4G. The connecting path of a well relay and its serving sink node may consist of multiple relays in between if the well relay is far away from the sink.Figure 1
The typical structure of IWSN for oil well monitoring.The traditional deployment issue in WSNs normally refers to the placement of sensors so as to minimize the number of deployed sensor nodes while ensuring the full coverage and connectivity of the monitoring area [6, 7]. However, contrary to generic WSNs, the node placement in IWSN for oil well monitoring is deterministic to meet the desired monitoring tasks. The deployment strategies of generic WSNs cannot be directly adopted for IWSNs for oil well monitoring due to the difference in requirements and work conditions such as critical requirement for reliability, the uniqueness of each sensor measurement [8, 9], and some application-specific features of IWSNs for oil well monitoring.With the prior information of each well relay’s position, the deployment problem in oil well monitoring mainly concerns the optimum placement of data delivery subsystem consisting of relays and sinks to ensure the full connectivity of the well relays while keeping the construction cost as low as possible and with least construction and maintenance complexity. In this paper, the deployment problem of IWSN for oil well monitoring is formulated as a constrained multiobjective optimization problem and solved with a novel multiobjective discrete binary particle swarm optimization algorithm. The major contributions of this paper are described as follows:(1)
A novel multiobjective model is provided to optimize the deployment problem of IWSN for oil well monitoring.
(2)
A discrete binary MOPSO framework with an efficient constraint handling method through the elitist seed repopulation and an infeasible solution archive is proposed.
(3)
A novel hybrid local attractor based position updating rule is introduced to improve the exploration and exploitation capabilities of the multiobjective particle swarm optimization (MOPSO) algorithm.The rest of this paper is organized as follows. Section2 discusses related work. Section 3 presents the mathematical formulation of the deployment problem including definition and objective functions. Section 4 describes the proposed approach for solving the multiobjective optimization of deployment problems as described. Section 5 illustrates the experimental results, followed by the conclusion in Section 6.


## 2. Related Work

The deployment problem of data delivery subsystem of industrial WSNs and nonindustrial WSNs has been studied well over the past years. The commonly used approach in early researches formulates the deployment problems as optimization problems in graph theory and solves them with polynomial time approximation algorithms in order to minimize the number of relays in order to meet certain connectivity or survivability requirements. Cheng et al. [10] proposed a scheme of placing the fewer amount of relay nodes in an area interconnecting all sensing nodes with the model formulated as a Steiner Minimum Tree with Minimum Number of Steiner Nodes problem. Pan et al. [11] presented the approaches to maximize the topological network lifetime for a two-tiered WSN by arranging the location of base stations and relay nodes under several definitions according to mission criticality. In [12], relay node placement was studied in a large-scale two-tiered WSN with the 2-connected network topology to provide fault-tolerance. Lloyd and Xue [13] studied the problems of minimizing the number of relay nodes by formulating the model as a Minimum Spanning Tree (MST) problem. Misra et al. [14] studied constrained version of the relay node placement problem. Some recent studies have focused on applying multiobjective optimization algorithms to relay nodes placement problems. A multiobjective evolutionary algorithm was proposed by Perez et al. [15] for the simultaneous optimization of the number of relays and the energy dissipation when deploying a WSN. Zhong and Zhang [16] presented a multiobjective memetic algorithm to optimize the relay node placement problem with consideration of the network lifetime and the number of relay nodes.Other than the researches on deployment of data delivery subsystem of WSNs, researches relevant to this study, such as the placement of base station in wireless networks, are described here. Chan et al. [17] proposed a multiobjective jumping genes paradigm for the base station placement optimization of factory WLAN network with the aim of optimization to obtain a good quality of service and good quality for communication traffics without suffering from the cost of overall system. Ting et al. [18] proposed a multiobjective variable-length genetic algorithm to search for the optimal number, types, and positions of heterogeneous base stations by considering coverage, cost, capacity, and overlap. An improved NSGA-II algorithm is proposed by Lakshminarasimman et al. [19] to determine optimal location of base station in cellular mobile communications network with consideration of coverage maximization and cost minimization. Wechtaisong et al. [20] presented a WiMAX network planning model to assign optimized amount and location of base stations in study area with the aims of minimizing installation cost and maximizing service coverage simultaneously; a weighted sum method was used to transform two objectives to single objective and the solution for the model was generated by CPLEX 5.2 optimization solver.All in all, the problem of data delivery subsystem deployment in oilfield monitoring has three features which could not be found in each of the approaches mentioned above:(1)
Energy efficiency is a key issue for the problem. The specific structure of IWSNs for oil well monitoring as illustrated in Figure1 is established to guarantee the maximum degree of energy efficiency. Based on this structural model, there is no need to consider the network lifetime or energy criterion in the optimization.
(2)
In different areas of oilfield, the distributions of oil wells are diverse; in particular, those in periphery are normally constructed in harsh environments with complicated landform. It imposes different challenges to the placement of relays and sinks.
(3)
Based on the structural model, the full coverage to the sensors is transformed to the full connection of the well relays.As in the structural model, the relay nodes and sink nodes constitute a data delivery subsystem of the IWSN; an optimized deployment of such a data delivery subsystem is economical for the construction of the network. Moreover, it is essential to ease the construction and operation maintenance of the wireless network for the generally complicated landform of oilfields and the diverse distributions of oil wells. Hence, it is necessary to minimize the cost of network while keeping the complexity of construction and maintenance of the IWSNs as simple as possible. However, these two objectives may not be reached simultaneously as they are contradictory.In this paper, a multiobjective discrete binary particle swarm optimization algorithm is studied in solving the deployment problem of IWSN for oilfield monitoring. The proposed scheme provides the nondominated tradeoff solutions so that the selective decision can be made among these solutions according to the design preferences.


## 3. Problem Description

Connectivity is the major concern in the deployment of WSNs [21]. Without loss of generality, in this paper, the communication ranges of different kinds of nodes are regarded as the same. Assume that any two nodes in the network, N
a and N
b, can directly communicate with each other if their Euclidian distance is within the communication range r
c; that is, |
N
a
N
b
|
≤
r
c. Due to complicated landform of oilfield, it may impose serious problems on the construction, maintenance, and communication reliability if there are too many relay nodes between the sink and the well relay in the link. In practice, the maximum hops in a single link should be limited to a brim, H
m
a
x. In other words, the maximum allowable number of relays between the sink and the well relay is confined to H
m
a
x
-
1.Suppose that the target area A is digitized intoL
×
W pixels with pixel size equal to 1. Although the oil wells are generally distributed around the target area dispersedly, several oil wells may be constructed closely for the requirement of operation. In this case, well relays which are able to communicate with at least another well relay can be grouped as one well relay cluster. The well relay cluster set on the target area can be defined as C
=
V
1
∪
V
2
⋯
∪
V
n with each cluster defined as V
i
=
{
c
i
1
,
c
i
2
,
…
,
c
i
l
}, where c
i
a is a well relay in the cluster, l is the number of well relays in the cluster, which varies in different clusters, and n is the number of the clusters, i
=
1,2
,
…
,
n. Each member c
i
a in the cluster V
i must satisfy the communication range requirement: (1)
c
i
a
c
i
b
≤
r
c
,
∃
c
i
b
∈
V
i
,
a
≠
b
.According to (1), the members of a cluster are interconnected with at least one neighbor; with at least one well relay connected to the sink, this will guarantee the connectivity of the whole cluster. The well relay of the cluster connected to the sink is called the entry well relay.Define the relay node set and sink node set on the target area asR and S, where R is of r
1
,
r
2
,
…
,
r
m and S is of s
1
,
s
2
,
…
,
s
k. The connectivity model of the IWSN system can be expressed as follows.A well relay clusterV
i is said to be connected to the data delivery subsystem, denoted as(2)
c
i
a
t
≤
r
c
,
∃
c
i
a
∈
V
i
,
∃
t
∈
R
∪
S
.Following the data path, where sensor data are delivered to the sink node, the model can be simplified as follows.A well relay clusterV
i is connected to the data delivery subsystem when it satisfies the following:(3)
c
i
a
t
≤
r
c
×
H
m
a
x
,
∃
c
i
a
∈
V
i
,
∃
t
∈
S
.Hence, the deployment problem can be regarded as a sink node placement problem.The event that the entry well relay of clusterV
i located at the pixel (
x
i
,
y
i
) connecting to the sink node s
j located at the pixel (
x
j
,
y
j
) is denoted as E
c
o
n
n
(
V
i
,
s
j
). It can be expressed by a two-valued function:(4)
E
c
o
n
n
V
i
,
s
j
=
1
,
i
f
x
j
-
x
i
2
+
y
j
-
y
i
2
≤
r
c
×
H
m
a
x
.
0
,
o
t
h
e
r
w
i
s
e
.If any entry well relay of cluster V
i is connected to any node in the sink node set, its connectivity is considered completed. As a result, the event that the well relay cluster V
i is connected to the sink node set can be denoted as(5)
E
c
o
n
n
V
i
,
S
=
1
-
∏
j
=
1
k
1
-
E
c
o
n
n
V
i
,
s
j
.Finally, the connectivity rate of the well relay cluster setR
c
(
C
) is defined as the ratio of the number of connected well relay clusters to the total number of the well relay clusters:(6)
R
c
C
=
∑
i
=
1
n
E
c
o
n
n
V
i
,
S
n
.The challenge of the deployment problem is to deploy the nodes of the data delivery subsystem in an optimal manner meeting the connectivity requirement and achieving a balance between the financial cost and the complexity of construction and maintenance.It can be illustrated by a simple example having six clusters of well relays as shown in Figure2.Figure 2
Scenario of solutions for the deployment problem.
(a)
(b)
(c)Figure2(a) illustrates one of the extreme solutions that the data delivery subsystem is comprised of only one sink node and multiple relay nodes. Unlike the normal relay simply having the capability of data forwarding, the sink node has higher operation capability and extra functions such as data storage and remote communication. It is much costly than a relay. Although this solution is less costly as only a single sink node is required, numerous relay nodes have to be installed on the specified locations. The complexity of construction and maintenance will be high and even not plausible due to the landform complexity. Besides, it is unreliable to have a single sink for the entire network since any failure may lead to a halt for the entire network and suspend the operation of the oil field.Figure2(b) shows the other extreme in which each cluster of the well relays is coupled with a nearby sink node. The complexity of construction and maintenance is rather low. However, it may lead to a very high cost due to the fact that a large number of sink nodes are needed.A compromised practical solution between the two extreme cases is shown in Figure2(c). The clusters of well relays can be grouped into superclusters, and each of the superclusters is connected to a single sink node. The financial cost of this solution is higher than the solution in Figure 2(a) but much lower than the solution in Figure 2(b). For the construction and maintenance complexity, on the contrary, it is higher than that of the solution in Figure 2(b) but lower than the solution in Figure 2(a). Obviously, the two criteria are contradictory such that a single optimum solution minimizing both objectives is implausible.Apparently, the cost of the data delivery subsystem is determined by the total cost of sink nodes and relay nodes of the communication link between the well cluster and the sink. It is desirable to keep the number of sinks as fewer as possible. On the other hand, more delivery nodes may lead to higher complexity for construction and maintenance. Because of the complicated landforms of oil fields and the dispersed distribution of oil wells, it is generally implausible to share the relay nodes amongst different data links. For such reason and the simplification of analysis as well, here we do not consider the sharing of relays between different links from well relays to sink nodes.Define the number of relay nodes in the links from well relay to sink node asQ
=
{
q
1
,
q
2
,
…
,
q
n
}, where n is the number of the links being the same as the number of clusters.Each element in the vector is referred to as(7)
q
i
=
c
i
e
s
j
r
c
-
1
,
i
f
c
i
e
s
j
≤
r
c
×
H
m
a
x
.
i
n
v
a
l
i
d
,
o
t
h
e
r
w
i
s
e
,where x is the ceiling function and s
j is the serving sink node for the well relay cluster V
i with the entry well relay c
i
e, i
=
1,2
,
…
,
n.As described above, the deployment problem can be formulated as a Constraint Multiobjective Optimization Problem (CMOP) for the following objectives.Objective
1. The financial cost of the data delivery subsystem is defined as(8)
Min
⁡
f
1
=
∑
i
=
1
n
q
i
×
F
r
+
M
×
F
s
,where F
r and F
s are the financial cost of relay node and sink node, respectively, and M is the number of the sink nodes in the data delivery subsystem.Objective
2. The complexity of the construction and maintenance is defined as(9)
Min
⁡
f
2
=
∑
i
=
1
n
q
i
+
Msubject to condition for the full connectivity of all well relay clusters defined as (10)
R
c
C
-
1
=
0
.To solve the deployment problem, a multiobjective particle swarm optimization is proposed. The algorithm is based on our recent work [22, 23] in solving multiobjective optimization problems.


## 4. The Proposed Algorithm

Particle swarm optimization is a population-based search algorithm inspired by observing the bird flocks and fish schools [24]. Kennedy and Eberhart [25] extended PSO to solve discrete problems in 1997. In their binary PSO (BPSO) algorithm, particles are bounded to include binary variables and the velocity serves as the probability to change a bit’s value. Afterward, various improved algorithms based on BPSO were proposed such as discrete multiphase PSO [26] and the angle modulated PSO [27].In order to provide discrete PSO in a more general framework, various approaches have been reported in the literatures. Clerc [28] defined the position of a particle as a permutation of numbers and the velocity as a set of swaps in discrete PSO (DPSO). Salman et al. [29] defined the position as a real vector and apply a space transformation technique to convert the position into its corresponding solution. Pang et al. [30] defined the position and velocity as a fuzzy matrix with the updating rules for the particles based on XOR operator. Chen et al. [31] proposed a set-based discrete PSO in which the candidate solution and velocity are defined as a crisp set and a set with possibilities, respectively.Nevertheless, most of these researches mainly concern with single objective optimization. In recent years, some attempts have been made to solve multiobjective optimization problems with DPSO. Palermo et al. [32] put forward a DPSO algorithm for multiobjective design space exploration based on an “aggregating” approach where the swarm is equally partitioned into subswarms, each of which uses a different cost-function determined by the product of the objectives combined with a set of randomly chosen exponents. Cordeiro et al. [33] adopted particles update equations in [32] and proposed an algorithm for the energy and performance optimization in IC design. Gong et al. [34] proposed a multiobjective discrete PSO algorithm based on decomposition.For the deployment problem, based on the binary coding scheme, we extend our recent works in [22, 23] with constraint handling through the elitist seed repopulation and an infeasible solution archive. In addition, a hybrid local attractor based position updating rule is employed to improve the exploration and exploitation capabilities of the algorithm. The pseudocode of the proposed elitist seed repopulation multiobjective discrete binary PSO with transposon (ES-MOBPSO-T) is given in Algorithm 1, and the details of the bolded and italic mechanisms in this algorithm will be elaborated later in this section.Algorithm 1:Pseudocode for the ES-MOBPSO-T algorithm.
Procedure of ES-MOBPSO-T
Initiate the position of each particle;
Evaluate each particle;
Store the useful individuals found into external archives;
Repeat
Elitist seed repopulation with constraint handling (Section4.2);
Select the leader from external archives (Section4.3);
Update the position of each particle (Section4.4);
Evaluate each particles;
Update pbests (Section4.5);
Update external archives (Section4.6);
Until maximum number of iterations is reached;
EndProcedure


### 4.1. Individual Representation

As the input, for each well relay, its location is defined by a two-dimensional coordinate and identified with a sequential number. The proposed deployment scheme aims at identifying the locations of the sink nodes for the system as the solution of the formulated multiobjective problem. According to the analysis in Section3, the maximum number of sink nodes is equal to the number of well relay clusters. In the proposed scheme, to serve n well relay clusters, each particle is represented by a binary bit-string a
=
a
1
,
a
2
,
…
,
a
n
∗
L as shown in Figure 3.Figure 3
Particle representation in ES-MOBPSO-T.It consists ofn tuples representing the total number of sink nodes. The length of the string depends on the size of the service area. Each tuple is the binary representation of the location of a sink node followed by a status bit describing the existence status of the sink node. If the status bit is “1,” it indicates that the sink exists; otherwise, the sink does not exist. The coordinates (
x
i
,
y
i
) of each sink node can be determined by (11)
x
i
=
∑
b
=
1
L
x
2
b
-
1
×
a
i
-
1
∗
L
+
b
y
i
=
∑
b
=
1
L
y
2
b
-
1
×
a
i
-
1
∗
L
+
L
x
+
b
,where L
x and L
y are the numbers of bit of the coordinates x
,
y and L is the bit length of a tuple.As described in Section3, a cluster is connected to the network if any member in the cluster is connected to a sink. To assign a serving sink node to a well relay cluster, the closest sink node to the well relay cluster will be selected as the serving sink provided the sink node exists according to the existence status of the sink and satisfies the multihop communication range constraint with any member of the cluster. If none of the members in a cluster having the distance to all the potential serving sink nodes satisfies the communication range constraint, the cluster is not connected to the delivery subsystem. In that case, the connectivity rate of the solution will be less than 1.


### 4.2. Elitist Seed Repopulation with Constraint Handling

As demonstrated in our previous work [22, 23], the combination of elitist seed repopulation and transposon operation improves both the exploitation and the exploration capabilities of MOPSO.The key of elitist seed repopulation procedure is the reinitialization of the swarm. Periodically, after certain normal search iteration, the whole swarm will be reinitialized with particles generated from the elitists selected from the external archives and processed by the diversification operations. The elitists are treated as the seeds to produce new subswarms for restarting a new cycle of swarm computation. Two parameters, the frequency of elitist seeding operation and the number of seeds, are used to specify the elitist seed repopulation.As described in Section3, the deployment problem is a constrained optimization problem. To handle the constraints, the commonly used approach is to make use of penalty functions. However, the penalty function approach has to deal with the difficulty of choosing the appropriate tuning parameters [35]. The proper fine-tuned penalty coefficients are capable of balancing the objective and penalty functions while inappropriate penalty coefficients will make the optimization problem intractable. In this paper, the works in [22, 23] are extended to deal with the constraints in obtaining solution through the elitist seed repopulation method. Instead of using a single external archive to maintain the nondominated feasible solutions, an extra external archive is used to keep the infeasible solutions with minor overall constraint violation. The two external archives are used in the elitist seed repopulation operation. The pseudocode of elitist seed repopulation with constraint handling is given in Algorithm 2. In the binary tournament selection, the comparison strategy similar to the method in [36] is adopted:(1)
If two solutions are feasible, the one in the less crowded area will be regarded as the better one.
(2)
If the two solutions are infeasible, the solution with smaller overall constraint violation will be regarded as the better one.
(3)
If one solution is feasible while the other is infeasible, the feasible solution will be regarded as the better one.In this paper, the diversification operations for elitist seed repopulation approach are transposon operations. A transposon is made of consecutive genes located in the randomly assigned position in each chromosome while the transposon operators are lateral movement operations that happen in one chromosome or between different ones. Generally speaking, there were two types of transposon operators [22, 37]: cut-and-paste and copy-and-paste as demonstrated in Figures 4 and 5, respectively. In PSO, a particle can be treated as a chromosome. The details of the transposon operation are shown in Algorithm 3.Algorithm 2:Pseudocode for elitist seed repopulation with constraint handling.
Procedure of Elitist seed repopulation with constraint handling
If operation frequency criterion is met
Combine the two external archive into a selection pool;
While the elitist pool is not full
Copy a solution to elitist pool based on binary tournament selection;
EndWhile
While the new swarm is not full
Randomly select two individual with indexi and j from elitist pool;
Repopulate new particles from elitist pool with transposon operation;
EndWhile
EndIf
EndProcedureAlgorithm 3:Pseudocode for the transposon operation.
Procedure of Transposon operation
If  i
=
j
For  k
=
1 to number of transposons
generate a random valuer
1 within the range of [
0,1
];
If  r
1 < jumping rate
generate a random valuer
2 within the range of [
0,1
];
If  r
2 > 0.5
perform cut and paste operation in the same individual;
Else
perform copy and paste operation in the same individual;
EndIf
EndIf
EndFor
Else
For  k
=
1 to number of transposons
generate a random valuer
1 within the range of [
0,1
];
If  r
1 < jumping rate
generate a random valuer
2 within the range of [
0,1
];
If  r
2 > 0.5
perform cut and paste operation in different individuals;
Else
perform copy and paste operation in different individuals;
EndIf
EndIf
EndFor
EndIf
Output the individual(s) after transposon operation
EndProcedureFigure 4
Cut-and-paste transposon operator. (a) Cut-and-paste in the same chromosome. (b) Cut-and-paste in different chromosomes.
(a)
(b)Figure 5
Copy-and-paste transposon operator. (a) Copy-and-paste in the same chromosome. (b) Copy-and-paste in different chromosomes.
(a)
(b)To control the transposon operation, three parameters, jumping rate, number of transposons, and length of transposon, are used.


### 4.3. Leader Selection

Similar to the method proposed in [38], the leader is randomly selected from a specified percentage of the external archive. As mentioned previously, there are two external archives which contain the feasible and infeasible superior solutions, respectively. Leader is selected from the feasible external archive sorted based on the degree of crowding if it is not empty; otherwise, it will be selected from the infeasible external archive sorted based on the overall constraint violation.


### 4.4. Position Updating Rule for Particles

From Clerc and Kennedy [39], the convergence of the PSO may be achieved if each particle converges to its local attractor p
i
=
(
p
i
,
1
,
p
i
,
2
,
…
,
p
i
,
D
). Each dimension of the local attractor is defined as (12)
p
i
,
d
t
=
φ
d
t
×
p
b
e
s
t
i
,
d
t
+
1
-
φ
d
t
×
g
b
e
s
t
d
t
,where p
b
e
s
t
i
,
d
t and g
b
e
s
t
d
t are the dth dimension of the personal best of particle i and the global best of the swarm in search iteration at t, respectively, and φ
d
(
t
) is a random number generated with the uniform probability distribution in the range of [
0,1
].To improve the exploration and exploitation capabilities of the algorithm, a hybrid position update approach is proposed. In each search cycle, the new position of each particle is generated by the transposon operation together with mutation on the position of its local attractor. The pseudocode of the proposed approach is given in Algorithm4, with transposon operation being the same as given in Algorithm 3.Algorithm 4:Pseudocode for the hybrid position update.
Procedure of Hybrid position update
Calculate the local attractor of each particle;
Set the position of each particle with the position of its local attractor;
Perform the transposon operation on the swarm;
Perform the mutation operation on each particle;
EndProcedureSince the position of a particle is binary encoded, the local attractor of each particle cannot be generated directly from (12). A different calculation method must be used to decide the local attractor of each particle. According to (12), the coordinate p
i
,
d of p
i lies between p
b
e
s
t
i
,
d and g
b
e
s
t
d; therefore, p
i distributes randomly in the hyperrectangle with p
b
e
s
t
i and g
b
e
s
t being its two diagonal ends. That is, (13)
p
i
-
p
b
e
s
t
i
≤
p
b
e
s
t
i
-
g
b
e
s
t
,
p
i
-
g
b
e
s
t
i
≤
p
b
e
s
t
i
-
g
b
e
s
t
.The discrete form of local attractorp
i can be generated through crossover operation [40]. p
i is randomly selected from the two offspring generated by crossover on p
b
e
s
t
i and g
b
e
s
t. Definitely, the position of p
i resulting from crossover satisfies (13) in terms of Hamming distance in binary space.For each particle, the local attractor is determined with its personal best and the global best through the crossover operation. After all the particles are updated with their local attractors, they are further processed with transposon and followed by mutation with one particle at a time. Apparently, the convergence is enhanced by setting the position of each particle with the position of its local attractor. Moreover, the transposon and mutation operations bring perturbation effects to the swarm and guide the particles to better positions beneficial for the promotion of solution diversity and coverage. In addition, the elitist seed repopulation mechanism will relocate the superior particles to some specified area which is beneficial to overcome the premature convergence.


### 4.5. pbest Update

The update procedure ofp
b
e
s
t of each particle is shown in Algorithm 5.Algorithm 5:Pseudocode for the pbest update.
Procedure ofpbest update
If  p
b
e
s
t
i and the new generated solution x
i both are feasible
If  x
i
≺
p
b
e
s
t
i
p
b
e
s
t
i = x
i;
Else
If the new generated solution x
i is feasible while p
b
e
s
t
i is infeasible
p
b
e
s
t
i = x
i;
Else
If  p
b
e
s
t
i and the new generated solution x
i both are infeasible
If constraint violation of x
i is smaller than constraint violation of p
b
e
s
t
i
p
b
e
s
t
i = x
i;
EndIf
EndIf
EndIf
EndIf
EndIf
EndProcedure


### 4.6. External Archives Update

New solutions will go through an acceptance test before entering the external archives set with a fixed maximal size. If the feasible external archive is fully filled with nondominated solutions, the crowding distances [36] of the solutions will be calculated and those solutions with larger crowding distance value will be kept while the others will be discarded. If the infeasible external archive is full, solutions with larger constraint violations will be discarded from the archive.


## 4.1. Individual Representation

As the input, for each well relay, its location is defined by a two-dimensional coordinate and identified with a sequential number. The proposed deployment scheme aims at identifying the locations of the sink nodes for the system as the solution of the formulated multiobjective problem. According to the analysis in Section3, the maximum number of sink nodes is equal to the number of well relay clusters. In the proposed scheme, to serve n well relay clusters, each particle is represented by a binary bit-string a
=
a
1
,
a
2
,
…
,
a
n
∗
L as shown in Figure 3.Figure 3
Particle representation in ES-MOBPSO-T.It consists ofn tuples representing the total number of sink nodes. The length of the string depends on the size of the service area. Each tuple is the binary representation of the location of a sink node followed by a status bit describing the existence status of the sink node. If the status bit is “1,” it indicates that the sink exists; otherwise, the sink does not exist. The coordinates (
x
i
,
y
i
) of each sink node can be determined by (11)
x
i
=
∑
b
=
1
L
x
2
b
-
1
×
a
i
-
1
∗
L
+
b
y
i
=
∑
b
=
1
L
y
2
b
-
1
×
a
i
-
1
∗
L
+
L
x
+
b
,where L
x and L
y are the numbers of bit of the coordinates x
,
y and L is the bit length of a tuple.As described in Section3, a cluster is connected to the network if any member in the cluster is connected to a sink. To assign a serving sink node to a well relay cluster, the closest sink node to the well relay cluster will be selected as the serving sink provided the sink node exists according to the existence status of the sink and satisfies the multihop communication range constraint with any member of the cluster. If none of the members in a cluster having the distance to all the potential serving sink nodes satisfies the communication range constraint, the cluster is not connected to the delivery subsystem. In that case, the connectivity rate of the solution will be less than 1.


## 4.2. Elitist Seed Repopulation with Constraint Handling

As demonstrated in our previous work [22, 23], the combination of elitist seed repopulation and transposon operation improves both the exploitation and the exploration capabilities of MOPSO.The key of elitist seed repopulation procedure is the reinitialization of the swarm. Periodically, after certain normal search iteration, the whole swarm will be reinitialized with particles generated from the elitists selected from the external archives and processed by the diversification operations. The elitists are treated as the seeds to produce new subswarms for restarting a new cycle of swarm computation. Two parameters, the frequency of elitist seeding operation and the number of seeds, are used to specify the elitist seed repopulation.As described in Section3, the deployment problem is a constrained optimization problem. To handle the constraints, the commonly used approach is to make use of penalty functions. However, the penalty function approach has to deal with the difficulty of choosing the appropriate tuning parameters [35]. The proper fine-tuned penalty coefficients are capable of balancing the objective and penalty functions while inappropriate penalty coefficients will make the optimization problem intractable. In this paper, the works in [22, 23] are extended to deal with the constraints in obtaining solution through the elitist seed repopulation method. Instead of using a single external archive to maintain the nondominated feasible solutions, an extra external archive is used to keep the infeasible solutions with minor overall constraint violation. The two external archives are used in the elitist seed repopulation operation. The pseudocode of elitist seed repopulation with constraint handling is given in Algorithm 2. In the binary tournament selection, the comparison strategy similar to the method in [36] is adopted:(1)
If two solutions are feasible, the one in the less crowded area will be regarded as the better one.
(2)
If the two solutions are infeasible, the solution with smaller overall constraint violation will be regarded as the better one.
(3)
If one solution is feasible while the other is infeasible, the feasible solution will be regarded as the better one.In this paper, the diversification operations for elitist seed repopulation approach are transposon operations. A transposon is made of consecutive genes located in the randomly assigned position in each chromosome while the transposon operators are lateral movement operations that happen in one chromosome or between different ones. Generally speaking, there were two types of transposon operators [22, 37]: cut-and-paste and copy-and-paste as demonstrated in Figures 4 and 5, respectively. In PSO, a particle can be treated as a chromosome. The details of the transposon operation are shown in Algorithm 3.Algorithm 2:Pseudocode for elitist seed repopulation with constraint handling.
Procedure of Elitist seed repopulation with constraint handling
If operation frequency criterion is met
Combine the two external archive into a selection pool;
While the elitist pool is not full
Copy a solution to elitist pool based on binary tournament selection;
EndWhile
While the new swarm is not full
Randomly select two individual with indexi and j from elitist pool;
Repopulate new particles from elitist pool with transposon operation;
EndWhile
EndIf
EndProcedureAlgorithm 3:Pseudocode for the transposon operation.
Procedure of Transposon operation
If  i
=
j
For  k
=
1 to number of transposons
generate a random valuer
1 within the range of [
0,1
];
If  r
1 < jumping rate
generate a random valuer
2 within the range of [
0,1
];
If  r
2 > 0.5
perform cut and paste operation in the same individual;
Else
perform copy and paste operation in the same individual;
EndIf
EndIf
EndFor
Else
For  k
=
1 to number of transposons
generate a random valuer
1 within the range of [
0,1
];
If  r
1 < jumping rate
generate a random valuer
2 within the range of [
0,1
];
If  r
2 > 0.5
perform cut and paste operation in different individuals;
Else
perform copy and paste operation in different individuals;
EndIf
EndIf
EndFor
EndIf
Output the individual(s) after transposon operation
EndProcedureFigure 4
Cut-and-paste transposon operator. (a) Cut-and-paste in the same chromosome. (b) Cut-and-paste in different chromosomes.
(a)
(b)Figure 5
Copy-and-paste transposon operator. (a) Copy-and-paste in the same chromosome. (b) Copy-and-paste in different chromosomes.
(a)
(b)To control the transposon operation, three parameters, jumping rate, number of transposons, and length of transposon, are used.


## 4.3. Leader Selection

Similar to the method proposed in [38], the leader is randomly selected from a specified percentage of the external archive. As mentioned previously, there are two external archives which contain the feasible and infeasible superior solutions, respectively. Leader is selected from the feasible external archive sorted based on the degree of crowding if it is not empty; otherwise, it will be selected from the infeasible external archive sorted based on the overall constraint violation.


## 4.4. Position Updating Rule for Particles

From Clerc and Kennedy [39], the convergence of the PSO may be achieved if each particle converges to its local attractor p
i
=
(
p
i
,
1
,
p
i
,
2
,
…
,
p
i
,
D
). Each dimension of the local attractor is defined as (12)
p
i
,
d
t
=
φ
d
t
×
p
b
e
s
t
i
,
d
t
+
1
-
φ
d
t
×
g
b
e
s
t
d
t
,where p
b
e
s
t
i
,
d
t and g
b
e
s
t
d
t are the dth dimension of the personal best of particle i and the global best of the swarm in search iteration at t, respectively, and φ
d
(
t
) is a random number generated with the uniform probability distribution in the range of [
0,1
].To improve the exploration and exploitation capabilities of the algorithm, a hybrid position update approach is proposed. In each search cycle, the new position of each particle is generated by the transposon operation together with mutation on the position of its local attractor. The pseudocode of the proposed approach is given in Algorithm4, with transposon operation being the same as given in Algorithm 3.Algorithm 4:Pseudocode for the hybrid position update.
Procedure of Hybrid position update
Calculate the local attractor of each particle;
Set the position of each particle with the position of its local attractor;
Perform the transposon operation on the swarm;
Perform the mutation operation on each particle;
EndProcedureSince the position of a particle is binary encoded, the local attractor of each particle cannot be generated directly from (12). A different calculation method must be used to decide the local attractor of each particle. According to (12), the coordinate p
i
,
d of p
i lies between p
b
e
s
t
i
,
d and g
b
e
s
t
d; therefore, p
i distributes randomly in the hyperrectangle with p
b
e
s
t
i and g
b
e
s
t being its two diagonal ends. That is, (13)
p
i
-
p
b
e
s
t
i
≤
p
b
e
s
t
i
-
g
b
e
s
t
,
p
i
-
g
b
e
s
t
i
≤
p
b
e
s
t
i
-
g
b
e
s
t
.The discrete form of local attractorp
i can be generated through crossover operation [40]. p
i is randomly selected from the two offspring generated by crossover on p
b
e
s
t
i and g
b
e
s
t. Definitely, the position of p
i resulting from crossover satisfies (13) in terms of Hamming distance in binary space.For each particle, the local attractor is determined with its personal best and the global best through the crossover operation. After all the particles are updated with their local attractors, they are further processed with transposon and followed by mutation with one particle at a time. Apparently, the convergence is enhanced by setting the position of each particle with the position of its local attractor. Moreover, the transposon and mutation operations bring perturbation effects to the swarm and guide the particles to better positions beneficial for the promotion of solution diversity and coverage. In addition, the elitist seed repopulation mechanism will relocate the superior particles to some specified area which is beneficial to overcome the premature convergence.


## 4.5. pbest Update

The update procedure ofp
b
e
s
t of each particle is shown in Algorithm 5.Algorithm 5:Pseudocode for the pbest update.
Procedure ofpbest update
If  p
b
e
s
t
i and the new generated solution x
i both are feasible
If  x
i
≺
p
b
e
s
t
i
p
b
e
s
t
i = x
i;
Else
If the new generated solution x
i is feasible while p
b
e
s
t
i is infeasible
p
b
e
s
t
i = x
i;
Else
If  p
b
e
s
t
i and the new generated solution x
i both are infeasible
If constraint violation of x
i is smaller than constraint violation of p
b
e
s
t
i
p
b
e
s
t
i = x
i;
EndIf
EndIf
EndIf
EndIf
EndIf
EndProcedure


## 4.6. External Archives Update

New solutions will go through an acceptance test before entering the external archives set with a fixed maximal size. If the feasible external archive is fully filled with nondominated solutions, the crowding distances [36] of the solutions will be calculated and those solutions with larger crowding distance value will be kept while the others will be discarded. If the infeasible external archive is full, solutions with larger constraint violations will be discarded from the archive.


## 5. Simulation Results

Simulation result for the verification of the proposed scheme is described in this section. The performance of the proposed multiobjective discrete binary particle swarm optimization is compared with three well-known multiobjective genetic algorithms, that is, NSGA-II [36], a well-perceived multiobjective optimization algorithm, JGGA which has demonstrated its superiority in solving a WLAN placement problem comparing to other well-known multiobjective optimization algorithms [17, 37], and SPEA2, a well-known algorithm in solving multiobjective optimization problems [41]. For the simulation experiment, 50 independent runs with random initial populations were conducted for each of the algorithms. All algorithms are implemented in MATLAB and all the simulations are conducted under the computing environment with a Lenovo notebook PC, 4 GB RAM, Core i3 2.13 GHz CPU clock speed, Microsoft Windows 7, and MATLAB 2010a.


### 5.1. Test Problems

Without loss of generality, all parameters represented in numerical value are without the measurement unit, and distances are measured in the unit of grid points for the simulations. All test problems are designed based on the reference of the oil wells distribution of Daqing and its periphery oilfields [42, 43]; the well densities and the scales of target areas were changed to construct different scenes for the empirical study. A typical distribution instance referred to in this study is shown in Figure 6. This area is a part of Y322 block at Yushulin oilfield in Daqing, China.Figure 6
A distribution instance of oil wells at Yushulin oilfield.Test Problem 1 (TP1). There are 100 well relays distributed in a 4096
∗
4096 target area. The total cluster number is calculated with (1) and turns out to be 25. As described in Section 4.1, there is a potential sink node for each cluster. Thus, the number of corresponding potential sink nodes is 25. According to the scale of the target area, the bit lengths of both x and y coordinates of each sink node’s location are 12. The maximum hop in each link is 4, and the communication range is 200. The financial costs of sink node and relay node are 8 and 1, respectively.Test Problem 2 (TP2). The well density is double as in test problem 1. There are 200 well relays clustered to 57 clusters with 57 potential sink nodes being distributed in a 4096
∗
4096 target area with other settings remaining the same as in test problem 1. This test problem is to test the search abilities of the algorithms in a more densely distributed oilfield.Test Problem 3 (TP3). The density is the same as in test problem 1 while the target area is four times larger. There are 400 well relays clustered to 124 clusters with 124 potential sink nodes being distributed in an 8192
∗
8192 target area. Accordingly, the bit lengths of both x and y coordinates of each sink node’s location are 13 with other settings remaining the same as in test problem 1. This test problem is to test the search abilities of the algorithms in a larger-scale oilfield.The deployment test problems were conducted with the proposed scheme. The potential solutions represented by the particles were reached through the searching process with the proposed ES-MOBPSO-T algorithm. Sink assignment for the potential solutions was carried out through the procedures as described in Section4 under the condition as described in Section 3. The Pareto-optimal front was located at the end of the searching process. The same procedures were performed on the compared algorithms for the test problems. It is worth noting that, following the idea in jMetal which is a famous framework for multiobjective optimization [44], the constraint handling method used in the experiment for SPEA2 is the same as in NSGA-II.


### 5.2. Comparison Metrics

According to the discussion in [45–47], no single unary metric can give a comprehensive measure on the performance of an MOEA. In order to evaluate the performance of the proposed deployment scheme comprehensively, two widely used unary metrics together with one binary indicator are adopted in this study. They are the maximum spread (MS) [46], the hypervolume (HV) indicator [48], and the binary ε-indicator [49].Maximum spread (MS) is a measure of how well the true Pareto-optimal front is covered by the generated Pareto-optimal front through hyperboxes formed by the extreme function values observed in the true Pareto-optimal front and the generated Pareto-optimal. The larger MS reflects a better coverage of the efficient frontier by the generated Pareto-optimal front. Since the true Pareto-optimal front of the problem is unknown, for each test instance, the maximum and minimum values in each of the optimization objectives are determined among all the Pareto-optimal fronts generated by all the algorithms in the simulation, and they are regarded as the extreme values of the true Pareto-optimal front in calculating the MS measure for each obtained solution.Hypervolume (HV) indicator measures the spread of the solutions along the Pareto fronts and the closeness of the solutions to the Pareto-optimal front. Since the calculation of HV can be performed without the information of the true Pareto front, it is suitable for evaluating the performance of the proposed ES-MOBPSO-T and the compared algorithms for the deployment problem. The higher the HV, the better the obtained solution set approximates the Pareto-optimal front and maintains the diversity. Based on the analysis of all of the obtained solutions, the reference points used for evaluating the HV of the obtained solutions in TP1, TP2, and TP3 are set at (173, 68), (244, 95), and (633, 222), respectively.Binaryε-indicator is a measure to compare the performance of two algorithms. Based on the values of I
ε, there are three cases of compared results of two algorithms: better, worse, and incomparable. For each algorithm, 50 simulation runs will produce 50 approximation sets. Hence, according to the definition of I
ε, there are 2500 pairs of comparison between the proposed algorithm and each of the compared algorithms.


### 5.3. Parameter Settings

To make a fair comparison, the positions of well relays were maintained the same for all algorithms and the stopping criteria were set at 600,000, the number of objective function evaluations for all algorithms in the simulation.For the proposed algorithm, an initial population of 50 was used and the sizes of both of the external archives are set to 100. Single-point crossover was used to generate the local attractor of each particle. Bitwise mutation operator was used in the mutation to each particle’s position and the probabilities for mutation were set as pm =1
/
l, where l is the binary string length. The jumping rate, number of transposons, and length of transposon were set at 0.1, 10, and 25 bits, respectively, in the position update procedure. In elitist seed repopulation, the jumping rate, number of transposons, and length of transposon were set at 1, 10, and 25 bits, respectively. The frequency of elitist seeding operation was set at 5 and the number of seeds was set at 4. An initial population of size of 100 was used for the compared algorithms. Single-point crossover and bitwise mutation operators were used in NSGA-II, JGGA, and SPEA2 with the probabilities for crossover and mutation set at 1 and 1
/
l, respectively. According to [37], the jumping rate, number of transposons, and length of transposon were 0.04, 1, and 3 bits, respectively, for JGGA.


### 5.4. Results and Analysis

For the MS and HV metrics results of the simulation experiment, the median (x
~) and interquartile range (IQR) of each test problem were recorded as the measure of location (or central tendency) and statistical dispersion. Since all the algorithms in the simulation experiment are stochastic, statistical analysis is essential to provide the confidential comparison. For MS and HV metrics on each test problem, the statistical analysis [50, 51] as shown in Figure 7 was conducted on the simulation results of the proposed algorithm and the algorithm with the best median value among the three compared algorithms. The Kolmogorov-Smirnov test is used to check whether the data sets of simulation result follow the Gaussian distribution or not. If the distribution is not normal, the nonparametric Kruskal-Wallis test is performed to compare the median results of the competing algorithms. For Gaussian distributions, the Levene test is used to check the homogeneity of the variances. And then if the variances are not equal, Welch test is invoked; otherwise, ANOVA test is performed to verify the confidence of comparisons. The results are shown in the last column of Tables 1 and 2, with the symbol “+” representing the performance difference between the proposed algorithm and the best among the three compared algorithms that is significant at the 95% confidence level, while the symbol “−” indicates that the difference is not significant. The simulation results in MS and HV metrics are shown in Tables 1 and 2, in which the best results for each problem are in bold. Table 3 shows the binary ε-indicator between ES-MOBPSO-T and the compared algorithms. TP1, TP2, and TP3 refer to the test problems 1, 2, and 3, respectively. The columns with titles “B,” “W,” and “I” in Table 3 refer to the number of times the ES-MOBPSO-T performs better than, worse than, and incomparable to the corresponding compared algorithms.Table 1
Median and interquartile range of the MS metric.
ES-MOBPSO-T
x
~ (IQR)
NSGA-II
x
~ (IQR)
JGGA
x
~ (IQR)
SPEA2
x
~ (IQR)
TP1
6.639
E
-
01
(
1.645
E
-
01
)
5
.
104
E
-
01
(
1
.
427
E
-
01
)
5.031
E
-
01
(
2.276
E
-
01
)
5.122
E
-
01
(
2.847
E
-
01
)
+
TP2
8.489
E
-
01
(
8.147
E
-
02
)
6.651
E
-
01
(
9
.
496
E
-
02
)
7.417
E
-
01
(
1.345
E
-
01
)
7
.
024
E
-
01
(
1
.
212
E
-
01
)
+
TP3
8.822
E
-
01
(
9.752
E
-
02
)
5
.
638
E
-
01
(
5
.
137
E
-
02
)
5.302
E
-
01
(
9.017
E
-
02
)
4.961
E
-
01
(
5
.
125
E
-
02
)
+Table 2
Median and interquartile range of the HV metric.
ES-MOBPSO-T
x
~ (IQR)
NSGA-II
x
~ (IQR)
JGGA
x
~ (IQR)
SPEA2
x
~ (IQR)
TP1
4.138
E
+
03
(
8.850
E
+
01
)
3.866
E
+
03
(
1.813
E
+
02
)
3.832
E
+
03
(
1.730
E
+
02
)
2.957
E
+
03
(
2.630
E
+
02
)
+
TP2
4.866
E
+
03
(
2.013
E
+
02
)
4.384
E
+
03
(
3.420
E
+
02
)
4.431
E
+
03
(
1.883
E
+
02
)
4.447
E
+
03
(
1.410
E
+
02
)
+
TP3
2.127
E
+
04
(
1.403
E
+
03
)
1.800
E
+
04
(
1.145
E
+
03
)
1.863
E
+
04
(
1.093
E
+
03
)
1.619
E
+
04
(
2.540
E
+
02
)
+Table 3
Comparison value on binaryε-indicator.
NSGA-II
JGGA
SPEA2
B
W
I
B
W
I
B
W
I
TP1
1912
24
564
1892
12
596
1755
42
703
TP2
1901
50
549
1853
135
512
1923
16
561
TP3
2100
0
400
1963
47
490
2023
0
477Figure 7
The general structure of the statistical analysis in this work.Table1 illustrates that, comparing to the other algorithms, the ES-MOBPSO-T is better in providing well-extended solutions for all test problems. On the other hand, Table 3 clearly reveals that the ES-MOBPSO-T algorithm exhibits very competitive performance to the compared algorithms in locating the Pareto-optimal front. In addition, from Table 2, we can see that the proposed ES-MOBPSO-T algorithm obtained the best values for the HV metric in all of the three test problems. It is evident that the ES-MOBPSO-T is competitive in both locating the Pareto-optimal front and maintaining the diversity.The Pareto front approximations obtained in the simulation run with the median MS value or with the value closest to the median MS value of ES-MOBPSO-T and the compared algorithms were shown in Figures8, 9, and 10. It should be noted that because solutions obtained in a single run may have the same objective values, there are some overlapped solutions in the target space. The algorithm with better capability to discover new and unknown areas in the search space will generate more scattered solutions. It is undoubted that ES-MOBPSO-T outperforms the compared algorithms in terms of solution coverage and diversity.Figure 8
The approximated Pareto-optimal front for test problem 1.Figure 9
The approximated Pareto-optimal front for test problem 2.Figure 10
The approximated Pareto-optimal front for test problem 3.In test problem 1, ES-MOBPSO-T obtains more candidate solutions than the compared algorithms, with solutions locating in the more competitive positions in the target space. In the test instances with larger density of wells or larger scale, as in the cases of test problem 2 and test problem 3, ES-MOBPSO-T outperforms the compared algorithms with solutions having more competitive objective function values and is able to obtain a much more well-extended Pareto-optimal front which almost totally covers the solutions generated by the compared algorithms.In the practical implementation of IWSNs for oil well monitoring system, one of the solutions obtained by ES-MOBPSO-T can be chosen and used to establish the fundamental positions of the relay nodes and sink nodes in data delivery subsystem based on the design preferences and the landform of oilfields.


## 5.1. Test Problems

Without loss of generality, all parameters represented in numerical value are without the measurement unit, and distances are measured in the unit of grid points for the simulations. All test problems are designed based on the reference of the oil wells distribution of Daqing and its periphery oilfields [42, 43]; the well densities and the scales of target areas were changed to construct different scenes for the empirical study. A typical distribution instance referred to in this study is shown in Figure 6. This area is a part of Y322 block at Yushulin oilfield in Daqing, China.Figure 6
A distribution instance of oil wells at Yushulin oilfield.Test Problem 1 (TP1). There are 100 well relays distributed in a 4096
∗
4096 target area. The total cluster number is calculated with (1) and turns out to be 25. As described in Section 4.1, there is a potential sink node for each cluster. Thus, the number of corresponding potential sink nodes is 25. According to the scale of the target area, the bit lengths of both x and y coordinates of each sink node’s location are 12. The maximum hop in each link is 4, and the communication range is 200. The financial costs of sink node and relay node are 8 and 1, respectively.Test Problem 2 (TP2). The well density is double as in test problem 1. There are 200 well relays clustered to 57 clusters with 57 potential sink nodes being distributed in a 4096
∗
4096 target area with other settings remaining the same as in test problem 1. This test problem is to test the search abilities of the algorithms in a more densely distributed oilfield.Test Problem 3 (TP3). The density is the same as in test problem 1 while the target area is four times larger. There are 400 well relays clustered to 124 clusters with 124 potential sink nodes being distributed in an 8192
∗
8192 target area. Accordingly, the bit lengths of both x and y coordinates of each sink node’s location are 13 with other settings remaining the same as in test problem 1. This test problem is to test the search abilities of the algorithms in a larger-scale oilfield.The deployment test problems were conducted with the proposed scheme. The potential solutions represented by the particles were reached through the searching process with the proposed ES-MOBPSO-T algorithm. Sink assignment for the potential solutions was carried out through the procedures as described in Section4 under the condition as described in Section 3. The Pareto-optimal front was located at the end of the searching process. The same procedures were performed on the compared algorithms for the test problems. It is worth noting that, following the idea in jMetal which is a famous framework for multiobjective optimization [44], the constraint handling method used in the experiment for SPEA2 is the same as in NSGA-II.


## 5.2. Comparison Metrics

According to the discussion in [45–47], no single unary metric can give a comprehensive measure on the performance of an MOEA. In order to evaluate the performance of the proposed deployment scheme comprehensively, two widely used unary metrics together with one binary indicator are adopted in this study. They are the maximum spread (MS) [46], the hypervolume (HV) indicator [48], and the binary ε-indicator [49].Maximum spread (MS) is a measure of how well the true Pareto-optimal front is covered by the generated Pareto-optimal front through hyperboxes formed by the extreme function values observed in the true Pareto-optimal front and the generated Pareto-optimal. The larger MS reflects a better coverage of the efficient frontier by the generated Pareto-optimal front. Since the true Pareto-optimal front of the problem is unknown, for each test instance, the maximum and minimum values in each of the optimization objectives are determined among all the Pareto-optimal fronts generated by all the algorithms in the simulation, and they are regarded as the extreme values of the true Pareto-optimal front in calculating the MS measure for each obtained solution.Hypervolume (HV) indicator measures the spread of the solutions along the Pareto fronts and the closeness of the solutions to the Pareto-optimal front. Since the calculation of HV can be performed without the information of the true Pareto front, it is suitable for evaluating the performance of the proposed ES-MOBPSO-T and the compared algorithms for the deployment problem. The higher the HV, the better the obtained solution set approximates the Pareto-optimal front and maintains the diversity. Based on the analysis of all of the obtained solutions, the reference points used for evaluating the HV of the obtained solutions in TP1, TP2, and TP3 are set at (173, 68), (244, 95), and (633, 222), respectively.Binaryε-indicator is a measure to compare the performance of two algorithms. Based on the values of I
ε, there are three cases of compared results of two algorithms: better, worse, and incomparable. For each algorithm, 50 simulation runs will produce 50 approximation sets. Hence, according to the definition of I
ε, there are 2500 pairs of comparison between the proposed algorithm and each of the compared algorithms.


## 5.3. Parameter Settings

To make a fair comparison, the positions of well relays were maintained the same for all algorithms and the stopping criteria were set at 600,000, the number of objective function evaluations for all algorithms in the simulation.For the proposed algorithm, an initial population of 50 was used and the sizes of both of the external archives are set to 100. Single-point crossover was used to generate the local attractor of each particle. Bitwise mutation operator was used in the mutation to each particle’s position and the probabilities for mutation were set as pm =1
/
l, where l is the binary string length. The jumping rate, number of transposons, and length of transposon were set at 0.1, 10, and 25 bits, respectively, in the position update procedure. In elitist seed repopulation, the jumping rate, number of transposons, and length of transposon were set at 1, 10, and 25 bits, respectively. The frequency of elitist seeding operation was set at 5 and the number of seeds was set at 4. An initial population of size of 100 was used for the compared algorithms. Single-point crossover and bitwise mutation operators were used in NSGA-II, JGGA, and SPEA2 with the probabilities for crossover and mutation set at 1 and 1
/
l, respectively. According to [37], the jumping rate, number of transposons, and length of transposon were 0.04, 1, and 3 bits, respectively, for JGGA.


## 5.4. Results and Analysis

For the MS and HV metrics results of the simulation experiment, the median (x
~) and interquartile range (IQR) of each test problem were recorded as the measure of location (or central tendency) and statistical dispersion. Since all the algorithms in the simulation experiment are stochastic, statistical analysis is essential to provide the confidential comparison. For MS and HV metrics on each test problem, the statistical analysis [50, 51] as shown in Figure 7 was conducted on the simulation results of the proposed algorithm and the algorithm with the best median value among the three compared algorithms. The Kolmogorov-Smirnov test is used to check whether the data sets of simulation result follow the Gaussian distribution or not. If the distribution is not normal, the nonparametric Kruskal-Wallis test is performed to compare the median results of the competing algorithms. For Gaussian distributions, the Levene test is used to check the homogeneity of the variances. And then if the variances are not equal, Welch test is invoked; otherwise, ANOVA test is performed to verify the confidence of comparisons. The results are shown in the last column of Tables 1 and 2, with the symbol “+” representing the performance difference between the proposed algorithm and the best among the three compared algorithms that is significant at the 95% confidence level, while the symbol “−” indicates that the difference is not significant. The simulation results in MS and HV metrics are shown in Tables 1 and 2, in which the best results for each problem are in bold. Table 3 shows the binary ε-indicator between ES-MOBPSO-T and the compared algorithms. TP1, TP2, and TP3 refer to the test problems 1, 2, and 3, respectively. The columns with titles “B,” “W,” and “I” in Table 3 refer to the number of times the ES-MOBPSO-T performs better than, worse than, and incomparable to the corresponding compared algorithms.Table 1
Median and interquartile range of the MS metric.
ES-MOBPSO-T
x
~ (IQR)
NSGA-II
x
~ (IQR)
JGGA
x
~ (IQR)
SPEA2
x
~ (IQR)
TP1
6.639
E
-
01
(
1.645
E
-
01
)
5
.
104
E
-
01
(
1
.
427
E
-
01
)
5.031
E
-
01
(
2.276
E
-
01
)
5.122
E
-
01
(
2.847
E
-
01
)
+
TP2
8.489
E
-
01
(
8.147
E
-
02
)
6.651
E
-
01
(
9
.
496
E
-
02
)
7.417
E
-
01
(
1.345
E
-
01
)
7
.
024
E
-
01
(
1
.
212
E
-
01
)
+
TP3
8.822
E
-
01
(
9.752
E
-
02
)
5
.
638
E
-
01
(
5
.
137
E
-
02
)
5.302
E
-
01
(
9.017
E
-
02
)
4.961
E
-
01
(
5
.
125
E
-
02
)
+Table 2
Median and interquartile range of the HV metric.
ES-MOBPSO-T
x
~ (IQR)
NSGA-II
x
~ (IQR)
JGGA
x
~ (IQR)
SPEA2
x
~ (IQR)
TP1
4.138
E
+
03
(
8.850
E
+
01
)
3.866
E
+
03
(
1.813
E
+
02
)
3.832
E
+
03
(
1.730
E
+
02
)
2.957
E
+
03
(
2.630
E
+
02
)
+
TP2
4.866
E
+
03
(
2.013
E
+
02
)
4.384
E
+
03
(
3.420
E
+
02
)
4.431
E
+
03
(
1.883
E
+
02
)
4.447
E
+
03
(
1.410
E
+
02
)
+
TP3
2.127
E
+
04
(
1.403
E
+
03
)
1.800
E
+
04
(
1.145
E
+
03
)
1.863
E
+
04
(
1.093
E
+
03
)
1.619
E
+
04
(
2.540
E
+
02
)
+Table 3
Comparison value on binaryε-indicator.
NSGA-II
JGGA
SPEA2
B
W
I
B
W
I
B
W
I
TP1
1912
24
564
1892
12
596
1755
42
703
TP2
1901
50
549
1853
135
512
1923
16
561
TP3
2100
0
400
1963
47
490
2023
0
477Figure 7
The general structure of the statistical analysis in this work.Table1 illustrates that, comparing to the other algorithms, the ES-MOBPSO-T is better in providing well-extended solutions for all test problems. On the other hand, Table 3 clearly reveals that the ES-MOBPSO-T algorithm exhibits very competitive performance to the compared algorithms in locating the Pareto-optimal front. In addition, from Table 2, we can see that the proposed ES-MOBPSO-T algorithm obtained the best values for the HV metric in all of the three test problems. It is evident that the ES-MOBPSO-T is competitive in both locating the Pareto-optimal front and maintaining the diversity.The Pareto front approximations obtained in the simulation run with the median MS value or with the value closest to the median MS value of ES-MOBPSO-T and the compared algorithms were shown in Figures8, 9, and 10. It should be noted that because solutions obtained in a single run may have the same objective values, there are some overlapped solutions in the target space. The algorithm with better capability to discover new and unknown areas in the search space will generate more scattered solutions. It is undoubted that ES-MOBPSO-T outperforms the compared algorithms in terms of solution coverage and diversity.Figure 8
The approximated Pareto-optimal front for test problem 1.Figure 9
The approximated Pareto-optimal front for test problem 2.Figure 10
The approximated Pareto-optimal front for test problem 3.In test problem 1, ES-MOBPSO-T obtains more candidate solutions than the compared algorithms, with solutions locating in the more competitive positions in the target space. In the test instances with larger density of wells or larger scale, as in the cases of test problem 2 and test problem 3, ES-MOBPSO-T outperforms the compared algorithms with solutions having more competitive objective function values and is able to obtain a much more well-extended Pareto-optimal front which almost totally covers the solutions generated by the compared algorithms.In the practical implementation of IWSNs for oil well monitoring system, one of the solutions obtained by ES-MOBPSO-T can be chosen and used to establish the fundamental positions of the relay nodes and sink nodes in data delivery subsystem based on the design preferences and the landform of oilfields.


## 6. Conclusion

A deployment scheme of industrial wireless sensor networks for oil well monitoring based on MOPSO was proposed in this paper. With the knowledge of well relay locations, an optimum sink node location and relay nodes number were designed to ensure the well relays' connectivity by considering the financial cost and the complexity of construction and maintenance. The simulation results show that the proposed ES-MOBPSO-T can produce diverse optimal solutions from the minimum financial cost to the minimum complexity of construction and maintenance and achieve a better performance than the other three state-of-the-art algorithms.Future research will be focused on more empirical studies and the implementation of the oil well monitoring system based on this work.



---

*Source: 9358358-2016-04-24.xml*