A column generation heuristic for optimal wireless sensor network design with mobile sinks

ARTICLE IN PRESS

JID: EOR

[m5G;December 23, 2016;8:34]

European Journal of Operational Research 0 0 0 (2016) 1–14

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Innovative Applications of O.R.

A column generation heuristic for optimal wireless sensor network design with mobile sinks Muhammed Emre Keskin Department of Industrial Engineering, Atatürk University, Erzurum, Turkey

a r t i c l e

i n f o

Article history: Received 27 April 2015 Accepted 2 December 2016 Available online xxx Keywords: Computing science Integer programing Column generation Lagrangian relaxation OR in telecommunications

a b s t r a c t Wireless Sensor Networks (WSNs) consist of a high number of tiny, multi-functional, electronic devices called sensors. They collectively provide a distributed environment that is capable of monitoring remote areas. Collected information is transmitted in a direct or multi-hop fashion to the gateway nodes called sinks. An even distribution of energy loads among the sensors is critical for elongating network lifetime. There are four main WSN design issues that substantially affect the distribution of the energy: locations of the sensors, schedule of the active and standby periods of the sensors, trajectory of the mobile sink(s) and routes for data flows. As a result, many studies try to make energy usage more efficient by optimal determination of these design issues. However, only a few of them provide a unified frame in which all four design issues are integrated. In this work, we follow this line of research and propose a column generation heuristic for a Mixed Integer Linear Programing (MILP) model that integrates all design issues. Based on the extensive numerical experiments, we can say that the heuristic outperforms its competitors in the literature. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Wireless Sensor Networks (WSNs) consist of high number of tiny, multi-functional, electronic devices called sensors which are deployed over a region of interest called the sensor field. A sensor can collect data from the neighboring area lying within its sensing range and the collected information is transmitted in a direct or multi-hop fashion to the gateway nodes called sinks. Communication between sensors or between a sensor and a sink is possible only if the receiver sensor or the sink lays within the communication range of the transmitter sensor. The collaborative effort of the sensors provides a distributed environment that is capable of monitoring remote areas. This is why WSNs have a wide range of application (Yick, Mukherjee, & Ghosal, 2008). Sensors can be in active or standby modes. An active sensor performs sensing, data transmission and data receiving duties and consumes energy. On the other hand, a standby sensor consumes negligible energy. Low energy capacities of the sensors is a prominent limitation of WSNs. As there are huge number of sensors in a typical WSN application, replacing depleted batteries is generally considered out of option. Moreover, an even distribution of energy loads among the sensors is critical for elongating network lifetime. Therefore, a

Fax: +90 442 231 6015. E-mail addresses: [email protected], [email protected]

precise management of energy affecting factors is critical in order to obtain network lifetimes which are long enough. There are four main WSN design issues that substantially affect the distribution of the energy among the sensors: locations of the sensors, schedule of the active and standby periods of the sensors, locations or trajectory of the sinks and routes for data flows. Locations and numbers of the sensors should be determined so that the coverage requirements of the sensor field are satisfied. This problem is called as the Coverage Problem (CP). It is generally assumed in WSN design studies that coverage requirements of the field are characterized by a discrete set of points called coverage points. As some parts of the sensor field can be more critical than the others, requirements of the coverage points can be differentiated implying a heterogeneous CP. Another concern is that some of the deployed sensors can be taken into standby mode in some periods. Hence, there should be enough number of sensors in order to provide some sort of flexibility for the activity schedules of the sensors. In addition, total sensor deployment cost cannot exceed the allocated budget. Next important WSN design issue is to schedule the activities of the deployed sensors which is called Activity Scheduling Problem (ASP). Sensors having relatively lower residual energies can be put into standby mode while some standby sensors are made active at the same time. Scheduling the set of active sensors results in a balanced distribution of the sensors’ residual energies. However, it should be noted that the coverage requirements of the field should be satisfied throughout

http://dx.doi.org/10.1016/j.ejor.2016.12.006 0377-2217/© 2016 Elsevier B.V. All rights reserved.

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

JID: EOR 2

ARTICLE IN PRESS

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

the network lifetime. Hence, the set of active sensors at any period should be able to meet the coverage requirements. Third energy affecting WSN design issue is the deployment of the sinks which is called Sink Location Problem (SLP). Locations of the sinks play an important role on the energy load distribution. Another form of this problem is named Sink Routing Problem (SRP) if sinks are mobile. It is a known phenomenon that sensors near to the sinks, which are called relay sensors, spend more energy than others. This is due to the fact that the collected data of the whole network are transmitted to the sinks through the relay sensors. This may cause the depletion of the relay sensor batteries at relatively early phases of the network lifetime. After the death of the relay sensors, sinks become disconnected from the network implying that data collection duty is over. This problem is called differently in several sources such as “crowded center effect” (Popa, Rostamizadeh, Karp, Papadimitriou, & Stoica, 2007), “energy hole problem” (Li & Mohapatra, 2007; Wu, Chen, & Das, 2008), and “sink neighborhood problem” (Basagni, Carosi, Melachrinoudis, Petrioli, & Wang, 2008). Changing the set of relay sensors by controlled sink mobility is offered as a remedy of this phenomenon. Hence, controlled mobility of the sinks plays a regulatory role on the energy distribution. Most of the mobile sink studies from the literature assume that the sink(s) has limitless energy and it instantaneously jumps from one point to another. Hence, sink travel times are usually taken as zero in the literature and the data collected during the sink travel times are also neglected. One of the rare studies that puts some limitations on the mobility of the sink is due to Liang, Luo, and Xu (2010) which considers the sink as an energy limited device that is powered by petrol or electricity. This implies that the sink has to go to the petrol stations or to the electric charge stations periodically to renew its energy. Keskin, Altınel, Aras, and Ersoy (2011) extends the idea of the limited sink mobility by assuming nonzero sink travel time for a single mobile sink. Sink travel time is considered as a part of the network lifetime and the data accumulated during the sink travel time is also taken into account. Authors manage to show that considering nonzero sink travel time produces networks with significantly longer lifetimes if the sink is to repeat its tour many times and the network field is relatively large. The case with multiple mobile sinks is analyzed by Keskin (2014) and it is shown in the study that considering nonzero sink travel times is important for the situations where the sink speeds are slower than 1 kilometer per hour, i.e., for the networks under extreme conditions. Hence, there is no harm in assuming zero sink travel times for WSN applications with multiple mobile sinks if they travel faster than 1 kilometer per hour. We also assume in this study that sink movements take negligible amount of time. In addition, the relocation costs of the sinks are also taken as zero. It should be noted that sinks travel from one point to another only between the periods and they stand on their locations within the periods implying that each sink travels a limited distance. Therefore, the relocation cost that occurs during sink travels is negligible compared to the revenue gained by having a network with longer lifetime. Final important WSN design issue affecting energy loads of the sensors is due to the routes for the data flows. The determination of the most energy efficient route is called Data Routing Problem (DRP). Although DRP is easy for given locations of the sinks and active sensors, problem complexity increases substantially when it is integrated with other problems. It should be noted that an energy affecting WSN design issue affects the quality of the other design issues. For instance, it is not possible to obtain a good activity schedule of the sensors if sensor deployment issue is handled miserably. Similarly, controlled sink mobility may not be able to provide balanced energy loads if the activity schedule of the sensors is poor. Finally, one cannot cre-

ate good data flow strategies for a bad set of active sensors and bad sink mobility structures. Therefore, these WSN design issues should be handled together within a unified frame implying integration of CP, ASP, SLP or SRP with DRP. However, most of the studies concentrate only on a subset of these problems and assume that the results of the other problems are ready at hand a priori. This approach produces suboptimal solutions since the optimality of the assumed results is not guaranteed. For instance, Altınel, Aras, Güney, and Ersoy (2008) concentrate only on the CP. Similarly, Wang, Basagni, Melachrinoudis, and Petrioli (2005), Basagni et al. (2008), Basagni, Carosi, Petrioli, and Phillips (2011), and Keskin et al. (2011) try to solve SRP for given sensor location, activity schedules and data flow protocol. In all these studies, activity schedules are taken exogenously and deployed sensors are kept active throughout the network lifetime. Only the studies that deliberately focus on ASP provide varying schedules. There are quite a number of studies which integrate SRP and DRP for given sensor locations and activity schedules. Some of the earliest studies of this kind are due to Gandham, Dawande, Prakash, and Venkatesan (2003), Azad and Chockalingam (2006) and Alsalih, Akl, and Hassanein (2007) in which energy efficient sink and data routes are sought. They all divide the time into equal-length periods but handle each period independently. This deficiency is resolved by Luo and Hubaux (2005) in which a Mixed Integer Linear Programing (MILP) model for optimum sink and data routes is found for multiple periods simultaneously. Xia and Shihada (2015) try to jointly minimize the energy consumption of the energy-critical sensors and the data transmission delay throughout the network. These five papers provide mathematical programing models but they concentrate on energy usage characteristics of the sensors rather than direct maximization of the WSN lifetime. On the contrary, Linear Programing (LP) model of Papadimitriou and Georgiadis (2005) combines DRP with SRP and maximizes the lifetime in the objective function. Gatzianas and Georgiadis (2008) revisit the model of Papadimitriou and Georgiadis (2005) to provide a distributed solution strategy which makes use of Lagrangian decomposition that is first offered by Madan and Lall (2006) for a model with static sinks. Yun and Xia (2010) extend the model of Papadimitriou and Georgiadis (2005) into two new models so that delay tolerant applications are also handled within. Yun, Xia, Behdani, and Smith (2010) and Behdani, Yun, Cole Smith, and Xia (2012) come up with decomposition strategies for one of the models of Yun and Xia (2010) which provide data and sink routes depending only on the local sensor characteristics. Finally, Güney, Aras, Altınel, and Ersoy (2010) extend the model of Papadimitriou and Georgiadis (2005) for multiple but static sinks while Luo and Hubaux (2010) provide a model including multiple mobile sinks. An interesting underwater application of WSNs combining SRP with DRP is due to Basagni et al. (2014) in which an autonomous underwater vehicle visits the sensors to collect their data. The data produced by each sensor is associated with a value indicating the importance of the knowledge captured within the data and this value is assumed to be decreasing by time. Therefore, the analysis is more suited for the event driven applications in which immediate intervention is required during event occurrences. The authors try to find a route for the autonomous underwater vehicle that maximizes the value of the collected information. They also offer a heuristic that delivers more than 80% of the theoretical maximum value of the information. Our analysis is different since we assume that each sensor is able to send its data to the sink either directly (if the sink is nearby) or through other sensors (if the sink is far away) immediately after producing it. Moreover, we also assume that each sensor periodically generates a scalar data. Maximization of the lifetime implies the maximization of the collected data in such a framework.

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

JID: EOR

ARTICLE IN PRESS

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

