Sensor Distributed Location through Linear Programming with Triangle Inequality Constraints
National Institute of Standards and Technology Wireless Communication Technologies Group Email: [e-mail]
Camillo Gentile
Abstract Interest in dense sensor networks due to falling price and reduced size has motivated research in sensor location in recent years. To our knowledge, the algorithm which achieves the best performance in sensor location solves an optimization program by relaxing the quadratic geometrical constraints of the network to render the program convex. In recent work we proposed solving the same program, however by applying convex geometrical constraints directly, necessitating no relaxation of the constraints and in turn ensuring a tighter solution. This paper proposes a distributed version of our algorithm which achieves the same globally optimal objective function as the decentralized version. We conduct extensive experimentation to substantiate the robustness of our algorithm even in the presence of high levels of noise, and report the messaging overhead for convergence. Index Terms Simplex method, Primal-dual method
This renders the problem automatically convex, necessitating no relaxation of the constraints and so guaranteeing a tighter solution, as conrmed through simulation. This paper proposes a distributed algorithm for sensor location which converges to the optimal solution to the problem dened in [8]. Section II states this problem formally. The primal-dual method explained in Section III claims a key advantage over the conventional simplex method to solving the problem in a decentralized fashion. Exploiting this advantage leads to our distributed version of the primal-dual method described in Section IV. We demonstrate the steps of this distributed algorithm through a simple example network in Section V. An extensive number of challenging tests conditions are reported in Section VI to substantiate the robustness of our algorithm to high levels of noise in comparison to the algorithm proposed by Biswas. We also report the messaging overhead of the algorithm. The last section provides conclusions and directions for further research.
I. I NTRODUCTION The falling price and reduced size of sensors in recent years have fueled the deployability of dense networks to monitor and relay environmental properties such as temperature, moisture, and light [1]. The ability to self-organize and nd their locations autonomously and with high accuracy proves particularly useful in military and public safety operations. In dense networks, multilateration can render good location accuracy despite signicant errors in range estimates between sensors. This has launched a research area known as sensor location which seeks to process potentially enormous quantities of data collectively to achieve optimal results. Most practical systems require local distributed processing to cope with dynamic links or nodes in motion to maintain a network updated; alternatively relaying information across a large network sanctions the centralized processing of obsolete data, limiting scalability. A recent paper on sensor location [2] provides an exhaustive survey of the available techniques for sensor location [3], [4], [6], [5]. To our knowledge, the two algorithms achieving the best performance in sensor location formulate a program with quadratic constraints to minimize a linear objective function [2], [7]. Since some of the constraints are non-convex, the papers differ primarily in their relaxation approaches to render the problem convex. The solution provided by Biswas et al. has greater applicability and yields better results than the one by Doherty et al. In recent work [8], we formulated a novel problem following their same approach, maintaining the efciency of convex optimization, however by applying linear triangle inequality constraints as opposed to quadratic ones.
II. P RELIMINARIES
Consider a network with two types of nodes: nA anchor nodes (or anchors) with known location and nS sensor nodes (or sensors) with unknown location, for a total of n = nA + nS nodes. For simplicity, let the nodes lie on a plane such that node i has location xi R2 indexed through i, i = 1 . . . nA for the anchors and i = nA + 1 . . . n for the sensors. The set N contains all pairs of nodes between which a link exists: (i, j), i < j; ||xi xj || < R, where || || denotes the Euclidean distance and the network parameter R is known as the radio range. The set M contains all triplets of nodes which form a triangle in the network: (i, j, k), (i, j) N ; (j, k) N ; (i, k) N . Neighboring nodes i and j measure the link distance dij between them through received-signal-strength or time-ofarrival techniques [9]. Given the locations of the anchor nodes and the measured distances between neighboring nodes in the network, this paper considers the following problem to solve for the locations of the sensors:
min
(i,j)N
|ij | (1)
min
(i,j)N
+ + ij ij
s.t. a1 (+ - ) + a1 (+ - ) + a1 (+ - ) - 1 = b1 ij ij ij,k jk,i ik,j ijk ijk jk jk ik ik a2 (+ - ) + a2 (+ - ) + a2 (+ - ) - 2 = b2 , ij,k ij ij jk,i ik,j ijk ijk jk jk ik ik + + + 3 3 3 3 3 aij,k (ij -ij ) + ajk,i (jk - jk ) + aik,j (ik - ik ) - ijk = bijk (i, j, k) M where dij = dij + ij . + 0 ij , (i, j) N The problem minimizes the sum of the absolute residuals 0 ij ij between the measured distances dij and the estimated (4) distances dij such that the latter conform to requisite geowhere metrical constraints. Using triangle inequality constraints as 1 opposed to quadratic constraints [2], [7] ensures the convexity 1 aij,k a1 a1 bijk 1 1 -1 -1 -1 1 dij jk,i ik,j of the problem without relaxing any of the original geometrical a2 a2 a2 = 1 - 1 1 , b2 = -1 1 - 1 djk , ij,k jk,i ik,j ijk constraints. Rewriting the problem in standard form removes 3 3 3 3 -1 1 1 1 -1 - 1 dik aij,k ajk,i aik,j bijk the absolute values and introduces bounding constraints: min
(i,j)N
s.t. dij + djk dik dij + dik djk , (i, j, k) M djk + dik dij
+ + ij ij
s.t. dij + djk dik dij + dik djk , (i, j, k) M djk + dik dij + 0 ij 0 ij , (i, j) N
(2)
where ij = + . The solution to the problem above does ij ij not directly yield the sensor locations xi , i = nA + 1...n, but only the values of the residuals of the link distances. Hence the complete algorithm requires aposteriori location propagation described in [8] to furnish the locations of the sensors from the residuals. Theorem: If any one of the three inequality constraints of a triangle is bound, then the other two are feasible. Proof: Assume without loss of generality that the rst inequality constraint is bound: dij + dik = djk , dij + djk = dik djk + dik = dij dij + dik dij + dik since dij , djk 0, but djk + dik djk + dik dij + dik djk .2 so djk + dik dij (3)
+ , are the 2|N | primal decision variables, and ij ij 1 , 2 , 3 are the 3|M | primal slack variables. ijk ijk ijk A basic solution to the primal problem contains exactly 3|M | nonzero (basic) variables and 2|N | zero (nonbasic) variables for the nongenerate case which assumes linear independence of all constraints in (4). Since the optimal solution is basic [10], the simplex method pivots between basic feasible solutions of the system to nd it. A pivot consists of raising an entering variable in the nonbasic set from zero that will improve the objective function. The entering variable can rise only a certain amount until a blocking variable in the basic set reduces to zero; hence the entering variable becomes basic and the blocking variable nonbasic at another basic feasible solution of the system. Determining this amount necessitates global knowledge of the values of all the basic variables, such that no variable loses feasibility and exactly 2|N | of them equal zero throughout the pivots. Hence the simplex method does not lend to distributed processing. Interior-point methods are suited for centralized computing even more than the simplex method, requiring the inversion of large sparse symmetric matrices [11]. B. The dual problem Each primal linear program has a unique dual linear program. The dual problem to (4) appears as [12]: max s.t. + : ij : ij 1 : ijk 2 : ijk 3 : ijk
(i,j,k)M 3 2 1 b1 ijk +b2 ijk + b3 ijk ijk ijk ijk
k, (i,j,k)M k, (i,j,k)M
III. T HE A. The primal problem
PRIMAL - DUAL METHOD
+ 1 2 3 a1 ijk + a2 ijk + a3 ijk +ij = 1 ij,k ij,k ij,k
1 2 3 - a1 ijk - a2 ijk -a3 ijk + ij = 1 ij,k ij,k ij,k
,
We denote the linear program (2) as the primal problem. Rewriting the primal in canonical form appears as
1 ijk 0 2 ijk 0 , (i, j, k) M 3 ijk 0
(i, j) N
(5)
1 2 3 where ijk , ijk , ijk are the dual decision variables and + ij , ij are the dual slack variables. As indicated, each primal variable has a corresponding dual constraint. The Complementary Slackness Theorem [10] states that any feasible primal and dual solutions are optimal if the following complementary slackness conditions hold:
l jk
i i j d jk k
l jk
i
j d jk
j, d jk
k k
(6) Rather than solve the primal problem through the simplex method, we formulate a distributed version of the primal-dual method. The key advantage of the latter relaxes the condition that a primal solution be basic. Our algorithm proceeds in the follow manner: rst a link in the network nds a local feasible primal solution (not necessarily basic) such that all incident triangles meet the triangle inequality constraints. Then the link applies the complementary slackness conditions given through this primal solution, dening the local restricted dual problem. If the dual solution to the restricted problem is also feasible, then the local primal solution is optimal through the Complementary Slackness Theorem; otherwise the primal solution is modied. Once all the links attain compatible locally compatible optimal solutions, the network achieves the globally optimal solution. Dantzig treats a full discussion on the primal-dual method [11]. IV. D ISTRIBUTED A. Network organization The nodes in the network transmit asynchronously. If node i wakes up after node j, then ni manages the link ij . The link manager maintains the information on ij : the measured distance dij and the residual ij initialized to zero. Fig. 1a illustrates an example network with three nodes ni , nj , nk , where nk woke up rst and nj second, assigning to nj manager of jk . When ni wakes up subsequently, it broadcasts a HELLO message containing its ID#: i. In Fig. 1b, nodes nj and nk respond with their ID#s; since nj serves as link manager, it also broadcasts the information on the links it manages. In receiving the messages from nj and nk in Fig. 1c, ni becomes manager of ij and ik and estimates dij and ik ; the two-way message exchange allows ni to measure these d distances asynchronously [13]. As manager, ni then broadcasts the information on the links it manages. Now both managers in Fig. 1d have access to information on all three links of ijk . While we refer to the links as the processing centers for the distributed algorithm in the sequel, the actual processing of course takes place at the link managers. B. A feasible primal solution Denote a triangle ijk as feasible if it meets all three of its triangle inequality constraints: u 0, u {1, 2, 3}. ijk
LOCATION
(a) (b) (c) (d)
u ijk u ijk v ij v ij
=0 >0 =0 >0
u ijk 0 u ijk = 0 v ij 0 v ij = 0
, ,
(i, j, k) M u {1, 2, 3} (i, j) N v { +, - }
(a)
l ij l ik
(b)
d ij d jk d ik
l ij
i j, d ij ,d ik
l ik
i
d ij d jk d ik
l jk
j d jk k
l jk
j d ij d jk d ik k
(c)
Fig. 1. Transmitted messages in network organization.
(d)
Suppose that a link ij changes value, rendering one of its incident triangles infeasible: u < 0, for u, k. Proof (3) ijk shows that if any one of the constraints is bound, then the other two are feasible; so by setting u = 0, ij restores feasibility ijk to ijk . Since the value of ij was just changed, it remains the same; rather the link selects one of the other two links on the triangle ( jk or ik ) to set the value of u to zero, say jk . ijk
jk Consider modifying v such that jk = jk + u u , ijk jk ijk u u where ijk = ijk represents the necessary change in u ijk to render the violated constraint feasible. The partial
v
v jk u ijk
ij ,ik
=
1 vau jk,i
(7)
while maintaining the value of the links ij and ik constant is computed by rewriting the primal constraint u in (4) as
vv jk
au (+ - ) + au (+ - ) + au (+ - ) - u = bu . jk,i ik,j ijk ijk ij,k ij ij jk jk ik ik (8) -v Note that v > 0 implies jk = 0. The program below jk summarizes the local pivot:
ij ik
Pivot I: Set u = 0 by modifying v ijk jk jk = jk +
v jk u ijk v jk u ijk
u ijk
(9)
u ijk
ij ,ik = u ijk
=
vau jk,i
1
In restoring feasibility to ijk , jk may in turn render a neighboring triangle jkl infeasible, analagous to the change in ij which rendered ijk infeasible in the previous step. Through this mechanism infeasibility propagates through the network between triangles, leaving those in its path feasible, and so obtaining a feasible primal solution locally at each step until termination at a certain triangle. In the worst case, propagation terminates at the edge of the network where jk has no neighboring jkl . Section IV C describes how to select jk optimally.
6 6
(6d), while unbounding the remaining dual constraints (6c). As any solution to a nondegenerate linear system contains at least the same number of unknowns as equations, any solution to the nondegenerate primal system in (4) contains at least the same number of nonzero primal decision variables as the number of bound primal constraints. So due to complementary slackness, the restricted dual contains no more unknowns than equations, allowing to solve for the unknowns through simple substitution.
6
1
2
12
2 123
23
3
1 234
2 245
5
4
24
1
2
3
3
4
2
5
5
9 8
(a)
4
6
4 2
1
2
3
3
1 1
3
5
7 8
(b)
4
6
4
Fig. 3.
Propagation in solving the restricted dual problem.
1
Fig. 2.
Propagation in nding a feasible primal solution.
The portion of a network in Fig. 2a consists of four feasible triangles (shaded) with the estimated distances d = d displayed on each link. Let d12 decrease to 1, rendering 123 infeasible. The infeasibility then propagates along the path indicated by the dashed arrow in Fig. 2b: d23 decreases to 3, restoring feasibility to 123 , but rendering 234 infeasible; d24 decreases to 7, restoring feasibility to 234 while maintaining 245 feasible. The propagation terminates at 245 with all four triangles newly feasible. C. The restricted dual problem 1) Dening and solving the restricted dual problem: We say that a link ij has a local feasible primal solution if all its incident triangles ijk , k, (i, j, k) M are feasible. Once a link obtains a local feasible primal solution, it applies the complementary slackness conditions which dene the local restricted dual: every bound primal constraint u = 0 admits ijk u one dual decision variable (unknown) ijk 0 to the restricted problem (6a), while setting the remaining dual decision variables to zero (6b); every nonzero primal decision variable v v > 0 admits one bound dual constraint (equation) ij = 0 ij
Fig. 3 graphically represents the restricted dual problem corresponding to the feasible primal solution in Fig. 2b. Each of the three darkly shaded triangles has one bound primal constraint, admitting one unknown per triangle to the dual: 2 1 123 (d12 + d13 = d23 ) at 123 , 234 (d23 + d34 = d24 ) 2 at 234 , and 245 (d24 + d25 = d45 ) at 245 ; each of the three boldfaced links with nonzero residuals ( = 3, = 12 23 2, = 2) admits one equation per link to the dual from 24 2 2 1 (5): (123 = 1) at 12 , (123 234 = 1) at 23 , and 1 2 (234 245 = 1) at 24 . Although 23 and 24 cannot solve for their two unknowns, 12 can solve for its single unknown 1 123 = 1, and then propagates the value to the other two 2 links of 123 ; now knowing the value of 123 , l23 solves 1 for its single remaining unknown 234 = 2 and propagates the value to the other two links of 234 ; now knowing the 1 value of 234 , 24 solves for its single remaining unknown 2 245 = 1 and propagates the value to the other two links of 245 . Paralleling the primal solution, the dual solution is also found through propagation. The dashed arrow in Fig. 3 indicates the direction of propagation. Modifying 2) the primal solution towards optimality: If the solution to the restricted dual is feasible (i.e. all the decision and slack variables are greater than or equal to zero), then the primal solution is optimal. Otherwise it can be shown [11] that setting an infeasible dual variable to zero improves the primal objective, provided that the nonzero feasible dual variables remain admitted to the restricted dual problem. If the restricted dual solution includes an infeasible decision u variable ijk < 0, then raising its corresponding primal slack u variable u from zero sets ijk = 0 through complementary ijk slackness. However entering u lowers a local blocking ijk
ijk variable u > 0, u = u, u = 1; conversely, reducing ijk
u
u ijk
to zero through Pivot I (9) raises u ijk
u ijk ij ,ik
ijk
u . ijk
The partial au jk,i
au jk,i
=
v jk
u ijk
u ijk v jk
=
(10)
is computed through (7). Note that Pivot I results in v > 0, jk v and in turn sets jk = 0 through complementary slackness; v so if jk > 0 before the pivot, the pivot removes a nonzero feasible dual variable from the restricted dual. Select v such jk v that jk 0. If the restricted dual solution includes an infeasible slack v variable ij < 0, then raising its corresponding primal v decision variable v from zero sets ij = 0 through comij plementary slackness. However entering v lowers a local ij
ij blocking variable v > 0, v = 1; conversely, reducing jk ij
the sensor location is known, it serves with another known sensor (or anchor) to determine the location of an unknown sensor neighboring the two. Fig. 4 displays two anchor nodes (shaded) and four sensor nodes. Anchors n1 and n2 propagate their locations to unknown sensor n3 ; anchor n2 and now known sensor n3 propagate their locations to unknown sensor n4 , anchor n2 and now known sensor n4 propagate their locations to unknown sensor n5 , and known sensor n3 and known sensor n4 propagate their locations to unknown sensor n6 . The location propagation is embedded with the primaldual pivots in our distributed location algorithm. The details as described in [8] and are omitted for space reasons.
6
v
1 3 4 2
v to zero through Pivot II below raises v . The partial ij jk is computed through (8). Pivot II: Set v = 0 by raising v ij jk
v ij v jk
ij = ij +
v ij v ,u =0 jk ik ijk v = v jk jk
v ij v jk
5
v , jk v au
ij,k
(11)
Fig. 4. Propagation in determining the sensor locations.
= vajk,i , u
Note that Pivot II affects the value of u , u {1, 2, 3}; ijk u so if ijk > 0 before the pivot, maintain u = 0 such that ijk the nonzero feasible dual variable remains in the restricted problem. Recall in Section IV-B that the initial primal feasible solution is found in the absence of any dual restrictions, and so the arbitrary selection of the link to modify in Pivot I may lead to a sub-optimal solution. When modifying the primal solution towards optimality, the local pivot may affect the primal feasibility of a non-incident triangle on the link. And so another link on that triangle maneuvers to restore its feasibility through Pivot I, but now using the local dual restrictions to decide which one. The new primal feasible solution of that triangle in turn generates other dual restrictions local to the selected link. As the pivots continue, the global primal problem becomes more and more restricted by the local duals until achieving global optimality. An example provided in the following section substantiates these ideas. D. Location propagation The estimated distances computed through (4) yield the desired sensor locations through location propagation [14]. The anchor nodes propagate their locations to the sensor nodes in a distributed fashion: if two anchor nodes share a neighboring sensor node, the sensor location can be determined from the two anchor locations coupled with the two estimated distances between the anchors and the sensor. Once
V. E XAMPLE NETWORK Consider the example network in Fig. 0 of Table I with four nodes and ve links. The measured distances d appear on each link. Table I displays the coefcients of the primal problem corresponding to this example network in A and b. The blank slots indicate zeros and are left so to reduce clutter. As the nodes wake up asynchronously, assume that links 12 and 34 are the last two established, completing the triangles 123 and 234 respectively. Note that 123 and 1 3 234 are initially infeasible with 123 = 2 and 134 = 1. In solving the initial primal problem and in the absence of any dual restrictions, 12 raises its value to 7(+ = 6) to 12 restore feasibility to 123 , and distributes this value to the other two links; likewise 34 raises its value to 2(+ = 1) 34 to restore feasibility to 234 , and distributes this value to the other two links. Now that each link ij in the network holds the primal decision variables of links jk and ik for all incident triangles ijk , it can compute the primal slack variables 1 , 2 , 3 . The initial (1) primal solution with ijk ijk ijk objective 7 appears in row v (1) and column u (1) of Table ij ijk I. The solution is indexed according to global solutions for the sake of clarity, but as just explained, the solutions are computed locally, distributedly, and asynchronously. With all the triangles incident on 12 now feasible, 12 applies the complementary slackness conditions from its local primal solution to dene the corresponding local restricted dual problem:
2
1 5
4
1
2
+ 7 12
1
2
+ 3 12
1
2
+ 2 12
1
2
3
3
2 3 123
3
5
3 134
3
1 2 123
5
3 134
3
1 2 123
4
3 134
13
3
4
1
+ 12 1 1 1
+ 2 34
4
3
+ 2 34
4
3
1
(0)
A 1 123 2 123 3 123 1 234 2 234 3 234 v (1) ij v ij (1) v (2) ij v ij (2) v (3) ij v ij (3) 12 1 1 1 + 23 1 1 1 23 1 1 1 + 13 1 1 1 1 1 1 3 2 0 2 2 2 3 13 1 1 1 1 1 1 1
(1)
+ 34 34 + 14 14 b 2 4 6 7 3 1
(2)
u (1) ijk 4 10 8 2 1
u ijk (1)
(3)
u (2) ijk 4 4 8 2 1
u ijk (2) 1
u (3) ijk 4 4 6 2
u ijk (3) 1
-1
6 2 2 2 1 2 0
0
1 1 1 1 1
1 1 1 2
1 1 1 0 0
1 1 1 2 2
0
1 1 1
2 1 1
TABLE I
1
T HE PRIMAL - DUAL TABLE
3 3 = 0 123 0 123 1 2 {1 , 2 } > 0 {123 , 123 } = 0 123 123 + + { , + , , + , } = 0 {12 , 23 , 23 , 13 , 13 } 0 12 23 23 13 13 + a+ > 0 12 = 0 12 3 So 12 processes a single unknown (123 0) and a single 3 equation (123 = 1) to solve for its local restricted dual. 3 It solves for 123 = 1 and distributes this value to the other two links of 123 . Through the same process, 34 3 solves for 234 = 1 in its restricted dual and distributes this value to the other two links of the triangle. Now that each link ij in the network holds the dual decision variables 1 2 3 ijk , ijk , ijk of all incident triangles ijk , it can compute + its slack dual variables ij , ij . The initial (1) dual solution u v appears in column ijk (1) and row ij (1) of Table I. The table evidences the complimentary slackness structure of the primal-dual solution, where the dual slack (decision) variable is zero if a primal decision (slack) variable is nonzero. The boxes indicate the admitted unknowns and equations in each restricted dual, and Fig. 1 graphically represents the restricted dual problem. The initial dual solution reveals two infeasible variables 3 123 = 1 and 13 = 1. Raising 3 or will improve 123 13
3 that 134 = 1 remains admitted to the restricted dual problem. This local pivot however raises 1 from zero, violating a 123 1 non-local dual restriction by removing 123 = 1 from the 1 same restricted dual. Link 13 sets 123 back to zero through Pivot I by modifying the value of 12 to 2 (+ = 1). The 12 third primal-dual solution with objective 2 appears in Table I, where all the feasible dual variables indicate the optimality of this primal solution.
the primal objective, provided that the nonzero feasible variables remain admitted to the restricted dual problem. Link 12 raises 3 through Pivot I, setting the local blocking variable 123 3 123 1 = 4, 1 = 1 to zero by modifying the value of 12 to 123 123 + 3 (+ = 2). Selecting to remove 12 = 0 from the restricted 12 dual problem ensures that this pivot improves the primal objective. The second primal solution remains feasible after Pivot I, necessitating no additional pivots to restore feasibility. Note that the primal objective equal to 3 has indeed improved. The second primal-dual solution appears in Table I. This dual solution still reveals the infeasible variable 13 = 1. Link 13 raises 13 through Pivot II, setting the local blocking variable 13 a+ = 1, + = 1 to zero while maintaining 3 = 0 such 134 34
34
VI. E XPERIMENTAL SETUP AND RESULTS In our previous work [8], we compared our centralized algorithm to Biswas by conducting experiments on a network with the same structure. The network contains 50 sensor nodes uniformly distributed throughout a one by one unit area. The three varying parameters are the number of anchor nodes, the radio range, and the noisy factor of the link distances. As Biswas, the ground-truth link distances dij between neighboring nodes i and j are perturbed with zero-mean unitvariance Gaussian noise N (0,1) and the varying parameter noise. So the link managers measure the noisy link distances dij = dij (1 + N (0, 1) noise). Fig. 5(a) illustrates a test network with three anchors, R = 0.25, and noise = 0.1. The anchors and sensors appear as dark and light asterisks respectively, and the links as dark lines between neighboring nodes. The network contains 225 links for an average node connectivity of 7.9623. The distributed location algorithm yields the estimated locations of the sensor nodes upon convergence. The true and estimated locations appear in Fig.5(b) as dark and light asterisks connected by an error line. The average location error is 0.0597. Biswas reports the results of a single trial network for the six test condition...