Studies that combine more than two WSN design issues are rare. One of the studies of this kind is due to Güney, Aras, Altınel, and Ersoy (2012) in which the model of Güney et al. (2010) is integrated with CP. Hence, optimal sensor deployment is also a concern in addition to sink location and data routing issues. This is the first attempt to the integration of CP, SLP and DRP. Authors propose a nested solution method in which sensor locations satisfying budget and coverage requirements are spanned by a Tabu search algorithm outside while best sink locations and data routes for given sensor locations are sought inside. Türkog˘ ulları, Aras, and Altınel (2010a) and Türkog˘ ulları, Aras, Altınel, and Ersoy (2010b); 2010c) extend the model of Güney et al. (2012) by integrating ASP on top of CP, SLP and DRP. These studies share the same model but provide different solution strategies. Türkog˘ ulları et al. (2010a) develop a Lagrangian heuristic first and then propose a two stage method in which locations and activity schedules of the sensors are determined at the first stage whereas the sink locations and data flows are calculated for given locations and activity schedules of the sensors at the second stage. Türkog˘ ulları et al. (2010b) first construct disjoint connected active sensor sets satisfying budget and coverage requirements and then try to find the best possible combination of them leading to the maximum lifetime. Finally, Türkog˘ ulları, Aras, Altınel, and Ersoy (2010c) implement a column generation algorithm which makes use of a reformulation of the original MILP model. The study by Castaño, Bourreau, Velasco, Rossi, and Sevaux (2015) also implements the column generation to find the feasible active set of sensors. They try to maximize the network lifetime by a MILP model combining CP, ASP and DRP. On the contrary, ASP is handled within a different context in the study by Lersteau, Rossi, and Sevaux (2016). The authors try to find the activity schedules of the sensors for the WSNs constructed for target tracking. Set of active sensors should be able to track the targets under uncertainty in such a framework. Besides, model of Keskin, Altınel, Aras, and Ersoy (2014) combines CP, ASP and SRP with DRP. They propose two solution strategies which they call Period Iteration Heuristic (PIH) and Sequential Assignment Heuristic (SAH). They both depend on the idea of reducing the number of binary variables of the MILP model by fixing a subset of them and solving the rest accordingly. PIH initially assumes that the number of periods is one, and keeps increasing it by one until no improvement in the network lifetime is observed. SAH also depends on the idea of period iteration but instead of using the solver at each iteration, it employs a sequential assignment procedure. Both methods have relatively simplistic structures as they depend on some intuitional rule of thumb procedures. Although it is shown in the paper that PIH and SAH have better performances than that of the commercial solver Gurobi (Gurobi Optimization, 2016), it is possible to obtain better WSN designs by making use of more theoretical decomposition strategies like column generation and Lagrangian relaxation. Alternatively, Keskin, Altınel, and Aras (2015) introduce a hybrid heuristic which combines simulated annealing with Lagrangian heuristic for a model with static sinks. Finally, Keskin, Altınel, Aras, and Ersoy (2016) provide four MILP models to analyze the effects of integration of the WSN design issues on the lifetime by gradually increasing the level of integration. The rest of the paper is organized as follows. In the next section, the MILP model which integrates CP, ASP, SRP and DRP is given as a reminder. Next, Lagrangian heuristic solution strategy is explained. We expose the success of the method in numerical results section. Finally, we conclude the paper and designate possible future research directions in the last section. 2. Mathematical model In this section, we first describe the sets, parameters, and variables that are used in the formulation. Later, the formulation of the

3

MILP model is provided. At the end, we provide valid constraints to strengthen the formulation. 2.1. Sets, parameters and variables We assume that the coverage requirements of the field are characterized by a set of discrete coverage points, which we call K. The set of candidate sensor locations on which deployment of the sensors is possible and the set of sink visit points which are potential sink locations are called S and N , respectively. We assume that there are more than one type of sensors and the set of sensor types is called R. A sensor of type r ∈ R deployed at candidate sensor location i ∈ S is called sensor (i, r) in the sequel. Data collected by sensor (i, r) can be transmitted directly to one of the sinks if there is a sink in the transmission range of sensor (i, r), or indirectly via other sensors lying in the transmission range. The set of potential sensor locations and sink visit points lying in the transmission range of sensor (i, r) are respectively named as Sir and Nir . In a similar manner, the set of coverage points lying in the sensing range of sensor (i, r) is represented by Kir . Finally, the set of periods is given by T . Parameters mostly depend on the above described sets. For instance, dk stands for the amount of coverage requirement of the coverage point k ∈ K. Hence, the number of active sensors that are capable of sensing point k should be at least as much as dk at any time throughout the network lifetime. er represents the initial battery energy of a type r sensor while hr is its data generation rate. cr and cs are the energy segments spent respectively for receiving a bit of data and sensing per unit time. Similarly, cti j represents the energy used by a sensor located at point i ∈ S for transmission of a single bit of information to the sensor or sink located at point j. fir is the cost of locating sensor (i, r). Finally, P represents the total number of actual sinks which are to be located at each period and B denotes the total allocated budget for sensor deployment. There are also big M parameters used in the formulation. The selection of their values is explained in Section 2.3.1. There are six set of decision variables of the model. The continuous variables are wt , xirjst and yirt . wt simply represents the length of period t. It should be noted that the period lengths are not predetermined as constants such as 12 hours or 24 hours as done in some previous studies (Türkog˘ ulları et al., 2010a; Türkog˘ ulları et al., 2010b; 2010c). In contrary, period lengths are optimally determined by the model which drastically reduces the number of decision variables and significantly increases the efficiency as shown in Keskin et al. (2014); Keskin et al. (2016). xirjst denotes the amount of data transmitted from sensor (i, r) to sensor (j, s) in period t while yirt stands for the amount of transmitted data from sensor (i, r) to a sink located at point  in period t. The binary variables of the model are named as pir , qirt and zt . pir indicates whether or not sensor (i, r) is deployed. It takes value 1 if sensor (i, r) is deployed and it is zero otherwise. Similarly, qirt represents the activity of sensor (i, r) in period t. It is set to 1 if sensor (i, r) is active in period t, and to zero otherwise. Finally, zt indicates whether or not a sink is located at point  in period t. Hence, it is equal to 1 if a sink is located at  in period t, and it is zero otherwise. A summary of the set, parameter and variable definitions is provided in Table 1 for the sake of convenience. 2.2. Mathematical model MILP model named CASD that integrates CP, ASP, SRP and DRP is given below.

CASD: max



wt

(1)

t∈T

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR 4

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

Table 1 Sets, parameters, and decision variables. Sets

Definition

K S N R T Kir Sir

Set of coverage points Set of potential sensor locations Set of potential sink locations Set of sensor types Set of time periods Set of coverage points in the sensing range of sensor (i, r) Set of potential sensor locations in the transmission range of sensor (i, r) Set of potential sink locations in the transmission range of sensor (i, r)

Nir Parameters B cr cs

Definition Available budget for sensor deployment Energy spent by sensors per bit of data received Sensing and processing energy spent by active sensors per unit time Energy spent by sensors per bit of data transmitted from location i to j Requirement of coverage point k Battery energy of a type-r sensor Deployment cost for sensor (i, r) Amount of data generated by an active type-r sensor per unit time Number of sinks to be located

cti j dk er fir hr P

Decision var. Definition Binary variable indicating whether sensor (i, r) is located or not pir Binary variable indicating whether sensor (i, r) is active or not in qirt period t Length of period t wt Amount of data flow from sensor (i, r) to sensor ( j, s) in period t xirjst Amount of data flow from sensor (i, r) to a sink located at point yirt  in period t Binary variable indicating whether there is a sink at point  in zt period t or not

s.t.  

x jsirt + hr wt qirt =

s∈R j:i∈S js



yirt +

∈Nir



xir jst

s∈R j∈Sir

i ∈ S, r ∈ R, t ∈ T 

 cs wt qirt + cr

(2)  

x jsirt +

s∈R j:i∈S js

t∈T



cti j xir jst +

s∈R j∈Sir

i ∈ S, r ∈ R  

1 yirt ≤ Mt zt

 ∈ N,t ∈ T



 cti yirt ≤ er

∈Nir

Lifetime is defined as the sum of the period lengths in the objective function (1) and it is maximized. Constraint (2) ensures that amount of data that comes to sensor (i, r) from neighboring sensors together with the data generated by sensor (i, r) itself is equal to amount of data sent out of sensor (i, r) to neighboring sensors or sinks for each period t ∈ T . Hence, constraint (2) stands for data flow balance for each sensor for each period. In the left hand side of constraint (3), we sum up over each period the energy segments used by sensor (i, r) respectively for sensing and processing data, receiving data from neighboring sensors and transmitting data to neighboring sensors and to neighboring sinks. The total energy used by sensor (i, r) is then forced to be less than or equal to the battery energy of the sensor. Next constraint (4) simply avoids data transmission to a sink visit point if there is no sink located at that point. Constraint (5) makes it sure that all P sinks are located in each period. We ensure satisfaction of coverage requirements for each k ∈ K by constraint (6). They simply state that the number of active sensors that are able to cover point k must be at least as much as its coverage demand dk for each period. Constraint (7) puts a budget restriction for the total deployment cost. Constraint (8) imposes that sensor (i, r) can be active in period t only if it is deployed at the first hand. Constraints (9) and (10) respectively state that there can be no outflow from and no inflow to sensor (i, r) in period t unless it is active in that period. Finally, constraints (11) and (12) stand for the usual nonnegativity and binary restrictions. As can be seen, the length of active time of sensor (i, r) is taken into consideration for calculation of data generated by sensor (i, r) in period t in constraint (2) and for calculation of energy used by sensor (i, r) for sensing and processing duties in constraint (3). This active time is simply taken as wt qirt since it will be zero if sensor (i, r) is in standby mode in period t, i.e., qirt = 0, and it will be equal to the length of period t if sensor (i, r) is active in period t, i.e., qirt = 1. However, by definition, wt qirt is the multiplication of a binary variable qirt and a continuous variable wt which introduces a nonlinearity to the model. In order to linearize the multiplication terms, we define a new continuous variable airt which is responsible for the active time of sensor (i, r) in period t, i.e., airt = wt qirt . We also introduce four new constraint sets, which are linear and force airt to be equal to wt qirt . These constraints are given in the following.

(3)

airt ≤ wt

i ∈ S, r ∈ R, t ∈ T

(4)

4 airt ≤ Mirt qirt

i ∈ S, r ∈ R, t ∈ T

(13) (14)

r∈R i:∈Nir



zt = P

t∈T

(5)

∈N

 

qirt ≥ dk

k ∈ K, t ∈ T

(6)

r∈R i:k∈Kir



fir pir ≤ B

(7)

i∈S r∈R

qirt ≤ pir 

i ∈ S, r ∈ R, t ∈ T

2 xir jst ≤ Mirt qirt

i ∈ S, r ∈ R, t ∈ T

(8) (9)

s∈R j∈Sir

 

x jsirt ≤

3 Mirt qirt

i ∈ S, r ∈ R, t ∈ T

(10)

i, j ∈ S, r, s ∈ R, n ∈ N , t ∈ T

(11)

 ∈ N , i ∈ S, r ∈ R, t ∈ T

(12)

s∈R j:i∈S js

wt , yirt , xir jst ≥ 0 zt , pir , qirt ∈ {0, 1}

4 airt ≥ wt − Mirt (1 − qirt )

airt ≥ 0

i ∈ S, r ∈ R, t ∈ T

i ∈ S, r ∈ R, t ∈ T

(15) (16)

As can be understood, if qirt takes value 1, then airt is equated to wt by means of constraints (13) and (15). Similarly, if qirt is zero, airt is set to zero by means of constraints (14) and (16). This implies that nonlinear multiplication wt qirt can be replaced with airt after adding constraints (13)–(16) to the model. 2.3. Strengthening the formulation Improving the upper bounds obtained from the linear relaxation of CASD is important in order to use it as a cut off point in an attempt to eliminate some nodes from the branch and bound tree. This can be achieved by determining the values of big M’s as small as possible. Similarly, some constraints can be made tighter or some additional valid constraints can be added to the formulation. Selection of big M values and efforts to have stronger and valid constraints are respectively explained in the following sections.

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

2.3.1. Value of big M’s 1, In order to strengthen the formulation CASD, the values of Mt 2 3 4 Mirt , Mirt and Mirt should be chosen appropriately so that these constraints are as tight as possible. In our formulation, all five constraints would be tighter as lower values are assigned to the big M’s since all constraints but (15) are in the form of less than or 4 is negative in constraint (15). equal to while the coefficient of Mirt Therefore, the lowest possible assignments lead to the strongest CASD formulation providing lowest linear relaxation upper bounds. 4 puts an upper bound on the amount of active time First of all, Mirt of sensor (i, r) in period t. An active sensor spends cs amount of energy per unit time for sensing and processing data. Moreover, assuming that sensor (i, r) does not receive any data from the neighboring sensors and transmits only its own data, the amount of energy it spends for data transmission per unit time is at least hr cti j

while j is the closest possible sensor or sink location to point i. Hence, since the minimum per unit time energy used by sensor (i, r) is cs + hr cti j , the amount of active time of sensor (i, r) cannot exceed

er cs +hr ct  ij

. We calculate this upper bound values for each

4 = sensor and assign Mirt

er cs +hr ct  ij

4 to calculate . We make use of Mirt

1 which limits the amount of total data transmisthe value of Mt sion to sink visit point  in period t. Total data transmitted by sen4 since M 4 is an upper sor (i, r) at any period cannot exceed hr Mirt irt bound for the active time of sensor for each period. It should be noted that although it is unlikely, all sensors can be active in a given period in theory. Moreover, hypothetically, all active sensors can send its data to a specific sink visit point. Hence, the upper bound for the total transmitted data to a sink visit point  is set   1 = 4 to Mt i∈S r∈R hr Mirt . A similar selection process is pursued 2 3 for Mirt and Mirt as well. It can be the case in theory that all sensors are active and send their data through sensor (i, r). Therefore,   2 = M3 = 4 Mirt j∈S s∈R hs M jst follows. For notational convenience, irt 4 and M to denote M 1 , M2 and M 3 in the we let M to denote Mirt t irt irt sequel.

5

s ∈ R, t ∈ T and similarly constraint (10) would be cr x jsirt ≤ er qirt i ∈ S, r ∈ R, j ∈ S, s ∈ R, t ∈ T . Observe that we get rid of big M’s of constraint (9) and constraint (10) as well. 4. A new constraint wt+1 ≤ wt t ∈ T \{|T |} can be added to the formulation in order to eliminate the identical solutions which are generated due to the symmetry. For instance, there may be solutions which are the same except values of the variables corresponding to specific two periods are interchanged. Such solutions can be eliminated by imposing a nondecreasing structure on the period lengths. 5. Lengths of some periods are found to be zero for some instances. Aggregation of such periods may help to eliminate the identical solutions. A new binary variable bt indicating whether or not period t has a positive length can be added to the formulation. Next, new set of constraints wt ≤ Mbt t ∈ T and bt+1 ≤ bt t ∈ T \{|T |} can be incorporated within formulation CASD. 6. A new small integer model can be formed to find the set of sensors covering the requirements of the field with minimum possible cost. We give its formulation below.

z1 = min



fir pir

(17)

k∈K

(18)

i ∈ S, r ∈ R

(19)

i∈S r∈R

s.t.  

pir ≥ dk

r∈R i:k∈Kir

pir ∈ {0, 1}

Here, z1 represents the optimal objective value of the model. After running the solver for a given time, a lower bound for z1 , say z1 , can be easily found. Note that the total cost of active sensors in each period must be larger than z1 . Hence, the con  straint i∈S r∈R fir qirt ≥ z1 t ∈ T is valid for CASD.  7. Let j represent the closest possible sensor or sink location to point i. Sensor (i, r) spends at least cs + cti j amount of energy per unit time during which it is active implying that the total amount of active time of sensor (i, r) cannot exceed s er t . c +c  ij

2.3.2. Stronger inequalities It is possible to strengthen some of the existing constraints of CASD. Moreover, some set of valid inequalities can be generated and added to the formulation. As a result, an improvement in the linear relaxation upper bounds at each node of the branch and bound tree can be observed. Therefore, it may be possible to cut off more nodes of the branch and bound tree as a result of the improved upper bounds. In addition, the time spent to find good feasible solutions would also be lowered. However, adding more constraints and in some cases more variables may negatively affect the time spent for the solution of the linear relaxation of the model in each node. Therefore, the combined effect of stronger constraints and/or valid constraints on the final feasible solution quality is ambiguous. One should test and observe the effect of them on a relatively small set of test instances before applying them. We itemize some of the possible strengthening efforts below. 1. Right hand side of constraint (3) can be replaced with er pir . It is easy to see that a sensor that is not located cannot consume any energy and hence the effective energy of sensor (i, r) is indeed er pir . 2. It should be possible to strengthen the linear relaxation by disaggregating constraint (4). For instance, constraint (4) could be rewritten as cti yirt ≤ er zt i ∈ S, r ∈ R,  ∈ N , t ∈ T . Note that by this transformation we get rid of the big M’s existing on the right hand side of the constraint. 3. A similar disaggregation logic can be imposed upon constraint (9) and constraint (10) as well. As a result, the new form of constraint (9) would be cti j xir jst ≤ er qirt i ∈ S, r ∈ R, j ∈ S,

Moreover, the active time of sensor (i, r) is practically zero  if it is not deployed at all. Therefore, constraint t∈T airt ≤ er t pir i ∈ S, r ∈ R can be added to CASD without cutting off s c +c  ij

any feasible point. 8. A new mathematical model minimizing the number of sensors covering the requirements of the field can be formed. This model is similar to the one given in item 6, but instead of the cost of the sensors, we minimize the total number of sensors in the objective function. Let z2 be a lower bound for the optimal objective value of the model found by running the solver for a     given time. Constraints i∈S r∈R pir ≥ z2 and i∈S r∈R qirt ≥ z2 t ∈ T are valid for CASD and can be added to the formulation. 9. Observe that the budget constraint (7) is a knapsack constraint. Hence, it is possible to find the cover inequalities that are violated by the optimal solution of the linear relaxation of CASD. Violated cover inequalities can be lifted by the standard lifting procedure and the lifted inequalities can be added to the formulation. Moreover, this procedure can be repeated until all the cover inequalities are satisfied by the optimal linear relaxation solution. We refer to each one of these items by their number, i.e., item 1, item 2, etc., in the sequel. 3. Solution method Computational time that the exact solution of CASD requires can be prohibitively large for even small and moderate size

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR 6

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

instances. One reason behind this difficulty is the large number of binary decision variables included in the formulation of CASD. This observation suggests referring to a heuristic method that decomposes the formulation into easy to solve smaller subproblems. We propose a Lagrangian heuristic for the base CASD formulation given in Section 2.2 for that purpose.



 cs airt + cr

 

x jsirt

s∈R j:i∈S js

t∈T

+



cti j xir jst

s∈R j∈Sir

airt ≤ wt





+

cti yirt

≤ er i ∈ S, r ∈ R

(23)

∈Nir

i ∈ S, r ∈ R, t ∈ T

(24)

3.1. Lagrangian relaxation scheme We relax constraints (4), (6), (9), (10), (14) and (15) and carry them to the objective function after multiplying with Lagrange multipliers θ , μ, λ, β, γ 1 and γ 2 . We let CASD(θ , μ, λ, β, γ 1 , γ 2 ) denote the remaining Lagrangian subproblem and give its formulation below.

CASD(θ , μ, λ, β, γ 1 , γ 2 ): max



wt +



 μkt

k∈K t∈T

+



θt M zt −

∈N t∈T

t∈T

+





 



M q

−

irt



+

+

yirt

irt

−



 

s.t. 





x jsirt



(20)

γ 1,

Note that the objective function of CASD(θ , μ, λ, β, is a function of original CASD variables w, a, x, y, p, z and q as well as a function of Lagrange multipliers θ , μ, λ, β, γ 1 and γ 2 . We let the optimal objective value of CASD(θ , μ, λ, β, γ 1 , γ 2 ) for given multipliers be denoted by Z (θ , μ, λ, β, γ 1 , γ 2 ). Lagrangian subproblem CASD(θ , μ, λ, β, γ 1 , γ 2 ) can be decomposed into three independent subproblems. The first subproblem, which is denoted by CASD1 (θ , λ, β, γ 1 , γ 2 ), consists of only the continuous decision variables and hence it is a Linear Program (LP). We also let Z 1 (θ , λ, β, γ 1 , γ 2 ) denote the optimal objective value of CASD1 (θ , λ, β, γ 1 , γ 2 ). The formulation of CASD1 (θ , λ, β, γ 1 , γ 2 ) is given in the following.

max

1−

t∈T

−





i∈S r∈R t∈T

+ βirt

 



γirt2 wt −

λirt

 ∈N t∈T

i∈S r∈R



θt

 

yirt

r∈R i:∈Nir



1 2 x jsirt + (γirt − γirt )airt

x jsirt + hr airt =

+

 s∈R j∈Sir

(29)

(21)

γ 2)

(30)

s.t. 

fir pir ≤ B

(31)

i∈S r∈R

pir ∈ {0, 1}

i ∈ S, r ∈ R 

max{0, M (λ

(32)

2 ) − γirt



where αir = t∈T + k∈K irt + βirt ) + ir μkt }. We make use of COMBO solver of Martello, Pisinger, and Toth (1999) for the solution of the binary knapsack problems. We observe in our numerical examples that the COMBO solver is very efficient so that it is able to solve all instances including the large ones in very short amount of times. Finally, the third subproblem contains only the binary variable set z related to the sink relocation. Parts of (20) containing z variables form the objective function of the third subproblem while the constraints of CASD(θ , μ, λ, β, γ 1 , γ 2 ) containing sink relocation variable consist of the constraints of the third subproblem. We name it as CASD3 (θ ) while Z 3 (θ ) denotes its optimal objective value for given θ values. Its formulation is provided below.

CASD3 (θ ):  max M θt zt

1 M (γirt

(33)

∈N t∈T

zt = P

t∈T

(34)

∈N

yirt

zt ∈ {0, 1}

∈Nir

xir jst i ∈ S, r ∈ R, t ∈ T

γ 1,

i∈S r∈R

s.t. 

s∈R j:i∈S js

i ∈ S, r ∈ R, t ∈ T .

CASD2 (μ, λ, β, γ 1 , γ 2 ):  max αir pir

γ 2)

xir jst

s∈R j∈Sir



(28)

A careful look to the formulation of λ, β, reveals that it can be transformed into a pure binary knapsack problem as

s∈R j:i∈S js

s.t.  

(27)

CASD2 (μ,

s.t. (2 ), (3 ), (5 ), (7 ), (8 ), (11 ) − (13 ), (16 ).



(26)

i ∈ S, r ∈ R, t ∈ T

pir , qirt ∈ {0, 1}

γirt1 (Mqirt − airt ) + γirt2 (airt − wt + M (1 − qirt ))





1 2 M λirt + M βirt + Mγirt − Mγirt qirt

fir pir ≤ B

qirt ≤ pir

xir jst

i∈S r∈R t∈T



r∈R i:k∈Kir

i∈S r∈R

s∈R j:i∈S js

CASD1 (θ , λ, β, γ 1 , γ 2 ):



i∈S r∈R t∈T

s∈R j∈Sir

i∈S r∈R t∈T

M q



(25)

CASD2 (μ, λ, β, γ 1 , γ 2 ):    max μkt qirt k∈K t∈T

qirt − dk



λirt



i, j ∈ S, r, s ∈ R,  ∈ N , t ∈ T ,

The second subproblem and its optimal objective value are respectively denoted by CASD2 (μ, λ, β, γ 1 , γ 2 ) and Z 2 (μ, λ, β, γ 1 , γ 2 ). We provide its formulation below.

r∈R i:∈Nir

r∈R i:k∈Kir



+ βirt

 

wt , airt , xir jst , yirt ≥ 0

(22)

 ∈ N,t ∈ T

(35)

CASD3 (θ )

It is easy to see that can be further decomposed into T smaller subproblems for each period t ∈ T . The formulation of

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

7

subproblem 3 for a given t is named as CASDt3 (θ ) and can be given as

Zt3 (θ ) ≤ u3 (t ) t ∈ T

(42)

CASDt3 (θ ):  max M θt zt

Z 4 ( μ, γ 2 ) ≤ u 4

(43)

u1 , u2 , u3 (t ), u4 ∈  t ∈ T

(44)

(36)

∈N

s.t. 

zt = P for t only

(37)

θ , μ, λ, β, γ 1 , γ 2 ≥ 0.



∈N

zt ∈ {0, 1}

∈N

(38)

Smaller subproblems CASDt3 (θ ) for each t ∈ T are easily solved by inspection for given θ values. If we let Zt3 (θ ) denote the optimal objective value of CASDt3 (θ ), then we can  write Z 3 (θ ) = t∈T Zt3 (θ ). Moreover, we let Z 4 (μ, γ 2 ) represent the constant part of the objective function, which is equal to      2 − M i∈S r∈R t∈T γirt k∈K t∈T μkt dk . Then, the optimal value of CASD(θ , μ, λ, β, γ 1 , γ 2 ) can be found by summing up the optimal values of the subproblems CASD1 (θ , λ, β, γ 1 , γ 2 ), CASD2 (μ, λ, β, γ 1 , γ 2 ) and CASD3 (θ ) with the constant part. Therefore, Z (θ , μ, λ, β, γ 1 , γ 2 ) = Z 1 (θ , λ, β, γ 1 , γ 2 )+Z 2 (μ, λ, β, γ 1 , γ 2 )  + t∈T Zt3 (θ ) + Z 4 (μ, γ 2 ). Note that all fractions of Z (θ , μ, λ, β, 1 γ , γ 2 ) can be calculated easily for given multiplier values. As a final remark, one can easily see that Z (θ , μ, λ, β, γ 1 , γ 2 ) always forms an upper bound for the optimum value of CASD for any given values of the multipliers θ , μ, λ, β, γ 1 and γ 2 . The optimal Lagrangian multipliers leading to the minimum possible Lagrangian upper bound is found by the Dantzig–Wolfe column generation method as described in the next section. We believe that subproblem solutions obtained at each iteration of the algorithm span a considerable part of the solution space. Therefore, subproblem solutions have the potential to produce good feasible solutions for the original CASD formulation. We computationally verify this idea in numerical results section. In other words, we exploit the column generation, that is originally employed to find the Lagrangian dual as described in the next section, as a tool to produce feasible solutions from the subproblem solutions and we show that this approach overperforms the heuristic solution procedures from the literature that are previously employed for the solution of CASD.

Finally, we define wth1 , airth1 , xir jsth1 , yirth1

(45)

h1 ∈H1

as the set of

extreme points of the feasible region of CASD1 (θ , λ, β, γ 1 , γ 2 ) defined by constraints (22)–(25). Here,  H1 is the index set labeling these extreme points. We also let pirh2 , qirth2 be the set of h2 ∈H2

2 points satisfying  (27)–(29) (i.e., the feasible region of CASD (μ, λ, β, γ 1 , γ 2 )), and zth3 (t ) be the set of points satisfying h3 (t )∈H3 (t )

(37) and (38) (i.e., the feasible region of CASDt3 (θ )) for each t ∈ T , where H2 and H3 (t ) respectively denote their index sets. The LD problem can now be written as

LD : ZLD = min u1 + u2 +



u3 (t ) + u4

(46)

t∈T

s.t. 

 1−



 γ

2 irt

wth1 −

−



 λirt

+βirt



 

yirth1

r∈R i:∈Nir

xir jsth1

s∈R j∈Sir

i∈S r∈R t∈T

 

θt

∈N t∈T

i∈S r∈R

t∈T



1 2 x jsirth1 + (γirt − γirt )airth1

 ≤ u1

h1 ∈ H 1

(47)

s∈R j:i∈S js



αir pirh2 ≤ u2 h2 ∈ H2

(48)

i∈S r∈R



M θt zth3 (t ) ≤ u3 (t ) t ∈ T , h3 (t ) ∈ H3 (t )

∈N

M

 i∈S r∈R t∈T

γirt2 −



μkt dk ≤ u4

(49) (50)

k∈K t∈T

3.2. Solving the Lagrangian dual problem

u1 , u2 , u3 (t ), u4 ∈  t ∈ T

(51)

The Lagrangian Dual (LD) problem searches for the optimal multiplier values leading to the minimum possible Z (θ , μ, λ, β, γ 1 , γ 2 ) value. Therefore, LD problem can be written as

θ , μ, λ, β, γ 1 , γ 2 ≥ 0.

(52)

LD:

min

θ ,μ,λ,β,γ 1 ,γ 2 ≥0

Z ( θ , μ , λ, β , γ 1 , γ 2 )

which is equivalent to

LD: +

min

θ ,μ,λ,β,γ 1 ,γ 2 ≥0



Z 1 (θ , λ, β, γ 1 , γ 2 ) + Z 2 (μ, λ, β, γ 1 , γ 2 )

Zt3 (θ ) + Z 4 (μ, γ 2 ).

t∈T

Defining u1 , u2 , u3 (t) for each t ∈ T and u4 as upper bounds respectively for Z 1 (θ , λ, β, γ 1 , γ 2 ), Z 2 (μ, λ, β, γ 1 , γ 2 ), Zt3 (θ ) for each t ∈ T and Z 4 (μ, γ 2 ), we have the following form of the LD problem:

LD : min u1 + u2 +



t∈T

u3 (t ) + u4

(39)

s.t. Z 1 ( θ , λ, β , γ 1 , γ 2 ) ≤ u 1 Z 2 ( μ , λ, β , γ 1 , γ 2 ) ≤ u 2

(40) (41)

Here, ZLD denotes its optimum value. There can be an enormous number of members in sets H1 , H2 and H3 (t ) implying that the number of constraints of type (47)–(49) can be prohibitively large. Moreover, it can be impossible to comprise all solutions from H1 , H2 and H3 (t ) for each t ∈ T in a reasonable amount of computation time. Therefore, instead of taking all constraints of type (47)– (49), one can start LD with a subset of them. Hence, a mathematical optimization formulation that is of the same form with LD but with less constraints can be formed, which is named as Restricted Lagrangian Dual Problem (RLD). It should be noted that the optimum value of RLD, say ZRLD , is a lower bound for ZLD since RLD is a relaxation of LD. On the other hand, Z (θ , μ, λ, β, γ 1 , γ 2 ) is an upper bound for ZLD for any value of the multipliers θ , μ, λ, β, γ 1 , γ 2 . Hence, lower and upper bounds can be generated in this manner for the minimum possible Lagrangian dual value, ZLD . Subgradient algorithm and multiplier adjustment method (Fisher, 1981) are two commonly employed ways used to find optimum Lagrange multiplier values leading to the minimum possible Lagrangian dual value. However, we observed in our numerical experiments that both of them are ineffective for CASD instances with respect to the quality of the feasible solutions. In

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

JID: EOR 8

ARTICLE IN PRESS

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

Algorithm 1 Dantzig–Wolfe column generation algorithm. Initiate multipliers θ , μ, λ, β, γ 1 , γ 2 , initiate  Let ZRLD = −∞, ZBEST = ∞, m = 1, n = 1, τ = 1/ min{50, 0.5 × (m + n )} Solve subproblems to obtain Z (θ , μ, λ, β, γ 1 , γ 2 ) Let θ BEST , μBEST , λBEST , βBEST , γ 1BEST , γ 2BEST be respectively equal to θ , μ, λ, β, γ 1 , γ 2 and ZBEST = Z (θ , μ, λ, β, γ 1 , γ 2 ) while (ZRLD +  < ZBEST ) AND (time limit is not exceeded) do Add new constraints to RLD in the form of (47)–(49) Solve RLD to obtain θ , μ, λ, β, γ 1 , γ 2 Update the multipliers as follows; θ = τ × θ + (1 − τ ) × θ BEST , μ = τ × μ + (1 − τ ) × μBEST , λ = τ × λ + (1 − τ ) × λBEST , β = τ × β + (1 − τ ) × βBEST , γ 1 = τ × γ 1 + (1 − τ ) × γ 1BEST , γ 2 = τ × γ 2 + (1 − τ ) × γ 2BEST Solve subproblems to obtain Z (θ , μ, λ, β, γ 1 , γ 2 ) if (ZBEST > Z (θ , μ, λ, β, γ 1 , γ 2 )) then ZBEST = Z (θ , μ, λ, β, γ 1 , γ 2 ) Let θ BEST , μBEST , λBEST , βBEST , γ 1BEST , γ 2BEST be respectively equal to θ , μ, λ, β, γ 1 , γ 2 n=n+1 end if τ = 1/ min{50, 0.5 × (m + n )} m←m+1 end while

other words, although both methods efficiently converge to the optimal Lagrangian dual multipliers, the feasible solutions obtained using subproblem solutions are relatively bad. Hence, we adopt the so-called Dantzig–Wolfe column generation method that is known to asymptotically converge to the Lagrangian dual. Applying column generation for the purpose of solving the Lagrangian dual problem is initially introduced by Wentges (1997) and successfully employed by Klose and Görtz (2007) for the capacitated facility location problem and by Keskin et al. (2015) for the Wireless Sensor Network design with static sinks. Study of Keskin et al. (2015) can be considered similar to the analysis provided by this study to some extent. However, our model assumes that the sinks are mobile. As indicated in Keskin et al., 2016), mobility of the sinks have significant effect on the lifetime of the networks especially for the large instances. For example, networks provided by our model have around 25% larger lifetimes than the lifetimes of the networks provided by the models of Keskin et al. (2015) for the instances with 200 candidate sensor locations on the average (the comparison is provided in the study by Keskin et al., 2016). Improvement made by considering sink mobility is even greater for larger networks. In addition, incorporating the mobility of the sinks into the formulation is achieved easily. Another important point is that the sink locations are determined by a simulated annealing heuristic in the study by Keskin et al. (2015) while the sink movement decisions are decided within the Lagrangian scheme in our analysis. Since sink movement has considerable effect on the lifetime, determining the sink movement decisions in a more accurate manner is equally important. We also report the results for the case in which the sink movements decisions are generated by the simulated annealing heuristic as it is done in the study of Keskin et al. (2015) in numerical results section for comparison purposes. It can be seen that determining the sink movement decisions within the Lagrangian heuristic framework produces networks with much larger lifetimes. Besides, Keskin et al. (2015) do not try to improve their MILP models by stronger and/or valid constraints as we do in this study. Finally, solutions coming from two subproblems are combined and a single constraint is added to the master problem in the work by Keskin et al. (2015) while we add separate constraints to the master problem for the solution of each subproblem. Hence, our master problem is stronger due to this disaggregation. Each iteration of the Dantzig–Wolfe column generation algorithm is comprised of two phases. In the first phase, CASD1 (θ , λ, β,

γ 1 , γ 2 ), CASD2 (μ, λ, β, γ 1 , γ 2 ) and CASD3 (θ ) are solved and successive constraints in form of (47)–(49) corresponding to the optimal solutions of the subproblems are added to RLD. Then, RLD is solved to optimality to find the current Lagrange multiplier values. Instead of directly using them, we compute the convex combination of the current multipliers with the best Lagrange multiplier values (those leading to the best upper bounds) and pass them to the Lagrangian subproblems. Note that the lower bound values, ZRLD , never decreases throughout the iterations since the optimum value of a minimization problem does not decrease by the increase in the number of constraints. Besides, the upper bounds, Z (θ , μ, λ, β, γ 1 , γ 2 ), fluctuate in general and the best upper bound is the one with the lowest value, which is named as ZBEST . Lagrange multipliers leading to ZBEST are recorded as the best multipliers. Weights of the current and best multipliers are chosen as τ = 1/ min{50, 0.5 × (m + n )} and 1 − τ (as suggested by Wentges, 1997) where m is the number of current iteration of the Dantzig– Wolfe decomposition algorithm and n stands for the number of improvements in the value of ZBEST . Convergence of the algorithm is ensured as the weights of the best multipliers is increased with the number of iterations and the number of updates in the best multiplier values. However, to keep the effect of the current multipliers above a threshold, we choose 1/50 as the coefficient of the current multipliers if 0.5 × (m + n ) becomes more than 50. Initial multiplier values are chosen as the optimal dual values of the corresponding constraints of the linear relaxation of CASD. Iterations of the algorithm are continued until the difference between ZBEST and ZRLD becomes less than a positive fraction, say  > 0. The procedure is outlined in Algorithm 1. Thus far, the methodology employed for solving LD problem is explained. In the next section, we describe how we generate feasible solutions for CASD from the optimal solutions of the subproblems in each step of Algorithm 1. 3.3. Comprising feasible solutions for CASD At each iteration of the Dantzig–Wolfe column generation algorithm we obtain optimal solutions of subproblems CASD1 (θ , λ, β, γ 1 , γ 2 ), CASD2 (μ, λ, β, γ 1 , γ 2 ) and CASD3 (θ ). Period lengths (w), active times of the sensors (a), flows between the sensors (x) and flows between sensors and sinks (y) come from CASD1 (θ , λ, β, γ 1 , γ 2 ). Sensor places (p) and activity schedules of the sensors (q) come from CASD2 (μ, λ, β, γ 1 , γ 2 ). Finally, sink relocation

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

decisions (z) come from CASD3 (θ ). These values are unlikely to form a feasible solution for CASD since they are not supposed to satisfy the relaxed six set of constraints. In this case, the easiest way of obtaining a feasible solution is to start from the binary p and q variables and to make them feasible by the minimum number of changes with respect to coverage constraint (6) and subsequently obtain the values of the continuous variables that satisfy the relaxed constraints. In order to search for the new sensor locations ( p) and activity schedule of the sensors ( q) that are closest to p and q and feasible with respect to the coverage constraint, we develop a new mathematical formulation called Coverage and Budget Model (CBM). To this end, we define δ ir and irt to denote the change in the value of pir and qirt . Formulation of CBM is given below. CBM:

min



δir +

i∈S r∈R

s.t.  

qirt ≥ dk



irt

(53)

i∈S r∈R t∈T

k ∈ K, t ∈ T

(54)

r∈R i:k∈Kir



fir pir ≤ B

(55)

i∈S r∈R

qirt ≤ pir

i ∈ S, r ∈ R, t ∈ T

(56)

δir ≥ pir − pir i ∈ S, r ∈ R

(57)

δir ≥ pir − pir i ∈ S, r ∈ R

(58)

irt ≥ qirt − qirt i ∈ S, r ∈ R, t ∈ T

(59)

irt ≥ qirt − qirt i ∈ S, r ∈ R, t ∈ T

(60)

δir , irt ≥ 0 i ∈ S, r ∈ R, t ∈ T

(61)

pir , qirt ∈ {0, 1}

(62)

i ∈ S, r ∈ R, t ∈ T

Note that pir and qirt values are calculated by solving the subproblems and considered as input parameters in the formulation CBM. In order p and q to be feasible with respect to coverage and budget constraints, constraints (54) and (55) which are similar to constraints (6) and (7) are added to CBM. Constraint (56) makes sure that no sensor can be active in a period without being located at the first hand. Constraints (57) and (58) equate δ ir to 1 if values of pir and pir are different. If they are different, then either the difference pir − pir or the difference pir − pir is equal to 1, and hence δ ir is at least 1. Due to the minimization objective, its value equals to 1 in an optimal solution. On the other hand, if the values of pir and pir are the same, then δ ir ≥ 0 is obtained implying that δ ir will be 0 in an optimal solution since we have a minimization objective. In a similar manner, irt is equated to 1 by constraints (59) and (60) if values of qirt and qirt are different, and to 0 otherwise. p and q values obtained from the optimal solution of CBM and z values obtained from the optimal solution of CASD3 (θ ) are then fixed in CASD, and CASD is solved as an LP only for the continuous w, x and y decision variables. As a result, we obtain and record a feasible solution and its objective value. The feasible solution having the maximum objective function value, which is denoted by ZLH , is reported as the result of the Lagrangian heuristic at the end. It is still possible to end up with an infeasible LP after fixing p, q and z values in CASD which imply that given active sensor and

9

sink locations form a disconnected network for all periods. In such cases, we provide no feasible solution for CASD for the current iteration. However, we see in our computational experiments that these iterations are very rare. 4. Numerical results In this section, we first explain the selection of parameters used in the formulation of CASD and then we indicate the superiority of Lagrangian heuristic over the commercial solver Gurobi and competitors from the literature using a large set of test instances. 4.1. Test problem generation and parameter selection 16 different problem sets with respectively 20, 30, 40, 50, 60, 70, 80, 90, 100, 125, 150, 175, 200, 300, 400 and 500 candidate sensor locations are randomly generated for numerical comparison purposes. Each of the sets include five different instances. Candidate sensor locations are assumed to possess a grid structure. Dimensions of the grid are chosen so that shape of the grid is as close as possible to a square. For instance, grid dimensions for 40 and 90 are respectively 5 × 8 and 9 × 10. Distance between neighboring grid points is chosen as 15 meters. This grid structure formed for the candidate sensor locations is called sensor grid in the sequel. Similar to the sensor grid, candidate sink locations are also assumed to have a grid shape which we call as the sink grid in the rest of the paper. We assume that sensor and sink grids have a nested structure in the sense that they both lay in the same sensor field. The number of candidate sink locations is either the half of the number of candidate sensor locations (if the number of candidate sensor locations is even) or it is half of the five more of the S number of candidate sensor locations otherwise, i.e., |N | = | 2 | if

|S | is even, or |N | =

|S |+5

otherwise. For instance, the number of 2 candidate sink locations is respectively 35 and 65 for 70 and 125 candidate sensor locations. The shape of the sink grid is also chosen so that it is as close as possible to a square. Outermost four corner points of the sink grid are chosen so that they are at the centroid of the outermost four squares of the sensor grid. As a consequence of this selection, vertical and horizontal distances between neighboring candidate sink locations are respectively given by (15 ∗ (n1 − 2 ))/(m1 − 1 ) and (15 ∗ (n2 − 2 ))/(m2 − 1 ) where n1 and n2 are the dimensions of the sensor grid and m1 and m2 are the dimensions of the sink grid. On the other hand, the actual number of sinks to be placed on the candidate sink locations are selected as 3, 5 and 7 for different settings. We assume that there are two types of sensors having respectively 15 and 22 meters sensing, and 50 and 80 meters communication ranges. Their initial battery energies are selected as e1 = 10, 0 0 0 Joules and e2 = 20, 0 0 0 Joules, respectively. We assume that both type of sensors produce the same amount of data per unit time. Hence, hr is selected as 4096 bits per hour for both sensor types. The deployment cost of a type 1 sensor is assumed to be determined randomly according to a uniform distribution between 1 and 10. Namely, fi1 ∼ U(1, 10). On the other hand, deploying a type 2 sensor is assumed to be more expensive and the amount of extra cost is also uniformly distributed between 0 and 5, i.e. fi2 = fi1 + U (0, 5 ). Moreover, total sensor deployment budget B is calculated as the weighted sum of the individual sensor deployment costs over candidate sensor locations. Weights are set  to 0.5 for all candidate sensor locations, i.e., B = 0.5 i∈S ( fi1 + fi2 ). On the other hand, we represent the sensor field by the help of a discrete set of points which we call coverage points and we assume that the coverage points coincide with the candidate sensor locations. They may have been determined in a different manner as well. For instance, the distribution of the coverage points does not have to possess a grid structure, they can be selected randomly.

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

JID: EOR 10

ARTICLE IN PRESS

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

The way chosen for the determination of the coverage points does not have an effect on our analysis. However, in our numerical experiments we take the coverage points and candidate sensor locations the same. For each of the coverage points we randomly assign 1 or 2 as their coverage demands in order to have a differentiated coverage structure. Finally, we assume that both types of sensors share the same energy usage characteristics. As suggested in Heinzelman, Chandrakasan, and Balakrishnan (20 0 0), the energy usage model implemented for the numerical experiments assumes that the energy used for the transmission is a function of the distance between the transmitting and receiving sensors. Hence, the power used by the transmitting sensor increases with the transmission distance. Total energy used during transmission is the sum of two expressions, the first one is a distance independent constant and the other part is linearly proportional to the square of the transmission distance when the path loss index is 2 as assumed in this study. Mathematically, the amount of energy used by a sensor located at point i for a transmission to a sensor located at j is equal to cti j = κ1 + κ2 di2j , where κ 1 , κ 2 and dij are the distance independent constant energy segment, the weight of the distance related part and the Euclidean distance between points i and j, respectively. In our numerical experiments, we select κ 1 as 50 microjoules per bit, and κ 2 as 100 nanojoules per bit square meters. Besides, the amount of energy used for receiving a bit of data is assumed to be independent from the transmission energy. It is implemented as 50 microjoules in our experiments, i.e., cr = 50 microjoules per bit. Finally, the amount of energy used by a sensor for sensing and processing one unit data (e.g. one bit) during a time unit (e.g. an hour) is selected as 50 nanojoules per hour bit, i.e. cs = 50 nanojoules per hour bit. It should be noted that this energy fraction is also responsible for timely and precise coordination between the sensors during the transmission. Hence, we implicitly assume that some data propagation protocols enabling timely and precise data transmissions are already employed for data transmissions between sensors and we represent the load of the employed protocol on each active sensor by cs . We make use of studies by Heinzelman et al. (20 0 0) and Türkog˘ ulları et al. (2010a) for the selection of the parameters values. 4.2. Performance analysis In this section, we assess the performance of the Lagrangian heuristic by comparing with the performance of the state-of-theart MILP solver Gurobi (Gurobi Optimization, 2016) on the described test bed. We code CASD and Lagrangian heuristic in Visual Studio environment by C# language and we carry out all experiments on a single core of a Dell computer having four Intel Xeon 5460 dual core and 28 gigabytes of RAM operating within Windows Server Edition 2003. We first analyze the effect of the strengthening items (given in Section 2.3.2) on the network lifetime for several problem sets. Later, we compare the results found by Lagrangian heuristic with the results found by Gurobi and some competitors from the literature. 4.2.1. Incorporation of stronger and valid constraints In order to assess the effect of the nine possible strengthening items, we run Gurobi for three hours for the base CASD formulation (that is the one given in Section 2.2), and for each of the models obtained by incorporating the listed items one by one for instances with 20, 50 and 100 candidate sensor locations. We report the average network lifetimes found by Gurobi for each set in Table 2. The first column of Table 2 indicates the version of CASD, for instance 0 stands for the base model while 1 represents the case in which the first item listed in Section 2.3.2 is imposed on CASD,

Table 2 Effect of stronger formulations. #

20

50

100

0 1 2 3 4 5 6 7 8 9

95419.89 92903.46 93164.63 91491.19 93415.58 97742.44 98086.43 94544.50 95257.47 95546.40

54054.28 60464.70 59228.43 59759.41 52767.69 52780.47 55239.77 62921.57 60081.23 58801.25

21509.12 23145.49 18015.95 12360.45 23511.6 16702.87 21273.26 36936.57 24631.72 28231.79

so on. Second, third and fourth columns keep the average network lifetimes found by Gurobi in three hours for 20, 50 and 100 candidate sensor locations, respectively. It should be noted that the effect of the items is not the same for each of the test sets. For instance, after imposing item one, average lifetime decreases from 95419.89 to 92903.46 for the instances with 20 candidate sensor locations while it increases from 54054.28 to 60464.70 and from 21509.12 to 23145.49 respectively for the instances with 50 and 100 candidate sensor locations. Similar variability can be observed for other items as well. Since our aim is to improve the performance of Gurobi on CASD formulation especially for relatively large instances, we concentrate on the effect of the items on the instances with 100 candidate sensor locations. In particular, we choose to implement the items that positively affect the lifetime for both of the instances with 50 and 100 candidate sensor locations, i.e., items 1, 7, 8, and 9. If positive effect is observed only for the instances with 100 candidate sensor locations and the negative effect on the lifetime for the instances with 50 candidate sensor locations is less than 5%, then we again choose to implement it, i.e., the effort 4. Therefore, CASD formulation is strengthened by incorporating items 1, 4, 7, 8, and 9 before it is run by Gurobi. 4.2.2. Performance of the Lagrangian heuristic We let Gurobi and Lagrangian heuristic run for at most 3 hours for each of the test instances and the maximum network lifetimes they find in the allotted running time are reported for P = 3 number of sinks in Table 3. We also report the network lifetimes found in three hours if the sink movement decisions are found by simulated annealing heuristic outside and a Lagrangian heuristic is run on the remaining model as done in the study by Keskin et al. (2015). This method is called SA in the sequel. Hence, we adapt the solution approach of Keskin et al. (2015) for CASD that is originally employed for the model with static sinks. The first and seventh columns of the table include the number of candidate sensor and sink locations. The second and eighth columns keep the maximum network lifetimes found by Gurobi for the corresponding test instances. Similarly, the third and ninth columns are responsible from the Lagrangian heuristic results while the fourth and tenth columns report the lifetimes found by SA. Fifth and eleventh columns keep the percent deviation between the Gurobi and LH results while sixth ant twelfth columns are responsible from the percent deviation between the Gurobi and SA results. Z −Z Mathematical formulation for percent deviations are 100 × LHZ G and 100 ×

ZSA −ZG ZSA ,

LH

respectively. Here, ZG , ZLH and ZSA respectively

represent the maximum network lifetimes found by Gurobi, Lagrangian heuristic and SA. As can be understood from Table 3, Gurobi is more successful than Lagrangian heuristic respectively for three and four instances of the problems with 20 and 40 candidate sensor locations which are among the smallest test instances and for one instance

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

11

Table 3 Maximum lifetimes and percent deviations for 3 sinks.

(|K|, |N | )

ZG

ZLH

ZSA

PD1

PD2

(|K|, |N | )

ZG

ZLH

ZSA

PD1

PD2

(20, (20, (20, (20, (20, (30, (30, (30, (30, (30, (40, (40, (40, (40, (40, (50, (50, (50, (50, (50, (60, (60, (60, (60, (60, (70, (70, (70, (70, (70, (80, (80, (80, (80, (80, (90, (90, (90, (90, (90,

93205.79 77547.56 89943.38 100509.65 86955.32 67765.16 77924.09 77803.88 79134.06 89150.70 82356.08 67312.33 63901.23 70123.69 65366.53 54049.07 60218.33 55289.40 54146.34 50356.20 48178.07 52300.90 49516.49 47085.65 52308.09 45192.37 43729.64 50460.68 50289.56 43957.68 40384.42 44007.21 46523.64 40825.14 40318.39 33613.17 42885.50 42403.48 32991.55 39007.04

89510.89 87780.91 88387.90 89510.89 89510.89 76029.88 79264.93 87730.84 76939.63 95204.30 81297.26 65109.67 62409.92 68784.39 65935.16 61808.04 65540.57 58878.16 61053.50 53703.92 49806.47 54459.75 52020.68 52784.39 58335.20 51715.69 51760.80 53137.54 51989.69 46984.31 46514.46 45368.07 43737.28 49247.79 44029.08 45382.49 44986.73 47981.01 42121.04 41752.37

59458.80 53014.27 53204.36 62299.91 53856.67 47764.76 44812.81 47625.76 40467.45 42489.22 43836.76 37933.33 37044.19 39294.25 40601.78 38641.83 38048.32 36268.03 37288.14 33925.30 31699.51 34358.91 33029.38 32945.41 35142.45 27951.89 30229.85 30027.36 30402.00 26378.69 26397.04 23425.29 26286.22 25690.19 24805.48 22176.54 23590.43 24380.65 23016.40 24427.78

−4.13 11.66 −1.76 −12.29 2.86 10.87 1.69 11.32 −2.85 6.36 −1.30 −3.38 −2.39 −1.95 0.86 12.55 8.12 6.10 11.31 6.23 3.27 3.96 4.81 10.80 10.33 12.61 15.52 5.04 3.27 6.44 13.18 3.00 −6.37 17.10 8.43 25.93 4.67 11.62 21.67 6.58

−56.76 −46.28 −69.05 −61.33 −61.46 −41.87 −73.89 −63.37 −95.55 −109.82 −87.87 −77.45 −72.50 −78.46 −60.99 −39.87 −58.27 −52.45 −45.21 −48.43 −51.98 −52.22 −49.92 −42.92 −48.85 −61.68 −44.66 −68.05 −65.42 −66.64 −52.99 −87.86 −76.99 −58.91 −62.54 −51.57 −81.79 −73.92 −43.34 −59.68

(100, 50) (100, 50) (100, 50) (100, 50) (100, 50) (125, 65) (125, 65) (125, 65) (125, 65) (125, 65) (150, 75) (150, 75) (150, 75) (150, 75) (150, 75) (175, 90) (175, 90) (175, 90) (175, 90) (175, 90) (20 0, 10 0) (20 0, 10 0) (20 0, 10 0) (20 0, 10 0) (20 0, 10 0) (300, 150) (300, 150) (300, 150) (300, 150) (300, 150) (40 0, 20 0) (40 0, 20 0) (40 0, 20 0) (40 0, 20 0) (40 0, 20 0) (500, 250) (500, 250) (500, 250) (500, 250) (500, 250)

37709.33 35534.75 37279.61 35901.49 35868.44 24722.43 19753.58 22234.37 18020.52 27662.14 19927.74 28079.55 24407.11 18266.96 27235.63 20243.11 0.00 26700.91 26638.60 0.00 26447.79 24563.88 24374.62 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

40462.08 37781.44 40130.08 39053.64 42077.87 29306.74 33487.02 32727.04 33781.73 32688.58 29915.96 28224.70 32601.14 30934.24 26007.08 26511.03 24922.03 29161.07 25772.54 29502.78 20065.50 25755.57 24840.60 21030.73 22017.86 14812.86 16650.77 17934.76 15145.28 14786.83 14151.12 11909.17 14778.77 12627.39 12774.04 11479.67 9479.30 10239.45 11511.24 10613.36

21226.17 18659.75 21647.43 19929.85 23332.48 21160.71 20948.44 21668.54 19929.85 20651.31 16533.27 17891.69 18538.72 17869.03 18085.48 11171.15 14001.01 13191.84 12392.62 12588.98 12203.50 13654.69 12483.56 12131.21 13757.76 11001.34 11297.39 10733.71 10993.38 11517.22 7811.57 8022.94 7874.26 8009.83 7840.83 5204.43 5308.53 4976.13 5237.18 5250.14

6.80 5.95 7.10 8.07 14.76 15.64 41.01 32.06 46.66 15.38 33.39 0.51 25.13 40.95 −4.72 23.64 100.00 8.44 −3.36 100.00 −31.81 4.63 −0.14 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

−77.65 −90.44 −72.21 −80.14 −53.73 −16.83 5.70 −2.61 9.58 −33.95 −20.53 −56.94 −31.65 −2.23 −50.59 −81.21 100.00 −102.40 −114.96 100.00 −116.72 −79.89 −99.26 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

10) 10) 10) 10) 10) 15) 15) 15) 15) 15) 20) 20) 20) 20) 20) 25) 25) 25) 25) 25) 30) 30) 30) 30) 30) 35) 35) 35) 35) 35) 40) 40) 40) 40) 40) 45) 45) 45) 45) 45)

from the sets with 30, 80, 150, 175 and 200 candidate sensor locations. For all other instances, Lagrangian heuristic produces better results than Gurobi. Therefore, Lagrangian heuristic definitely performs better than Gurobi does for CASD instances. Furthermore, the success of the Lagrangian heuristic becomes even clearer for larger instances as the percent deviation between Gurobi and Lagrangian heuristic increases as the number of candidate sensor locations increases. Moreover, Gurobi is unable to find a network with positive network lifetime for two instances of the networks having 175 and 200 candidate sensor locations and for all instances of the networks having more than 200 candidate sensor locations implying that Gurobi is not a good WSN design alternative for real large size applications. On the other hand, Gurobi results are better than SA results for all instances with less than or equal to 100 candidate sensor locations. However, the gap between the Gurobi and SA results tends to decrease with the increasing problem size. Therefore, SA can be considered as a solution alternative for large instances. Nevertheless, LH is better than SA for all instances of each set. Hence, one should resort to LH if available. A further elaboration reveals the reasons of the clear dominance of LH over SA. First of all, heuristics that search the solution space starting from a local neighborhood like the simulated annealing should be able to assess the quality of the neighboring solutions within a very small amount of time to be able to span a satisfactory number of solutions. However, we run a Lagrangian heuristic for the reduced model (the model left after sink movement decisions are set in CASD) after jumping to a neighboring solution by moving one of the sinks in each iteration

of SA. Hence, assessing time of each neighboring solution is equal to the time that inner Lagrangian heuristic spends to converge a solution, which is too much for a successful simulated annealing application. Moreover, the solution space for the sink movement decisions is |T | times larger than the solution space of the sink locations from the study of Keskin et al. (2015). Since the total computation time is taken as constant, we either search only for a very limited part of the solution space for the sink movement decisions or we cut the inner Lagrangian heuristic before its convergence. Both alternatives give poor sink movement decisions and relatively lower network lifetimes at the end. Besides, feasibility and the quality of the different variable sets of the complicated MILP models like CASD depend on each other. For instance, a sensor cannot be active without being located at the first hand. Similarly, sensor locations and sink movement decisions should be chosen in such a way that each sensor is able to send its data to a sink location. Dependence between the different decisions are regularized by strict constraints. A random change in the value of a variable in an effort to jump to a neighboring solution would probably result in infeasibility in such a framework. Therefore, it takes considerable amount of time to restore feasibility at each jump or to search for a feasible neighbor at each iteration. Thus, heuristic procedures employing random jump strategies are not expected to be successful for the applications where feasible solutions are mostly surrounded by infeasible ones as in the case for CASD. In order to have a clearer picture of the success of the Lagrangian heuristic we compute the average network lifetime of the five instances for each of the test problems with different

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR 12

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

9

×10

4

Gurobi LH PIH SAH SA

8

Network Lifetime (Hours)

7 6 5 4 3 2 1 0

0

2030405060708090100 125 150 175 200 300 Number of Candidate Sensor Locations

400

500

Fig. 1. Maximum network lifetimes by Gurobi, LH, PIH, SAH and SA for 3 sinks.

100 Gurobi vs LH Gurobi vs PIH Gurobi vs SAH Gurobi vs SA

80

Percent Deviation (%)

60 40 20 0 -20 -40 -60 -80

0

2030405060708090100 125 150 175 200

300

400

500

Number of Candidate Sensor Locations Fig. 2. LH, PIH, SAH, and SA vs Gurobi: percent deviations for 3 sinks.

candidate sensor locations for Gurobi, Lagrangian heuristic and SA results. Then, we plot them against the number of candidate sensor locations in Fig. 1. In addition, Fig. 1 contains the average of the PIH and SAH results given in the paper by Keskin et al. (2014) for comparison purposes. Moreover, we also give the plot of percent deviations between Gurobi and PIH, Gurobi and SAH, Gurobi and SA, and Gurobi and LH in Fig. 2. As can be understood from Figs. 1 and 2, LH produces networks with longer lifetimes than its competitors Gurobi, PIH, SAH and SA for most test instances on the basis of the average network lifetimes. SAH produces networks having slightly longer lifetimes for the sets with 200 and 300 candidate sensor locations but its performance rapidly decreases afterward. Moreover, LH is more robust than PIH and SAH in the sense that it produces networks with satisfactorily long lifetimes for very large instances such as the networks with 400 and 500 candidate sensor locations for which PIH and SAH are not able to produce networks with positive lifetimes. On the other hand, SA produces modest results for small sized

instances but it is able to provide similar network lifetimes for very large instances as well. Another important point that should be underlined is that Gurobi is able to keep up with PIH and SAH and able to produce better networks than SA for small and medium sized networks. The improvement between Gurobi results reported in this study and the ones given in the study by Keskin et al. (2014) is mainly due to the strengthening items explained in Section 2.3.2. We also observe from Fig. 1 that the average network lifetime decreases with the increasing number of candidate sensor locations regardless of the employed method. There may be two reasons behind this phenomenon. First of all, as the size of the network gets larger, the amount of data flowing through the relay sensors increases which causes faster depletion of the batteries and consequently, lower lifetimes. Secondly, it is possible to relocate mobile sinks near the critical sensors having relatively lower residual energies in order to keep their energy usage rates at their minimum. However, as the network size increases, the number of

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

ARTICLE IN PRESS

JID: EOR

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14 Table 4 Effect of increase in the number of sinks. P/|S |

3 5 7

150

175

200

Gurobi

LH

Gurobi

LH

Gurobi

LH

23583.40 28433.87 25478.78

29536.62 39345.26 46404.27

14716.52 23274.14 22133.77

27173.89 34462.98 43146.21

15177.26 12477.43 21867.67

22742.05 31111.21 35248.56

critical sensors becomes too high so that it would be impossible to help all of them with the limited number of sinks at hand. In order to justify this claim and to observe the effect of the number of sinks on the network lifetime, we also run Gurobi and LH for 5 and 7 sinks and for 150, 175 and 200 number of candidate sensor locations. We tabulate the average of the five instances for each experimental shift in Table 4 below. We also give results of the 3 sink cases for corresponding test instances for comparison purposes. The results of LH from Table 4 reveals that the number of sinks has an increasing effect on the average network lifetime since the average lifetime found by LH always increases by the increase in the number of sinks. The increase in the number of sinks positively affects the lifetime for most of the Gurobi results as well. Average lifetime that Gurobi finds for the instances respectively with 150 and 175 candidate sensor locations decreases from 28433.87 to 25478.78 and from 23274.14 to 22133.77 as response to an increase in the number of sinks from 5 to 7. In addition, the average lifetime that Gurobi finds for the networks with 200 candidate sensor locations decreases from 15177.26 to 12477.43 after the number of sinks is increased from 3 to 5. Nevertheless, this unexpected results are due to the fact that Gurobi is very sensitive to the increase in the problem size. In other words, although it is expected that an increase in the number of sinks positively affects the average network lifetime, this effect is not always observed in Gurobi results since Gurobi is also negatively affected by the increase in the problem size that occurs as a result of increase in the number of sinks. Hence, LH is more robust than Gurobi in this sense as well since LH results do not heavily depend on the problem size. 5. Conclusions and future research directions This paper presents a column generation heuristic for CASD formulation which is first introduced by Keskin et al. (2014) for integrated modeling of sensor deployment, activity scheduling of the sensors, data and mobile sink routing. Some of the constraints of CASD that avoid decomposition of the model are relaxed and carried to the objective function. This relaxation decomposes CASD into three Lagrangian subproblems. The sum of the optimal objective values of the subproblems for any Lagrangian multiplier values provides an upper bound for the Lagrangian dual value. On the other hand, a Lagrangian Dual formulation (LD) with excessive number of constraints is formulated. A different version of the LD called the Restricted Lagrangian Dual (RLD) includes only a subset of the constraints of LD and the solution of RLD gives a lower bound for the Lagrangian dual value. Upper and lower bounds of the Lagrangian dual value are improved by Dantzig–Wolfe column generation idea which is initially used by Wentges (1997) for the solution of the Lagrangian dual problem. The same method is successfully employed in Klose and Görtz (2007) for the capacitated facility location problem and in Keskin et al. (2015) for WSN design problem with multiple static sinks. A feasible solution of CASD at each iteration of the algorithm is constructed by making use of the solutions of Lagrangian subproblems. A new mathematical model called the Coverage and Budget Model (CBM), which finds the nearest feasible solution of CASD to the subproblem solutions, is constructed. Finally, the success and robustness of the LH results are demonstrated by comparing them with the lifetimes found by

[m5G;December 23, 2016;8:34] 13

the competitor solution methods Gurobi, PIH, SAH and SA on a large set of test instances. This work can be extended by testing the networks found by LH with some real life WSN simulators such as NS-2 (Institute, 2012), TOSSIM (Levis, Lee, Welsh, & Culler, 2003) and OMNeT++ (Varga & Hornig, 2008). Such an attempt would reveal the performance of LH with respect to some other WSN performance criteria like message latency, message loss rate, throughput rates. As another consequence, more realistic parameter values for CASD would be determined. Moreover, other Lagrange multiplier update methods like stabilized column generation of Du Merle, Villeneuve, Desrosiers, and Hansen (1999), the bundle code due to Frangioni (2002) or proximal analytic center cutting plane method of Nesterov and Vial (1999) can be incorporated instead of weighted Dantzig–Wolfe decomposition. References Alsalih, W., Akl, S., & Hassanein, H. (2007). Placement of multiple mobile base stations in wireless sensor networks. In Proceedings of the IEEE international symposium on signal processing and information technology (pp. 229–233). Altınel, I.˙ K., Aras, N., Güney, E., & Ersoy, C. (2008). Binary integer programming formulation and heuristics for differentiated coverage in heterogeneous sensor networks. Computer Networks, 52(12), 2419–2431. Azad, A., & Chockalingam, A. (2006). Mobile base stations placement and energy aware routing in wireless sensor networks. In Proceedings of the IEEE wireless communications and networking conference: 1 (pp. 264–269). Basagni, S., Boloni, L., Gjanci, P., Petrioli, C., Phillips, C. A., & Turgut, D. (2014). Maximizing the value of sensed information in underwater wireless sensor networks via an autonomous underwater vehicle. In Proceedings of the IEEE conference on computer communications, INFOCOM 2014 (pp. 988–996). IEEE. Basagni, S., Carosi, A., Melachrinoudis, E., Petrioli, C., & Wang, Z. M. (2008). Controlled sink mobility for prolonging wireless sensor networks lifetime. Wireless Networks, 14(6), 831–858. Basagni, S., Carosi, A., Petrioli, C., & Phillips, C. A. (2011). Coordinated and controlled mobility of multiple sinks for maximizing the lifetime of wireless sensor networks. Wireless Networks, 17(3), 759–778. Behdani, B., Yun, Y. S., Cole Smith, J., & Xia, Y. (2012). Decomposition algorithms for maximizing the lifetime of wireless sensor networks with mobile sinks. Computers and Operations Research, 39(5), 1054–1061. Castaño, F., Bourreau, E., Velasco, N., Rossi, A., & Sevaux, M. (2015). Exact approaches for lifetime maximization in connectivity constrained wireless multi-role sensor networks. European Journal of Operational Research, 241(1), 28–38. Du Merle, O., Villeneuve, D., Desrosiers, J., & Hansen, P. (1999). Stabilized column generation. Discrete Mathematics, 194(1), 229–237. Fisher, M. L. (1981). The Lagrangian relaxation method for solving integer programming problems. Management Science, 27(1), 1–18. Frangioni, A. (2002). Generalized bundle methods. SIAM Journal on Optimization, 13(1), 117–156. Gandham, S. R., Dawande, M., Prakash, R., & Venkatesan, S. (2003). Energy efficient schemes for wireless sensor networks with multiple mobile base stations. In Proceedings of the IEEE global telecommunications conference: 1 (pp. 377–381). Gatzianas, M., & Georgiadis, L. (2008). A distributed algorithm for maximum lifetime routing in sensor networks with mobile sink. IEEE Transactions on Wireless Communications, 7(3), 984–994. Güney, E., Aras, N., Altınel, I.˙ K., & Ersoy, C. (2010). Efficient integer programming formulations for optimum sink location and routing in heterogeneous wireless sensor networks. Computer Networks, 54(11), 1805–1822. Güney, E., Aras, N., Altınel, I.˙ K., & Ersoy, C. (2012). Efficient solution techniques for the integrated coverage, sink location and routing problem in wireless sensor networks. Computers and Operations Research, 39(7), 1530–1539. Gurobi Optimization, I., 2016. Gurobi optimizer reference manual, http://www. gurobi.com. (Accessed September 2016) Heinzelman, W. R., Chandrakasan, A., & Balakrishnan, H. (20 0 0). Energy-efficient communication protocol for wireless microsensor networks. In Proceedings of the thirty third annual Hawaii international conference on system sciences (pp. 10–19). Institute, I. S. (2012). The network simulator. University of Southern California. http: //www.isi.edu/nsam/ns/. Keskin, M. E. (2014). Controlled sink mobility and wireless sensor network lifetime maximization. University of Bog˘ aziçi. Ph.D. thesis. Keskin, M. E., Altınel, I.˙ K., & Aras, N. (2015). Combining simulated annealing with Lagrangian relaxation and weighted Dantzig–Wolfe decomposition for integrated design decisions in wireless sensor networks. Computers and Operations Research, 59, 132–143. Keskin, M. E., Altınel, I.˙ K., Aras, N., & Ersoy, C. (2011). Lifetime maximization in wireless sensor networks using a mobile sink with nonzero traveling time. The Computer Journal, 54(12), 1987–1999. Keskin, M. E., Altınel, I.˙ K., Aras, N., & Ersoy, C. (2014). Wireless sensor network lifetime maximization by optimal sensor deployment, activity scheduling, data routing and sink mobility. Ad Hoc Networks, 17, 18–36.

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006

JID: EOR 14

ARTICLE IN PRESS

[m5G;December 23, 2016;8:34]

M.E. Keskin / European Journal of Operational Research 000 (2016) 1–14

Keskin, M. E., Altınel, I. K., Aras, N., & Ersoy, C. (2016). Wireless sensor network design by lifetime maximisation: An empirical evaluation of integrating major design issues and sink mobility. International Journal of Sensor Networks, 20(3), 131–146. Klose, A., & Görtz, S. (2007). A branch-and-price algorithm for the capacitated facility location problem. European Journal of Operational Research, 179(3), 1109–1125. Lersteau, C., Rossi, A., & Sevaux, M. (2016). Robust scheduling of wireless sensor networks for target tracking under uncertainty. European Journal of Operational Research, 252(2), 407–417. Levis, P., Lee, N., Welsh, M., & Culler, D. (2003). Tossim: Accurate and scalable simulation of entire tiny OS applications. In Proceedings of the first international conference on embedded networked sensor systems (pp. 126–137). Li, J., & Mohapatra, P. (2007). Analytical modeling and mitigation techniques for the energy hole problem in sensor networks. Pervasive and Mobile Computing, 3(3), 233–254. Liang, W., Luo, J., & Xu, X. (2010). Prolonging network lifetime via a controlled mobile sink in wireless sensor networks. In Proceedings of the IEEE global telecommunications conference (pp. 1–6). Luo, J., & Hubaux, J.-P. (2005). Joint mobility and routing for lifetime elongation in wireless sensor networks. In Proceedings of the twenty fourth annual joint conference of the IEEE computer and communications societies: 3 (pp. 1735–1746). Luo, J., & Hubaux, J.-P. (2010). Joint sink mobility and routing to maximize the lifetime of wireless sensor networks: the case of constrained mobility. IEEE/ACM Transactions on Networking, 18(3), 871–884. Madan, R., & Lall, S. (2006). Distributed algorithms for maximum lifetime routing in wireless sensor networks. IEEE Transactions on Wireless Communications, 5(8), 2185–2193. Martello, S., Pisinger, D., & Toth, P. (1999). Dynamic programming and strong bounds for the 0–1 knapsack problem. Management Science, 45(3), 414–424. Nesterov, Y., & Vial, J. P. (1999). Homogeneous analytic center cutting plane methods for convex problems and variational inequalities. SIAM Journal on Optimization, 9(3), 707–728. Papadimitriou, I., & Georgiadis, L. (2005). Maximum lifetime routing to mobile sink in wireless sensor networks. In Proceedings of the thirteenth IEEE international conference on software in telecommunications and computer networks, SoftCom.

Popa, L., Rostamizadeh, A., Karp, R., Papadimitriou, C., & Stoica, I. (2007). Balancing traffic load in wireless networks with curveball routing. In Proceedings of the international symposium on mobile ad hoc networking and computing: 9 (pp. 170–179). Türkog˘ ulları, Y. B., Aras, N., & Altınel, I.˙ K. (2010a). Optimal placement, scheduling, and routing to maximize lifetime in sensor networks. Journal of the Operational Research Society, 61(6), 10 0 0–1012. Türkog˘ ulları, Y. B., Aras, N., Altınel, I.˙ K., & Ersoy, C. (2010b). A column generation based heuristic for sensor placement, activity scheduling and data routing in wireless sensor networks. European Journal of Operational Research, 207(2), 1014–1026. Türkog˘ ulları, Y. B., Aras, N., Altınel, I.˙ K., & Ersoy, C. (2010c). An efficient heuristic for placement, scheduling and routing in wireless sensor networks. Ad Hoc Networks, 8(6), 654–667. Varga, A., & Hornig, R. (2008). An overview of the OMNeT++ simulation environment. In Proceedings of the first international conference on simulation tools and techniques for communications, networks and systems and workshops (pp. 60–69). Wang, Z. M., Basagni, S., Melachrinoudis, E., & Petrioli, C. (2005). Exploiting sink mobility for maximizing sensor networks lifetime. In Proceedings of the thirty eighth annual Hawaii international conference on system sciences. 287a–287a. Wentges, P. (1997). Weighted Dantzig–Wolfe decomposition for linear mixed-integer programming. International Transactions in Operational Research, 4(2), 151–162. Wu, X., Chen, G., & Das, S. K. (2008). Avoiding energy holes in wireless sensor networks with nonuniform node distribution. IEEE Transactions on Parallel and Distributed Systems, 19(5), 710–720. Xia, L., & Shihada, B. (2015). A jackson network model and threshold policy for joint optimization of energy and delay in multi-hop wireless networks. European Journal of Operational Research, 242(3), 778–787. Yick, J., Mukherjee, B., & Ghosal, D. (2008). Wireless sensor network survey. Computer Networks, 52(12), 2292–2330. Yun, Y., & Xia, Y. (2010). Maximizing the lifetime of wireless sensor networks with mobile sink in delay-tolerant applications. IEEE Transactions on Mobile Computing, 9(9), 1308–1318. Yun, Y., Xia, Y., Behdani, B., & Smith, J. C. (2010). Distributed algorithm for lifetime maximization in delay-tolerant wireless sensor network with mobile sink. In Proceedings of the forty ninth IEEE conference on decision and control (pp. 370–375).

Please cite this article as: M.E. Keskin, A column generation heuristic for optimal wireless sensor network design with mobile sinks, European Journal of Operational Research (2016), http://dx.doi.org/10.1016/j.ejor.2016.12.006